src/gromacs/gmxpreprocess/solvate.cpp

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