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