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