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