src/gromacs/gmxpreprocess/solvate.cpp

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