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