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, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
39 #include "insert-molecules.h"
47 #include "gromacs/commandline/cmdlineoptionsmodule.h"
48 #include "gromacs/fileio/confio.h"
49 #include "gromacs/fileio/filetypes.h"
50 #include "gromacs/fileio/xvgr.h"
51 #include "gromacs/gmxlib/conformation-utilities.h"
52 #include "gromacs/gmxpreprocess/read-conformation.h"
53 #include "gromacs/math/functions.h"
54 #include "gromacs/math/utilities.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/options/basicoptions.h"
57 #include "gromacs/options/filenameoption.h"
58 #include "gromacs/options/ioptionscontainer.h"
59 #include "gromacs/pbcutil/pbc.h"
60 #include "gromacs/random/threefry.h"
61 #include "gromacs/random/uniformrealdistribution.h"
62 #include "gromacs/selection/nbsearch.h"
63 #include "gromacs/selection/selection.h"
64 #include "gromacs/selection/selectioncollection.h"
65 #include "gromacs/selection/selectionoption.h"
66 #include "gromacs/selection/selectionoptionbehavior.h"
67 #include "gromacs/topology/atomprop.h"
68 #include "gromacs/topology/atoms.h"
69 #include "gromacs/topology/atomsbuilder.h"
70 #include "gromacs/topology/mtop_util.h"
71 #include "gromacs/topology/topology.h"
72 #include "gromacs/trajectory/trajectoryframe.h"
73 #include "gromacs/utility/cstringutil.h"
74 #include "gromacs/utility/exceptions.h"
75 #include "gromacs/utility/fatalerror.h"
76 #include "gromacs/utility/smalloc.h"
80 /* enum for random rotations of inserted solutes */
82 en_rotXYZ, en_rotZ, en_rotNone
84 const char *const cRotationEnum[] = {"xyz", "z", "none"};
86 static void center_molecule(std::vector<RVec> *x)
88 const size_t atomCount = x->size();
91 for (size_t i = 0; i < atomCount; ++i)
93 rvec_inc(center, (*x)[i]);
95 svmul(1.0/atomCount, center, center);
96 for (size_t i = 0; i < atomCount; ++i)
98 rvec_dec((*x)[i], center);
102 static void generate_trial_conf(const std::vector<RVec> &xin,
103 const rvec offset, RotationType enum_rot,
104 gmx::DefaultRandomEngine * rng,
105 std::vector<RVec> *xout)
107 gmx::UniformRealDistribution<real> dist(0, 2.0*M_PI);
110 real alfa = 0.0, beta = 0.0, gamma = 0.0;
123 alfa = beta = gamma = 0.;
126 if (enum_rot == en_rotXYZ || enum_rot == en_rotZ)
128 rotate_conf(xout->size(), as_rvec_array(xout->data()), nullptr, alfa, beta, gamma);
130 for (size_t i = 0; i < xout->size(); ++i)
132 rvec_inc((*xout)[i], offset);
136 static bool isInsertionAllowed(gmx::AnalysisNeighborhoodSearch *search,
137 const std::vector<real> &exclusionDistances,
138 const std::vector<RVec> &x,
139 const std::vector<real> &exclusionDistances_insrt,
140 const t_atoms &atoms,
141 const std::set<int> &removableAtoms,
142 gmx::AtomsRemover *remover)
144 gmx::AnalysisNeighborhoodPositions pos(x);
145 gmx::AnalysisNeighborhoodPairSearch pairSearch = search->startPairSearch(pos);
146 gmx::AnalysisNeighborhoodPair pair;
147 while (pairSearch.findNextPair(&pair))
149 const real r1 = exclusionDistances[pair.refIndex()];
150 const real r2 = exclusionDistances_insrt[pair.testIndex()];
151 if (pair.distance2() < gmx::square(r1 + r2))
153 if (removableAtoms.count(pair.refIndex()) == 0)
157 // TODO: If molecule information is available, this should ideally
158 // use it to remove whole molecules.
159 remover->markResidue(atoms, pair.refIndex(), true);
165 static void insert_mols(int nmol_insrt, int ntry, int seed,
166 real defaultDistance, real scaleFactor,
167 t_atoms *atoms, t_symtab *symtab, std::vector<RVec> *x,
168 const std::set<int> &removableAtoms,
169 const t_atoms &atoms_insrt, const std::vector<RVec> &x_insrt,
170 int ePBC, matrix box,
171 const std::string &posfn, const rvec deltaR,
172 RotationType enum_rot)
174 fprintf(stderr, "Initialising inter-atomic distances...\n");
175 gmx_atomprop_t aps = gmx_atomprop_init();
176 std::vector<real> exclusionDistances(
177 makeExclusionDistances(atoms, aps, defaultDistance, scaleFactor));
178 const std::vector<real> exclusionDistances_insrt(
179 makeExclusionDistances(&atoms_insrt, aps, defaultDistance, scaleFactor));
180 gmx_atomprop_destroy(aps);
182 const real maxInsertRadius
183 = *std::max_element(exclusionDistances_insrt.begin(),
184 exclusionDistances_insrt.end());
185 real maxRadius = maxInsertRadius;
186 if (!exclusionDistances.empty())
188 const real maxExistingRadius
189 = *std::max_element(exclusionDistances.begin(),
190 exclusionDistances.end());
191 maxRadius = std::max(maxInsertRadius, maxExistingRadius);
194 // TODO: Make all of this exception-safe.
195 gmx::AnalysisNeighborhood nb;
196 nb.setCutoff(maxInsertRadius + maxRadius);
201 seed = static_cast<int>(gmx::makeRandomSeed());
203 fprintf(stderr, "Using random seed %d\n", seed);
205 gmx::DefaultRandomEngine rng(seed);
208 set_pbc(&pbc, ePBC, box);
210 /* With -ip, take nmol_insrt from file posfn */
211 double **rpos = nullptr;
212 const bool insertAtPositions = !posfn.empty();
213 if (insertAtPositions)
216 nmol_insrt = read_xvg(posfn.c_str(), &rpos, &ncol);
219 gmx_fatal(FARGS, "Expected 3 columns (x/y/z coordinates) in file %s\n",
222 fprintf(stderr, "Read %d positions from file %s\n\n",
223 nmol_insrt, posfn.c_str());
226 gmx::AtomsBuilder builder(atoms, symtab);
227 gmx::AtomsRemover remover(*atoms);
229 const int finalAtomCount = atoms->nr + nmol_insrt * atoms_insrt.nr;
230 const int finalResidueCount = atoms->nres + nmol_insrt * atoms_insrt.nres;
231 builder.reserve(finalAtomCount, finalResidueCount);
232 x->reserve(finalAtomCount);
233 exclusionDistances.reserve(finalAtomCount);
236 std::vector<RVec> x_n(x_insrt.size());
242 gmx::UniformRealDistribution<real> dist;
244 while (mol < nmol_insrt && trial < ntry*nmol_insrt)
246 // cppcheck 1.72 complains about uninitialized variables in the
247 // assignments below otherwise...
249 if (!insertAtPositions)
251 // Insert at random positions.
252 offset_x[XX] = box[XX][XX] * dist(rng);
253 offset_x[YY] = box[YY][YY] * dist(rng);
254 offset_x[ZZ] = box[ZZ][ZZ] * dist(rng);
258 // Skip a position if ntry trials were not successful.
259 if (trial >= firstTrial + ntry)
261 fprintf(stderr, " skipped position (%.3f, %.3f, %.3f)\n",
262 rpos[XX][mol], rpos[YY][mol], rpos[ZZ][mol]);
268 // Insert at positions taken from option -ip file.
269 offset_x[XX] = rpos[XX][mol] + deltaR[XX]*(2 * dist(rng)-1);
270 offset_x[YY] = rpos[YY][mol] + deltaR[YY]*(2 * dist(rng)-1);
271 offset_x[ZZ] = rpos[ZZ][mol] + deltaR[ZZ]*(2 * dist(rng)-1);
273 fprintf(stderr, "\rTry %d", ++trial);
276 generate_trial_conf(x_insrt, offset_x, enum_rot, &rng, &x_n);
277 gmx::AnalysisNeighborhoodPositions pos(*x);
278 gmx::AnalysisNeighborhoodSearch search = nb.initSearch(&pbc, pos);
279 if (isInsertionAllowed(&search, exclusionDistances, x_n, exclusionDistances_insrt,
280 *atoms, removableAtoms, &remover))
282 x->insert(x->end(), x_n.begin(), x_n.end());
283 exclusionDistances.insert(exclusionDistances.end(),
284 exclusionDistances_insrt.begin(),
285 exclusionDistances_insrt.end());
286 builder.mergeAtoms(atoms_insrt);
289 fprintf(stderr, " success (now %d atoms)!\n", builder.currentAtomCount());
293 fprintf(stderr, "\n");
294 /* print number of molecules added */
295 fprintf(stderr, "Added %d molecules (out of %d requested)\n",
296 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",
306 originalResidueCount - atoms->nres,
307 originalAtomCount - atoms->nr);
312 for (int i = 0; i < DIM; ++i)
326 class InsertMolecules : public ICommandLineOptionsModule, public ITopologyProvider
330 : bBox_(false), nmolIns_(0), nmolTry_(10), seed_(0),
331 defaultDistance_(0.105), scaleFactor_(0.57), enumRot_(en_rotXYZ),
332 top_(nullptr), ePBC_(-1)
338 virtual ~InsertMolecules()
347 // From ITopologyProvider
348 virtual gmx_mtop_t *getTopology(bool /*required*/) { return top_; }
349 virtual int getAtomCount() { return 0; }
351 // From ICommandLineOptionsModule
352 virtual void init(CommandLineModuleSettings * /*settings*/)
355 virtual void initOptions(IOptionsContainer *options,
356 ICommandLineOptionsModuleSettings *settings);
357 virtual void optionsFinished();
363 SelectionCollection selections_;
365 std::string inputConfFile_;
366 std::string insertConfFile_;
367 std::string positionFile_;
368 std::string outputConfFile_;
374 real defaultDistance_;
377 RotationType enumRot_;
378 Selection replaceSel_;
381 std::vector<RVec> x_;
386 void InsertMolecules::initOptions(IOptionsContainer *options,
387 ICommandLineOptionsModuleSettings *settings)
389 const char *const desc[] = {
390 "[THISMODULE] inserts [TT]-nmol[tt] copies of the system specified in",
391 "the [TT]-ci[tt] input file. The insertions take place either into",
392 "vacant space in the solute conformation given with [TT]-f[tt], or",
393 "into an empty box given by [TT]-box[tt]. Specifying both [TT]-f[tt]",
394 "and [TT]-box[tt] behaves like [TT]-f[tt], but places a new box",
395 "around the solute before insertions. Any velocities present are",
398 "It is possible to also insert into a solvated configuration and",
399 "replace solvent atoms with the inserted atoms. To do this, use",
400 "[TT]-replace[tt] to specify a selection that identifies the atoms",
401 "that can be replaced. The tool assumes that all molecules in this",
402 "selection consist of single residues: each residue from this",
403 "selection that overlaps with the inserted molecules will be removed",
404 "instead of preventing insertion.",
406 "By default, the insertion positions are random (with initial seed",
407 "specified by [TT]-seed[tt]). The program iterates until [TT]-nmol[tt]",
408 "molecules have been inserted in the box. Molecules are not inserted",
409 "where the distance between any existing atom and any atom of the",
410 "inserted molecule is less than the sum based on the van der Waals",
411 "radii of both atoms. A database ([TT]vdwradii.dat[tt]) of van der",
412 "Waals radii is read by the program, and the resulting radii scaled",
413 "by [TT]-scale[tt]. If radii are not found in the database, those"
414 "atoms are assigned the (pre-scaled) distance [TT]-radius[tt].",
416 "A total of [TT]-nmol[tt] * [TT]-try[tt] insertion attempts are made",
417 "before giving up. Increase [TT]-try[tt] if you have several small",
418 "holes to fill. Option [TT]-rot[tt] specifies whether the insertion",
419 "molecules are randomly oriented before insertion attempts.",
421 "Alternatively, the molecules can be inserted only at positions defined in",
422 "positions.dat ([TT]-ip[tt]). That file should have 3 columns (x,y,z),",
423 "that give the displacements compared to the input molecule position",
424 "([TT]-ci[tt]). Hence, if that file should contain the absolute",
425 "positions, the molecule must be centered on (0,0,0) before using",
426 "[THISMODULE] (e.g. from [gmx-editconf] [TT]-center[tt]).",
427 "Comments in that file starting with # are ignored. Option [TT]-dr[tt]",
428 "defines the maximally allowed displacements during insertial trials.",
429 "[TT]-try[tt] and [TT]-rot[tt] work as in the default mode (see above)."
432 settings->setHelpText(desc);
434 std::shared_ptr<SelectionOptionBehavior> selectionOptionBehavior(
435 new SelectionOptionBehavior(&selections_, this));
436 settings->addOptionsBehavior(selectionOptionBehavior);
438 // TODO: Replace use of legacyType.
439 options->addOption(FileNameOption("f")
440 .legacyType(efSTX).inputFile()
441 .store(&inputConfFile_)
442 .defaultBasename("protein")
443 .description("Existing configuration to insert into"));
444 options->addOption(FileNameOption("ci")
445 .legacyType(efSTX).inputFile().required()
446 .store(&insertConfFile_)
447 .defaultBasename("insert")
448 .description("Configuration to insert"));
449 options->addOption(FileNameOption("ip")
450 .filetype(eftGenericData).inputFile()
451 .store(&positionFile_)
452 .defaultBasename("positions")
453 .description("Predefined insertion trial positions"));
454 options->addOption(FileNameOption("o")
455 .legacyType(efSTO).outputFile().required()
456 .store(&outputConfFile_)
457 .defaultBasename("out")
458 .description("Output configuration after insertion"));
460 options->addOption(SelectionOption("replace").onlyAtoms()
462 .description("Atoms that can be removed if overlapping"));
463 selectionOptionBehavior->initOptions(options);
465 options->addOption(RealOption("box").vector()
466 .store(newBox_).storeIsSet(&bBox_)
467 .description("Box size (in nm)"));
468 options->addOption(IntegerOption("nmol")
470 .description("Number of extra molecules to insert"));
471 options->addOption(IntegerOption("try")
473 .description("Try inserting [TT]-nmol[tt] times [TT]-try[tt] times"));
474 options->addOption(IntegerOption("seed")
476 .description("Random generator seed (0 means generate)"));
477 options->addOption(RealOption("radius")
478 .store(&defaultDistance_)
479 .description("Default van der Waals distance"));
480 options->addOption(RealOption("scale")
481 .store(&scaleFactor_)
482 .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."));
483 options->addOption(RealOption("dr").vector()
485 .description("Allowed displacement in x/y/z from positions in [TT]-ip[tt] file"));
486 options->addOption(EnumOption<RotationType>("rot").enumValue(cRotationEnum)
488 .description("Rotate inserted molecules randomly"));
491 void InsertMolecules::optionsFinished()
493 if (nmolIns_ <= 0 && positionFile_.empty())
495 GMX_THROW(InconsistentInputError("Either -nmol must be larger than 0, "
496 "or positions must be given with -ip."));
498 if (inputConfFile_.empty() && !bBox_)
500 GMX_THROW(InconsistentInputError("When no solute (-f) is specified, "
501 "a box size (-box) must be specified."));
503 if (replaceSel_.isValid() && inputConfFile_.empty())
505 GMX_THROW(InconsistentInputError("Replacement (-replace) only makes sense "
506 "together with an existing configuration (-f)."));
510 if (!inputConfFile_.empty())
512 readConformation(inputConfFile_.c_str(), top_, &x_, nullptr,
513 &ePBC_, box_, "solute");
514 if (top_->natoms == 0)
516 fprintf(stderr, "Note: no atoms in %s\n", inputConfFile_.c_str());
521 int InsertMolecules::run()
523 std::set<int> removableAtoms;
524 if (replaceSel_.isValid())
527 set_pbc(&pbc, ePBC_, box_);
530 fr->natoms = x_.size();
532 fr->x = as_rvec_array(x_.data());
533 selections_.evaluate(fr, &pbc);
535 removableAtoms.insert(replaceSel_.atomIndices().begin(),
536 replaceSel_.atomIndices().end());
537 // TODO: It could be nice to check that removableAtoms contains full
538 // residues, since we anyways remove whole residues instead of
546 box_[XX][XX] = newBox_[XX];
547 box_[YY][YY] = newBox_[YY];
548 box_[ZZ][ZZ] = newBox_[ZZ];
552 gmx_fatal(FARGS, "Undefined solute box.\nCreate one with gmx editconf "
553 "or give explicit -box command line option");
556 t_topology *top_insrt;
557 std::vector<RVec> x_insrt;
562 readConformation(insertConfFile_.c_str(), top_insrt, &x_insrt,
563 nullptr, &ePBC_dummy, box_dummy, "molecule");
564 if (top_insrt->atoms.nr == 0)
566 gmx_fatal(FARGS, "No molecule in %s, please check your input",
567 insertConfFile_.c_str());
569 if (top_->name == nullptr)
571 top_->name = top_insrt->name;
573 if (positionFile_.empty())
575 center_molecule(&x_insrt);
579 // TODO: Adapt to use mtop throughout.
580 t_atoms atoms = gmx_mtop_global_atoms(top_);
582 /* add nmol_ins molecules of atoms_ins
583 in random orientation at random place */
584 insert_mols(nmolIns_, nmolTry_, seed_, defaultDistance_, scaleFactor_,
585 &atoms, &top_->symtab, &x_, removableAtoms, top_insrt->atoms, x_insrt,
586 ePBC_, box_, positionFile_, deltaR_, enumRot_);
588 /* write new configuration to file confout */
589 fprintf(stderr, "Writing generated configuration to %s\n",
590 outputConfFile_.c_str());
591 write_sto_conf(outputConfFile_.c_str(), *top_->name, &atoms,
592 as_rvec_array(x_.data()), nullptr, ePBC_, box_);
594 /* print size of generated configuration */
595 fprintf(stderr, "\nOutput configuration contains %d atoms in %d residues\n",
596 atoms.nr, atoms.nres);
607 const char InsertMoleculesInfo::name[] = "insert-molecules";
608 const char InsertMoleculesInfo::shortDescription[] =
609 "Insert molecules into existing vacancies";
610 ICommandLineOptionsModulePointer InsertMoleculesInfo::create()
612 return ICommandLineOptionsModulePointer(new InsertMolecules());