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,2021, 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.
48 #include "gromacs/commandline/pargs.h"
49 #include "gromacs/fileio/confio.h"
50 #include "gromacs/fileio/pdbio.h"
51 #include "gromacs/gmxlib/conformation_utilities.h"
52 #include "gromacs/gmxpreprocess/makeexclusiondistances.h"
53 #include "gromacs/math/functions.h"
54 #include "gromacs/math/units.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/pbcutil/boxutilities.h"
57 #include "gromacs/pbcutil/pbc.h"
58 #include "gromacs/selection/nbsearch.h"
59 #include "gromacs/topology/atomprop.h"
60 #include "gromacs/topology/atoms.h"
61 #include "gromacs/topology/atomsbuilder.h"
62 #include "gromacs/topology/mtop_util.h"
63 #include "gromacs/topology/topology.h"
64 #include "gromacs/utility/arraysize.h"
65 #include "gromacs/utility/cstringutil.h"
66 #include "gromacs/utility/fatalerror.h"
67 #include "gromacs/utility/futil.h"
68 #include "gromacs/utility/gmxassert.h"
69 #include "gromacs/utility/smalloc.h"
73 /*! \brief Describes a molecule type, and keeps track of the number of these molecules
75 * Used for sorting coordinate file data after solvation
81 //! number of atoms in the molecule
83 //! number of occurences of molecule
87 static void sort_molecule(t_atoms** atoms_solvt, t_atoms** newatoms, std::vector<RVec>* x, std::vector<RVec>* v)
90 fprintf(stderr, "Sorting configuration\n");
91 t_atoms* atoms = *atoms_solvt;
93 /* copy each residue from *atoms to a molecule in *molecule */
94 std::vector<MoleculeType> molTypes;
95 for (int i = 0; i < atoms->nr; i++)
97 if ((i == 0) || (atoms->atom[i].resind != atoms->atom[i - 1].resind))
99 /* see if this was a molecule type we haven't had yet: */
100 auto matchingMolType = std::find_if(
101 molTypes.begin(), molTypes.end(), [atoms, i](const MoleculeType& molecule) {
102 return molecule.name == *atoms->resinfo[atoms->atom[i].resind].name;
104 if (matchingMolType == molTypes.end())
106 int numAtomsInMolType = 0;
107 while ((i + numAtomsInMolType < atoms->nr)
108 && (atoms->atom[i].resind == atoms->atom[i + numAtomsInMolType].resind))
112 molTypes.emplace_back(MoleculeType{
113 *atoms->resinfo[atoms->atom[i].resind].name, numAtomsInMolType, 1 });
117 matchingMolType->numMolecules++;
123 "Found %zu%s molecule type%s:\n",
125 molTypes.size() == 1 ? "" : " different",
126 molTypes.size() == 1 ? "" : "s");
127 for (const auto& molType : molTypes)
129 fprintf(stderr, "%7s (%4d atoms): %5d residues\n", molType.name.c_str(), molType.numAtoms, molType.numMolecules);
132 /* if we have only 1 moleculetype, we don't have to sort */
133 if (molTypes.size() > 1)
135 /* now put them there: */
137 init_t_atoms(*newatoms, atoms->nr, FALSE);
138 (*newatoms)->nres = atoms->nres;
139 srenew((*newatoms)->resinfo, atoms->nres);
140 std::vector<RVec> newx(x->size());
141 std::vector<RVec> newv(v->size());
145 for (const auto& moleculeType : molTypes)
148 while (i < atoms->nr)
150 int residOfCurrAtom = atoms->atom[i].resind;
151 if (moleculeType.name == *atoms->resinfo[residOfCurrAtom].name)
153 /* Copy the residue info */
154 (*newatoms)->resinfo[residIndex] = atoms->resinfo[residOfCurrAtom];
155 // Residue numbering starts from 1, so +1 from the index
156 (*newatoms)->resinfo[residIndex].nr = residIndex + 1;
157 /* Copy the atom info */
160 (*newatoms)->atom[atomIndex] = atoms->atom[i];
161 (*newatoms)->atomname[atomIndex] = atoms->atomname[i];
162 (*newatoms)->atom[atomIndex].resind = residIndex;
163 copy_rvec((*x)[i], newx[atomIndex]);
166 copy_rvec((*v)[i], newv[atomIndex]);
170 } while (i < atoms->nr && atoms->atom[i].resind == residOfCurrAtom);
171 /* Increase the new residue counters */
176 /* Skip this residue */
180 } while (i < atoms->nr && atoms->atom[i].resind == residOfCurrAtom);
185 /* put them back into the original arrays and throw away temporary arrays */
187 *atoms_solvt = (*newatoms);
193 static void rm_res_pbc(const t_atoms* atoms, std::vector<RVec>* x, matrix box)
201 for (int n = 0; n < atoms->nr; n++)
203 if (!is_hydrogen(*(atoms->atomname[n])))
206 rvec_inc(xcg, (*x)[n]);
208 if ((n + 1 == atoms->nr) || (atoms->atom[n + 1].resind != atoms->atom[n].resind))
210 /* if nat==0 we have only hydrogens in the solvent,
211 we take last coordinate as cg */
215 copy_rvec((*x)[n], xcg);
217 svmul(1.0 / nat, xcg, xcg);
218 for (int d = 0; d < DIM; d++)
222 for (int i = start; i <= n; i++)
224 (*x)[i][d] += box[d][d];
228 while (xcg[d] >= box[d][d])
230 for (int i = start; i <= n; i++)
232 (*x)[i][d] -= box[d][d];
245 * Generates a solvent configuration of desired size by stacking solvent boxes.
247 * \param[in,out] atoms Solvent atoms.
248 * \param[in,out] x Solvent positions.
249 * \param[in,out] v Solvent velocities (`*v` can be NULL).
250 * \param[in,out] r Solvent exclusion radii.
251 * \param[in] box Initial solvent box.
252 * \param[in] boxTarget Target box size.
254 * The solvent box of desired size is created by stacking the initial box in
255 * the smallest k*l*m array that covers the box, and then removing any residue
256 * where all atoms are outside the target box (with a small margin).
257 * This function does not remove overlap between solvent atoms across the
260 * Note that the input configuration should be in the rectangular unit cell and
261 * have whole residues.
263 static void replicateSolventBox(t_atoms* atoms,
264 std::vector<RVec>* x,
265 std::vector<RVec>* v,
266 std::vector<real>* r,
268 const matrix boxTarget)
270 // Calculate the box multiplication factors.
273 for (int i = 0; i < DIM; ++i)
276 while (n_box[i] * box[i][i] < boxTarget[i][i])
282 fprintf(stderr, "Will generate new solvent configuration of %dx%dx%d boxes\n", n_box[XX], n_box[YY], n_box[ZZ]);
284 // Create arrays for storing the generated system (cannot be done in-place
285 // in case the target box is smaller than the original in one dimension,
288 init_t_atoms(&newAtoms, 0, FALSE);
289 gmx::AtomsBuilder builder(&newAtoms, nullptr);
290 builder.reserve(atoms->nr * nmol, atoms->nres * nmol);
291 std::vector<RVec> newX(atoms->nr * nmol);
292 std::vector<RVec> newV(!v->empty() ? atoms->nr * nmol : 0);
293 std::vector<real> newR(atoms->nr * nmol);
295 const real maxRadius = *std::max_element(r->begin(), r->end());
297 for (int i = 0; i < DIM; ++i)
299 // The code below is only interested about the box diagonal.
300 boxWithMargin[i] = boxTarget[i][i] + 3 * maxRadius;
303 for (int ix = 0; ix < n_box[XX]; ++ix)
306 delta[XX] = ix * box[XX][XX];
307 for (int iy = 0; iy < n_box[YY]; ++iy)
309 delta[YY] = iy * box[YY][YY];
310 for (int iz = 0; iz < n_box[ZZ]; ++iz)
312 delta[ZZ] = iz * box[ZZ][ZZ];
313 bool bKeepResidue = false;
314 for (int i = 0; i < atoms->nr; ++i)
316 const int newIndex = builder.currentAtomCount();
317 bool bKeepAtom = true;
318 for (int m = 0; m < DIM; ++m)
320 const real newCoord = delta[m] + (*x)[i][m];
321 bKeepAtom = bKeepAtom && (newCoord < boxWithMargin[m]);
322 newX[newIndex][m] = newCoord;
324 bKeepResidue = bKeepResidue || bKeepAtom;
327 copy_rvec((*v)[i], newV[newIndex]);
329 newR[newIndex] = (*r)[i];
330 builder.addAtom(*atoms, i);
331 if (i == atoms->nr - 1 || atoms->atom[i + 1].resind != atoms->atom[i].resind)
335 builder.finishResidue(atoms->resinfo[atoms->atom[i].resind]);
339 builder.discardCurrentResidue();
341 // Reset state for the next residue.
342 bKeepResidue = false;
349 sfree(atoms->atomname);
350 sfree(atoms->resinfo);
351 atoms->nr = newAtoms.nr;
352 atoms->nres = newAtoms.nres;
353 atoms->atom = newAtoms.atom;
354 atoms->atomname = newAtoms.atomname;
355 atoms->resinfo = newAtoms.resinfo;
357 newX.resize(atoms->nr);
361 newV.resize(atoms->nr);
364 newR.resize(atoms->nr);
367 fprintf(stderr, "Solvent box contains %d atoms in %d residues\n", atoms->nr, atoms->nres);
371 * Removes overlap of solvent atoms across the edges.
373 * \param[in,out] atoms Solvent atoms.
374 * \param[in,out] x Solvent positions.
375 * \param[in,out] v Solvent velocities (can be empty).
376 * \param[in,out] r Solvent exclusion radii.
377 * \param[in] pbc PBC information.
379 * Solvent residues that lay on the edges that do not touch the origin are
380 * removed if they overlap with other solvent atoms across the PBC.
381 * This is done in this way as the assumption is that the input solvent
382 * configuration is already equilibrated, and so does not contain any
383 * undesirable overlap. The only overlap that should be removed is caused by
384 * cutting the box in half in replicateSolventBox() and leaving a margin of
385 * solvent outside those box edges; these atoms can then overlap with those on
386 * the opposite box edge in a way that is not part of the pre-equilibrated
389 static void removeSolventBoxOverlap(t_atoms* atoms,
390 std::vector<RVec>* x,
391 std::vector<RVec>* v,
392 std::vector<real>* r,
395 gmx::AtomsRemover remover(*atoms);
397 // TODO: We could limit the amount of pairs searched significantly,
398 // since we are only interested in pairs where the positions are on
400 const real maxRadius = *std::max_element(r->begin(), r->end());
401 gmx::AnalysisNeighborhood nb;
402 nb.setCutoff(2 * maxRadius);
403 gmx::AnalysisNeighborhoodPositions pos(*x);
404 gmx::AnalysisNeighborhoodSearch search = nb.initSearch(&pbc, pos);
405 gmx::AnalysisNeighborhoodPairSearch pairSearch = search.startPairSearch(pos);
406 gmx::AnalysisNeighborhoodPair pair;
407 while (pairSearch.findNextPair(&pair))
409 const int i1 = pair.refIndex();
410 const int i2 = pair.testIndex();
411 if (remover.isMarked(i2))
413 pairSearch.skipRemainingPairsForTestPosition();
416 if (remover.isMarked(i1) || atoms->atom[i1].resind == atoms->atom[i2].resind)
420 if (pair.distance2() < gmx::square((*r)[i1] + (*r)[i2]))
423 rvec_sub((*x)[i2], (*x)[i1], dx);
424 bool bCandidate1 = false, bCandidate2 = false;
425 // To satisfy Clang static analyzer.
426 GMX_ASSERT(pbc.ndim_ePBC <= DIM, "Too many periodic dimensions");
427 for (int d = 0; d < pbc.ndim_ePBC; ++d)
429 // If the distance in some dimension is larger than the
430 // cutoff, then it means that the distance has been computed
431 // over the PBC. Mark the position with a larger coordinate
432 // for potential removal.
433 if (dx[d] > maxRadius)
437 else if (dx[d] < -maxRadius)
442 // Only mark one of the positions for removal if both were
444 if (bCandidate2 && (!bCandidate1 || i2 > i1))
446 remover.markResidue(*atoms, i2, true);
447 pairSearch.skipRemainingPairsForTestPosition();
449 else if (bCandidate1)
451 remover.markResidue(*atoms, i1, true);
456 remover.removeMarkedElements(x);
459 remover.removeMarkedElements(v);
461 remover.removeMarkedElements(r);
462 const int originalAtomCount = atoms->nr;
463 remover.removeMarkedAtoms(atoms);
464 fprintf(stderr, "Removed %d solvent atoms due to solvent-solvent overlap\n", originalAtomCount - atoms->nr);
468 * Remove all solvent molecules outside a give radius from the solute.
470 * \param[in,out] atoms Solvent atoms.
471 * \param[in,out] x_solvent Solvent positions.
472 * \param[in,out] v_solvent Solvent velocities.
473 * \param[in,out] r Atomic exclusion radii.
474 * \param[in] pbc PBC information.
475 * \param[in] x_solute Solute positions.
476 * \param[in] rshell The radius outside the solute molecule.
478 static void removeSolventOutsideShell(t_atoms* atoms,
479 std::vector<RVec>* x_solvent,
480 std::vector<RVec>* v_solvent,
481 std::vector<real>* r,
483 const std::vector<RVec>& x_solute,
486 gmx::AtomsRemover remover(*atoms);
487 gmx::AnalysisNeighborhood nb;
488 nb.setCutoff(rshell);
489 gmx::AnalysisNeighborhoodPositions posSolute(x_solute);
490 gmx::AnalysisNeighborhoodSearch search = nb.initSearch(&pbc, posSolute);
491 gmx::AnalysisNeighborhoodPositions pos(*x_solvent);
492 gmx::AnalysisNeighborhoodPairSearch pairSearch = search.startPairSearch(pos);
493 gmx::AnalysisNeighborhoodPair pair;
497 // Now put back those within the shell without checking for overlap
498 while (pairSearch.findNextPair(&pair))
500 remover.markResidue(*atoms, pair.testIndex(), false);
501 pairSearch.skipRemainingPairsForTestPosition();
503 remover.removeMarkedElements(x_solvent);
504 if (!v_solvent->empty())
506 remover.removeMarkedElements(v_solvent);
508 remover.removeMarkedElements(r);
509 const int originalAtomCount = atoms->nr;
510 remover.removeMarkedAtoms(atoms);
512 "Removed %d solvent atoms more than %f nm from solute.\n",
513 originalAtomCount - atoms->nr,
518 * Removes solvent molecules that overlap with the solute.
520 * \param[in,out] atoms Solvent atoms.
521 * \param[in,out] x Solvent positions.
522 * \param[in,out] v Solvent velocities (can be empty).
523 * \param[in,out] r Solvent exclusion radii.
524 * \param[in] pbc PBC information.
525 * \param[in] x_solute Solute positions.
526 * \param[in] r_solute Solute exclusion radii.
528 static void removeSolventOverlappingWithSolute(t_atoms* atoms,
529 std::vector<RVec>* x,
530 std::vector<RVec>* v,
531 std::vector<real>* r,
533 const std::vector<RVec>& x_solute,
534 const std::vector<real>& r_solute)
536 gmx::AtomsRemover remover(*atoms);
537 const real maxRadius1 = *std::max_element(r->begin(), r->end());
538 const real maxRadius2 = *std::max_element(r_solute.begin(), r_solute.end());
540 // Now check for overlap.
541 gmx::AnalysisNeighborhood nb;
542 gmx::AnalysisNeighborhoodPair pair;
543 nb.setCutoff(maxRadius1 + maxRadius2);
544 gmx::AnalysisNeighborhoodPositions posSolute(x_solute);
545 gmx::AnalysisNeighborhoodSearch search = nb.initSearch(&pbc, posSolute);
546 gmx::AnalysisNeighborhoodPositions pos(*x);
547 gmx::AnalysisNeighborhoodPairSearch pairSearch = search.startPairSearch(pos);
548 while (pairSearch.findNextPair(&pair))
550 if (remover.isMarked(pair.testIndex()))
552 pairSearch.skipRemainingPairsForTestPosition();
555 const real r1 = r_solute[pair.refIndex()];
556 const real r2 = (*r)[pair.testIndex()];
557 const bool bRemove = (pair.distance2() < gmx::square(r1 + r2));
558 remover.markResidue(*atoms, pair.testIndex(), bRemove);
561 remover.removeMarkedElements(x);
564 remover.removeMarkedElements(v);
566 remover.removeMarkedElements(r);
567 const int originalAtomCount = atoms->nr;
568 remover.removeMarkedAtoms(atoms);
569 fprintf(stderr, "Removed %d solvent atoms due to solute-solvent overlap\n", originalAtomCount - atoms->nr);
573 * Removes a given number of solvent residues.
575 * \param[in,out] atoms Solvent atoms.
576 * \param[in,out] x Solvent positions.
577 * \param[in,out] v Solvent velocities (can be empty).
578 * \param[in] numberToRemove Number of residues to remove.
580 * This function is called last in the process of creating the solvent box,
581 * so it does not operate on the exclusion radii, as no code after this needs
584 static void removeExtraSolventMolecules(t_atoms* atoms, std::vector<RVec>* x, std::vector<RVec>* v, int numberToRemove)
586 gmx::AtomsRemover remover(*atoms);
587 std::random_device rd;
588 std::mt19937 randomNumberGenerator(rd());
589 std::uniform_int_distribution<> randomDistribution(0, atoms->nr - 1);
590 while (numberToRemove > 0)
592 int atomIndex = randomDistribution(randomNumberGenerator);
593 if (!remover.isMarked(atomIndex))
595 remover.markResidue(*atoms, atomIndex, true);
599 remover.removeMarkedElements(x);
602 remover.removeMarkedElements(v);
604 remover.removeMarkedAtoms(atoms);
607 static void add_solv(const char* filename,
610 std::vector<RVec>* x,
611 std::vector<RVec>* v,
615 real defaultDistance,
620 gmx_mtop_t topSolvent;
621 std::vector<RVec> xSolvent, vSolvent;
622 matrix boxSolvent = { { 0 } };
623 PbcType pbcTypeSolvent;
625 fprintf(stderr, "Reading solvent configuration\n");
626 bool bTprFileWasRead;
627 rvec *temporaryX = nullptr, *temporaryV = nullptr;
628 readConfAndTopology(gmx::findLibraryFile(filename).c_str(),
635 t_atoms* atomsSolvent;
636 snew(atomsSolvent, 1);
637 *atomsSolvent = gmx_mtop_global_atoms(topSolvent);
638 xSolvent.assign(temporaryX, temporaryX + topSolvent.natoms);
640 vSolvent.assign(temporaryV, temporaryV + topSolvent.natoms);
642 if (gmx::boxIsZero(boxSolvent))
645 "No box information for solvent in %s, please use a properly formatted file\n",
648 if (0 == atomsSolvent->nr)
650 gmx_fatal(FARGS, "No solvent in %s, please check your input\n", filename);
652 fprintf(stderr, "\n");
654 /* initialise distance arrays for solvent configuration */
655 fprintf(stderr, "Initialising inter-atomic distances...\n");
656 const std::vector<real> exclusionDistances(
657 makeExclusionDistances(atoms, aps, defaultDistance, scaleFactor));
658 std::vector<real> exclusionDistances_solvt(
659 makeExclusionDistances(atomsSolvent, aps, defaultDistance, scaleFactor));
661 /* generate a new solvent configuration */
662 fprintf(stderr, "Generating solvent configuration\n");
664 set_pbc(&pbc, pbcType, box);
665 if (!gmx::boxesAreEqual(boxSolvent, box))
667 if (TRICLINIC(boxSolvent))
670 "Generating from non-rectangular solvent boxes is currently not supported.\n"
671 "You can try to pass the same box for -cp and -cs.");
673 /* apply pbc for solvent configuration for whole molecules */
674 rm_res_pbc(atomsSolvent, &xSolvent, boxSolvent);
675 replicateSolventBox(atomsSolvent, &xSolvent, &vSolvent, &exclusionDistances_solvt, boxSolvent, box);
676 if (pbcType != PbcType::No)
678 removeSolventBoxOverlap(atomsSolvent, &xSolvent, &vSolvent, &exclusionDistances_solvt, pbc);
685 removeSolventOutsideShell(
686 atomsSolvent, &xSolvent, &vSolvent, &exclusionDistances_solvt, pbc, *x, rshell);
688 removeSolventOverlappingWithSolute(
689 atomsSolvent, &xSolvent, &vSolvent, &exclusionDistances_solvt, pbc, *x, exclusionDistances);
692 if (max_sol > 0 && atomsSolvent->nres > max_sol)
694 const int numberToRemove = atomsSolvent->nres - max_sol;
695 removeExtraSolventMolecules(atomsSolvent, &xSolvent, &vSolvent, numberToRemove);
698 /* Sort the solvent mixture, not the protein... */
699 t_atoms* newatoms = nullptr;
700 // The sort_molecule function does something creative with the
701 // t_atoms pointers. We need to make sure we neither leak, nor
702 // double-free, so make a shallow pointer that is fine for it to
704 t_atoms* sortedAtomsSolvent = atomsSolvent;
705 sort_molecule(&sortedAtomsSolvent, &newatoms, &xSolvent, &vSolvent);
707 // Merge the two configurations.
708 x->insert(x->end(), xSolvent.begin(), xSolvent.end());
711 v->insert(v->end(), vSolvent.begin(), vSolvent.end());
714 gmx::AtomsBuilder builder(atoms, symtab);
715 builder.mergeAtoms(*sortedAtomsSolvent);
718 "Generated solvent containing %d atoms in %d residues\n",
729 done_atom(atomsSolvent);
734 static void update_top(t_atoms* atoms,
735 int firstSolventResidueIndex,
742 char buf[STRLEN * 2], buf2[STRLEN], *temp;
743 const char* topinout;
750 int nsol = atoms->nres - firstSolventResidueIndex;
753 for (i = 0; (i < atoms->nr); i++)
755 aps->setAtomProperty(epropMass,
756 std::string(*atoms->resinfo[atoms->atom[i].resind].name),
757 std::string(*atoms->atomname[i]),
764 fprintf(stderr, "Volume : %10g (nm^3)\n", vol);
765 fprintf(stderr, "Density : %10g (g/l)\n", (mtot * 1e24) / (gmx::c_avogadro * vol));
766 fprintf(stderr, "Number of solvent molecules: %5d \n\n", nsol);
768 /* open topology file and append sol molecules */
769 topinout = ftp2fn(efTOP, NFILE, fnm);
770 if (ftp2bSet(efTOP, NFILE, fnm))
772 char temporary_filename[STRLEN];
773 strncpy(temporary_filename, "temp.topXXXXXX", STRLEN);
775 fprintf(stderr, "Processing topology\n");
776 fpin = gmx_ffopen(topinout, "r");
777 fpout = gmx_fopen_temporary(temporary_filename);
780 while (fgets(buf, STRLEN, fpin))
784 if ((temp = strchr(buf2, '\n')) != nullptr)
792 if ((temp = strchr(buf2, '\n')) != nullptr)
797 if (buf2[strlen(buf2) - 1] == ']')
799 buf2[strlen(buf2) - 1] = '\0';
802 bSystem = (gmx_strcasecmp(buf2, "system") == 0);
805 else if (bSystem && nsol && (buf[0] != ';'))
807 /* if sol present, append "in water" to system name */
809 if (buf2[0] && (!strstr(buf2, " water")))
811 sprintf(buf, "%s in water\n", buf2);
815 fprintf(fpout, "%s", buf);
819 // Add new solvent molecules to the topology
822 std::string currRes = *atoms->resinfo[firstSolventResidueIndex].name;
825 // Iterate through solvent molecules and increment a count until new resname found
826 for (int i = firstSolventResidueIndex; i < atoms->nres; i++)
828 if ((currRes == *atoms->resinfo[i].name))
834 // Change topology and restart count
836 "Adding line for %d solvent molecules with resname (%s) to "
837 "topology file (%s)\n",
841 fprintf(fpout, "%-15s %5d\n", currRes.c_str(), resCount);
842 currRes = *atoms->resinfo[i].name;
846 // One more print needed for last residue type
848 "Adding line for %d solvent molecules with resname (%s) to "
849 "topology file (%s)\n",
853 fprintf(fpout, "%-15s %5d\n", currRes.c_str(), resCount);
856 make_backup(topinout);
857 gmx_file_rename(temporary_filename, topinout);
861 int gmx_solvate(int argc, char* argv[])
863 const char* desc[] = {
864 "[THISMODULE] can do one of 2 things:[PAR]",
866 "1) Generate a box of solvent. Specify [TT]-cs[tt] and [TT]-box[tt].",
867 "Or specify [TT]-cs[tt] and [TT]-cp[tt] with a structure file with",
868 "a box, but without atoms.[PAR]",
870 "2) Solvate a solute configuration, e.g. a protein, in a bath of solvent ",
871 "molecules. Specify [TT]-cp[tt] (solute) and [TT]-cs[tt] (solvent). ",
872 "The box specified in the solute coordinate file ([TT]-cp[tt]) is used,",
873 "unless [TT]-box[tt] is set.",
874 "If you want the solute to be centered in the box,",
875 "the program [gmx-editconf] has sophisticated options",
876 "to change the box dimensions and center the solute.",
877 "Solvent molecules are removed from the box where the ",
878 "distance between any atom of the solute molecule(s) and any atom of ",
879 "the solvent molecule is less than the sum of the scaled van der Waals",
880 "radii of both atoms. A database ([TT]vdwradii.dat[tt]) of van der",
881 "Waals radii is read by the program, and the resulting radii scaled",
882 "by [TT]-scale[tt]. If radii are not found in the database, those",
883 "atoms are assigned the (pre-scaled) distance [TT]-radius[tt].",
884 "Note that the usefulness of those radii depends on the atom names,",
885 "and thus varies widely with force field.",
887 "The default solvent is Simple Point Charge water (SPC), with coordinates ",
888 "from [TT]$GMXLIB/spc216.gro[tt]. These coordinates can also be used",
889 "for other 3-site water models, since a short equibilibration will remove",
890 "the small differences between the models.",
891 "Other solvents are also supported, as well as mixed solvents. The",
892 "only restriction to solvent types is that a solvent molecule consists",
893 "of exactly one residue. The residue information in the coordinate",
894 "files is used, and should therefore be more or less consistent.",
895 "In practice this means that two subsequent solvent molecules in the ",
896 "solvent coordinate file should have different residue number.",
897 "The box of solute is built by stacking the coordinates read from",
898 "the coordinate file. This means that these coordinates should be ",
899 "equlibrated in periodic boundary conditions to ensure a good",
900 "alignment of molecules on the stacking interfaces.",
901 "The [TT]-maxsol[tt] option simply adds only the first [TT]-maxsol[tt]",
902 "solvent molecules and leaves out the rest that would have fitted",
903 "into the box. This can create a void that can cause problems later.",
904 "Choose your volume wisely.[PAR]",
906 "Setting [TT]-shell[tt] larger than zero will place a layer of water of",
907 "the specified thickness (nm) around the solute. Hint: it is a good",
908 "idea to put the protein in the center of a box first (using [gmx-editconf]).",
911 "Finally, [THISMODULE] will optionally remove lines from your topology file in ",
912 "which a number of solvent molecules is already added, and adds a ",
913 "line with the total number of solvent molecules in your coordinate file."
916 const char* bugs[] = {
917 "Molecules must be whole in the initial configurations.",
921 gmx_bool bProt, bBox;
922 const char *conf_prot, *confout;
925 { efSTX, "-cp", "protein", ffOPTRD },
926 { efSTX, "-cs", "spc216", ffLIBRD },
927 { efSTO, nullptr, nullptr, ffWRITE },
928 { efTOP, nullptr, nullptr, ffOPTRW },
930 #define NFILE asize(fnm)
932 real defaultDistance = 0.105, r_shell = 0, scaleFactor = 0.57;
933 rvec new_box = { 0.0, 0.0, 0.0 };
934 gmx_bool bReadV = FALSE;
936 int firstSolventResidueIndex = 0;
937 gmx_output_env_t* oenv;
939 { "-box", FALSE, etRVEC, { new_box }, "Box size (in nm)" },
940 { "-radius", FALSE, etREAL, { &defaultDistance }, "Default van der Waals distance" },
945 "Scale factor to multiply Van der Waals radii from the database in "
946 "share/gromacs/top/vdwradii.dat. The default value of 0.57 yields density close to 1000 "
947 "g/l for proteins in water." },
948 { "-shell", FALSE, etREAL, { &r_shell }, "Thickness of optional water layer around solute" },
953 "Maximum number of solvent molecules to add if they fit in the box. If zero (default) "
955 { "-vel", FALSE, etBOOL, { &bReadV }, "Keep velocities from input solute and solvent" },
958 if (!parse_common_args(
959 &argc, argv, 0, NFILE, fnm, asize(pa), pa, asize(desc), desc, asize(bugs), bugs, &oenv))
964 const char* solventFileName = opt2fn("-cs", NFILE, fnm);
965 bProt = opt2bSet("-cp", NFILE, fnm);
966 bBox = opt2parg_bSet("-box", asize(pa), pa);
972 "When no solute (-cp) is specified, "
973 "a box size (-box) must be specified");
978 /* solute configuration data */
980 std::vector<RVec> x, v;
981 matrix box = { { 0 } };
982 PbcType pbcType = PbcType::Unset;
987 /* Generate a solute configuration */
988 conf_prot = opt2fn("-cp", NFILE, fnm);
989 fprintf(stderr, "Reading solute configuration%s\n", bReadV ? " and velocities" : "");
990 bool bTprFileWasRead;
991 rvec *temporaryX = nullptr, *temporaryV = nullptr;
993 conf_prot, &bTprFileWasRead, &top, &pbcType, &temporaryX, bReadV ? &temporaryV : nullptr, box);
994 *atoms = gmx_mtop_global_atoms(top);
995 x.assign(temporaryX, temporaryX + top.natoms);
999 v.assign(temporaryV, temporaryV + top.natoms);
1004 fprintf(stderr, "Note: no velocities found\n");
1008 fprintf(stderr, "Note: no atoms in %s\n", conf_prot);
1013 firstSolventResidueIndex = atoms->nres;
1016 PbcType pbcTypeForOutput = pbcType;
1019 pbcTypeForOutput = PbcType::Xyz;
1021 box[XX][XX] = new_box[XX];
1022 box[YY][YY] = new_box[YY];
1023 box[ZZ][ZZ] = new_box[ZZ];
1028 "Undefined solute box.\nCreate one with gmx editconf "
1029 "or give explicit -box command line option");
1032 add_solv(solventFileName, atoms, &top.symtab, &x, &v, pbcTypeForOutput, box, &aps, defaultDistance, scaleFactor, r_shell, max_sol);
1034 /* write new configuration 1 to file confout */
1035 confout = ftp2fn(efSTO, NFILE, fnm);
1036 fprintf(stderr, "Writing generated configuration to %s\n", confout);
1037 const char* outputTitle = (bProt ? *top.name : "Generated by gmx solvate");
1038 write_sto_conf(confout,
1041 as_rvec_array(x.data()),
1042 !v.empty() ? as_rvec_array(v.data()) : nullptr,
1046 /* print size of generated configuration */
1047 fprintf(stderr, "\nOutput configuration contains %d atoms in %d residues\n", atoms->nr, atoms->nres);
1048 update_top(atoms, firstSolventResidueIndex, box, NFILE, fnm, &aps);
1052 output_env_done(oenv);