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, 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"
44 #include "gromacs/commandline/cmdlineoptionsmodule.h"
45 #include "gromacs/fileio/confio.h"
46 #include "gromacs/fileio/filenm.h"
47 #include "gromacs/fileio/xvgr.h"
48 #include "gromacs/gmxlib/conformation-utilities.h"
49 #include "gromacs/gmxpreprocess/read-conformation.h"
50 #include "gromacs/math/utilities.h"
51 #include "gromacs/math/vec.h"
52 #include "gromacs/options/basicoptions.h"
53 #include "gromacs/options/filenameoption.h"
54 #include "gromacs/options/options.h"
55 #include "gromacs/pbcutil/pbc.h"
56 #include "gromacs/random/random.h"
57 #include "gromacs/selection/nbsearch.h"
58 #include "gromacs/topology/atomprop.h"
59 #include "gromacs/topology/atoms.h"
60 #include "gromacs/utility/cstringutil.h"
61 #include "gromacs/utility/fatalerror.h"
62 #include "gromacs/utility/smalloc.h"
64 /* enum for random rotations of inserted solutes */
66 en_rotXYZ, en_rotZ, en_rotNone
69 static void center_molecule(int atomCount, rvec x[])
73 for (int i = 0; i < atomCount; ++i)
75 rvec_inc(center, x[i]);
77 svmul(1.0/atomCount, center, center);
78 for (int i = 0; i < atomCount; ++i)
80 rvec_dec(x[i], center);
84 static void generate_trial_conf(int atomCount, const rvec xin[],
85 const rvec offset, int enum_rot, gmx_rng_t rng,
88 for (int i = 0; i < atomCount; ++i)
90 copy_rvec(xin[i], xout[i]);
92 real alfa = 0.0, beta = 0.0, gamma = 0.0;
96 alfa = 2*M_PI * gmx_rng_uniform_real(rng);
97 beta = 2*M_PI * gmx_rng_uniform_real(rng);
98 gamma = 2*M_PI * gmx_rng_uniform_real(rng);
102 gamma = 2*M_PI * gmx_rng_uniform_real(rng);
105 alfa = beta = gamma = 0.;
108 if (enum_rot == en_rotXYZ || (enum_rot == en_rotZ))
110 rotate_conf(atomCount, xout, NULL, alfa, beta, gamma);
112 for (int i = 0; i < atomCount; ++i)
114 rvec_inc(xout[i], offset);
118 static bool is_insertion_allowed(gmx::AnalysisNeighborhoodSearch *search,
119 const real *exclusionDistances,
120 int atomCount, const rvec *x,
121 const real *exclusionDistances_insrt)
123 gmx::AnalysisNeighborhoodPositions pos(x, atomCount);
124 gmx::AnalysisNeighborhoodPairSearch pairSearch = search->startPairSearch(pos);
125 gmx::AnalysisNeighborhoodPair pair;
126 while (pairSearch.findNextPair(&pair))
128 const real r1 = exclusionDistances[pair.refIndex()];
129 const real r2 = exclusionDistances_insrt[pair.testIndex()];
130 if (pair.distance2() < sqr(r1 + r2))
138 static void merge_atoms_noalloc(t_atoms *atoms, const t_atoms *atoms_add)
143 resnr = atoms->resinfo[atoms->atom[atoms->nr-1].resind].nr;
146 for (int i = 0; i < atoms_add->nr; ++i)
148 if (atoms_add->atom[i].resind != prevResInd)
150 prevResInd = atoms_add->atom[i].resind;
152 atoms->resinfo[atoms->nres] = atoms_add->resinfo[prevResInd];
153 atoms->resinfo[atoms->nres].nr = resnr;
156 atoms->atom[atoms->nr] = atoms_add->atom[i];
157 atoms->atomname[atoms->nr] = atoms_add->atomname[i];
158 atoms->atom[atoms->nr].resind = atoms->nres-1;
163 static void insert_mols(int nmol_insrt, int ntry, int seed,
164 real defaultDistance, real scaleFactor,
165 t_atoms *atoms, rvec **x,
166 const t_atoms *atoms_insrt, const rvec *x_insrt,
167 int ePBC, matrix box,
168 const std::string &posfn, const rvec deltaR, int enum_rot)
173 fprintf(stderr, "Initialising inter-atomic distances...\n");
174 gmx_atomprop_t aps = gmx_atomprop_init();
175 real *exclusionDistances
176 = makeExclusionDistances(atoms, aps, defaultDistance, scaleFactor);
177 real *exclusionDistances_insrt
178 = makeExclusionDistances(atoms_insrt, aps, defaultDistance, scaleFactor);
179 gmx_atomprop_destroy(aps);
181 const real maxInsertRadius
182 = *std::max_element(exclusionDistances_insrt,
183 exclusionDistances_insrt + atoms_insrt->nr);
184 real maxRadius = maxInsertRadius;
187 const real maxExistingRadius
188 = *std::max_element(exclusionDistances,
189 exclusionDistances + atoms->nr);
190 maxRadius = std::max(maxInsertRadius, maxExistingRadius);
193 // TODO: Make all of this exception-safe.
194 gmx::AnalysisNeighborhood nb;
195 nb.setCutoff(maxInsertRadius + maxRadius);
197 gmx_rng_t rng = gmx_rng_init(seed);
198 set_pbc(&pbc, ePBC, box);
200 snew(x_n, atoms_insrt->nr);
202 /* With -ip, take nmol_insrt from file posfn */
203 double **rpos = NULL;
207 nmol_insrt = read_xvg(posfn.c_str(), &rpos, &ncol);
210 gmx_fatal(FARGS, "Expected 3 columns (x/y/z coordinates) in file %s\n",
213 fprintf(stderr, "Read %d positions from file %s\n\n",
214 nmol_insrt, posfn.c_str());
218 const int finalAtomCount = atoms->nr + nmol_insrt * atoms_insrt->nr;
219 const int finalResidueCount = atoms->nres + nmol_insrt * atoms_insrt->nres;
220 srenew(atoms->resinfo, finalResidueCount);
221 srenew(atoms->atomname, finalAtomCount);
222 srenew(atoms->atom, finalAtomCount);
223 srenew(*x, finalAtomCount);
224 srenew(exclusionDistances, finalAtomCount);
231 while ((mol < nmol_insrt) && (trial < ntry*nmol_insrt))
236 // Insert at random positions.
237 offset_x[XX] = box[XX][XX] * gmx_rng_uniform_real(rng);
238 offset_x[YY] = box[YY][YY] * gmx_rng_uniform_real(rng);
239 offset_x[ZZ] = box[ZZ][ZZ] * gmx_rng_uniform_real(rng);
243 // Skip a position if ntry trials were not successful.
244 if (trial >= firstTrial + ntry)
246 fprintf(stderr, " skipped position (%.3f, %.3f, %.3f)\n",
247 rpos[XX][mol], rpos[YY][mol], rpos[ZZ][mol]);
251 // Insert at positions taken from option -ip file.
252 offset_x[XX] = rpos[XX][mol] + deltaR[XX]*(2 * gmx_rng_uniform_real(rng)-1);
253 offset_x[YY] = rpos[YY][mol] + deltaR[YY]*(2 * gmx_rng_uniform_real(rng)-1);
254 offset_x[ZZ] = rpos[ZZ][mol] + deltaR[ZZ]*(2 * gmx_rng_uniform_real(rng)-1);
256 fprintf(stderr, "\rTry %d", ++trial);
257 generate_trial_conf(atoms_insrt->nr, x_insrt, offset_x, enum_rot, rng,
259 gmx::AnalysisNeighborhoodPositions pos(*x, atoms->nr);
260 gmx::AnalysisNeighborhoodSearch search = nb.initSearch(&pbc, pos);
261 if (is_insertion_allowed(&search, exclusionDistances, atoms_insrt->nr,
262 x_n, exclusionDistances_insrt))
264 const int firstIndex = atoms->nr;
265 for (int i = 0; i < atoms_insrt->nr; ++i)
267 copy_rvec(x_n[i], (*x)[firstIndex + i]);
268 exclusionDistances[firstIndex + i] = exclusionDistances_insrt[i];
270 merge_atoms_noalloc(atoms, atoms_insrt);
273 fprintf(stderr, " success (now %d atoms)!\n", atoms->nr);
276 gmx_rng_destroy(rng);
277 srenew(atoms->resinfo, atoms->nres);
278 srenew(atoms->atomname, atoms->nr);
279 srenew(atoms->atom, atoms->nr);
280 srenew(*x, atoms->nr);
282 fprintf(stderr, "\n");
283 /* print number of molecules added */
284 fprintf(stderr, "Added %d molecules (out of %d requested)\n",
285 mol - failed, nmol_insrt);
288 sfree(exclusionDistances);
289 sfree(exclusionDistances_insrt);
292 for (int i = 0; i < DIM; ++i)
306 class InsertMolecules : public ICommandLineOptionsModule
310 : bBox_(false), nmolIns_(0), nmolTry_(10), seed_(1997),
311 defaultDistance_(0.105), scaleFactor_(0.57), enumRot_(en_rotXYZ)
317 virtual void init(CommandLineModuleSettings * /*settings*/)
321 virtual void initOptions(Options *options,
322 ICommandLineOptionsModuleSettings *settings);
323 virtual void optionsFinished(Options *options);
328 std::string inputConfFile_;
329 std::string insertConfFile_;
330 std::string positionFile_;
331 std::string outputConfFile_;
337 real defaultDistance_;
343 void InsertMolecules::initOptions(Options *options,
344 ICommandLineOptionsModuleSettings *settings)
346 const char *const desc[] = {
347 "[THISMODULE] inserts [TT]-nmol[tt] copies of the system specified in",
348 "the [TT]-ci[tt] input file. The insertions take place either into",
349 "vacant space in the solute conformation given with [TT]-f[tt], or",
350 "into an empty box given by [TT]-box[tt]. Specifying both [TT]-f[tt]",
351 "and [TT]-box[tt] behaves like [TT]-f[tt], but places a new box",
352 "around the solute before insertions. Any velocities present are",
355 "By default, the insertion positions are random (with initial seed",
356 "specified by [TT]-seed[tt]). The program iterates until [TT]-nmol[tt]",
357 "molecules have been inserted in the box. Molecules are not inserted",
358 "where the distance between any existing atom and any atom of the",
359 "inserted molecule is less than the sum based on the van der Waals",
360 "radii of both atoms. A database ([TT]vdwradii.dat[tt]) of van der",
361 "Waals radii is read by the program, and the resulting radii scaled",
362 "by [TT]-scale[tt]. If radii are not found in the database, those"
363 "atoms are assigned the (pre-scaled) distance [TT]-radius[tt].[PAR]",
365 "A total of [TT]-nmol[tt] * [TT]-try[tt] insertion attempts are made",
366 "before giving up. Increase [TT]-try[tt] if you have several small",
367 "holes to fill. Option [TT]-rot[tt] specifies whether the insertion",
368 "molecules are randomly oriented before insertion attempts.[PAR]",
370 "Alternatively, the molecules can be inserted only at positions defined in",
371 "positions.dat ([TT]-ip[tt]). That file should have 3 columns (x,y,z),",
372 "that give the displacements compared to the input molecule position",
373 "([TT]-ci[tt]). Hence, if that file should contain the absolute",
374 "positions, the molecule must be centered on (0,0,0) before using",
375 "[THISMODULE] (e.g. from [gmx-editconf] [TT]-center[tt]).",
376 "Comments in that file starting with # are ignored. Option [TT]-dr[tt]",
377 "defines the maximally allowed displacements during insertial trials.",
378 "[TT]-try[tt] and [TT]-rot[tt] work as in the default mode (see above)."
381 settings->setHelpText(desc);
383 // TODO: Replace use of legacyType.
384 options->addOption(FileNameOption("f")
385 .legacyType(efSTX).inputFile()
386 .store(&inputConfFile_)
387 .defaultBasename("protein")
388 .description("Existing configuration to insert into"));
389 options->addOption(FileNameOption("ci")
390 .legacyType(efSTX).inputFile().required()
391 .store(&insertConfFile_)
392 .defaultBasename("insert")
393 .description("Configuration to insert"));
394 options->addOption(FileNameOption("ip")
395 .filetype(eftGenericData).inputFile()
396 .store(&positionFile_)
397 .defaultBasename("positions")
398 .description("Predefined insertion trial positions"));
399 options->addOption(FileNameOption("o")
400 .legacyType(efSTO).outputFile().required()
401 .store(&outputConfFile_)
402 .defaultBasename("out")
403 .description("Output configuration after insertion"));
405 options->addOption(RealOption("box").vector()
407 .description("Box size (in nm)"));
408 options->addOption(IntegerOption("nmol")
410 .description("Number of extra molecules to insert"));
411 options->addOption(IntegerOption("try")
413 .description("Try inserting [TT]-nmol[tt] times [TT]-try[tt] times"));
414 options->addOption(IntegerOption("seed")
416 .description("Random generator seed"));
417 options->addOption(RealOption("radius")
418 .store(&defaultDistance_)
419 .description("Default van der Waals distance"));
420 options->addOption(RealOption("scale")
421 .store(&scaleFactor_)
422 .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."));
423 options->addOption(RealOption("dr").vector()
425 .description("Allowed displacement in x/y/z from positions in [TT]-ip[tt] file"));
426 const char *const cRotationEnum[] = {"xyz", "z", "none"};
427 options->addOption(StringOption("rot").enumValue(cRotationEnum)
428 .storeEnumIndex(&enumRot_)
429 .description("Rotate inserted molecules randomly"));
432 void InsertMolecules::optionsFinished(Options *options)
434 bBox_ = options->isSet("box");
437 int InsertMolecules::run()
439 const bool bProt = !inputConfFile_.empty();
442 if (nmolIns_ <= 0 && positionFile_.empty())
444 gmx_fatal(FARGS, "Either -nmol must be larger than 0, "
445 "or positions must be given with -ip");
447 if (!bProt && !bBox_)
449 gmx_fatal(FARGS, "When no solute (-f) is specified, "
450 "a box size (-box) must be specified");
459 init_t_atoms(atoms, 0, FALSE);
462 /* Generate a solute configuration */
463 title = readConformation(inputConfFile_.c_str(), atoms, &x, NULL,
464 &ePBC, box, "solute");
467 fprintf(stderr, "Note: no atoms in %s\n", inputConfFile_.c_str());
476 box[XX][XX] = newBox_[XX];
477 box[YY][YY] = newBox_[YY];
478 box[ZZ][ZZ] = newBox_[ZZ];
482 gmx_fatal(FARGS, "Undefined solute box.\nCreate one with gmx editconf "
483 "or give explicit -box command line option");
486 t_atoms *atoms_insrt;
487 rvec *x_insrt = NULL;
488 snew(atoms_insrt, 1);
489 init_t_atoms(atoms_insrt, 0, FALSE);
494 = readConformation(insertConfFile_.c_str(), atoms_insrt, &x_insrt,
495 NULL, &ePBC_dummy, box_dummy, "molecule");
496 if (atoms_insrt->nr == 0)
498 gmx_fatal(FARGS, "No molecule in %s, please check your input",
499 insertConfFile_.c_str());
509 if (positionFile_.empty())
511 center_molecule(atoms_insrt->nr, x_insrt);
515 /* add nmol_ins molecules of atoms_ins
516 in random orientation at random place */
517 insert_mols(nmolIns_, nmolTry_, seed_, defaultDistance_, scaleFactor_,
518 atoms, &x, atoms_insrt, x_insrt,
519 ePBC, box, positionFile_, deltaR_, enumRot_);
521 /* write new configuration to file confout */
522 fprintf(stderr, "Writing generated configuration to %s\n",
523 outputConfFile_.c_str());
524 write_sto_conf(outputConfFile_.c_str(), title, atoms, x, NULL, ePBC, box);
526 /* print size of generated configuration */
527 fprintf(stderr, "\nOutput configuration contains %d atoms in %d residues\n",
528 atoms->nr, atoms->nres);
533 done_atom(atoms_insrt);
543 const char InsertMoleculesInfo::name[] = "insert-molecules";
544 const char InsertMoleculesInfo::shortDescription[] =
545 "Insert molecules into existing vacancies";
546 ICommandLineOptionsModule *InsertMoleculesInfo::create()
548 return new InsertMolecules();