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, AtomProperties *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,
 755                        int             firstSolventResidueIndex,
 756                        matrix          box,
 757                        int             NFILE,
 758                        t_filenm        fnm[],
 759                        AtomProperties *aps)
 760 {
 761     FILE        *fpin, *fpout;
 762     char         buf[STRLEN*2], buf2[STRLEN], *temp;
 763     const char  *topinout;
 764     int          line;
 765     bool         bSystem;
 766     int          i;
 767     double       mtot;
 768     real         vol, mm;
 769
 770     int          nsol = atoms->nres - firstSolventResidueIndex;
 771
 772     mtot = 0;
 773     for (i = 0; (i < atoms->nr); i++)
 774     {
 775         aps->setAtomProperty(epropMass,
 776                              std::string(*atoms->resinfo[atoms->atom[i].resind].name),
 777                              std::string(*atoms->atomname[i]), &mm);
 778         mtot += mm;
 779     }
 780
 781     vol = det(box);
 782
 783     fprintf(stderr, "Volume                 :  %10g (nm^3)\n", vol);
 784     fprintf(stderr, "Density                :  %10g (g/l)\n",
 785             (mtot*1e24)/(AVOGADRO*vol));
 786     fprintf(stderr, "Number of solvent molecules:  %5d   \n\n", nsol);
 787
 788     /* open topology file and append sol molecules */
 789     topinout  = ftp2fn(efTOP, NFILE, fnm);
 790     if (ftp2bSet(efTOP, NFILE, fnm) )
 791     {
 792         char temporary_filename[STRLEN];
 793         strncpy(temporary_filename, "temp.topXXXXXX", STRLEN);
 794
 795         fprintf(stderr, "Processing topology\n");
 796         fpin    = gmx_ffopen(topinout, "r");
 797         fpout   = gmx_fopen_temporary(temporary_filename);
 798         line    = 0;
 799         bSystem = false;
 800         while (fgets(buf, STRLEN, fpin))
 801         {
 802             line++;
 803             strcpy(buf2, buf);
 804             if ((temp = strchr(buf2, '\n')) != nullptr)
 805             {
 806                 temp[0] = '\0';
 807             }
 808             ltrim(buf2);
 809             if (buf2[0] == '[')
 810             {
 811                 buf2[0] = ' ';
 812                 if ((temp = strchr(buf2, '\n')) != nullptr)
 813                 {
 814                     temp[0] = '\0';
 815                 }
 816                 rtrim(buf2);
 817                 if (buf2[strlen(buf2)-1] == ']')
 818                 {
 819                     buf2[strlen(buf2)-1] = '\0';
 820                     ltrim(buf2);
 821                     rtrim(buf2);
 822                     bSystem    = (gmx_strcasecmp(buf2, "system") == 0);
 823                 }
 824             }
 825             else if (bSystem && nsol && (buf[0] != ';') )
 826             {
 827                 /* if sol present, append "in water" to system name */
 828                 rtrim(buf2);
 829                 if (buf2[0] && (!strstr(buf2, " water")) )
 830                 {
 831                     sprintf(buf, "%s in water\n", buf2);
 832                     bSystem = false;
 833                 }
 834             }
 835             fprintf(fpout, "%s", buf);
 836         }
 837         gmx_ffclose(fpin);
 838
 839         // Add new solvent molecules to the topology
 840         if (nsol > 0)
 841         {
 842             std::string currRes  = *atoms->resinfo[firstSolventResidueIndex].name;
 843             int         resCount = 0;
 844
 845             // Iterate through solvent molecules and increment a count until new resname found
 846             for (int i = firstSolventResidueIndex; i < atoms->nres; i++)
 847             {
 848                 if ((currRes == *atoms->resinfo[i].name))
 849                 {
 850                     resCount += 1;
 851                 }
 852                 else
 853                 {
 854                     // Change topology and restart count
 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                     currRes  = *atoms->resinfo[i].name;
 859                     resCount = 1;
 860                 }
 861             }
 862             // One more print needed for last residue type
 863             fprintf(stdout, "Adding line for %d solvent molecules with resname (%s) to "
 864                     "topology file (%s)\n", resCount, currRes.c_str(), topinout);
 865             fprintf(fpout, "%-15s %5d\n", currRes.c_str(), resCount);
 866         }
 867         gmx_ffclose(fpout);
 868         make_backup(topinout);
 869         gmx_file_rename(temporary_filename, topinout);
 870     }
 871 }
 872
 873 int gmx_solvate(int argc, char *argv[])
 874 {
 875     const char *desc[] = {
 876         "[THISMODULE] can do one of 2 things:[PAR]",
 877
 878         "1) Generate a box of solvent. Specify [TT]-cs[tt] and [TT]-box[tt].",
 879         "Or specify [TT]-cs[tt] and [TT]-cp[tt] with a structure file with",
 880         "a box, but without atoms.[PAR]",
 881
 882         "2) Solvate a solute configuration, e.g. a protein, in a bath of solvent ",
 883         "molecules. Specify [TT]-cp[tt] (solute) and [TT]-cs[tt] (solvent). ",
 884         "The box specified in the solute coordinate file ([TT]-cp[tt]) is used,",
 885         "unless [TT]-box[tt] is set.",
 886         "If you want the solute to be centered in the box,",
 887         "the program [gmx-editconf] has sophisticated options",
 888         "to change the box dimensions and center the solute.",
 889         "Solvent molecules are removed from the box where the ",
 890         "distance between any atom of the solute molecule(s) and any atom of ",
 891         "the solvent molecule is less than the sum of the scaled van der Waals",
 892         "radii of both atoms. A database ([TT]vdwradii.dat[tt]) of van der",
 893         "Waals radii is read by the program, and the resulting radii scaled",
 894         "by [TT]-scale[tt]. If radii are not found in the database, those",
 895         "atoms are assigned the (pre-scaled) distance [TT]-radius[tt].",
 896         "Note that the usefulness of those radii depends on the atom names,",
 897         "and thus varies widely with force field.",
 898         "",
 899         "The default solvent is Simple Point Charge water (SPC), with coordinates ",
 900         "from [TT]$GMXLIB/spc216.gro[tt]. These coordinates can also be used",
 901         "for other 3-site water models, since a short equibilibration will remove",
 902         "the small differences between the models.",
 903         "Other solvents are also supported, as well as mixed solvents. The",
 904         "only restriction to solvent types is that a solvent molecule consists",
 905         "of exactly one residue. The residue information in the coordinate",
 906         "files is used, and should therefore be more or less consistent.",
 907         "In practice this means that two subsequent solvent molecules in the ",
 908         "solvent coordinate file should have different residue number.",
 909         "The box of solute is built by stacking the coordinates read from",
 910         "the coordinate file. This means that these coordinates should be ",
 911         "equlibrated in periodic boundary conditions to ensure a good",
 912         "alignment of molecules on the stacking interfaces.",
 913         "The [TT]-maxsol[tt] option simply adds only the first [TT]-maxsol[tt]",
 914         "solvent molecules and leaves out the rest that would have fitted",
 915         "into the box. This can create a void that can cause problems later.",
 916         "Choose your volume wisely.[PAR]",
 917
 918         "Setting [TT]-shell[tt] larger than zero will place a layer of water of",
 919         "the specified thickness (nm) around the solute. Hint: it is a good",
 920         "idea to put the protein in the center of a box first (using [gmx-editconf]).",
 921         "[PAR]",
 922
 923         "Finally, [THISMODULE] will optionally remove lines from your topology file in ",
 924         "which a number of solvent molecules is already added, and adds a ",
 925         "line with the total number of solvent molecules in your coordinate file."
 926     };
 927
 928     const char *bugs[] = {
 929         "Molecules must be whole in the initial configurations.",
 930     };
 931
 932     /* parameter data */
 933     gmx_bool       bProt, bBox;
 934     const char    *conf_prot, *confout;
 935
 936     t_filenm       fnm[] = {
 937         { efSTX, "-cp", "protein", ffOPTRD },
 938         { efSTX, "-cs", "spc216",  ffLIBRD},
 939         { efSTO, nullptr,  nullptr,      ffWRITE},
 940         { efTOP, nullptr,  nullptr,      ffOPTRW},
 941     };
 942 #define NFILE asize(fnm)
 943
 944     real              defaultDistance          = 0.105, r_shell = 0, scaleFactor = 0.57;
 945     rvec              new_box                  = {0.0, 0.0, 0.0};
 946     gmx_bool          bReadV                   = FALSE;
 947     int               max_sol                  = 0;
 948     int               firstSolventResidueIndex = 0;
 949     gmx_output_env_t *oenv;
 950     t_pargs           pa[]              = {
 951         { "-box",    FALSE, etRVEC, {new_box},
 952           "Box size (in nm)" },
 953         { "-radius",   FALSE, etREAL, {&defaultDistance},
 954           "Default van der Waals distance"},
 955         { "-scale", FALSE, etREAL, {&scaleFactor},
 956           "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." },
 957         { "-shell",  FALSE, etREAL, {&r_shell},
 958           "Thickness of optional water layer around solute" },
 959         { "-maxsol", FALSE, etINT,  {&max_sol},
 960           "Maximum number of solvent molecules to add if they fit in the box. If zero (default) this is ignored" },
 961         { "-vel",    FALSE, etBOOL, {&bReadV},
 962           "Keep velocities from input solute and solvent" },
 963     };
 964
 965     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
 966                            asize(desc), desc, asize(bugs), bugs, &oenv))
 967     {
 968         return 0;
 969     }
 970
 971     const char *solventFileName = opt2fn("-cs", NFILE, fnm);
 972     bProt     = opt2bSet("-cp", NFILE, fnm);
 973     bBox      = opt2parg_bSet("-box", asize(pa), pa);
 974
 975     /* check input */
 976     if (!bProt && !bBox)
 977     {
 978         gmx_fatal(FARGS, "When no solute (-cp) is specified, "
 979                   "a box size (-box) must be specified");
 980     }
 981
 982     AtomProperties           aps;
 983
 984     /* solute configuration data */
 985     gmx_mtop_t        top;
 986     std::vector<RVec> x, v;
 987     matrix            box  = {{ 0 }};
 988     int               ePBC = -1;
 989     t_atoms          *atoms;
 990     snew(atoms, 1);
 991     if (bProt)
 992     {
 993         /* Generate a solute configuration */
 994         conf_prot = opt2fn("-cp", NFILE, fnm);
 995         fprintf(stderr, "Reading solute configuration%s\n",
 996                 bReadV ? " and velocities" : "");
 997         bool  bTprFileWasRead;
 998         rvec *temporaryX = nullptr, *temporaryV = nullptr;
 999         readConfAndTopology(conf_prot, &bTprFileWasRead, &top,
1000                             &ePBC, &temporaryX, bReadV ? &temporaryV : nullptr, box);
1001         *atoms = gmx_mtop_global_atoms(&top);
1002         x.assign(temporaryX, temporaryX + top.natoms);
1003         sfree(temporaryX);
1004         if (temporaryV)
1005         {
1006             v.assign(temporaryV, temporaryV + top.natoms);
1007             sfree(temporaryV);
1008         }
1009         else if (bReadV)
1010         {
1011             fprintf(stderr, "Note: no velocities found\n");
1012         }
1013         if (atoms->nr == 0)
1014         {
1015             fprintf(stderr, "Note: no atoms in %s\n", conf_prot);
1016             bProt = FALSE;
1017         }
1018         else
1019         {
1020             firstSolventResidueIndex = atoms->nres;
1021         }
1022     }
1023     int ePBCForOutput = ePBC;
1024     if (bBox)
1025     {
1026         ePBCForOutput = epbcXYZ;
1027         clear_mat(box);
1028         box[XX][XX] = new_box[XX];
1029         box[YY][YY] = new_box[YY];
1030         box[ZZ][ZZ] = new_box[ZZ];
1031     }
1032     if (det(box) == 0)
1033     {
1034         gmx_fatal(FARGS, "Undefined solute box.\nCreate one with gmx editconf "
1035                   "or give explicit -box command line option");
1036     }
1037
1038     add_solv(solventFileName, atoms, &top.symtab, &x, &v, ePBCForOutput, box,
1039              &aps, defaultDistance, scaleFactor, r_shell, max_sol);
1040
1041     /* write new configuration 1 to file confout */
1042     confout = ftp2fn(efSTO, NFILE, fnm);
1043     fprintf(stderr, "Writing generated configuration to %s\n", confout);
1044     const char *outputTitle = (bProt ? *top.name : "Generated by gmx solvate");
1045     write_sto_conf(confout, outputTitle, atoms, as_rvec_array(x.data()),
1046                    !v.empty() ? as_rvec_array(v.data()) : nullptr, ePBCForOutput, box);
1047
1048     /* print size of generated configuration */
1049     fprintf(stderr, "\nOutput configuration contains %d atoms in %d residues\n",
1050             atoms->nr, atoms->nres);
1051     update_top(atoms, firstSolventResidueIndex, box, NFILE, fnm, &aps);
1052
1053     done_atom(atoms);
1054     sfree(atoms);
1055     output_env_done(oenv);
1056
1057     return 0;
1058 }