2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
40 #include "insert_molecules.h"
48 #include "gromacs/commandline/cmdlineoptionsmodule.h"
49 #include "gromacs/fileio/confio.h"
50 #include "gromacs/fileio/filetypes.h"
51 #include "gromacs/fileio/xvgr.h"
52 #include "gromacs/gmxlib/conformation_utilities.h"
53 #include "gromacs/gmxpreprocess/makeexclusiondistances.h"
54 #include "gromacs/math/functions.h"
55 #include "gromacs/math/utilities.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/options/basicoptions.h"
58 #include "gromacs/options/filenameoption.h"
59 #include "gromacs/options/ioptionscontainer.h"
60 #include "gromacs/pbcutil/pbc.h"
61 #include "gromacs/random/threefry.h"
62 #include "gromacs/random/uniformrealdistribution.h"
63 #include "gromacs/selection/nbsearch.h"
64 #include "gromacs/selection/selection.h"
65 #include "gromacs/selection/selectioncollection.h"
66 #include "gromacs/selection/selectionoption.h"
67 #include "gromacs/selection/selectionoptionbehavior.h"
68 #include "gromacs/topology/atomprop.h"
69 #include "gromacs/topology/atoms.h"
70 #include "gromacs/topology/atomsbuilder.h"
71 #include "gromacs/topology/mtop_util.h"
72 #include "gromacs/topology/symtab.h"
73 #include "gromacs/topology/topology.h"
74 #include "gromacs/trajectory/trajectoryframe.h"
75 #include "gromacs/utility/cstringutil.h"
76 #include "gromacs/utility/exceptions.h"
77 #include "gromacs/utility/fatalerror.h"
78 #include "gromacs/utility/smalloc.h"
82 /* enum for random rotations of inserted solutes */
89 const char* const cRotationEnum
[] = { "xyz", "z", "none" };
91 static void center_molecule(gmx::ArrayRef
<RVec
> x
)
98 svmul(1.0 / x
.size(), center
, center
);
101 rvec_dec(xi
, center
);
105 static void generate_trial_conf(gmx::ArrayRef
<RVec
> xin
,
107 RotationType enum_rot
,
108 gmx::DefaultRandomEngine
* rng
,
109 std::vector
<RVec
>* xout
)
111 gmx::UniformRealDistribution
<real
> dist(0, 2.0 * M_PI
);
112 xout
->assign(xin
.begin(), xin
.end());
114 real alfa
= 0.0, beta
= 0.0, gamma
= 0.0;
126 case en_rotNone
: alfa
= beta
= gamma
= 0.; break;
128 if (enum_rot
== en_rotXYZ
|| enum_rot
== en_rotZ
)
130 rotate_conf(xout
->size(), as_rvec_array(xout
->data()), nullptr, alfa
, beta
, gamma
);
132 for (size_t i
= 0; i
< xout
->size(); ++i
)
134 rvec_inc((*xout
)[i
], offset
);
138 static bool isInsertionAllowed(gmx::AnalysisNeighborhoodSearch
* search
,
139 const std::vector
<real
>& exclusionDistances
,
140 const std::vector
<RVec
>& x
,
141 const std::vector
<real
>& exclusionDistances_insrt
,
142 const t_atoms
& atoms
,
143 const std::set
<int>& removableAtoms
,
144 gmx::AtomsRemover
* remover
)
146 gmx::AnalysisNeighborhoodPositions
pos(x
);
147 gmx::AnalysisNeighborhoodPairSearch pairSearch
= search
->startPairSearch(pos
);
148 gmx::AnalysisNeighborhoodPair pair
;
149 while (pairSearch
.findNextPair(&pair
))
151 const real r1
= exclusionDistances
[pair
.refIndex()];
152 const real r2
= exclusionDistances_insrt
[pair
.testIndex()];
153 if (pair
.distance2() < gmx::square(r1
+ r2
))
155 if (removableAtoms
.count(pair
.refIndex()) == 0)
159 // TODO: If molecule information is available, this should ideally
160 // use it to remove whole molecules.
161 remover
->markResidue(atoms
, pair
.refIndex(), true);
167 static void insert_mols(int nmol_insrt
,
170 real defaultDistance
,
174 std::vector
<RVec
>* x
,
175 const std::set
<int>& removableAtoms
,
176 const t_atoms
& atoms_insrt
,
177 gmx::ArrayRef
<RVec
> x_insrt
,
180 const std::string
& posfn
,
182 RotationType enum_rot
)
184 fprintf(stderr
, "Initialising inter-atomic distances...\n");
186 std::vector
<real
> exclusionDistances(makeExclusionDistances(atoms
, &aps
, defaultDistance
, scaleFactor
));
187 const std::vector
<real
> exclusionDistances_insrt(
188 makeExclusionDistances(&atoms_insrt
, &aps
, defaultDistance
, scaleFactor
));
190 const real maxInsertRadius
=
191 *std::max_element(exclusionDistances_insrt
.begin(), exclusionDistances_insrt
.end());
192 real maxRadius
= maxInsertRadius
;
193 if (!exclusionDistances
.empty())
195 const real maxExistingRadius
=
196 *std::max_element(exclusionDistances
.begin(), exclusionDistances
.end());
197 maxRadius
= std::max(maxInsertRadius
, maxExistingRadius
);
200 // TODO: Make all of this exception-safe.
201 gmx::AnalysisNeighborhood nb
;
202 nb
.setCutoff(maxInsertRadius
+ maxRadius
);
207 seed
= static_cast<int>(gmx::makeRandomSeed());
209 fprintf(stderr
, "Using random seed %d\n", seed
);
211 gmx::DefaultRandomEngine
rng(seed
);
214 set_pbc(&pbc
, pbcType
, box
);
216 /* With -ip, take nmol_insrt from file posfn */
217 double** rpos
= nullptr;
218 const bool insertAtPositions
= !posfn
.empty();
219 if (insertAtPositions
)
222 nmol_insrt
= read_xvg(posfn
.c_str(), &rpos
, &ncol
);
225 gmx_fatal(FARGS
, "Expected 3 columns (x/y/z coordinates) in file %s\n", posfn
.c_str());
227 fprintf(stderr
, "Read %d positions from file %s\n\n", nmol_insrt
, posfn
.c_str());
230 gmx::AtomsBuilder
builder(atoms
, symtab
);
231 gmx::AtomsRemover
remover(*atoms
);
233 const int finalAtomCount
= atoms
->nr
+ nmol_insrt
* atoms_insrt
.nr
;
234 const int finalResidueCount
= atoms
->nres
+ nmol_insrt
* atoms_insrt
.nres
;
235 builder
.reserve(finalAtomCount
, finalResidueCount
);
236 x
->reserve(finalAtomCount
);
237 exclusionDistances
.reserve(finalAtomCount
);
240 std::vector
<RVec
> x_n(x_insrt
.size());
246 gmx::UniformRealDistribution
<real
> dist
;
248 while (mol
< nmol_insrt
&& trial
< ntry
* nmol_insrt
)
251 if (!insertAtPositions
)
253 // Insert at random positions.
254 offset_x
[XX
] = box
[XX
][XX
] * dist(rng
);
255 offset_x
[YY
] = box
[YY
][YY
] * dist(rng
);
256 offset_x
[ZZ
] = box
[ZZ
][ZZ
] * dist(rng
);
260 // Skip a position if ntry trials were not successful.
261 if (trial
>= firstTrial
+ ntry
)
263 fprintf(stderr
, " skipped position (%.3f, %.3f, %.3f)\n", rpos
[XX
][mol
],
264 rpos
[YY
][mol
], rpos
[ZZ
][mol
]);
270 // Insert at positions taken from option -ip file.
271 offset_x
[XX
] = rpos
[XX
][mol
] + deltaR
[XX
] * (2 * dist(rng
) - 1);
272 offset_x
[YY
] = rpos
[YY
][mol
] + deltaR
[YY
] * (2 * dist(rng
) - 1);
273 offset_x
[ZZ
] = rpos
[ZZ
][mol
] + deltaR
[ZZ
] * (2 * dist(rng
) - 1);
275 fprintf(stderr
, "\rTry %d", ++trial
);
278 generate_trial_conf(x_insrt
, offset_x
, enum_rot
, &rng
, &x_n
);
279 gmx::AnalysisNeighborhoodPositions
pos(*x
);
280 gmx::AnalysisNeighborhoodSearch search
= nb
.initSearch(&pbc
, pos
);
281 if (isInsertionAllowed(&search
, exclusionDistances
, x_n
, exclusionDistances_insrt
, *atoms
,
282 removableAtoms
, &remover
))
284 x
->insert(x
->end(), x_n
.begin(), x_n
.end());
285 exclusionDistances
.insert(exclusionDistances
.end(), exclusionDistances_insrt
.begin(),
286 exclusionDistances_insrt
.end());
287 builder
.mergeAtoms(atoms_insrt
);
290 fprintf(stderr
, " success (now %d atoms)!\n", builder
.currentAtomCount());
294 fprintf(stderr
, "\n");
295 /* print number of molecules added */
296 fprintf(stderr
, "Added %d molecules (out of %d requested)\n", mol
- failed
, nmol_insrt
);
298 const int originalAtomCount
= atoms
->nr
;
299 const int originalResidueCount
= atoms
->nres
;
300 remover
.refreshAtomCount(*atoms
);
301 remover
.removeMarkedElements(x
);
302 remover
.removeMarkedAtoms(atoms
);
303 if (atoms
->nr
< originalAtomCount
)
305 fprintf(stderr
, "Replaced %d residues (%d atoms)\n", originalResidueCount
- atoms
->nres
,
306 originalAtomCount
- atoms
->nr
);
311 for (int i
= 0; i
< DIM
; ++i
)
325 class InsertMolecules
: public ICommandLineOptionsModule
, public ITopologyProvider
333 defaultDistance_(0.105),
342 // From ITopologyProvider
343 gmx_mtop_t
* getTopology(bool /*required*/) override
{ return &top_
; }
344 int getAtomCount() override
{ return 0; }
346 // From ICommandLineOptionsModule
347 void init(CommandLineModuleSettings
* /*settings*/) override
{}
348 void initOptions(IOptionsContainer
* options
, ICommandLineOptionsModuleSettings
* settings
) override
;
349 void optionsFinished() override
;
353 SelectionCollection selections_
;
355 std::string inputConfFile_
;
356 std::string insertConfFile_
;
357 std::string positionFile_
;
358 std::string outputConfFile_
;
364 real defaultDistance_
;
367 RotationType enumRot_
;
368 Selection replaceSel_
;
371 std::vector
<RVec
> x_
;
376 void InsertMolecules::initOptions(IOptionsContainer
* options
, ICommandLineOptionsModuleSettings
* settings
)
378 const char* const desc
[] = {
379 "[THISMODULE] inserts [TT]-nmol[tt] copies of the system specified in",
380 "the [TT]-ci[tt] input file. The insertions take place either into",
381 "vacant space in the solute conformation given with [TT]-f[tt], or",
382 "into an empty box given by [TT]-box[tt]. Specifying both [TT]-f[tt]",
383 "and [TT]-box[tt] behaves like [TT]-f[tt], but places a new box",
384 "around the solute before insertions. Any velocities present are",
387 "It is possible to also insert into a solvated configuration and",
388 "replace solvent atoms with the inserted atoms. To do this, use",
389 "[TT]-replace[tt] to specify a selection that identifies the atoms",
390 "that can be replaced. The tool assumes that all molecules in this",
391 "selection consist of single residues: each residue from this",
392 "selection that overlaps with the inserted molecules will be removed",
393 "instead of preventing insertion.",
395 "By default, the insertion positions are random (with initial seed",
396 "specified by [TT]-seed[tt]). The program iterates until [TT]-nmol[tt]",
397 "molecules have been inserted in the box. Molecules are not inserted",
398 "where the distance between any existing atom and any atom of the",
399 "inserted molecule is less than the sum based on the van der Waals",
400 "radii of both atoms. A database ([TT]vdwradii.dat[tt]) of van der",
401 "Waals radii is read by the program, and the resulting radii scaled",
402 "by [TT]-scale[tt]. If radii are not found in the database, those",
403 "atoms are assigned the (pre-scaled) distance [TT]-radius[tt].",
404 "Note that the usefulness of those radii depends on the atom names,",
405 "and thus varies widely with force field.",
407 "A total of [TT]-nmol[tt] * [TT]-try[tt] insertion attempts are made",
408 "before giving up. Increase [TT]-try[tt] if you have several small",
409 "holes to fill. Option [TT]-rot[tt] specifies whether the insertion",
410 "molecules are randomly oriented before insertion attempts.",
412 "Alternatively, the molecules can be inserted only at positions defined in",
413 "positions.dat ([TT]-ip[tt]). That file should have 3 columns (x,y,z),",
414 "that give the displacements compared to the input molecule position",
415 "([TT]-ci[tt]). Hence, if that file should contain the absolute",
416 "positions, the molecule must be centered on (0,0,0) before using",
417 "[THISMODULE] (e.g. from [gmx-editconf] [TT]-center[tt]).",
418 "Comments in that file starting with # are ignored. Option [TT]-dr[tt]",
419 "defines the maximally allowed displacements during insertial trials.",
420 "[TT]-try[tt] and [TT]-rot[tt] work as in the default mode (see above)."
423 settings
->setHelpText(desc
);
425 std::shared_ptr
<SelectionOptionBehavior
> selectionOptionBehavior(
426 new SelectionOptionBehavior(&selections_
, this));
427 settings
->addOptionsBehavior(selectionOptionBehavior
);
429 // TODO: Replace use of legacyType.
430 options
->addOption(FileNameOption("f")
433 .store(&inputConfFile_
)
434 .defaultBasename("protein")
435 .description("Existing configuration to insert into"));
436 options
->addOption(FileNameOption("ci")
440 .store(&insertConfFile_
)
441 .defaultBasename("insert")
442 .description("Configuration to insert"));
443 options
->addOption(FileNameOption("ip")
444 .filetype(eftGenericData
)
446 .store(&positionFile_
)
447 .defaultBasename("positions")
448 .description("Predefined insertion trial positions"));
449 options
->addOption(FileNameOption("o")
453 .store(&outputConfFile_
)
454 .defaultBasename("out")
455 .description("Output configuration after insertion"));
458 SelectionOption("replace").onlyAtoms().store(&replaceSel_
).description("Atoms that can be removed if overlapping"));
459 selectionOptionBehavior
->initOptions(options
);
461 options
->addOption(RealOption("box").vector().store(newBox_
).storeIsSet(&bBox_
).description(
462 "Box size (in nm)"));
463 options
->addOption(IntegerOption("nmol").store(&nmolIns_
).description(
464 "Number of extra molecules to insert"));
465 options
->addOption(IntegerOption("try").store(&nmolTry_
).description(
466 "Try inserting [TT]-nmol[tt] times [TT]-try[tt] times"));
467 options
->addOption(IntegerOption("seed").store(&seed_
).description(
468 "Random generator seed (0 means generate)"));
470 RealOption("radius").store(&defaultDistance_
).description("Default van der Waals distance"));
472 RealOption("scale").store(&scaleFactor_
).description("Scale factor to multiply Van der Waals radii from the database in share/gromacs/top/vdwradii.dat. The default value of 0.57 yields density close to 1000 g/l for proteins in water."));
473 options
->addOption(RealOption("dr").vector().store(deltaR_
).description(
474 "Allowed displacement in x/y/z from positions in [TT]-ip[tt] file"));
475 options
->addOption(EnumOption
<RotationType
>("rot")
476 .enumValue(cRotationEnum
)
478 .description("Rotate inserted molecules randomly"));
481 void InsertMolecules::optionsFinished()
483 if (nmolIns_
<= 0 && positionFile_
.empty())
486 InconsistentInputError("Either -nmol must be larger than 0, "
487 "or positions must be given with -ip."));
489 if (inputConfFile_
.empty() && !bBox_
)
492 InconsistentInputError("When no solute (-f) is specified, "
493 "a box size (-box) must be specified."));
495 if (replaceSel_
.isValid() && inputConfFile_
.empty())
498 InconsistentInputError("Replacement (-replace) only makes sense "
499 "together with an existing configuration (-f)."));
502 if (!inputConfFile_
.empty())
504 bool bTprFileWasRead
;
505 rvec
* temporaryX
= nullptr;
506 fprintf(stderr
, "Reading solute configuration\n");
507 readConfAndTopology(inputConfFile_
.c_str(), &bTprFileWasRead
, &top_
, &pbcType_
, &temporaryX
,
509 x_
.assign(temporaryX
, temporaryX
+ top_
.natoms
);
511 if (top_
.natoms
== 0)
513 fprintf(stderr
, "Note: no atoms in %s\n", inputConfFile_
.c_str());
518 int InsertMolecules::run()
520 std::set
<int> removableAtoms
;
521 if (replaceSel_
.isValid())
524 set_pbc(&pbc
, pbcType_
, box_
);
527 fr
->natoms
= x_
.size();
529 fr
->x
= as_rvec_array(x_
.data());
530 selections_
.evaluate(fr
, &pbc
);
532 removableAtoms
.insert(replaceSel_
.atomIndices().begin(), replaceSel_
.atomIndices().end());
533 // TODO: It could be nice to check that removableAtoms contains full
534 // residues, since we anyways remove whole residues instead of
538 PbcType pbcTypeForOutput
= pbcType_
;
541 pbcTypeForOutput
= PbcType::Xyz
;
543 box_
[XX
][XX
] = newBox_
[XX
];
544 box_
[YY
][YY
] = newBox_
[YY
];
545 box_
[ZZ
][ZZ
] = newBox_
[ZZ
];
550 "Undefined solute box.\nCreate one with gmx editconf "
551 "or give explicit -box command line option");
554 gmx_mtop_t topInserted
;
555 t_atoms atomsInserted
;
556 std::vector
<RVec
> xInserted
;
558 bool bTprFileWasRead
;
559 PbcType pbcType_dummy
;
562 readConfAndTopology(insertConfFile_
.c_str(), &bTprFileWasRead
, &topInserted
, &pbcType_dummy
,
563 &temporaryX
, nullptr, box_dummy
);
564 xInserted
.assign(temporaryX
, temporaryX
+ topInserted
.natoms
);
566 atomsInserted
= gmx_mtop_global_atoms(&topInserted
);
567 if (atomsInserted
.nr
== 0)
569 gmx_fatal(FARGS
, "No molecule in %s, please check your input", insertConfFile_
.c_str());
571 if (top_
.name
== nullptr)
573 top_
.name
= topInserted
.name
;
575 if (positionFile_
.empty())
577 center_molecule(xInserted
);
581 t_atoms atoms
= gmx_mtop_global_atoms(&top_
);
583 /* add nmol_ins molecules of atoms_ins
584 in random orientation at random place */
585 insert_mols(nmolIns_
, nmolTry_
, seed_
, defaultDistance_
, scaleFactor_
, &atoms
, &top_
.symtab
,
586 &x_
, removableAtoms
, atomsInserted
, xInserted
, pbcTypeForOutput
, box_
,
587 positionFile_
, deltaR_
, enumRot_
);
589 /* write new configuration to file confout */
590 fprintf(stderr
, "Writing generated configuration to %s\n", outputConfFile_
.c_str());
591 write_sto_conf(outputConfFile_
.c_str(), *top_
.name
, &atoms
, as_rvec_array(x_
.data()), nullptr,
592 pbcTypeForOutput
, box_
);
594 /* print size of generated configuration */
595 fprintf(stderr
, "\nOutput configuration contains %d atoms in %d residues\n", atoms
.nr
, atoms
.nres
);
598 done_atom(&atomsInserted
);
606 const char* InsertMoleculesInfo::name()
608 static const char* name
= "insert-molecules";
612 const char* InsertMoleculesInfo::shortDescription()
614 static const char* shortDescription
= "Insert molecules into existing vacancies";
615 return shortDescription
;
618 ICommandLineOptionsModulePointer
InsertMoleculesInfo::create()
620 return ICommandLineOptionsModulePointer(new InsertMolecules());