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/pbc.h"
55 #include "gromacs/selection/nbsearch.h"
56 #include "gromacs/topology/atomprop.h"
57 #include "gromacs/topology/atoms.h"
58 #include "gromacs/topology/atomsbuilder.h"
59 #include "gromacs/topology/topology.h"
60 #include "gromacs/utility/arraysize.h"
61 #include "gromacs/utility/cstringutil.h"
62 #include "gromacs/utility/fatalerror.h"
63 #include "gromacs/utility/futil.h"
64 #include "gromacs/utility/gmxassert.h"
65 #include "gromacs/utility/smalloc.h"
77 static void sort_molecule(t_atoms
**atoms_solvt
, std::vector
<RVec
> *x
,
80 int atnr
, i
, j
, moltp
= 0, nrmoltypes
, resi_o
, resi_n
, resnr
;
82 t_atoms
*atoms
, *newatoms
;
84 fprintf(stderr
, "Sorting configuration\n");
88 /* copy each residue from *atoms to a molecule in *molecule */
91 for (i
= 0; i
< atoms
->nr
; i
++)
93 if ( (i
== 0) || (atoms
->atom
[i
].resind
!= atoms
->atom
[i
-1].resind
) )
95 /* see if this was a molecule type we haven't had yet: */
97 for (j
= 0; (j
< nrmoltypes
) && (moltp
== -1); j
++)
99 /* cppcheck-suppress nullPointer
100 * moltypes is guaranteed to be allocated because otherwise
101 * nrmoltypes is 0. */
102 if (strcmp(*(atoms
->resinfo
[atoms
->atom
[i
].resind
].name
), moltypes
[j
].name
) == 0)
111 srenew(moltypes
, nrmoltypes
);
112 moltypes
[moltp
].name
= *(atoms
->resinfo
[atoms
->atom
[i
].resind
].name
);
114 while ((i
+atnr
< atoms
->nr
) &&
115 (atoms
->atom
[i
].resind
== atoms
->atom
[i
+atnr
].resind
))
119 moltypes
[moltp
].natoms
= atnr
;
120 moltypes
[moltp
].nmol
= 0;
122 moltypes
[moltp
].nmol
++;
126 fprintf(stderr
, "Found %d%s molecule type%s:\n",
127 nrmoltypes
, nrmoltypes
== 1 ? "" : " different", nrmoltypes
== 1 ? "" : "s");
128 for (j
= 0; j
< nrmoltypes
; j
++)
132 moltypes
[j
].res0
= 0;
136 moltypes
[j
].res0
= moltypes
[j
-1].res0
+moltypes
[j
-1].nmol
;
138 fprintf(stderr
, "%7s (%4d atoms): %5d residues\n",
139 moltypes
[j
].name
, moltypes
[j
].natoms
, moltypes
[j
].nmol
);
142 /* if we have only 1 moleculetype, we don't have to sort */
145 /* find out which molecules should go where: */
146 moltypes
[0].i
= moltypes
[0].i0
= 0;
147 for (j
= 1; j
< nrmoltypes
; j
++)
151 moltypes
[j
-1].i0
+moltypes
[j
-1].natoms
*moltypes
[j
-1].nmol
;
154 /* now put them there: */
156 init_t_atoms(newatoms
, atoms
->nr
, FALSE
);
157 newatoms
->nres
= atoms
->nres
;
158 snew(newatoms
->resinfo
, atoms
->nres
);
159 std::vector
<RVec
> newx(x
->size());
160 std::vector
<RVec
> newv(v
->size());
165 for (moltp
= 0; moltp
< nrmoltypes
; moltp
++)
168 while (i
< atoms
->nr
)
170 resi_o
= atoms
->atom
[i
].resind
;
171 if (strcmp(*atoms
->resinfo
[resi_o
].name
, moltypes
[moltp
].name
) == 0)
173 /* Copy the residue info */
174 newatoms
->resinfo
[resi_n
] = atoms
->resinfo
[resi_o
];
175 newatoms
->resinfo
[resi_n
].nr
= resnr
;
176 /* Copy the atom info */
179 newatoms
->atom
[j
] = atoms
->atom
[i
];
180 newatoms
->atomname
[j
] = atoms
->atomname
[i
];
181 newatoms
->atom
[j
].resind
= resi_n
;
182 copy_rvec((*x
)[i
], newx
[j
]);
185 copy_rvec((*v
)[i
], newv
[j
]);
190 while (i
< atoms
->nr
&& atoms
->atom
[i
].resind
== resi_o
);
191 /* Increase the new residue counters */
197 /* Skip this residue */
202 while (i
< atoms
->nr
&& atoms
->atom
[i
].resind
== resi_o
);
207 /* put them back into the original arrays and throw away temporary arrays */
208 sfree(atoms
->atomname
);
209 sfree(atoms
->resinfo
);
212 *atoms_solvt
= newatoms
;
219 static void rm_res_pbc(const t_atoms
*atoms
, std::vector
<RVec
> *x
, matrix box
)
227 for (int n
= 0; n
< atoms
->nr
; n
++)
229 if (!is_hydrogen(*(atoms
->atomname
[n
])))
232 rvec_inc(xcg
, (*x
)[n
]);
234 if ( (n
+1 == atoms
->nr
) ||
235 (atoms
->atom
[n
+1].resind
!= atoms
->atom
[n
].resind
) )
237 /* if nat==0 we have only hydrogens in the solvent,
238 we take last coordinate as cg */
242 copy_rvec((*x
)[n
], xcg
);
244 svmul(1.0/nat
, xcg
, xcg
);
245 for (int d
= 0; d
< DIM
; d
++)
249 for (int i
= start
; i
<= n
; i
++)
251 (*x
)[i
][d
] += box
[d
][d
];
255 while (xcg
[d
] >= box
[d
][d
])
257 for (int i
= start
; i
<= n
; i
++)
259 (*x
)[i
][d
] -= box
[d
][d
];
272 * Generates a solvent configuration of desired size by stacking solvent boxes.
274 * \param[in,out] atoms Solvent atoms.
275 * \param[in,out] x Solvent positions.
276 * \param[in,out] v Solvent velocities (`*v` can be NULL).
277 * \param[in,out] r Solvent exclusion radii.
278 * \param[in] box Initial solvent box.
279 * \param[in] boxTarget Target box size.
281 * The solvent box of desired size is created by stacking the initial box in
282 * the smallest k*l*m array that covers the box, and then removing any residue
283 * where all atoms are outside the target box (with a small margin).
284 * This function does not remove overlap between solvent atoms across the
287 * Note that the input configuration should be in the rectangular unit cell and
288 * have whole residues.
290 static void replicateSolventBox(t_atoms
*atoms
, std::vector
<RVec
> *x
,
291 std::vector
<RVec
> *v
, std::vector
<real
> *r
,
292 const matrix box
, const matrix boxTarget
)
294 // Calculate the box multiplication factors.
297 for (int i
= 0; i
< DIM
; ++i
)
300 while (n_box
[i
] * box
[i
][i
] < boxTarget
[i
][i
])
306 fprintf(stderr
, "Will generate new solvent configuration of %dx%dx%d boxes\n",
307 n_box
[XX
], n_box
[YY
], n_box
[ZZ
]);
309 // Create arrays for storing the generated system (cannot be done in-place
310 // in case the target box is smaller than the original in one dimension,
313 init_t_atoms(&newAtoms
, 0, FALSE
);
314 gmx::AtomsBuilder
builder(&newAtoms
, nullptr);
315 builder
.reserve(atoms
->nr
* nmol
, atoms
->nres
* nmol
);
316 std::vector
<RVec
> newX(atoms
->nr
* nmol
);
317 std::vector
<RVec
> newV(!v
->empty() ? atoms
->nr
* nmol
: 0);
318 std::vector
<real
> newR(atoms
->nr
* nmol
);
320 const real maxRadius
= *std::max_element(r
->begin(), r
->end());
322 for (int i
= 0; i
< DIM
; ++i
)
324 // The code below is only interested about the box diagonal.
325 boxWithMargin
[i
] = boxTarget
[i
][i
] + 3*maxRadius
;
328 for (int ix
= 0; ix
< n_box
[XX
]; ++ix
)
331 delta
[XX
] = ix
*box
[XX
][XX
];
332 for (int iy
= 0; iy
< n_box
[YY
]; ++iy
)
334 delta
[YY
] = iy
*box
[YY
][YY
];
335 for (int iz
= 0; iz
< n_box
[ZZ
]; ++iz
)
337 delta
[ZZ
] = iz
*box
[ZZ
][ZZ
];
338 bool bKeepResidue
= false;
339 for (int i
= 0; i
< atoms
->nr
; ++i
)
341 const int newIndex
= builder
.currentAtomCount();
342 bool bKeepAtom
= true;
343 for (int m
= 0; m
< DIM
; ++m
)
345 const real newCoord
= delta
[m
] + (*x
)[i
][m
];
346 bKeepAtom
= bKeepAtom
&& (newCoord
< boxWithMargin
[m
]);
347 newX
[newIndex
][m
] = newCoord
;
349 bKeepResidue
= bKeepResidue
|| bKeepAtom
;
352 copy_rvec((*v
)[i
], newV
[newIndex
]);
354 newR
[newIndex
] = (*r
)[i
];
355 builder
.addAtom(*atoms
, i
);
356 if (i
== atoms
->nr
- 1
357 || atoms
->atom
[i
+1].resind
!= atoms
->atom
[i
].resind
)
361 builder
.finishResidue(atoms
->resinfo
[atoms
->atom
[i
].resind
]);
365 builder
.discardCurrentResidue();
367 // Reset state for the next residue.
368 bKeepResidue
= false;
375 sfree(atoms
->atomname
);
376 sfree(atoms
->resinfo
);
377 atoms
->nr
= newAtoms
.nr
;
378 atoms
->nres
= newAtoms
.nres
;
379 atoms
->atom
= newAtoms
.atom
;
380 atoms
->atomname
= newAtoms
.atomname
;
381 atoms
->resinfo
= newAtoms
.resinfo
;
383 newX
.resize(atoms
->nr
);
387 newV
.resize(atoms
->nr
);
390 newR
.resize(atoms
->nr
);
393 fprintf(stderr
, "Solvent box contains %d atoms in %d residues\n",
394 atoms
->nr
, atoms
->nres
);
398 * Removes overlap of solvent atoms across the edges.
400 * \param[in,out] atoms Solvent atoms.
401 * \param[in,out] x Solvent positions.
402 * \param[in,out] v Solvent velocities (can be empty).
403 * \param[in,out] r Solvent exclusion radii.
404 * \param[in] pbc PBC information.
406 * Solvent residues that lay on the edges that do not touch the origin are
407 * removed if they overlap with other solvent atoms across the PBC.
408 * This is done in this way as the assumption is that the input solvent
409 * configuration is already equilibrated, and so does not contain any
410 * undesirable overlap. The only overlap that should be removed is caused by
411 * cutting the box in half in replicateSolventBox() and leaving a margin of
412 * solvent outside those box edges; these atoms can then overlap with those on
413 * the opposite box edge in a way that is not part of the pre-equilibrated
416 static void removeSolventBoxOverlap(t_atoms
*atoms
, std::vector
<RVec
> *x
,
417 std::vector
<RVec
> *v
, std::vector
<real
> *r
,
420 gmx::AtomsRemover
remover(*atoms
);
422 // TODO: We could limit the amount of pairs searched significantly,
423 // since we are only interested in pairs where the positions are on
425 const real maxRadius
= *std::max_element(r
->begin(), r
->end());
426 gmx::AnalysisNeighborhood nb
;
427 nb
.setCutoff(2*maxRadius
);
428 gmx::AnalysisNeighborhoodPositions
pos(*x
);
429 gmx::AnalysisNeighborhoodSearch search
= nb
.initSearch(&pbc
, pos
);
430 gmx::AnalysisNeighborhoodPairSearch pairSearch
= search
.startPairSearch(pos
);
431 gmx::AnalysisNeighborhoodPair pair
;
432 while (pairSearch
.findNextPair(&pair
))
434 const int i1
= pair
.refIndex();
435 const int i2
= pair
.testIndex();
436 if (remover
.isMarked(i2
))
438 pairSearch
.skipRemainingPairsForTestPosition();
441 if (remover
.isMarked(i1
) || atoms
->atom
[i1
].resind
== atoms
->atom
[i2
].resind
)
445 if (pair
.distance2() < gmx::square((*r
)[i1
] + (*r
)[i2
]))
448 rvec_sub((*x
)[i2
], (*x
)[i1
], dx
);
449 bool bCandidate1
= false, bCandidate2
= false;
450 // To satisfy Clang static analyzer.
451 GMX_ASSERT(pbc
.ndim_ePBC
<= DIM
, "Too many periodic dimensions");
452 for (int d
= 0; d
< pbc
.ndim_ePBC
; ++d
)
454 // If the distance in some dimension is larger than the
455 // cutoff, then it means that the distance has been computed
456 // over the PBC. Mark the position with a larger coordinate
457 // for potential removal.
458 if (dx
[d
] > maxRadius
)
462 else if (dx
[d
] < -maxRadius
)
467 // Only mark one of the positions for removal if both were
469 if (bCandidate2
&& (!bCandidate1
|| i2
> i1
))
471 remover
.markResidue(*atoms
, i2
, true);
472 pairSearch
.skipRemainingPairsForTestPosition();
474 else if (bCandidate1
)
476 remover
.markResidue(*atoms
, i1
, true);
481 remover
.removeMarkedElements(x
);
484 remover
.removeMarkedElements(v
);
486 remover
.removeMarkedElements(r
);
487 const int originalAtomCount
= atoms
->nr
;
488 remover
.removeMarkedAtoms(atoms
);
489 fprintf(stderr
, "Removed %d solvent atoms due to solvent-solvent overlap\n",
490 originalAtomCount
- atoms
->nr
);
494 * Remove all solvent molecules outside a give radius from the solute.
496 * \param[in,out] atoms Solvent atoms.
497 * \param[in,out] x_solvent Solvent positions.
498 * \param[in,out] v_solvent Solvent velocities.
499 * \param[in,out] r Atomic exclusion radii.
500 * \param[in] pbc PBC information.
501 * \param[in] x_solute Solute positions.
502 * \param[in] rshell The radius outside the solute molecule.
504 static void removeSolventOutsideShell(t_atoms
*atoms
,
505 std::vector
<RVec
> *x_solvent
,
506 std::vector
<RVec
> *v_solvent
,
507 std::vector
<real
> *r
,
509 const std::vector
<RVec
> &x_solute
,
512 gmx::AtomsRemover
remover(*atoms
);
513 gmx::AnalysisNeighborhood nb
;
514 nb
.setCutoff(rshell
);
515 gmx::AnalysisNeighborhoodPositions
posSolute(x_solute
);
516 gmx::AnalysisNeighborhoodSearch search
= nb
.initSearch(&pbc
, posSolute
);
517 gmx::AnalysisNeighborhoodPositions
pos(*x_solvent
);
518 gmx::AnalysisNeighborhoodPairSearch pairSearch
= search
.startPairSearch(pos
);
519 gmx::AnalysisNeighborhoodPair pair
;
523 // Now put back those within the shell without checking for overlap
524 while (pairSearch
.findNextPair(&pair
))
526 remover
.markResidue(*atoms
, pair
.testIndex(), false);
527 pairSearch
.skipRemainingPairsForTestPosition();
529 remover
.removeMarkedElements(x_solvent
);
530 if (!v_solvent
->empty())
532 remover
.removeMarkedElements(v_solvent
);
534 remover
.removeMarkedElements(r
);
535 const int originalAtomCount
= atoms
->nr
;
536 remover
.removeMarkedAtoms(atoms
);
537 fprintf(stderr
, "Removed %d solvent atoms more than %f nm from solute.\n",
538 originalAtomCount
- atoms
->nr
, rshell
);
542 * Removes solvent molecules that overlap with the solute.
544 * \param[in,out] atoms Solvent atoms.
545 * \param[in,out] x Solvent positions.
546 * \param[in,out] v Solvent velocities (can be empty).
547 * \param[in,out] r Solvent exclusion radii.
548 * \param[in] pbc PBC information.
549 * \param[in] x_solute Solute positions.
550 * \param[in] r_solute Solute exclusion radii.
552 static void removeSolventOverlappingWithSolute(t_atoms
*atoms
,
553 std::vector
<RVec
> *x
,
554 std::vector
<RVec
> *v
,
555 std::vector
<real
> *r
,
557 const std::vector
<RVec
> &x_solute
,
558 const std::vector
<real
> &r_solute
)
560 gmx::AtomsRemover
remover(*atoms
);
561 const real maxRadius1
562 = *std::max_element(r
->begin(), r
->end());
563 const real maxRadius2
564 = *std::max_element(r_solute
.begin(), r_solute
.end());
566 // Now check for overlap.
567 gmx::AnalysisNeighborhood nb
;
568 gmx::AnalysisNeighborhoodPair pair
;
569 nb
.setCutoff(maxRadius1
+ maxRadius2
);
570 gmx::AnalysisNeighborhoodPositions
posSolute(x_solute
);
571 gmx::AnalysisNeighborhoodSearch search
= nb
.initSearch(&pbc
, posSolute
);
572 gmx::AnalysisNeighborhoodPositions
pos(*x
);
573 gmx::AnalysisNeighborhoodPairSearch pairSearch
= search
.startPairSearch(pos
);
574 while (pairSearch
.findNextPair(&pair
))
576 if (remover
.isMarked(pair
.testIndex()))
578 pairSearch
.skipRemainingPairsForTestPosition();
581 const real r1
= r_solute
[pair
.refIndex()];
582 const real r2
= (*r
)[pair
.testIndex()];
583 const bool bRemove
= (pair
.distance2() < gmx::square(r1
+ r2
));
584 remover
.markResidue(*atoms
, pair
.testIndex(), bRemove
);
587 remover
.removeMarkedElements(x
);
590 remover
.removeMarkedElements(v
);
592 remover
.removeMarkedElements(r
);
593 const int originalAtomCount
= atoms
->nr
;
594 remover
.removeMarkedAtoms(atoms
);
595 fprintf(stderr
, "Removed %d solvent atoms due to solute-solvent overlap\n",
596 originalAtomCount
- atoms
->nr
);
600 * Removes a given number of solvent residues.
602 * \param[in,out] atoms Solvent atoms.
603 * \param[in,out] x Solvent positions.
604 * \param[in,out] v Solvent velocities (can be empty).
605 * \param[in] numberToRemove Number of residues to remove.
607 * This function is called last in the process of creating the solvent box,
608 * so it does not operate on the exclusion radii, as no code after this needs
611 static void removeExtraSolventMolecules(t_atoms
*atoms
, std::vector
<RVec
> *x
,
612 std::vector
<RVec
> *v
,
615 gmx::AtomsRemover
remover(*atoms
);
616 // TODO: It might be nicer to remove a random set of residues, but
617 // in practice this should give a roughly uniform spatial distribution.
618 const int stride
= atoms
->nr
/ numberToRemove
;
619 for (int i
= 0; i
< numberToRemove
; ++i
)
621 int atomIndex
= (i
+1)*stride
- 1;
622 while (remover
.isMarked(atomIndex
))
625 if (atomIndex
== atoms
->nr
)
630 remover
.markResidue(*atoms
, atomIndex
, true);
632 remover
.removeMarkedElements(x
);
635 remover
.removeMarkedElements(v
);
637 remover
.removeMarkedAtoms(atoms
);
640 static void add_solv(const char *fn
, t_topology
*top
,
641 std::vector
<RVec
> *x
, std::vector
<RVec
> *v
,
642 int ePBC
, matrix box
, gmx_atomprop_t aps
,
643 real defaultDistance
, real scaleFactor
,
644 real rshell
, int max_sol
)
646 t_topology
*top_solvt
;
647 std::vector
<RVec
> x_solvt
;
648 std::vector
<RVec
> v_solvt
;
652 char *filename
= gmxlibfn(fn
);
654 readConformation(filename
, top_solvt
, &x_solvt
, !v
->empty() ? &v_solvt
: nullptr,
655 &ePBC_solvt
, box_solvt
, "solvent");
656 t_atoms
*atoms_solvt
= &top_solvt
->atoms
;
657 if (0 == atoms_solvt
->nr
)
659 gmx_fatal(FARGS
, "No solvent in %s, please check your input\n", filename
);
662 fprintf(stderr
, "\n");
664 /* apply pbc for solvent configuration for whole molecules */
665 rm_res_pbc(atoms_solvt
, &x_solvt
, box_solvt
);
667 /* initialise distance arrays for solvent configuration */
668 fprintf(stderr
, "Initialising inter-atomic distances...\n");
669 const std::vector
<real
> exclusionDistances(
670 makeExclusionDistances(&top
->atoms
, aps
, defaultDistance
, scaleFactor
));
671 std::vector
<real
> exclusionDistances_solvt(
672 makeExclusionDistances(atoms_solvt
, aps
, defaultDistance
, scaleFactor
));
674 /* generate a new solvent configuration */
675 fprintf(stderr
, "Generating solvent configuration\n");
677 set_pbc(&pbc
, ePBC
, box
);
678 replicateSolventBox(atoms_solvt
, &x_solvt
, &v_solvt
, &exclusionDistances_solvt
,
680 if (ePBC
!= epbcNONE
)
682 removeSolventBoxOverlap(atoms_solvt
, &x_solvt
, &v_solvt
, &exclusionDistances_solvt
, pbc
);
684 if (top
->atoms
.nr
> 0)
688 removeSolventOutsideShell(atoms_solvt
, &x_solvt
, &v_solvt
,
689 &exclusionDistances_solvt
, pbc
, *x
, rshell
);
691 removeSolventOverlappingWithSolute(atoms_solvt
, &x_solvt
, &v_solvt
,
692 &exclusionDistances_solvt
, pbc
, *x
,
696 if (max_sol
> 0 && atoms_solvt
->nres
> max_sol
)
698 const int numberToRemove
= atoms_solvt
->nres
- max_sol
;
699 removeExtraSolventMolecules(atoms_solvt
, &x_solvt
, &v_solvt
, numberToRemove
);
702 /* Sort the solvent mixture, not the protein... */
703 sort_molecule(&atoms_solvt
, &x_solvt
, &v_solvt
);
705 // Merge the two configurations.
706 x
->insert(x
->end(), x_solvt
.begin(), x_solvt
.end());
709 v
->insert(v
->end(), v_solvt
.begin(), v_solvt
.end());
712 gmx::AtomsBuilder
builder(&top
->atoms
, &top
->symtab
);
713 builder
.mergeAtoms(*atoms_solvt
);
715 fprintf(stderr
, "Generated solvent containing %d atoms in %d residues\n",
716 atoms_solvt
->nr
, atoms_solvt
->nres
);
722 static void update_top(t_atoms
*atoms
, matrix box
, int NFILE
, t_filenm fnm
[],
726 char buf
[STRLEN
], buf2
[STRLEN
], *temp
;
727 const char *topinout
;
729 gmx_bool bSystem
, bMolecules
, bSkip
;
734 for (i
= 0; (i
< atoms
->nres
); i
++)
736 /* calculate number of SOLvent molecules */
737 if ( (strcmp(*atoms
->resinfo
[i
].name
, "SOL") == 0) ||
738 (strcmp(*atoms
->resinfo
[i
].name
, "WAT") == 0) ||
739 (strcmp(*atoms
->resinfo
[i
].name
, "HOH") == 0) )
745 for (i
= 0; (i
< atoms
->nr
); i
++)
747 gmx_atomprop_query(aps
, epropMass
,
748 *atoms
->resinfo
[atoms
->atom
[i
].resind
].name
,
749 *atoms
->atomname
[i
], &mm
);
755 fprintf(stderr
, "Volume : %10g (nm^3)\n", vol
);
756 fprintf(stderr
, "Density : %10g (g/l)\n",
757 (mtot
*1e24
)/(AVOGADRO
*vol
));
758 fprintf(stderr
, "Number of SOL molecules: %5d \n\n", nsol
);
760 /* open topology file and append sol molecules */
761 topinout
= ftp2fn(efTOP
, NFILE
, fnm
);
762 if (ftp2bSet(efTOP
, NFILE
, fnm
) )
764 char temporary_filename
[STRLEN
];
765 strncpy(temporary_filename
, "temp.topXXXXXX", STRLEN
);
767 fprintf(stderr
, "Processing topology\n");
768 fpin
= gmx_ffopen(topinout
, "r");
769 fpout
= gmx_fopen_temporary(temporary_filename
);
771 bSystem
= bMolecules
= FALSE
;
772 while (fgets(buf
, STRLEN
, fpin
))
777 if ((temp
= strchr(buf2
, '\n')) != nullptr)
785 if ((temp
= strchr(buf2
, '\n')) != nullptr)
790 if (buf2
[strlen(buf2
)-1] == ']')
792 buf2
[strlen(buf2
)-1] = '\0';
795 bSystem
= (gmx_strcasecmp(buf2
, "system") == 0);
796 bMolecules
= (gmx_strcasecmp(buf2
, "molecules") == 0);
799 else if (bSystem
&& nsol
&& (buf
[0] != ';') )
801 /* if sol present, append "in water" to system name */
803 if (buf2
[0] && (!strstr(buf2
, " water")) )
805 sprintf(buf
, "%s in water\n", buf2
);
811 /* check if this is a line with solvent molecules */
812 sscanf(buf
, "%4095s", buf2
);
813 if (strcmp(buf2
, "SOL") == 0)
815 sscanf(buf
, "%*4095s %20d", &i
);
826 if ((temp
= strchr(buf
, '\n')) != nullptr)
830 fprintf(stdout
, "Removing line #%d '%s' from topology file (%s)\n",
831 line
, buf
, topinout
);
835 fprintf(fpout
, "%s", buf
);
841 fprintf(stdout
, "Adding line for %d solvent molecules to "
842 "topology file (%s)\n", nsol
, topinout
);
843 fprintf(fpout
, "%-15s %5d\n", "SOL", nsol
);
846 make_backup(topinout
);
847 gmx_file_rename(temporary_filename
, topinout
);
851 int gmx_solvate(int argc
, char *argv
[])
853 const char *desc
[] = {
854 "[THISMODULE] can do one of 2 things:[PAR]",
856 "1) Generate a box of solvent. Specify [TT]-cs[tt] and [TT]-box[tt].",
857 "Or specify [TT]-cs[tt] and [TT]-cp[tt] with a structure file with",
858 "a box, but without atoms.[PAR]",
860 "2) Solvate a solute configuration, e.g. a protein, in a bath of solvent ",
861 "molecules. Specify [TT]-cp[tt] (solute) and [TT]-cs[tt] (solvent). ",
862 "The box specified in the solute coordinate file ([TT]-cp[tt]) is used,",
863 "unless [TT]-box[tt] is set.",
864 "If you want the solute to be centered in the box,",
865 "the program [gmx-editconf] has sophisticated options",
866 "to change the box dimensions and center the solute.",
867 "Solvent molecules are removed from the box where the ",
868 "distance between any atom of the solute molecule(s) and any atom of ",
869 "the solvent molecule is less than the sum of the scaled van der Waals",
870 "radii of both atoms. A database ([TT]vdwradii.dat[tt]) of van der",
871 "Waals radii is read by the program, and the resulting radii scaled",
872 "by [TT]-scale[tt]. If radii are not found in the database, those",
873 "atoms are assigned the (pre-scaled) distance [TT]-radius[tt].",
874 "Note that the usefulness of those radii depends on the atom names,",
875 "and thus varies widely with force field.",
877 "The default solvent is Simple Point Charge water (SPC), with coordinates ",
878 "from [TT]$GMXLIB/spc216.gro[tt]. These coordinates can also be used",
879 "for other 3-site water models, since a short equibilibration will remove",
880 "the small differences between the models.",
881 "Other solvents are also supported, as well as mixed solvents. The",
882 "only restriction to solvent types is that a solvent molecule consists",
883 "of exactly one residue. The residue information in the coordinate",
884 "files is used, and should therefore be more or less consistent.",
885 "In practice this means that two subsequent solvent molecules in the ",
886 "solvent coordinate file should have different residue number.",
887 "The box of solute is built by stacking the coordinates read from",
888 "the coordinate file. This means that these coordinates should be ",
889 "equlibrated in periodic boundary conditions to ensure a good",
890 "alignment of molecules on the stacking interfaces.",
891 "The [TT]-maxsol[tt] option simply adds only the first [TT]-maxsol[tt]",
892 "solvent molecules and leaves out the rest that would have fitted",
893 "into the box. This can create a void that can cause problems later.",
894 "Choose your volume wisely.[PAR]",
896 "Setting [TT]-shell[tt] larger than zero will place a layer of water of",
897 "the specified thickness (nm) around the solute. Hint: it is a good",
898 "idea to put the protein in the center of a box first (using [gmx-editconf]).",
901 "Finally, [THISMODULE] will optionally remove lines from your topology file in ",
902 "which a number of solvent molecules is already added, and adds a ",
903 "line with the total number of solvent molecules in your coordinate file."
906 const char *bugs
[] = {
907 "Molecules must be whole in the initial configurations.",
911 gmx_bool bProt
, bBox
;
912 const char *conf_prot
, *confout
;
915 /* solute configuration data */
921 { efSTX
, "-cp", "protein", ffOPTRD
},
922 { efSTX
, "-cs", "spc216", ffLIBRD
},
923 { efSTO
, nullptr, nullptr, ffWRITE
},
924 { efTOP
, nullptr, nullptr, ffOPTRW
},
926 #define NFILE asize(fnm)
928 static real defaultDistance
= 0.105, r_shell
= 0, scaleFactor
= 0.57;
929 static rvec new_box
= {0.0, 0.0, 0.0};
930 static gmx_bool bReadV
= FALSE
;
931 static int max_sol
= 0;
932 gmx_output_env_t
*oenv
;
934 { "-box", FALSE
, etRVEC
, {new_box
},
935 "Box size (in nm)" },
936 { "-radius", FALSE
, etREAL
, {&defaultDistance
},
937 "Default van der Waals distance"},
938 { "-scale", FALSE
, etREAL
, {&scaleFactor
},
939 "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." },
940 { "-shell", FALSE
, etREAL
, {&r_shell
},
941 "Thickness of optional water layer around solute" },
942 { "-maxsol", FALSE
, etINT
, {&max_sol
},
943 "Maximum number of solvent molecules to add if they fit in the box. If zero (default) this is ignored" },
944 { "-vel", FALSE
, etBOOL
, {&bReadV
},
945 "Keep velocities from input solute and solvent" },
948 if (!parse_common_args(&argc
, argv
, 0, NFILE
, fnm
, asize(pa
), pa
,
949 asize(desc
), desc
, asize(bugs
), bugs
, &oenv
))
954 const char *solventFileName
= opt2fn("-cs", NFILE
, fnm
);
955 bProt
= opt2bSet("-cp", NFILE
, fnm
);
956 bBox
= opt2parg_bSet("-box", asize(pa
), pa
);
961 gmx_fatal(FARGS
, "When no solute (-cp) is specified, "
962 "a box size (-box) must be specified");
965 aps
= gmx_atomprop_init();
972 /* Generate a solute configuration */
973 conf_prot
= opt2fn("-cp", NFILE
, fnm
);
974 readConformation(conf_prot
, top
, &x
,
975 bReadV
? &v
: nullptr, &ePBC
, box
, "solute");
976 if (bReadV
&& v
.empty())
978 fprintf(stderr
, "Note: no velocities found\n");
980 if (top
->atoms
.nr
== 0)
982 fprintf(stderr
, "Note: no atoms in %s\n", conf_prot
);
990 box
[XX
][XX
] = new_box
[XX
];
991 box
[YY
][YY
] = new_box
[YY
];
992 box
[ZZ
][ZZ
] = new_box
[ZZ
];
996 gmx_fatal(FARGS
, "Undefined solute box.\nCreate one with gmx editconf "
997 "or give explicit -box command line option");
1000 add_solv(solventFileName
, top
, &x
, &v
, ePBC
, box
,
1001 aps
, defaultDistance
, scaleFactor
, r_shell
, max_sol
);
1003 /* write new configuration 1 to file confout */
1004 confout
= ftp2fn(efSTO
, NFILE
, fnm
);
1005 fprintf(stderr
, "Writing generated configuration to %s\n", confout
);
1006 const char *outputTitle
= (bProt
? *top
->name
: "Generated by gmx solvate");
1007 write_sto_conf(confout
, outputTitle
, &top
->atoms
, as_rvec_array(x
.data()),
1008 !v
.empty() ? as_rvec_array(v
.data()) : nullptr, ePBC
, box
);
1010 /* print size of generated configuration */
1011 fprintf(stderr
, "\nOutput configuration contains %d atoms in %d residues\n",
1012 top
->atoms
.nr
, top
->atoms
.nres
);
1013 update_top(&top
->atoms
, box
, NFILE
, fnm
, aps
);
1015 gmx_atomprop_destroy(aps
);
1018 output_env_done(oenv
);
1019 done_filenms(NFILE
, fnm
);