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, 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     bool        bSystem, bMolecules, bSkip;
 740     int         i, nsol = 0;
 741     double      mtot;
 742     real        vol, mm;
 743
 744     for (i = 0; (i < atoms->nres); i++)
 745     {
 746         /* calculate number of SOLvent molecules */
 747         if ( (strcmp(*atoms->resinfo[i].name, "SOL") == 0) ||
 748              (strcmp(*atoms->resinfo[i].name, "WAT") == 0) ||
 749              (strcmp(*atoms->resinfo[i].name, "HOH") == 0) )
 750         {
 751             nsol++;
 752         }
 753     }
 754     mtot = 0;
 755     for (i = 0; (i < atoms->nr); i++)
 756     {
 757         gmx_atomprop_query(aps, epropMass,
 758                            *atoms->resinfo[atoms->atom[i].resind].name,
 759                            *atoms->atomname[i], &mm);
 760         mtot += mm;
 761     }
 762
 763     vol = det(box);
 764
 765     fprintf(stderr, "Volume                 :  %10g (nm^3)\n", vol);
 766     fprintf(stderr, "Density                :  %10g (g/l)\n",
 767             (mtot*1e24)/(AVOGADRO*vol));
 768     fprintf(stderr, "Number of SOL molecules:  %5d   \n\n", nsol);
 769
 770     /* open topology file and append sol molecules */
 771     topinout  = ftp2fn(efTOP, NFILE, fnm);
 772     if (ftp2bSet(efTOP, NFILE, fnm) )
 773     {
 774         char temporary_filename[STRLEN];
 775         strncpy(temporary_filename, "temp.topXXXXXX", STRLEN);
 776
 777         fprintf(stderr, "Processing topology\n");
 778         fpin    = gmx_ffopen(topinout, "r");
 779         fpout   = gmx_fopen_temporary(temporary_filename);
 780         line    = 0;
 781         bSystem = bMolecules = false;
 782         while (fgets(buf, STRLEN, fpin))
 783         {
 784             bSkip = false;
 785             line++;
 786             strcpy(buf2, buf);
 787             if ((temp = strchr(buf2, '\n')) != nullptr)
 788             {
 789                 temp[0] = '\0';
 790             }
 791             ltrim(buf2);
 792             if (buf2[0] == '[')
 793             {
 794                 buf2[0] = ' ';
 795                 if ((temp = strchr(buf2, '\n')) != nullptr)
 796                 {
 797                     temp[0] = '\0';
 798                 }
 799                 rtrim(buf2);
 800                 if (buf2[strlen(buf2)-1] == ']')
 801                 {
 802                     buf2[strlen(buf2)-1] = '\0';
 803                     ltrim(buf2);
 804                     rtrim(buf2);
 805                     bSystem    = (gmx_strcasecmp(buf2, "system") == 0);
 806                     bMolecules = (gmx_strcasecmp(buf2, "molecules") == 0);
 807                 }
 808             }
 809             else if (bSystem && nsol && (buf[0] != ';') )
 810             {
 811                 /* if sol present, append "in water" to system name */
 812                 rtrim(buf2);
 813                 if (buf2[0] && (!strstr(buf2, " water")) )
 814                 {
 815                     sprintf(buf, "%s in water\n", buf2);
 816                     bSystem = false;
 817                 }
 818             }
 819             else if (bMolecules)
 820             {
 821                 /* check if this is a line with solvent molecules */
 822                 sscanf(buf, "%4095s", buf2);
 823                 if (strcmp(buf2, "SOL") == 0)
 824                 {
 825                     sscanf(buf, "%*4095s %20d", &i);
 826                     nsol -= i;
 827                     if (nsol < 0)
 828                     {
 829                         bSkip = true;
 830                         nsol += i;
 831                     }
 832                 }
 833             }
 834             if (bSkip)
 835             {
 836                 if ((temp = strchr(buf, '\n')) != nullptr)
 837                 {
 838                     temp[0] = '\0';
 839                 }
 840                 fprintf(stdout, "Removing line #%d '%s' from topology file (%s)\n",
 841                         line, buf, topinout);
 842             }
 843             else
 844             {
 845                 fprintf(fpout, "%s", buf);
 846             }
 847         }
 848         gmx_ffclose(fpin);
 849         if (nsol)
 850         {
 851             fprintf(stdout, "Adding line for %d solvent molecules to "
 852                     "topology file (%s)\n", nsol, topinout);
 853             fprintf(fpout, "%-15s %5d\n", "SOL", nsol);
 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     gmx_atomprop_t aps;
 924
 925     /* solute configuration data */
 926     t_topology *top;
 927     int         ePBC = -1;
 928     matrix      box;
 929
 930     t_filenm    fnm[] = {
 931         { efSTX, "-cp", "protein", ffOPTRD },
 932         { efSTX, "-cs", "spc216",  ffLIBRD},
 933         { efSTO, nullptr,  nullptr,      ffWRITE},
 934         { efTOP, nullptr,  nullptr,      ffOPTRW},
 935     };
 936 #define NFILE asize(fnm)
 937
 938     real              defaultDistance = 0.105, r_shell = 0, scaleFactor = 0.57;
 939     rvec              new_box         = {0.0, 0.0, 0.0};
 940     gmx_bool          bReadV          = FALSE;
 941     int               max_sol         = 0;
 942     gmx_output_env_t *oenv;
 943     t_pargs           pa[]              = {
 944         { "-box",    FALSE, etRVEC, {new_box},
 945           "Box size (in nm)" },
 946         { "-radius",   FALSE, etREAL, {&defaultDistance},
 947           "Default van der Waals distance"},
 948         { "-scale", FALSE, etREAL, {&scaleFactor},
 949           "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." },
 950         { "-shell",  FALSE, etREAL, {&r_shell},
 951           "Thickness of optional water layer around solute" },
 952         { "-maxsol", FALSE, etINT,  {&max_sol},
 953           "Maximum number of solvent molecules to add if they fit in the box. If zero (default) this is ignored" },
 954         { "-vel",    FALSE, etBOOL, {&bReadV},
 955           "Keep velocities from input solute and solvent" },
 956     };
 957
 958     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
 959                            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, "When no solute (-cp) is specified, "
 972                   "a box size (-box) must be specified");
 973     }
 974
 975     aps = gmx_atomprop_init();
 976
 977     std::vector<RVec> x;
 978     std::vector<RVec> v;
 979     snew(top, 1);
 980     if (bProt)
 981     {
 982         /* Generate a solute configuration */
 983         conf_prot = opt2fn("-cp", NFILE, fnm);
 984         readConformation(conf_prot, top, &x,
 985                          bReadV ? &v : nullptr, &ePBC, box, "solute");
 986         if (bReadV && v.empty())
 987         {
 988             fprintf(stderr, "Note: no velocities found\n");
 989         }
 990         if (top->atoms.nr == 0)
 991         {
 992             fprintf(stderr, "Note: no atoms in %s\n", conf_prot);
 993             bProt = FALSE;
 994         }
 995     }
 996     if (bBox)
 997     {
 998         ePBC = epbcXYZ;
 999         clear_mat(box);
1000         box[XX][XX] = new_box[XX];
1001         box[YY][YY] = new_box[YY];
1002         box[ZZ][ZZ] = new_box[ZZ];
1003     }
1004     if (det(box) == 0)
1005     {
1006         gmx_fatal(FARGS, "Undefined solute box.\nCreate one with gmx editconf "
1007                   "or give explicit -box command line option");
1008     }
1009
1010     add_solv(solventFileName, top, &x, &v, ePBC, box,
1011              aps, defaultDistance, scaleFactor, r_shell, max_sol);
1012
1013     /* write new configuration 1 to file confout */
1014     confout = ftp2fn(efSTO, NFILE, fnm);
1015     fprintf(stderr, "Writing generated configuration to %s\n", confout);
1016     const char *outputTitle = (bProt ? *top->name : "Generated by gmx solvate");
1017     write_sto_conf(confout, outputTitle, &top->atoms, as_rvec_array(x.data()),
1018                    !v.empty() ? as_rvec_array(v.data()) : nullptr, ePBC, box);
1019
1020     /* print size of generated configuration */
1021     fprintf(stderr, "\nOutput configuration contains %d atoms in %d residues\n",
1022             top->atoms.nr, top->atoms.nres);
1023     update_top(&top->atoms, box, NFILE, fnm, aps);
1024
1025     gmx_atomprop_destroy(aps);
1026     done_top(top);
1027     sfree(top);
1028     output_env_done(oenv);
1029
1030     return 0;
1031 }