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.
46 #include "gromacs/commandline/pargs.h"
47 #include "gromacs/fileio/confio.h"
48 #include "gromacs/fileio/pdbio.h"
49 #include "gromacs/gmxlib/conformation-utilities.h"
50 #include "gromacs/gmxpreprocess/read-conformation.h"
51 #include "gromacs/math/functions.h"
52 #include "gromacs/math/units.h"
53 #include "gromacs/math/vec.h"
54 #include "gromacs/pbcutil/boxutilities.h"
55 #include "gromacs/pbcutil/pbc.h"
56 #include "gromacs/selection/nbsearch.h"
57 #include "gromacs/topology/atomprop.h"
58 #include "gromacs/topology/atoms.h"
59 #include "gromacs/topology/atomsbuilder.h"
60 #include "gromacs/topology/topology.h"
61 #include "gromacs/utility/arraysize.h"
62 #include "gromacs/utility/cstringutil.h"
63 #include "gromacs/utility/fatalerror.h"
64 #include "gromacs/utility/futil.h"
65 #include "gromacs/utility/gmxassert.h"
66 #include "gromacs/utility/smalloc.h"
78 static void sort_molecule(t_atoms
**atoms_solvt
, std::vector
<RVec
> *x
,
81 int atnr
, i
, j
, moltp
= 0, nrmoltypes
, resi_o
, resi_n
, resnr
;
83 t_atoms
*atoms
, *newatoms
;
85 fprintf(stderr
, "Sorting configuration\n");
89 /* copy each residue from *atoms to a molecule in *molecule */
92 for (i
= 0; i
< atoms
->nr
; i
++)
94 if ( (i
== 0) || (atoms
->atom
[i
].resind
!= atoms
->atom
[i
-1].resind
) )
96 /* see if this was a molecule type we haven't had yet: */
98 for (j
= 0; (j
< nrmoltypes
) && (moltp
== -1); j
++)
100 /* cppcheck-suppress nullPointer
101 * moltypes is guaranteed to be allocated because otherwise
102 * nrmoltypes is 0. */
103 if (strcmp(*(atoms
->resinfo
[atoms
->atom
[i
].resind
].name
), moltypes
[j
].name
) == 0)
112 srenew(moltypes
, nrmoltypes
);
113 moltypes
[moltp
].name
= *(atoms
->resinfo
[atoms
->atom
[i
].resind
].name
);
115 while ((i
+atnr
< atoms
->nr
) &&
116 (atoms
->atom
[i
].resind
== atoms
->atom
[i
+atnr
].resind
))
120 moltypes
[moltp
].natoms
= atnr
;
121 moltypes
[moltp
].nmol
= 0;
123 moltypes
[moltp
].nmol
++;
127 fprintf(stderr
, "Found %d%s molecule type%s:\n",
128 nrmoltypes
, nrmoltypes
== 1 ? "" : " different", nrmoltypes
== 1 ? "" : "s");
129 for (j
= 0; j
< nrmoltypes
; j
++)
133 moltypes
[j
].res0
= 0;
137 moltypes
[j
].res0
= moltypes
[j
-1].res0
+moltypes
[j
-1].nmol
;
139 fprintf(stderr
, "%7s (%4d atoms): %5d residues\n",
140 moltypes
[j
].name
, moltypes
[j
].natoms
, moltypes
[j
].nmol
);
143 /* if we have only 1 moleculetype, we don't have to sort */
146 /* find out which molecules should go where: */
147 moltypes
[0].i
= moltypes
[0].i0
= 0;
148 for (j
= 1; j
< nrmoltypes
; j
++)
152 moltypes
[j
-1].i0
+moltypes
[j
-1].natoms
*moltypes
[j
-1].nmol
;
155 /* now put them there: */
157 init_t_atoms(newatoms
, atoms
->nr
, FALSE
);
158 newatoms
->nres
= atoms
->nres
;
159 snew(newatoms
->resinfo
, atoms
->nres
);
160 std::vector
<RVec
> newx(x
->size());
161 std::vector
<RVec
> newv(v
->size());
166 for (moltp
= 0; moltp
< nrmoltypes
; moltp
++)
169 while (i
< atoms
->nr
)
171 resi_o
= atoms
->atom
[i
].resind
;
172 if (strcmp(*atoms
->resinfo
[resi_o
].name
, moltypes
[moltp
].name
) == 0)
174 /* Copy the residue info */
175 newatoms
->resinfo
[resi_n
] = atoms
->resinfo
[resi_o
];
176 newatoms
->resinfo
[resi_n
].nr
= resnr
;
177 /* Copy the atom info */
180 newatoms
->atom
[j
] = atoms
->atom
[i
];
181 newatoms
->atomname
[j
] = atoms
->atomname
[i
];
182 newatoms
->atom
[j
].resind
= resi_n
;
183 copy_rvec((*x
)[i
], newx
[j
]);
186 copy_rvec((*v
)[i
], newv
[j
]);
191 while (i
< atoms
->nr
&& atoms
->atom
[i
].resind
== resi_o
);
192 /* Increase the new residue counters */
198 /* Skip this residue */
203 while (i
< atoms
->nr
&& atoms
->atom
[i
].resind
== resi_o
);
208 /* put them back into the original arrays and throw away temporary arrays */
209 sfree(atoms
->atomname
);
210 sfree(atoms
->resinfo
);
213 *atoms_solvt
= newatoms
;
220 static void rm_res_pbc(const t_atoms
*atoms
, std::vector
<RVec
> *x
, matrix box
)
228 for (int n
= 0; n
< atoms
->nr
; n
++)
230 if (!is_hydrogen(*(atoms
->atomname
[n
])))
233 rvec_inc(xcg
, (*x
)[n
]);
235 if ( (n
+1 == atoms
->nr
) ||
236 (atoms
->atom
[n
+1].resind
!= atoms
->atom
[n
].resind
) )
238 /* if nat==0 we have only hydrogens in the solvent,
239 we take last coordinate as cg */
243 copy_rvec((*x
)[n
], xcg
);
245 svmul(1.0/nat
, xcg
, xcg
);
246 for (int d
= 0; d
< DIM
; d
++)
250 for (int i
= start
; i
<= n
; i
++)
252 (*x
)[i
][d
] += box
[d
][d
];
256 while (xcg
[d
] >= box
[d
][d
])
258 for (int i
= start
; i
<= n
; i
++)
260 (*x
)[i
][d
] -= box
[d
][d
];
273 * Generates a solvent configuration of desired size by stacking solvent boxes.
275 * \param[in,out] atoms Solvent atoms.
276 * \param[in,out] x Solvent positions.
277 * \param[in,out] v Solvent velocities (`*v` can be NULL).
278 * \param[in,out] r Solvent exclusion radii.
279 * \param[in] box Initial solvent box.
280 * \param[in] boxTarget Target box size.
282 * The solvent box of desired size is created by stacking the initial box in
283 * the smallest k*l*m array that covers the box, and then removing any residue
284 * where all atoms are outside the target box (with a small margin).
285 * This function does not remove overlap between solvent atoms across the
288 * Note that the input configuration should be in the rectangular unit cell and
289 * have whole residues.
291 static void replicateSolventBox(t_atoms
*atoms
, std::vector
<RVec
> *x
,
292 std::vector
<RVec
> *v
, std::vector
<real
> *r
,
293 const matrix box
, const matrix boxTarget
)
295 // Calculate the box multiplication factors.
298 for (int i
= 0; i
< DIM
; ++i
)
301 while (n_box
[i
] * box
[i
][i
] < boxTarget
[i
][i
])
307 fprintf(stderr
, "Will generate new solvent configuration of %dx%dx%d boxes\n",
308 n_box
[XX
], n_box
[YY
], n_box
[ZZ
]);
310 // Create arrays for storing the generated system (cannot be done in-place
311 // in case the target box is smaller than the original in one dimension,
314 init_t_atoms(&newAtoms
, 0, FALSE
);
315 gmx::AtomsBuilder
builder(&newAtoms
, nullptr);
316 builder
.reserve(atoms
->nr
* nmol
, atoms
->nres
* nmol
);
317 std::vector
<RVec
> newX(atoms
->nr
* nmol
);
318 std::vector
<RVec
> newV(!v
->empty() ? atoms
->nr
* nmol
: 0);
319 std::vector
<real
> newR(atoms
->nr
* nmol
);
321 const real maxRadius
= *std::max_element(r
->begin(), r
->end());
323 for (int i
= 0; i
< DIM
; ++i
)
325 // The code below is only interested about the box diagonal.
326 boxWithMargin
[i
] = boxTarget
[i
][i
] + 3*maxRadius
;
329 for (int ix
= 0; ix
< n_box
[XX
]; ++ix
)
332 delta
[XX
] = ix
*box
[XX
][XX
];
333 for (int iy
= 0; iy
< n_box
[YY
]; ++iy
)
335 delta
[YY
] = iy
*box
[YY
][YY
];
336 for (int iz
= 0; iz
< n_box
[ZZ
]; ++iz
)
338 delta
[ZZ
] = iz
*box
[ZZ
][ZZ
];
339 bool bKeepResidue
= false;
340 for (int i
= 0; i
< atoms
->nr
; ++i
)
342 const int newIndex
= builder
.currentAtomCount();
343 bool bKeepAtom
= true;
344 for (int m
= 0; m
< DIM
; ++m
)
346 const real newCoord
= delta
[m
] + (*x
)[i
][m
];
347 bKeepAtom
= bKeepAtom
&& (newCoord
< boxWithMargin
[m
]);
348 newX
[newIndex
][m
] = newCoord
;
350 bKeepResidue
= bKeepResidue
|| bKeepAtom
;
353 copy_rvec((*v
)[i
], newV
[newIndex
]);
355 newR
[newIndex
] = (*r
)[i
];
356 builder
.addAtom(*atoms
, i
);
357 if (i
== atoms
->nr
- 1
358 || atoms
->atom
[i
+1].resind
!= atoms
->atom
[i
].resind
)
362 builder
.finishResidue(atoms
->resinfo
[atoms
->atom
[i
].resind
]);
366 builder
.discardCurrentResidue();
368 // Reset state for the next residue.
369 bKeepResidue
= false;
376 sfree(atoms
->atomname
);
377 sfree(atoms
->resinfo
);
378 atoms
->nr
= newAtoms
.nr
;
379 atoms
->nres
= newAtoms
.nres
;
380 atoms
->atom
= newAtoms
.atom
;
381 atoms
->atomname
= newAtoms
.atomname
;
382 atoms
->resinfo
= newAtoms
.resinfo
;
384 newX
.resize(atoms
->nr
);
388 newV
.resize(atoms
->nr
);
391 newR
.resize(atoms
->nr
);
394 fprintf(stderr
, "Solvent box contains %d atoms in %d residues\n",
395 atoms
->nr
, atoms
->nres
);
399 * Removes overlap of solvent atoms across the edges.
401 * \param[in,out] atoms Solvent atoms.
402 * \param[in,out] x Solvent positions.
403 * \param[in,out] v Solvent velocities (can be empty).
404 * \param[in,out] r Solvent exclusion radii.
405 * \param[in] pbc PBC information.
407 * Solvent residues that lay on the edges that do not touch the origin are
408 * removed if they overlap with other solvent atoms across the PBC.
409 * This is done in this way as the assumption is that the input solvent
410 * configuration is already equilibrated, and so does not contain any
411 * undesirable overlap. The only overlap that should be removed is caused by
412 * cutting the box in half in replicateSolventBox() and leaving a margin of
413 * solvent outside those box edges; these atoms can then overlap with those on
414 * the opposite box edge in a way that is not part of the pre-equilibrated
417 static void removeSolventBoxOverlap(t_atoms
*atoms
, std::vector
<RVec
> *x
,
418 std::vector
<RVec
> *v
, std::vector
<real
> *r
,
421 gmx::AtomsRemover
remover(*atoms
);
423 // TODO: We could limit the amount of pairs searched significantly,
424 // since we are only interested in pairs where the positions are on
426 const real maxRadius
= *std::max_element(r
->begin(), r
->end());
427 gmx::AnalysisNeighborhood nb
;
428 nb
.setCutoff(2*maxRadius
);
429 gmx::AnalysisNeighborhoodPositions
pos(*x
);
430 gmx::AnalysisNeighborhoodSearch search
= nb
.initSearch(&pbc
, pos
);
431 gmx::AnalysisNeighborhoodPairSearch pairSearch
= search
.startPairSearch(pos
);
432 gmx::AnalysisNeighborhoodPair pair
;
433 while (pairSearch
.findNextPair(&pair
))
435 const int i1
= pair
.refIndex();
436 const int i2
= pair
.testIndex();
437 if (remover
.isMarked(i2
))
439 pairSearch
.skipRemainingPairsForTestPosition();
442 if (remover
.isMarked(i1
) || atoms
->atom
[i1
].resind
== atoms
->atom
[i2
].resind
)
446 if (pair
.distance2() < gmx::square((*r
)[i1
] + (*r
)[i2
]))
449 rvec_sub((*x
)[i2
], (*x
)[i1
], dx
);
450 bool bCandidate1
= false, bCandidate2
= false;
451 // To satisfy Clang static analyzer.
452 GMX_ASSERT(pbc
.ndim_ePBC
<= DIM
, "Too many periodic dimensions");
453 for (int d
= 0; d
< pbc
.ndim_ePBC
; ++d
)
455 // If the distance in some dimension is larger than the
456 // cutoff, then it means that the distance has been computed
457 // over the PBC. Mark the position with a larger coordinate
458 // for potential removal.
459 if (dx
[d
] > maxRadius
)
463 else if (dx
[d
] < -maxRadius
)
468 // Only mark one of the positions for removal if both were
470 if (bCandidate2
&& (!bCandidate1
|| i2
> i1
))
472 remover
.markResidue(*atoms
, i2
, true);
473 pairSearch
.skipRemainingPairsForTestPosition();
475 else if (bCandidate1
)
477 remover
.markResidue(*atoms
, i1
, true);
482 remover
.removeMarkedElements(x
);
485 remover
.removeMarkedElements(v
);
487 remover
.removeMarkedElements(r
);
488 const int originalAtomCount
= atoms
->nr
;
489 remover
.removeMarkedAtoms(atoms
);
490 fprintf(stderr
, "Removed %d solvent atoms due to solvent-solvent overlap\n",
491 originalAtomCount
- atoms
->nr
);
495 * Remove all solvent molecules outside a give radius from the solute.
497 * \param[in,out] atoms Solvent atoms.
498 * \param[in,out] x_solvent Solvent positions.
499 * \param[in,out] v_solvent Solvent velocities.
500 * \param[in,out] r Atomic exclusion radii.
501 * \param[in] pbc PBC information.
502 * \param[in] x_solute Solute positions.
503 * \param[in] rshell The radius outside the solute molecule.
505 static void removeSolventOutsideShell(t_atoms
*atoms
,
506 std::vector
<RVec
> *x_solvent
,
507 std::vector
<RVec
> *v_solvent
,
508 std::vector
<real
> *r
,
510 const std::vector
<RVec
> &x_solute
,
513 gmx::AtomsRemover
remover(*atoms
);
514 gmx::AnalysisNeighborhood nb
;
515 nb
.setCutoff(rshell
);
516 gmx::AnalysisNeighborhoodPositions
posSolute(x_solute
);
517 gmx::AnalysisNeighborhoodSearch search
= nb
.initSearch(&pbc
, posSolute
);
518 gmx::AnalysisNeighborhoodPositions
pos(*x_solvent
);
519 gmx::AnalysisNeighborhoodPairSearch pairSearch
= search
.startPairSearch(pos
);
520 gmx::AnalysisNeighborhoodPair pair
;
524 // Now put back those within the shell without checking for overlap
525 while (pairSearch
.findNextPair(&pair
))
527 remover
.markResidue(*atoms
, pair
.testIndex(), false);
528 pairSearch
.skipRemainingPairsForTestPosition();
530 remover
.removeMarkedElements(x_solvent
);
531 if (!v_solvent
->empty())
533 remover
.removeMarkedElements(v_solvent
);
535 remover
.removeMarkedElements(r
);
536 const int originalAtomCount
= atoms
->nr
;
537 remover
.removeMarkedAtoms(atoms
);
538 fprintf(stderr
, "Removed %d solvent atoms more than %f nm from solute.\n",
539 originalAtomCount
- atoms
->nr
, rshell
);
543 * Removes solvent molecules that overlap with the solute.
545 * \param[in,out] atoms Solvent atoms.
546 * \param[in,out] x Solvent positions.
547 * \param[in,out] v Solvent velocities (can be empty).
548 * \param[in,out] r Solvent exclusion radii.
549 * \param[in] pbc PBC information.
550 * \param[in] x_solute Solute positions.
551 * \param[in] r_solute Solute exclusion radii.
553 static void removeSolventOverlappingWithSolute(t_atoms
*atoms
,
554 std::vector
<RVec
> *x
,
555 std::vector
<RVec
> *v
,
556 std::vector
<real
> *r
,
558 const std::vector
<RVec
> &x_solute
,
559 const std::vector
<real
> &r_solute
)
561 gmx::AtomsRemover
remover(*atoms
);
562 const real maxRadius1
563 = *std::max_element(r
->begin(), r
->end());
564 const real maxRadius2
565 = *std::max_element(r_solute
.begin(), r_solute
.end());
567 // Now check for overlap.
568 gmx::AnalysisNeighborhood nb
;
569 gmx::AnalysisNeighborhoodPair pair
;
570 nb
.setCutoff(maxRadius1
+ maxRadius2
);
571 gmx::AnalysisNeighborhoodPositions
posSolute(x_solute
);
572 gmx::AnalysisNeighborhoodSearch search
= nb
.initSearch(&pbc
, posSolute
);
573 gmx::AnalysisNeighborhoodPositions
pos(*x
);
574 gmx::AnalysisNeighborhoodPairSearch pairSearch
= search
.startPairSearch(pos
);
575 while (pairSearch
.findNextPair(&pair
))
577 if (remover
.isMarked(pair
.testIndex()))
579 pairSearch
.skipRemainingPairsForTestPosition();
582 const real r1
= r_solute
[pair
.refIndex()];
583 const real r2
= (*r
)[pair
.testIndex()];
584 const bool bRemove
= (pair
.distance2() < gmx::square(r1
+ r2
));
585 remover
.markResidue(*atoms
, pair
.testIndex(), bRemove
);
588 remover
.removeMarkedElements(x
);
591 remover
.removeMarkedElements(v
);
593 remover
.removeMarkedElements(r
);
594 const int originalAtomCount
= atoms
->nr
;
595 remover
.removeMarkedAtoms(atoms
);
596 fprintf(stderr
, "Removed %d solvent atoms due to solute-solvent overlap\n",
597 originalAtomCount
- atoms
->nr
);
601 * Removes a given number of solvent residues.
603 * \param[in,out] atoms Solvent atoms.
604 * \param[in,out] x Solvent positions.
605 * \param[in,out] v Solvent velocities (can be empty).
606 * \param[in] numberToRemove Number of residues to remove.
608 * This function is called last in the process of creating the solvent box,
609 * so it does not operate on the exclusion radii, as no code after this needs
612 static void removeExtraSolventMolecules(t_atoms
*atoms
, std::vector
<RVec
> *x
,
613 std::vector
<RVec
> *v
,
616 gmx::AtomsRemover
remover(*atoms
);
617 // TODO: It might be nicer to remove a random set of residues, but
618 // in practice this should give a roughly uniform spatial distribution.
619 const int stride
= atoms
->nr
/ numberToRemove
;
620 for (int i
= 0; i
< numberToRemove
; ++i
)
622 int atomIndex
= (i
+1)*stride
- 1;
623 while (remover
.isMarked(atomIndex
))
626 if (atomIndex
== atoms
->nr
)
631 remover
.markResidue(*atoms
, atomIndex
, true);
633 remover
.removeMarkedElements(x
);
636 remover
.removeMarkedElements(v
);
638 remover
.removeMarkedAtoms(atoms
);
641 static void add_solv(const char *fn
, t_topology
*top
,
642 std::vector
<RVec
> *x
, std::vector
<RVec
> *v
,
643 int ePBC
, matrix box
, gmx_atomprop_t aps
,
644 real defaultDistance
, real scaleFactor
,
645 real rshell
, int max_sol
)
647 t_topology
*top_solvt
;
648 std::vector
<RVec
> x_solvt
;
649 std::vector
<RVec
> v_solvt
;
653 char *filename
= gmxlibfn(fn
);
655 readConformation(filename
, top_solvt
, &x_solvt
, !v
->empty() ? &v_solvt
: nullptr,
656 &ePBC_solvt
, box_solvt
, "solvent");
657 t_atoms
*atoms_solvt
= &top_solvt
->atoms
;
658 if (0 == atoms_solvt
->nr
)
660 gmx_fatal(FARGS
, "No solvent in %s, please check your input\n", filename
);
663 fprintf(stderr
, "\n");
665 /* initialise distance arrays for solvent configuration */
666 fprintf(stderr
, "Initialising inter-atomic distances...\n");
667 const std::vector
<real
> exclusionDistances(
668 makeExclusionDistances(&top
->atoms
, aps
, defaultDistance
, scaleFactor
));
669 std::vector
<real
> exclusionDistances_solvt(
670 makeExclusionDistances(atoms_solvt
, aps
, defaultDistance
, scaleFactor
));
672 /* generate a new solvent configuration */
673 fprintf(stderr
, "Generating solvent configuration\n");
675 set_pbc(&pbc
, ePBC
, box
);
676 if (!gmx::boxesAreEqual(box_solvt
, box
))
678 if (TRICLINIC(box_solvt
))
680 gmx_fatal(FARGS
, "Generating from non-rectangular solvent boxes is currently not supported.\n"
681 "You can try to pass the same box for -cp and -cs.");
683 /* apply pbc for solvent configuration for whole molecules */
684 rm_res_pbc(atoms_solvt
, &x_solvt
, box_solvt
);
685 replicateSolventBox(atoms_solvt
, &x_solvt
, &v_solvt
, &exclusionDistances_solvt
,
687 if (ePBC
!= epbcNONE
)
689 removeSolventBoxOverlap(atoms_solvt
, &x_solvt
, &v_solvt
, &exclusionDistances_solvt
, pbc
);
692 if (top
->atoms
.nr
> 0)
696 removeSolventOutsideShell(atoms_solvt
, &x_solvt
, &v_solvt
,
697 &exclusionDistances_solvt
, pbc
, *x
, rshell
);
699 removeSolventOverlappingWithSolute(atoms_solvt
, &x_solvt
, &v_solvt
,
700 &exclusionDistances_solvt
, pbc
, *x
,
704 if (max_sol
> 0 && atoms_solvt
->nres
> max_sol
)
706 const int numberToRemove
= atoms_solvt
->nres
- max_sol
;
707 removeExtraSolventMolecules(atoms_solvt
, &x_solvt
, &v_solvt
, numberToRemove
);
710 /* Sort the solvent mixture, not the protein... */
711 sort_molecule(&atoms_solvt
, &x_solvt
, &v_solvt
);
713 // Merge the two configurations.
714 x
->insert(x
->end(), x_solvt
.begin(), x_solvt
.end());
717 v
->insert(v
->end(), v_solvt
.begin(), v_solvt
.end());
720 gmx::AtomsBuilder
builder(&top
->atoms
, &top
->symtab
);
721 builder
.mergeAtoms(*atoms_solvt
);
723 fprintf(stderr
, "Generated solvent containing %d atoms in %d residues\n",
724 atoms_solvt
->nr
, atoms_solvt
->nres
);
730 static void update_top(t_atoms
*atoms
, matrix box
, int NFILE
, t_filenm fnm
[],
734 char buf
[STRLEN
], buf2
[STRLEN
], *temp
;
735 const char *topinout
;
737 gmx_bool bSystem
, bMolecules
, bSkip
;
742 for (i
= 0; (i
< atoms
->nres
); i
++)
744 /* calculate number of SOLvent molecules */
745 if ( (strcmp(*atoms
->resinfo
[i
].name
, "SOL") == 0) ||
746 (strcmp(*atoms
->resinfo
[i
].name
, "WAT") == 0) ||
747 (strcmp(*atoms
->resinfo
[i
].name
, "HOH") == 0) )
753 for (i
= 0; (i
< atoms
->nr
); i
++)
755 gmx_atomprop_query(aps
, epropMass
,
756 *atoms
->resinfo
[atoms
->atom
[i
].resind
].name
,
757 *atoms
->atomname
[i
], &mm
);
763 fprintf(stderr
, "Volume : %10g (nm^3)\n", vol
);
764 fprintf(stderr
, "Density : %10g (g/l)\n",
765 (mtot
*1e24
)/(AVOGADRO
*vol
));
766 fprintf(stderr
, "Number of SOL 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
);
779 bSystem
= bMolecules
= FALSE
;
780 while (fgets(buf
, STRLEN
, fpin
))
785 if ((temp
= strchr(buf2
, '\n')) != nullptr)
793 if ((temp
= strchr(buf2
, '\n')) != nullptr)
798 if (buf2
[strlen(buf2
)-1] == ']')
800 buf2
[strlen(buf2
)-1] = '\0';
803 bSystem
= (gmx_strcasecmp(buf2
, "system") == 0);
804 bMolecules
= (gmx_strcasecmp(buf2
, "molecules") == 0);
807 else if (bSystem
&& nsol
&& (buf
[0] != ';') )
809 /* if sol present, append "in water" to system name */
811 if (buf2
[0] && (!strstr(buf2
, " water")) )
813 sprintf(buf
, "%s in water\n", buf2
);
819 /* check if this is a line with solvent molecules */
820 sscanf(buf
, "%4095s", buf2
);
821 if (strcmp(buf2
, "SOL") == 0)
823 sscanf(buf
, "%*4095s %20d", &i
);
834 if ((temp
= strchr(buf
, '\n')) != nullptr)
838 fprintf(stdout
, "Removing line #%d '%s' from topology file (%s)\n",
839 line
, buf
, topinout
);
843 fprintf(fpout
, "%s", buf
);
849 fprintf(stdout
, "Adding line for %d solvent molecules to "
850 "topology file (%s)\n", nsol
, topinout
);
851 fprintf(fpout
, "%-15s %5d\n", "SOL", nsol
);
854 make_backup(topinout
);
855 gmx_file_rename(temporary_filename
, topinout
);
859 int gmx_solvate(int argc
, char *argv
[])
861 const char *desc
[] = {
862 "[THISMODULE] can do one of 2 things:[PAR]",
864 "1) Generate a box of solvent. Specify [TT]-cs[tt] and [TT]-box[tt].",
865 "Or specify [TT]-cs[tt] and [TT]-cp[tt] with a structure file with",
866 "a box, but without atoms.[PAR]",
868 "2) Solvate a solute configuration, e.g. a protein, in a bath of solvent ",
869 "molecules. Specify [TT]-cp[tt] (solute) and [TT]-cs[tt] (solvent). ",
870 "The box specified in the solute coordinate file ([TT]-cp[tt]) is used,",
871 "unless [TT]-box[tt] is set.",
872 "If you want the solute to be centered in the box,",
873 "the program [gmx-editconf] has sophisticated options",
874 "to change the box dimensions and center the solute.",
875 "Solvent molecules are removed from the box where the ",
876 "distance between any atom of the solute molecule(s) and any atom of ",
877 "the solvent molecule is less than the sum of the scaled van der Waals",
878 "radii of both atoms. A database ([TT]vdwradii.dat[tt]) of van der",
879 "Waals radii is read by the program, and the resulting radii scaled",
880 "by [TT]-scale[tt]. If radii are not found in the database, those",
881 "atoms are assigned the (pre-scaled) distance [TT]-radius[tt].",
882 "Note that the usefulness of those radii depends on the atom names,",
883 "and thus varies widely with force field.",
885 "The default solvent is Simple Point Charge water (SPC), with coordinates ",
886 "from [TT]$GMXLIB/spc216.gro[tt]. These coordinates can also be used",
887 "for other 3-site water models, since a short equibilibration will remove",
888 "the small differences between the models.",
889 "Other solvents are also supported, as well as mixed solvents. The",
890 "only restriction to solvent types is that a solvent molecule consists",
891 "of exactly one residue. The residue information in the coordinate",
892 "files is used, and should therefore be more or less consistent.",
893 "In practice this means that two subsequent solvent molecules in the ",
894 "solvent coordinate file should have different residue number.",
895 "The box of solute is built by stacking the coordinates read from",
896 "the coordinate file. This means that these coordinates should be ",
897 "equlibrated in periodic boundary conditions to ensure a good",
898 "alignment of molecules on the stacking interfaces.",
899 "The [TT]-maxsol[tt] option simply adds only the first [TT]-maxsol[tt]",
900 "solvent molecules and leaves out the rest that would have fitted",
901 "into the box. This can create a void that can cause problems later.",
902 "Choose your volume wisely.[PAR]",
904 "Setting [TT]-shell[tt] larger than zero will place a layer of water of",
905 "the specified thickness (nm) around the solute. Hint: it is a good",
906 "idea to put the protein in the center of a box first (using [gmx-editconf]).",
909 "Finally, [THISMODULE] will optionally remove lines from your topology file in ",
910 "which a number of solvent molecules is already added, and adds a ",
911 "line with the total number of solvent molecules in your coordinate file."
914 const char *bugs
[] = {
915 "Molecules must be whole in the initial configurations.",
919 gmx_bool bProt
, bBox
;
920 const char *conf_prot
, *confout
;
923 /* solute configuration data */
929 { efSTX
, "-cp", "protein", ffOPTRD
},
930 { efSTX
, "-cs", "spc216", ffLIBRD
},
931 { efSTO
, nullptr, nullptr, ffWRITE
},
932 { efTOP
, nullptr, nullptr, ffOPTRW
},
934 #define NFILE asize(fnm)
936 static real defaultDistance
= 0.105, r_shell
= 0, scaleFactor
= 0.57;
937 static rvec new_box
= {0.0, 0.0, 0.0};
938 static gmx_bool bReadV
= FALSE
;
939 static int max_sol
= 0;
940 gmx_output_env_t
*oenv
;
942 { "-box", FALSE
, etRVEC
, {new_box
},
943 "Box size (in nm)" },
944 { "-radius", FALSE
, etREAL
, {&defaultDistance
},
945 "Default van der Waals distance"},
946 { "-scale", FALSE
, etREAL
, {&scaleFactor
},
947 "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." },
948 { "-shell", FALSE
, etREAL
, {&r_shell
},
949 "Thickness of optional water layer around solute" },
950 { "-maxsol", FALSE
, etINT
, {&max_sol
},
951 "Maximum number of solvent molecules to add if they fit in the box. If zero (default) this is ignored" },
952 { "-vel", FALSE
, etBOOL
, {&bReadV
},
953 "Keep velocities from input solute and solvent" },
956 if (!parse_common_args(&argc
, argv
, 0, NFILE
, fnm
, asize(pa
), pa
,
957 asize(desc
), desc
, asize(bugs
), bugs
, &oenv
))
962 const char *solventFileName
= opt2fn("-cs", NFILE
, fnm
);
963 bProt
= opt2bSet("-cp", NFILE
, fnm
);
964 bBox
= opt2parg_bSet("-box", asize(pa
), pa
);
969 gmx_fatal(FARGS
, "When no solute (-cp) is specified, "
970 "a box size (-box) must be specified");
973 aps
= gmx_atomprop_init();
980 /* Generate a solute configuration */
981 conf_prot
= opt2fn("-cp", NFILE
, fnm
);
982 readConformation(conf_prot
, top
, &x
,
983 bReadV
? &v
: nullptr, &ePBC
, box
, "solute");
984 if (bReadV
&& v
.empty())
986 fprintf(stderr
, "Note: no velocities found\n");
988 if (top
->atoms
.nr
== 0)
990 fprintf(stderr
, "Note: no atoms in %s\n", conf_prot
);
998 box
[XX
][XX
] = new_box
[XX
];
999 box
[YY
][YY
] = new_box
[YY
];
1000 box
[ZZ
][ZZ
] = new_box
[ZZ
];
1004 gmx_fatal(FARGS
, "Undefined solute box.\nCreate one with gmx editconf "
1005 "or give explicit -box command line option");
1008 add_solv(solventFileName
, top
, &x
, &v
, ePBC
, box
,
1009 aps
, defaultDistance
, scaleFactor
, r_shell
, max_sol
);
1011 /* write new configuration 1 to file confout */
1012 confout
= ftp2fn(efSTO
, NFILE
, fnm
);
1013 fprintf(stderr
, "Writing generated configuration to %s\n", confout
);
1014 const char *outputTitle
= (bProt
? *top
->name
: "Generated by gmx solvate");
1015 write_sto_conf(confout
, outputTitle
, &top
->atoms
, as_rvec_array(x
.data()),
1016 !v
.empty() ? as_rvec_array(v
.data()) : nullptr, ePBC
, box
);
1018 /* print size of generated configuration */
1019 fprintf(stderr
, "\nOutput configuration contains %d atoms in %d residues\n",
1020 top
->atoms
.nr
, top
->atoms
.nres
);
1021 update_top(&top
->atoms
, box
, NFILE
, fnm
, aps
);
1023 gmx_atomprop_destroy(aps
);
1026 output_env_done(oenv
);
1027 done_filenms(NFILE
, fnm
);