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