src/gromacs/gmxpreprocess/solvate.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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, 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.
  11  *
  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.
  16  *
  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.
  21  *
  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.
  26  *
  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.
  34  *
  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.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "solvate.h"
  41
  42 #include <cstring>
  43
  44 #include <algorithm>
  45 #include <random>
  46 #include <vector>
  47
  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"
  70
  71 using gmx::RVec;
  72
  73 /*! \brief Describes a molecule type, and keeps track of the number of these molecules
  74  *
  75  *  Used for sorting coordinate file data after solvation
  76  */
  77 struct MoleculeType
  78 {
  79     //! molecule name
  80     std::string name;
  81     //! number of atoms in the molecule
  82     int numAtoms = 0;
  83     //! number of occurences of molecule
  84     int numMolecules = 0;
  85 };
  86
  87 static void sort_molecule(t_atoms** atoms_solvt, t_atoms** newatoms, std::vector<RVec>* x, std::vector<RVec>* v)
  88 {
  89
  90     fprintf(stderr, "Sorting configuration\n");
  91     t_atoms* atoms = *atoms_solvt;
  92
  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++)
  96     {
  97         if ((i == 0) || (atoms->atom[i].resind != atoms->atom[i - 1].resind))
  98         {
  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;
 103                     });
 104             if (matchingMolType == molTypes.end())
 105             {
 106                 int numAtomsInMolType = 0;
 107                 while ((i + numAtomsInMolType < atoms->nr)
 108                        && (atoms->atom[i].resind == atoms->atom[i + numAtomsInMolType].resind))
 109                 {
 110                     numAtomsInMolType++;
 111                 }
 112                 molTypes.emplace_back(MoleculeType{ *atoms->resinfo[atoms->atom[i].resind].name,
 113                                                     numAtomsInMolType, 1 });
 114             }
 115             else
 116             {
 117                 matchingMolType->numMolecules++;
 118             }
 119         }
 120     }
 121
 122     fprintf(stderr, "Found %zu%s molecule type%s:\n", molTypes.size(),
 123             molTypes.size() == 1 ? "" : " different", molTypes.size() == 1 ? "" : "s");
 124     for (const auto& molType : molTypes)
 125     {
 126         fprintf(stderr, "%7s (%4d atoms): %5d residues\n", molType.name.c_str(), molType.numAtoms,
 127                 molType.numMolecules);
 128     }
 129
 130     /* if we have only 1 moleculetype, we don't have to sort */
 131     if (molTypes.size() > 1)
 132     {
 133         /* now put them there: */
 134         snew(*newatoms, 1);
 135         init_t_atoms(*newatoms, atoms->nr, FALSE);
 136         (*newatoms)->nres = atoms->nres;
 137         srenew((*newatoms)->resinfo, atoms->nres);
 138         std::vector<RVec> newx(x->size());
 139         std::vector<RVec> newv(v->size());
 140
 141         int residIndex = 0;
 142         int atomIndex  = 0;
 143         for (const auto& moleculeType : molTypes)
 144         {
 145             int i = 0;
 146             while (i < atoms->nr)
 147             {
 148                 int residOfCurrAtom = atoms->atom[i].resind;
 149                 if (moleculeType.name == *atoms->resinfo[residOfCurrAtom].name)
 150                 {
 151                     /* Copy the residue info */
 152                     (*newatoms)->resinfo[residIndex] = atoms->resinfo[residOfCurrAtom];
 153                     // Residue numbering starts from 1, so +1 from the index
 154                     (*newatoms)->resinfo[residIndex].nr = residIndex + 1;
 155                     /* Copy the atom info */
 156                     do
 157                     {
 158                         (*newatoms)->atom[atomIndex]        = atoms->atom[i];
 159                         (*newatoms)->atomname[atomIndex]    = atoms->atomname[i];
 160                         (*newatoms)->atom[atomIndex].resind = residIndex;
 161                         copy_rvec((*x)[i], newx[atomIndex]);
 162                         if (!v->empty())
 163                         {
 164                             copy_rvec((*v)[i], newv[atomIndex]);
 165                         }
 166                         i++;
 167                         atomIndex++;
 168                     } while (i < atoms->nr && atoms->atom[i].resind == residOfCurrAtom);
 169                     /* Increase the new residue counters */
 170                     residIndex++;
 171                 }
 172                 else
 173                 {
 174                     /* Skip this residue */
 175                     do
 176                     {
 177                         i++;
 178                     } while (i < atoms->nr && atoms->atom[i].resind == residOfCurrAtom);
 179                 }
 180             }
 181         }
 182
 183         /* put them back into the original arrays and throw away temporary arrays */
 184         done_atom(atoms);
 185         *atoms_solvt = (*newatoms);
 186         std::swap(*x, newx);
 187         std::swap(*v, newv);
 188     }
 189 }
 190
 191 static void rm_res_pbc(const t_atoms* atoms, std::vector<RVec>* x, matrix box)
 192 {
 193     int  start, nat;
 194     rvec xcg;
 195
 196     start = 0;
 197     nat   = 0;
 198     clear_rvec(xcg);
 199     for (int n = 0; n < atoms->nr; n++)
 200     {
 201         if (!is_hydrogen(*(atoms->atomname[n])))
 202         {
 203             nat++;
 204             rvec_inc(xcg, (*x)[n]);
 205         }
 206         if ((n + 1 == atoms->nr) || (atoms->atom[n + 1].resind != atoms->atom[n].resind))
 207         {
 208             /* if nat==0 we have only hydrogens in the solvent,
 209                we take last coordinate as cg */
 210             if (nat == 0)
 211             {
 212                 nat = 1;
 213                 copy_rvec((*x)[n], xcg);
 214             }
 215             svmul(1.0 / nat, xcg, xcg);
 216             for (int d = 0; d < DIM; d++)
 217             {
 218                 while (xcg[d] < 0)
 219                 {
 220                     for (int i = start; i <= n; i++)
 221                     {
 222                         (*x)[i][d] += box[d][d];
 223                     }
 224                     xcg[d] += box[d][d];
 225                 }
 226                 while (xcg[d] >= box[d][d])
 227                 {
 228                     for (int i = start; i <= n; i++)
 229                     {
 230                         (*x)[i][d] -= box[d][d];
 231                     }
 232                     xcg[d] -= box[d][d];
 233                 }
 234             }
 235             start = n + 1;
 236             nat   = 0;
 237             clear_rvec(xcg);
 238         }
 239     }
 240 }
 241
 242 /*! \brief
 243  * Generates a solvent configuration of desired size by stacking solvent boxes.
 244  *
 245  * \param[in,out] atoms      Solvent atoms.
 246  * \param[in,out] x          Solvent positions.
 247  * \param[in,out] v          Solvent velocities (`*v` can be NULL).
 248  * \param[in,out] r          Solvent exclusion radii.
 249  * \param[in]     box        Initial solvent box.
 250  * \param[in]     boxTarget  Target box size.
 251  *
 252  * The solvent box of desired size is created by stacking the initial box in
 253  * the smallest k*l*m array that covers the box, and then removing any residue
 254  * where all atoms are outside the target box (with a small margin).
 255  * This function does not remove overlap between solvent atoms across the
 256  * edges.
 257  *
 258  * Note that the input configuration should be in the rectangular unit cell and
 259  * have whole residues.
 260  */
 261 static void replicateSolventBox(t_atoms*           atoms,
 262                                 std::vector<RVec>* x,
 263                                 std::vector<RVec>* v,
 264                                 std::vector<real>* r,
 265                                 const matrix       box,
 266                                 const matrix       boxTarget)
 267 {
 268     // Calculate the box multiplication factors.
 269     ivec n_box;
 270     int  nmol = 1;
 271     for (int i = 0; i < DIM; ++i)
 272     {
 273         n_box[i] = 1;
 274         while (n_box[i] * box[i][i] < boxTarget[i][i])
 275         {
 276             n_box[i]++;
 277         }
 278         nmol *= n_box[i];
 279     }
 280     fprintf(stderr, "Will generate new solvent configuration of %dx%dx%d boxes\n", n_box[XX],
 281             n_box[YY], n_box[ZZ]);
 282
 283     // Create arrays for storing the generated system (cannot be done in-place
 284     // in case the target box is smaller than the original in one dimension,
 285     // but not in all).
 286     t_atoms newAtoms;
 287     init_t_atoms(&newAtoms, 0, FALSE);
 288     gmx::AtomsBuilder builder(&newAtoms, nullptr);
 289     builder.reserve(atoms->nr * nmol, atoms->nres * nmol);
 290     std::vector<RVec> newX(atoms->nr * nmol);
 291     std::vector<RVec> newV(!v->empty() ? atoms->nr * nmol : 0);
 292     std::vector<real> newR(atoms->nr * nmol);
 293
 294     const real maxRadius = *std::max_element(r->begin(), r->end());
 295     rvec       boxWithMargin;
 296     for (int i = 0; i < DIM; ++i)
 297     {
 298         // The code below is only interested about the box diagonal.
 299         boxWithMargin[i] = boxTarget[i][i] + 3 * maxRadius;
 300     }
 301
 302     for (int ix = 0; ix < n_box[XX]; ++ix)
 303     {
 304         rvec delta;
 305         delta[XX] = ix * box[XX][XX];
 306         for (int iy = 0; iy < n_box[YY]; ++iy)
 307         {
 308             delta[YY] = iy * box[YY][YY];
 309             for (int iz = 0; iz < n_box[ZZ]; ++iz)
 310             {
 311                 delta[ZZ]         = iz * box[ZZ][ZZ];
 312                 bool bKeepResidue = false;
 313                 for (int i = 0; i < atoms->nr; ++i)
 314                 {
 315                     const int newIndex  = builder.currentAtomCount();
 316                     bool      bKeepAtom = true;
 317                     for (int m = 0; m < DIM; ++m)
 318                     {
 319                         const real newCoord = delta[m] + (*x)[i][m];
 320                         bKeepAtom           = bKeepAtom && (newCoord < boxWithMargin[m]);
 321                         newX[newIndex][m]   = newCoord;
 322                     }
 323                     bKeepResidue = bKeepResidue || bKeepAtom;
 324                     if (!v->empty())
 325                     {
 326                         copy_rvec((*v)[i], newV[newIndex]);
 327                     }
 328                     newR[newIndex] = (*r)[i];
 329                     builder.addAtom(*atoms, i);
 330                     if (i == atoms->nr - 1 || atoms->atom[i + 1].resind != atoms->atom[i].resind)
 331                     {
 332                         if (bKeepResidue)
 333                         {
 334                             builder.finishResidue(atoms->resinfo[atoms->atom[i].resind]);
 335                         }
 336                         else
 337                         {
 338                             builder.discardCurrentResidue();
 339                         }
 340                         // Reset state for the next residue.
 341                         bKeepResidue = false;
 342                     }
 343                 }
 344             }
 345         }
 346     }
 347     sfree(atoms->atom);
 348     sfree(atoms->atomname);
 349     sfree(atoms->resinfo);
 350     atoms->nr       = newAtoms.nr;
 351     atoms->nres     = newAtoms.nres;
 352     atoms->atom     = newAtoms.atom;
 353     atoms->atomname = newAtoms.atomname;
 354     atoms->resinfo  = newAtoms.resinfo;
 355
 356     newX.resize(atoms->nr);
 357     std::swap(*x, newX);
 358     if (!v->empty())
 359     {
 360         newV.resize(atoms->nr);
 361         std::swap(*v, newV);
 362     }
 363     newR.resize(atoms->nr);
 364     std::swap(*r, newR);
 365
 366     fprintf(stderr, "Solvent box contains %d atoms in %d residues\n", atoms->nr, atoms->nres);
 367 }
 368
 369 /*! \brief
 370  * Removes overlap of solvent atoms across the edges.
 371  *
 372  * \param[in,out] atoms      Solvent atoms.
 373  * \param[in,out] x          Solvent positions.
 374  * \param[in,out] v          Solvent velocities (can be empty).
 375  * \param[in,out] r          Solvent exclusion radii.
 376  * \param[in]     pbc        PBC information.
 377  *
 378  * Solvent residues that lay on the edges that do not touch the origin are
 379  * removed if they overlap with other solvent atoms across the PBC.
 380  * This is done in this way as the assumption is that the input solvent
 381  * configuration is already equilibrated, and so does not contain any
 382  * undesirable overlap.  The only overlap that should be removed is caused by
 383  * cutting the box in half in replicateSolventBox() and leaving a margin of
 384  * solvent outside those box edges; these atoms can then overlap with those on
 385  * the opposite box edge in a way that is not part of the pre-equilibrated
 386  * configuration.
 387  */
 388 static void removeSolventBoxOverlap(t_atoms*           atoms,
 389                                     std::vector<RVec>* x,
 390                                     std::vector<RVec>* v,
 391                                     std::vector<real>* r,
 392                                     const t_pbc&       pbc)
 393 {
 394     gmx::AtomsRemover remover(*atoms);
 395
 396     // TODO: We could limit the amount of pairs searched significantly,
 397     // since we are only interested in pairs where the positions are on
 398     // opposite edges.
 399     const real                maxRadius = *std::max_element(r->begin(), r->end());
 400     gmx::AnalysisNeighborhood nb;
 401     nb.setCutoff(2 * maxRadius);
 402     gmx::AnalysisNeighborhoodPositions  pos(*x);
 403     gmx::AnalysisNeighborhoodSearch     search     = nb.initSearch(&pbc, pos);
 404     gmx::AnalysisNeighborhoodPairSearch pairSearch = search.startPairSearch(pos);
 405     gmx::AnalysisNeighborhoodPair       pair;
 406     while (pairSearch.findNextPair(&pair))
 407     {
 408         const int i1 = pair.refIndex();
 409         const int i2 = pair.testIndex();
 410         if (remover.isMarked(i2))
 411         {
 412             pairSearch.skipRemainingPairsForTestPosition();
 413             continue;
 414         }
 415         if (remover.isMarked(i1) || atoms->atom[i1].resind == atoms->atom[i2].resind)
 416         {
 417             continue;
 418         }
 419         if (pair.distance2() < gmx::square((*r)[i1] + (*r)[i2]))
 420         {
 421             rvec dx;
 422             rvec_sub((*x)[i2], (*x)[i1], dx);
 423             bool bCandidate1 = false, bCandidate2 = false;
 424             // To satisfy Clang static analyzer.
 425             GMX_ASSERT(pbc.ndim_ePBC <= DIM, "Too many periodic dimensions");
 426             for (int d = 0; d < pbc.ndim_ePBC; ++d)
 427             {
 428                 // If the distance in some dimension is larger than the
 429                 // cutoff, then it means that the distance has been computed
 430                 // over the PBC.  Mark the position with a larger coordinate
 431                 // for potential removal.
 432                 if (dx[d] > maxRadius)
 433                 {
 434                     bCandidate2 = true;
 435                 }
 436                 else if (dx[d] < -maxRadius)
 437                 {
 438                     bCandidate1 = true;
 439                 }
 440             }
 441             // Only mark one of the positions for removal if both were
 442             // candidates.
 443             if (bCandidate2 && (!bCandidate1 || i2 > i1))
 444             {
 445                 remover.markResidue(*atoms, i2, true);
 446                 pairSearch.skipRemainingPairsForTestPosition();
 447             }
 448             else if (bCandidate1)
 449             {
 450                 remover.markResidue(*atoms, i1, true);
 451             }
 452         }
 453     }
 454
 455     remover.removeMarkedElements(x);
 456     if (!v->empty())
 457     {
 458         remover.removeMarkedElements(v);
 459     }
 460     remover.removeMarkedElements(r);
 461     const int originalAtomCount = atoms->nr;
 462     remover.removeMarkedAtoms(atoms);
 463     fprintf(stderr, "Removed %d solvent atoms due to solvent-solvent overlap\n",
 464             originalAtomCount - atoms->nr);
 465 }
 466
 467 /*! \brief
 468  * Remove all solvent molecules outside a give radius from the solute.
 469  *
 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.
 477  */
 478 static void removeSolventOutsideShell(t_atoms*                 atoms,
 479                                       std::vector<RVec>*       x_solvent,
 480                                       std::vector<RVec>*       v_solvent,
 481                                       std::vector<real>*       r,
 482                                       const t_pbc&             pbc,
 483                                       const std::vector<RVec>& x_solute,
 484                                       real                     rshell)
 485 {
 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;
 494
 495     // Remove everything
 496     remover.markAll();
 497     // Now put back those within the shell without checking for overlap
 498     while (pairSearch.findNextPair(&pair))
 499     {
 500         remover.markResidue(*atoms, pair.testIndex(), false);
 501         pairSearch.skipRemainingPairsForTestPosition();
 502     }
 503     remover.removeMarkedElements(x_solvent);
 504     if (!v_solvent->empty())
 505     {
 506         remover.removeMarkedElements(v_solvent);
 507     }
 508     remover.removeMarkedElements(r);
 509     const int originalAtomCount = atoms->nr;
 510     remover.removeMarkedAtoms(atoms);
 511     fprintf(stderr, "Removed %d solvent atoms more than %f nm from solute.\n",
 512             originalAtomCount - atoms->nr, rshell);
 513 }
 514
 515 /*! \brief
 516  * Removes solvent molecules that overlap with the solute.
 517  *
 518  * \param[in,out] atoms    Solvent atoms.
 519  * \param[in,out] x        Solvent positions.
 520  * \param[in,out] v        Solvent velocities (can be empty).
 521  * \param[in,out] r        Solvent exclusion radii.
 522  * \param[in]     pbc      PBC information.
 523  * \param[in]     x_solute Solute positions.
 524  * \param[in]     r_solute Solute exclusion radii.
 525  */
 526 static void removeSolventOverlappingWithSolute(t_atoms*                 atoms,
 527                                                std::vector<RVec>*       x,
 528                                                std::vector<RVec>*       v,
 529                                                std::vector<real>*       r,
 530                                                const t_pbc&             pbc,
 531                                                const std::vector<RVec>& x_solute,
 532                                                const std::vector<real>& r_solute)
 533 {
 534     gmx::AtomsRemover remover(*atoms);
 535     const real        maxRadius1 = *std::max_element(r->begin(), r->end());
 536     const real        maxRadius2 = *std::max_element(r_solute.begin(), r_solute.end());
 537
 538     // Now check for overlap.
 539     gmx::AnalysisNeighborhood     nb;
 540     gmx::AnalysisNeighborhoodPair pair;
 541     nb.setCutoff(maxRadius1 + maxRadius2);
 542     gmx::AnalysisNeighborhoodPositions  posSolute(x_solute);
 543     gmx::AnalysisNeighborhoodSearch     search = nb.initSearch(&pbc, posSolute);
 544     gmx::AnalysisNeighborhoodPositions  pos(*x);
 545     gmx::AnalysisNeighborhoodPairSearch pairSearch = search.startPairSearch(pos);
 546     while (pairSearch.findNextPair(&pair))
 547     {
 548         if (remover.isMarked(pair.testIndex()))
 549         {
 550             pairSearch.skipRemainingPairsForTestPosition();
 551             continue;
 552         }
 553         const real r1      = r_solute[pair.refIndex()];
 554         const real r2      = (*r)[pair.testIndex()];
 555         const bool bRemove = (pair.distance2() < gmx::square(r1 + r2));
 556         remover.markResidue(*atoms, pair.testIndex(), bRemove);
 557     }
 558
 559     remover.removeMarkedElements(x);
 560     if (!v->empty())
 561     {
 562         remover.removeMarkedElements(v);
 563     }
 564     remover.removeMarkedElements(r);
 565     const int originalAtomCount = atoms->nr;
 566     remover.removeMarkedAtoms(atoms);
 567     fprintf(stderr, "Removed %d solvent atoms due to solute-solvent overlap\n",
 568             originalAtomCount - atoms->nr);
 569 }
 570
 571 /*! \brief
 572  * Removes a given number of solvent residues.
 573  *
 574  * \param[in,out] atoms           Solvent atoms.
 575  * \param[in,out] x               Solvent positions.
 576  * \param[in,out] v               Solvent velocities (can be empty).
 577  * \param[in]     numberToRemove  Number of residues to remove.
 578  *
 579  * This function is called last in the process of creating the solvent box,
 580  * so it does not operate on the exclusion radii, as no code after this needs
 581  * them.
 582  */
 583 static void removeExtraSolventMolecules(t_atoms* atoms, std::vector<RVec>* x, std::vector<RVec>* v, int numberToRemove)
 584 {
 585     gmx::AtomsRemover               remover(*atoms);
 586     std::random_device              rd;
 587     std::mt19937                    randomNumberGenerator(rd());
 588     std::uniform_int_distribution<> randomDistribution(0, atoms->nr - 1);
 589     while (numberToRemove > 0)
 590     {
 591         int atomIndex = randomDistribution(randomNumberGenerator);
 592         if (!remover.isMarked(atomIndex))
 593         {
 594             remover.markResidue(*atoms, atomIndex, true);
 595             numberToRemove--;
 596         }
 597     }
 598     remover.removeMarkedElements(x);
 599     if (!v->empty())
 600     {
 601         remover.removeMarkedElements(v);
 602     }
 603     remover.removeMarkedAtoms(atoms);
 604 }
 605
 606 static void add_solv(const char*        filename,
 607                      t_atoms*           atoms,
 608                      t_symtab*          symtab,
 609                      std::vector<RVec>* x,
 610                      std::vector<RVec>* v,
 611                      PbcType            pbcType,
 612                      matrix             box,
 613                      AtomProperties*    aps,
 614                      real               defaultDistance,
 615                      real               scaleFactor,
 616                      real               rshell,
 617                      int                max_sol)
 618 {
 619     gmx_mtop_t        topSolvent;
 620     std::vector<RVec> xSolvent, vSolvent;
 621     matrix            boxSolvent = { { 0 } };
 622     PbcType           pbcTypeSolvent;
 623
 624     fprintf(stderr, "Reading solvent configuration\n");
 625     bool  bTprFileWasRead;
 626     rvec *temporaryX = nullptr, *temporaryV = nullptr;
 627     readConfAndTopology(gmx::findLibraryFile(filename).c_str(), &bTprFileWasRead, &topSolvent,
 628                         &pbcTypeSolvent, &temporaryX, &temporaryV, boxSolvent);
 629     t_atoms* atomsSolvent;
 630     snew(atomsSolvent, 1);
 631     *atomsSolvent = gmx_mtop_global_atoms(&topSolvent);
 632     xSolvent.assign(temporaryX, temporaryX + topSolvent.natoms);
 633     sfree(temporaryX);
 634     vSolvent.assign(temporaryV, temporaryV + topSolvent.natoms);
 635     sfree(temporaryV);
 636     if (gmx::boxIsZero(boxSolvent))
 637     {
 638         gmx_fatal(FARGS,
 639                   "No box information for solvent in %s, please use a properly formatted file\n",
 640                   filename);
 641     }
 642     if (0 == atomsSolvent->nr)
 643     {
 644         gmx_fatal(FARGS, "No solvent in %s, please check your input\n", filename);
 645     }
 646     fprintf(stderr, "\n");
 647
 648     /* initialise distance arrays for solvent configuration */
 649     fprintf(stderr, "Initialising inter-atomic distances...\n");
 650     const std::vector<real> exclusionDistances(
 651             makeExclusionDistances(atoms, aps, defaultDistance, scaleFactor));
 652     std::vector<real> exclusionDistances_solvt(
 653             makeExclusionDistances(atomsSolvent, aps, defaultDistance, scaleFactor));
 654
 655     /* generate a new solvent configuration */
 656     fprintf(stderr, "Generating solvent configuration\n");
 657     t_pbc pbc;
 658     set_pbc(&pbc, pbcType, box);
 659     if (!gmx::boxesAreEqual(boxSolvent, box))
 660     {
 661         if (TRICLINIC(boxSolvent))
 662         {
 663             gmx_fatal(FARGS,
 664                       "Generating from non-rectangular solvent boxes is currently not supported.\n"
 665                       "You can try to pass the same box for -cp and -cs.");
 666         }
 667         /* apply pbc for solvent configuration for whole molecules */
 668         rm_res_pbc(atomsSolvent, &xSolvent, boxSolvent);
 669         replicateSolventBox(atomsSolvent, &xSolvent, &vSolvent, &exclusionDistances_solvt, boxSolvent, box);
 670         if (pbcType != PbcType::No)
 671         {
 672             removeSolventBoxOverlap(atomsSolvent, &xSolvent, &vSolvent, &exclusionDistances_solvt, pbc);
 673         }
 674     }
 675     if (atoms->nr > 0)
 676     {
 677         if (rshell > 0.0)
 678         {
 679             removeSolventOutsideShell(atomsSolvent, &xSolvent, &vSolvent, &exclusionDistances_solvt,
 680                                       pbc, *x, rshell);
 681         }
 682         removeSolventOverlappingWithSolute(atomsSolvent, &xSolvent, &vSolvent,
 683                                            &exclusionDistances_solvt, pbc, *x, exclusionDistances);
 684     }
 685
 686     if (max_sol > 0 && atomsSolvent->nres > max_sol)
 687     {
 688         const int numberToRemove = atomsSolvent->nres - max_sol;
 689         removeExtraSolventMolecules(atomsSolvent, &xSolvent, &vSolvent, numberToRemove);
 690     }
 691
 692     /* Sort the solvent mixture, not the protein... */
 693     t_atoms* newatoms = nullptr;
 694     // The sort_molecule function does something creative with the
 695     // t_atoms pointers. We need to make sure we neither leak, nor
 696     // double-free, so make a shallow pointer that is fine for it to
 697     // change.
 698     t_atoms* sortedAtomsSolvent = atomsSolvent;
 699     sort_molecule(&sortedAtomsSolvent, &newatoms, &xSolvent, &vSolvent);
 700
 701     // Merge the two configurations.
 702     x->insert(x->end(), xSolvent.begin(), xSolvent.end());
 703     if (!v->empty())
 704     {
 705         v->insert(v->end(), vSolvent.begin(), vSolvent.end());
 706     }
 707     {
 708         gmx::AtomsBuilder builder(atoms, symtab);
 709         builder.mergeAtoms(*sortedAtomsSolvent);
 710     }
 711     fprintf(stderr, "Generated solvent containing %d atoms in %d residues\n", atomsSolvent->nr,
 712             atomsSolvent->nres);
 713
 714     if (newatoms)
 715     {
 716         done_atom(newatoms);
 717         sfree(newatoms);
 718     }
 719     if (atomsSolvent)
 720     {
 721         done_atom(atomsSolvent);
 722         sfree(atomsSolvent);
 723     }
 724 }
 725
 726 static void update_top(t_atoms*        atoms,
 727                        int             firstSolventResidueIndex,
 728                        matrix          box,
 729                        int             NFILE,
 730                        t_filenm        fnm[],
 731                        AtomProperties* aps)
 732 {
 733     FILE *      fpin, *fpout;
 734     char        buf[STRLEN * 2], buf2[STRLEN], *temp;
 735     const char* topinout;
 736     int         line;
 737     bool        bSystem;
 738     int         i;
 739     double      mtot;
 740     real        vol, mm;
 741
 742     int nsol = atoms->nres - firstSolventResidueIndex;
 743
 744     mtot = 0;
 745     for (i = 0; (i < atoms->nr); i++)
 746     {
 747         aps->setAtomProperty(epropMass, std::string(*atoms->resinfo[atoms->atom[i].resind].name),
 748                              std::string(*atoms->atomname[i]), &mm);
 749         mtot += mm;
 750     }
 751
 752     vol = det(box);
 753
 754     fprintf(stderr, "Volume                 :  %10g (nm^3)\n", vol);
 755     fprintf(stderr, "Density                :  %10g (g/l)\n", (mtot * 1e24) / (AVOGADRO * vol));
 756     fprintf(stderr, "Number of solvent molecules:  %5d   \n\n", nsol);
 757
 758     /* open topology file and append sol molecules */
 759     topinout = ftp2fn(efTOP, NFILE, fnm);
 760     if (ftp2bSet(efTOP, NFILE, fnm))
 761     {
 762         char temporary_filename[STRLEN];
 763         strncpy(temporary_filename, "temp.topXXXXXX", STRLEN);
 764
 765         fprintf(stderr, "Processing topology\n");
 766         fpin    = gmx_ffopen(topinout, "r");
 767         fpout   = gmx_fopen_temporary(temporary_filename);
 768         line    = 0;
 769         bSystem = false;
 770         while (fgets(buf, STRLEN, fpin))
 771         {
 772             line++;
 773             strcpy(buf2, buf);
 774             if ((temp = strchr(buf2, '\n')) != nullptr)
 775             {
 776                 temp[0] = '\0';
 777             }
 778             ltrim(buf2);
 779             if (buf2[0] == '[')
 780             {
 781                 buf2[0] = ' ';
 782                 if ((temp = strchr(buf2, '\n')) != nullptr)
 783                 {
 784                     temp[0] = '\0';
 785                 }
 786                 rtrim(buf2);
 787                 if (buf2[strlen(buf2) - 1] == ']')
 788                 {
 789                     buf2[strlen(buf2) - 1] = '\0';
 790                     ltrim(buf2);
 791                     rtrim(buf2);
 792                     bSystem = (gmx_strcasecmp(buf2, "system") == 0);
 793                 }
 794             }
 795             else if (bSystem && nsol && (buf[0] != ';'))
 796             {
 797                 /* if sol present, append "in water" to system name */
 798                 rtrim(buf2);
 799                 if (buf2[0] && (!strstr(buf2, " water")))
 800                 {
 801                     sprintf(buf, "%s in water\n", buf2);
 802                     bSystem = false;
 803                 }
 804             }
 805             fprintf(fpout, "%s", buf);
 806         }
 807         gmx_ffclose(fpin);
 808
 809         // Add new solvent molecules to the topology
 810         if (nsol > 0)
 811         {
 812             std::string currRes  = *atoms->resinfo[firstSolventResidueIndex].name;
 813             int         resCount = 0;
 814
 815             // Iterate through solvent molecules and increment a count until new resname found
 816             for (int i = firstSolventResidueIndex; i < atoms->nres; i++)
 817             {
 818                 if ((currRes == *atoms->resinfo[i].name))
 819                 {
 820                     resCount += 1;
 821                 }
 822                 else
 823                 {
 824                     // Change topology and restart count
 825                     fprintf(stdout,
 826                             "Adding line for %d solvent molecules with resname (%s) to "
 827                             "topology file (%s)\n",
 828                             resCount, currRes.c_str(), topinout);
 829                     fprintf(fpout, "%-15s %5d\n", currRes.c_str(), resCount);
 830                     currRes  = *atoms->resinfo[i].name;
 831                     resCount = 1;
 832                 }
 833             }
 834             // One more print needed for last residue type
 835             fprintf(stdout,
 836                     "Adding line for %d solvent molecules with resname (%s) to "
 837                     "topology file (%s)\n",
 838                     resCount, currRes.c_str(), topinout);
 839             fprintf(fpout, "%-15s %5d\n", currRes.c_str(), resCount);
 840         }
 841         gmx_ffclose(fpout);
 842         make_backup(topinout);
 843         gmx_file_rename(temporary_filename, topinout);
 844     }
 845 }
 846
 847 int gmx_solvate(int argc, char* argv[])
 848 {
 849     const char* desc[] = {
 850         "[THISMODULE] can do one of 2 things:[PAR]",
 851
 852         "1) Generate a box of solvent. Specify [TT]-cs[tt] and [TT]-box[tt].",
 853         "Or specify [TT]-cs[tt] and [TT]-cp[tt] with a structure file with",
 854         "a box, but without atoms.[PAR]",
 855
 856         "2) Solvate a solute configuration, e.g. a protein, in a bath of solvent ",
 857         "molecules. Specify [TT]-cp[tt] (solute) and [TT]-cs[tt] (solvent). ",
 858         "The box specified in the solute coordinate file ([TT]-cp[tt]) is used,",
 859         "unless [TT]-box[tt] is set.",
 860         "If you want the solute to be centered in the box,",
 861         "the program [gmx-editconf] has sophisticated options",
 862         "to change the box dimensions and center the solute.",
 863         "Solvent molecules are removed from the box where the ",
 864         "distance between any atom of the solute molecule(s) and any atom of ",
 865         "the solvent molecule is less than the sum of the scaled van der Waals",
 866         "radii of both atoms. A database ([TT]vdwradii.dat[tt]) of van der",
 867         "Waals radii is read by the program, and the resulting radii scaled",
 868         "by [TT]-scale[tt]. If radii are not found in the database, those",
 869         "atoms are assigned the (pre-scaled) distance [TT]-radius[tt].",
 870         "Note that the usefulness of those radii depends on the atom names,",
 871         "and thus varies widely with force field.",
 872         "",
 873         "The default solvent is Simple Point Charge water (SPC), with coordinates ",
 874         "from [TT]$GMXLIB/spc216.gro[tt]. These coordinates can also be used",
 875         "for other 3-site water models, since a short equibilibration will remove",
 876         "the small differences between the models.",
 877         "Other solvents are also supported, as well as mixed solvents. The",
 878         "only restriction to solvent types is that a solvent molecule consists",
 879         "of exactly one residue. The residue information in the coordinate",
 880         "files is used, and should therefore be more or less consistent.",
 881         "In practice this means that two subsequent solvent molecules in the ",
 882         "solvent coordinate file should have different residue number.",
 883         "The box of solute is built by stacking the coordinates read from",
 884         "the coordinate file. This means that these coordinates should be ",
 885         "equlibrated in periodic boundary conditions to ensure a good",
 886         "alignment of molecules on the stacking interfaces.",
 887         "The [TT]-maxsol[tt] option simply adds only the first [TT]-maxsol[tt]",
 888         "solvent molecules and leaves out the rest that would have fitted",
 889         "into the box. This can create a void that can cause problems later.",
 890         "Choose your volume wisely.[PAR]",
 891
 892         "Setting [TT]-shell[tt] larger than zero will place a layer of water of",
 893         "the specified thickness (nm) around the solute. Hint: it is a good",
 894         "idea to put the protein in the center of a box first (using [gmx-editconf]).",
 895         "[PAR]",
 896
 897         "Finally, [THISMODULE] will optionally remove lines from your topology file in ",
 898         "which a number of solvent molecules is already added, and adds a ",
 899         "line with the total number of solvent molecules in your coordinate file."
 900     };
 901
 902     const char* bugs[] = {
 903         "Molecules must be whole in the initial configurations.",
 904     };
 905
 906     /* parameter data */
 907     gmx_bool    bProt, bBox;
 908     const char *conf_prot, *confout;
 909
 910     t_filenm fnm[] = {
 911         { efSTX, "-cp", "protein", ffOPTRD },
 912         { efSTX, "-cs", "spc216", ffLIBRD },
 913         { efSTO, nullptr, nullptr, ffWRITE },
 914         { efTOP, nullptr, nullptr, ffOPTRW },
 915     };
 916 #define NFILE asize(fnm)
 917
 918     real              defaultDistance = 0.105, r_shell = 0, scaleFactor = 0.57;
 919     rvec              new_box                  = { 0.0, 0.0, 0.0 };
 920     gmx_bool          bReadV                   = FALSE;
 921     int               max_sol                  = 0;
 922     int               firstSolventResidueIndex = 0;
 923     gmx_output_env_t* oenv;
 924     t_pargs           pa[] = {
 925         { "-box", FALSE, etRVEC, { new_box }, "Box size (in nm)" },
 926         { "-radius", FALSE, etREAL, { &defaultDistance }, "Default van der Waals distance" },
 927         { "-scale",
 928           FALSE,
 929           etREAL,
 930           { &scaleFactor },
 931           "Scale factor to multiply Van der Waals radii from the database in "
 932           "share/gromacs/top/vdwradii.dat. The default value of 0.57 yields density close to 1000 "
 933           "g/l for proteins in water." },
 934         { "-shell", FALSE, etREAL, { &r_shell }, "Thickness of optional water layer around solute" },
 935         { "-maxsol",
 936           FALSE,
 937           etINT,
 938           { &max_sol },
 939           "Maximum number of solvent molecules to add if they fit in the box. If zero (default) "
 940           "this is ignored" },
 941         { "-vel", FALSE, etBOOL, { &bReadV }, "Keep velocities from input solute and solvent" },
 942     };
 943
 944     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa, asize(desc), desc,
 945                            asize(bugs), bugs, &oenv))
 946     {
 947         return 0;
 948     }
 949
 950     const char* solventFileName = opt2fn("-cs", NFILE, fnm);
 951     bProt                       = opt2bSet("-cp", NFILE, fnm);
 952     bBox                        = opt2parg_bSet("-box", asize(pa), pa);
 953
 954     /* check input */
 955     if (!bProt && !bBox)
 956     {
 957         gmx_fatal(FARGS,
 958                   "When no solute (-cp) is specified, "
 959                   "a box size (-box) must be specified");
 960     }
 961
 962     AtomProperties aps;
 963
 964     /* solute configuration data */
 965     gmx_mtop_t        top;
 966     std::vector<RVec> x, v;
 967     matrix            box     = { { 0 } };
 968     PbcType           pbcType = PbcType::Unset;
 969     t_atoms*          atoms;
 970     snew(atoms, 1);
 971     if (bProt)
 972     {
 973         /* Generate a solute configuration */
 974         conf_prot = opt2fn("-cp", NFILE, fnm);
 975         fprintf(stderr, "Reading solute configuration%s\n", bReadV ? " and velocities" : "");
 976         bool  bTprFileWasRead;
 977         rvec *temporaryX = nullptr, *temporaryV = nullptr;
 978         readConfAndTopology(conf_prot, &bTprFileWasRead, &top, &pbcType, &temporaryX,
 979                             bReadV ? &temporaryV : nullptr, box);
 980         *atoms = gmx_mtop_global_atoms(&top);
 981         x.assign(temporaryX, temporaryX + top.natoms);
 982         sfree(temporaryX);
 983         if (temporaryV)
 984         {
 985             v.assign(temporaryV, temporaryV + top.natoms);
 986             sfree(temporaryV);
 987         }
 988         else if (bReadV)
 989         {
 990             fprintf(stderr, "Note: no velocities found\n");
 991         }
 992         if (atoms->nr == 0)
 993         {
 994             fprintf(stderr, "Note: no atoms in %s\n", conf_prot);
 995             bProt = FALSE;
 996         }
 997         else
 998         {
 999             firstSolventResidueIndex = atoms->nres;
1000         }
1001     }
1002     PbcType pbcTypeForOutput = pbcType;
1003     if (bBox)
1004     {
1005         pbcTypeForOutput = PbcType::Xyz;
1006         clear_mat(box);
1007         box[XX][XX] = new_box[XX];
1008         box[YY][YY] = new_box[YY];
1009         box[ZZ][ZZ] = new_box[ZZ];
1010     }
1011     if (det(box) == 0)
1012     {
1013         gmx_fatal(FARGS,
1014                   "Undefined solute box.\nCreate one with gmx editconf "
1015                   "or give explicit -box command line option");
1016     }
1017
1018     add_solv(solventFileName, atoms, &top.symtab, &x, &v, pbcTypeForOutput, box, &aps,
1019              defaultDistance, scaleFactor, r_shell, max_sol);
1020
1021     /* write new configuration 1 to file confout */
1022     confout = ftp2fn(efSTO, NFILE, fnm);
1023     fprintf(stderr, "Writing generated configuration to %s\n", confout);
1024     const char* outputTitle = (bProt ? *top.name : "Generated by gmx solvate");
1025     write_sto_conf(confout, outputTitle, atoms, as_rvec_array(x.data()),
1026                    !v.empty() ? as_rvec_array(v.data()) : nullptr, pbcTypeForOutput, box);
1027
1028     /* print size of generated configuration */
1029     fprintf(stderr, "\nOutput configuration contains %d atoms in %d residues\n", atoms->nr, atoms->nres);
1030     update_top(atoms, firstSolventResidueIndex, box, NFILE, fnm, &aps);
1031
1032     done_atom(atoms);
1033     sfree(atoms);
1034     output_env_done(oenv);
1035
1036     return 0;
1037 }