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, 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 <string.h>
  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/pbc.h"
  55 #include "gromacs/selection/nbsearch.h"
  56 #include "gromacs/topology/atomprop.h"
  57 #include "gromacs/topology/atoms.h"
  58 #include "gromacs/topology/atomsbuilder.h"
  59 #include "gromacs/topology/topology.h"
  60 #include "gromacs/utility/arraysize.h"
  61 #include "gromacs/utility/cstringutil.h"
  62 #include "gromacs/utility/fatalerror.h"
  63 #include "gromacs/utility/futil.h"
  64 #include "gromacs/utility/gmxassert.h"
  65 #include "gromacs/utility/smalloc.h"
  66
  67 using gmx::RVec;
  68
  69 typedef struct {
  70     char *name;
  71     int   natoms;
  72     int   nmol;
  73     int   i, i0;
  74     int   res0;
  75 } t_moltypes;
  76
  77 static void sort_molecule(t_atoms **atoms_solvt, std::vector<RVec> *x,
  78                           std::vector<RVec> *v)
  79 {
  80     int         atnr, i, j, moltp = 0, nrmoltypes, resi_o, resi_n, resnr;
  81     t_moltypes *moltypes;
  82     t_atoms    *atoms, *newatoms;
  83
  84     fprintf(stderr, "Sorting configuration\n");
  85
  86     atoms = *atoms_solvt;
  87
  88     /* copy each residue from *atoms to a molecule in *molecule */
  89     moltypes   = NULL;
  90     nrmoltypes = 0;
  91     for (i = 0; i < atoms->nr; i++)
  92     {
  93         if ( (i == 0) || (atoms->atom[i].resind != atoms->atom[i-1].resind) )
  94         {
  95             /* see if this was a molecule type we haven't had yet: */
  96             moltp = -1;
  97             for (j = 0; (j < nrmoltypes) && (moltp == -1); j++)
  98             {
  99                 /* cppcheck-suppress nullPointer
 100                  * moltypes is guaranteed to be allocated because otherwise
 101                  * nrmoltypes is 0. */
 102                 if (strcmp(*(atoms->resinfo[atoms->atom[i].resind].name), moltypes[j].name) == 0)
 103                 {
 104                     moltp = j;
 105                 }
 106             }
 107             if (moltp == -1)
 108             {
 109                 moltp = nrmoltypes;
 110                 nrmoltypes++;
 111                 srenew(moltypes, nrmoltypes);
 112                 moltypes[moltp].name = *(atoms->resinfo[atoms->atom[i].resind].name);
 113                 atnr                 = 0;
 114                 while ((i+atnr < atoms->nr) &&
 115                        (atoms->atom[i].resind == atoms->atom[i+atnr].resind))
 116                 {
 117                     atnr++;
 118                 }
 119                 moltypes[moltp].natoms = atnr;
 120                 moltypes[moltp].nmol   = 0;
 121             }
 122             moltypes[moltp].nmol++;
 123         }
 124     }
 125
 126     fprintf(stderr, "Found %d%s molecule type%s:\n",
 127             nrmoltypes, nrmoltypes == 1 ? "" : " different", nrmoltypes == 1 ? "" : "s");
 128     for (j = 0; j < nrmoltypes; j++)
 129     {
 130         if (j == 0)
 131         {
 132             moltypes[j].res0 = 0;
 133         }
 134         else
 135         {
 136             moltypes[j].res0 = moltypes[j-1].res0+moltypes[j-1].nmol;
 137         }
 138         fprintf(stderr, "%7s (%4d atoms): %5d residues\n",
 139                 moltypes[j].name, moltypes[j].natoms, moltypes[j].nmol);
 140     }
 141
 142     /* if we have only 1 moleculetype, we don't have to sort */
 143     if (nrmoltypes > 1)
 144     {
 145         /* find out which molecules should go where: */
 146         moltypes[0].i = moltypes[0].i0 = 0;
 147         for (j = 1; j < nrmoltypes; j++)
 148         {
 149             moltypes[j].i      =
 150                 moltypes[j].i0 =
 151                     moltypes[j-1].i0+moltypes[j-1].natoms*moltypes[j-1].nmol;
 152         }
 153
 154         /* now put them there: */
 155         snew(newatoms, 1);
 156         init_t_atoms(newatoms, atoms->nr, FALSE);
 157         newatoms->nres = atoms->nres;
 158         snew(newatoms->resinfo, atoms->nres);
 159         std::vector<RVec> newx(x->size());
 160         std::vector<RVec> newv(v->size());
 161
 162         resi_n = 0;
 163         resnr  = 1;
 164         j      = 0;
 165         for (moltp = 0; moltp < nrmoltypes; moltp++)
 166         {
 167             i = 0;
 168             while (i < atoms->nr)
 169             {
 170                 resi_o = atoms->atom[i].resind;
 171                 if (strcmp(*atoms->resinfo[resi_o].name, moltypes[moltp].name) == 0)
 172                 {
 173                     /* Copy the residue info */
 174                     newatoms->resinfo[resi_n]    = atoms->resinfo[resi_o];
 175                     newatoms->resinfo[resi_n].nr = resnr;
 176                     /* Copy the atom info */
 177                     do
 178                     {
 179                         newatoms->atom[j]        = atoms->atom[i];
 180                         newatoms->atomname[j]    = atoms->atomname[i];
 181                         newatoms->atom[j].resind = resi_n;
 182                         copy_rvec((*x)[i], newx[j]);
 183                         if (!v->empty())
 184                         {
 185                             copy_rvec((*v)[i], newv[j]);
 186                         }
 187                         i++;
 188                         j++;
 189                     }
 190                     while (i < atoms->nr && atoms->atom[i].resind == resi_o);
 191                     /* Increase the new residue counters */
 192                     resi_n++;
 193                     resnr++;
 194                 }
 195                 else
 196                 {
 197                     /* Skip this residue */
 198                     do
 199                     {
 200                         i++;
 201                     }
 202                     while (i < atoms->nr && atoms->atom[i].resind == resi_o);
 203                 }
 204             }
 205         }
 206
 207         /* put them back into the original arrays and throw away temporary arrays */
 208         sfree(atoms->atomname);
 209         sfree(atoms->resinfo);
 210         sfree(atoms->atom);
 211         sfree(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, NULL);
 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  * Removes solvent molecules that overlap with the solute, and optionally also
 495  * those that are outside a given shell radius from the solute.
 496  *
 497  * \param[in,out] atoms      Solvent atoms.
 498  * \param[in,out] x          Solvent positions.
 499  * \param[in,out] v          Solvent velocities (can be empty).
 500  * \param[in,out] r          Solvent exclusion radii.
 501  * \param[in]     pbc        PBC information.
 502  * \param[in]     x_solute   Solute positions.
 503  * \param[in]     r_solute   Solute exclusion radii.
 504  * \param[in]     rshell     If >0, only keep solvent atoms within a shell of
 505  *     this size from the solute.
 506  */
 507 static void removeSoluteOverlap(t_atoms *atoms, std::vector<RVec> *x,
 508                                 std::vector<RVec> *v, std::vector<real> *r,
 509                                 const t_pbc &pbc,
 510                                 const std::vector<RVec> &x_solute,
 511                                 const std::vector<real> &r_solute,
 512                                 real rshell)
 513 {
 514     const real                          maxRadius1
 515         = *std::max_element(r->begin(), r->end());
 516     const real                          maxRadius2
 517         = *std::max_element(r_solute.begin(), r_solute.end());
 518
 519     gmx::AtomsRemover                   remover(*atoms);
 520     // If rshell is >0, the neighborhood search looks at all pairs
 521     // within rshell, and unmarks those that are within the cutoff.
 522     // This line marks everything, so that solvent outside rshell remains
 523     // marked after the loop.
 524     // Without rshell, the neighborhood search only marks the overlapping
 525     // solvent atoms, and all others are left alone.
 526     if (rshell > 0.0)
 527     {
 528         remover.markAll();
 529     }
 530
 531     gmx::AnalysisNeighborhood           nb;
 532     nb.setCutoff(std::max(maxRadius1 + maxRadius2, rshell));
 533     gmx::AnalysisNeighborhoodPositions  posSolute(x_solute);
 534     gmx::AnalysisNeighborhoodSearch     search     = nb.initSearch(&pbc, posSolute);
 535     gmx::AnalysisNeighborhoodPositions  pos(*x);
 536     gmx::AnalysisNeighborhoodPairSearch pairSearch = search.startPairSearch(pos);
 537     gmx::AnalysisNeighborhoodPair       pair;
 538     while (pairSearch.findNextPair(&pair))
 539     {
 540         if (remover.isMarked(pair.testIndex()))
 541         {
 542             pairSearch.skipRemainingPairsForTestPosition();
 543             continue;
 544         }
 545         const real r1      = r_solute[pair.refIndex()];
 546         const real r2      = (*r)[pair.testIndex()];
 547         const bool bRemove = (pair.distance2() < gmx::square(r1 + r2));
 548         remover.markResidue(*atoms, pair.testIndex(), bRemove);
 549     }
 550
 551     remover.removeMarkedElements(x);
 552     if (!v->empty())
 553     {
 554         remover.removeMarkedElements(v);
 555     }
 556     remover.removeMarkedElements(r);
 557     const int originalAtomCount = atoms->nr;
 558     remover.removeMarkedAtoms(atoms);
 559     fprintf(stderr, "Removed %d solvent atoms due to solute-solvent overlap\n",
 560             originalAtomCount - atoms->nr);
 561 }
 562
 563 /*! \brief
 564  * Removes a given number of solvent residues.
 565  *
 566  * \param[in,out] atoms           Solvent atoms.
 567  * \param[in,out] x               Solvent positions.
 568  * \param[in,out] v               Solvent velocities (can be empty).
 569  * \param[in]     numberToRemove  Number of residues to remove.
 570  *
 571  * This function is called last in the process of creating the solvent box,
 572  * so it does not operate on the exclusion radii, as no code after this needs
 573  * them.
 574  */
 575 static void removeExtraSolventMolecules(t_atoms *atoms, std::vector<RVec> *x,
 576                                         std::vector<RVec> *v,
 577                                         int numberToRemove)
 578 {
 579     gmx::AtomsRemover remover(*atoms);
 580     // TODO: It might be nicer to remove a random set of residues, but
 581     // in practice this should give a roughly uniform spatial distribution.
 582     const int stride = atoms->nr / numberToRemove;
 583     for (int i = 0; i < numberToRemove; ++i)
 584     {
 585         int atomIndex = (i+1)*stride - 1;
 586         while (remover.isMarked(atomIndex))
 587         {
 588             ++atomIndex;
 589             if (atomIndex == atoms->nr)
 590             {
 591                 atomIndex = 0;
 592             }
 593         }
 594         remover.markResidue(*atoms, atomIndex, true);
 595     }
 596     remover.removeMarkedElements(x);
 597     if (!v->empty())
 598     {
 599         remover.removeMarkedElements(v);
 600     }
 601     remover.removeMarkedAtoms(atoms);
 602 }
 603
 604 static void add_solv(const char *fn, t_topology *top,
 605                      std::vector<RVec> *x, std::vector<RVec> *v,
 606                      int ePBC, matrix box, gmx_atomprop_t aps,
 607                      real defaultDistance, real scaleFactor,
 608                      real rshell, int max_sol)
 609 {
 610     t_topology       *top_solvt;
 611     std::vector<RVec> x_solvt;
 612     std::vector<RVec> v_solvt;
 613     int               ePBC_solvt;
 614     matrix            box_solvt;
 615
 616     char             *filename = gmxlibfn(fn);
 617     snew(top_solvt, 1);
 618     readConformation(filename, top_solvt, &x_solvt, !v->empty() ? &v_solvt : NULL,
 619                      &ePBC_solvt, box_solvt, "solvent");
 620     t_atoms *atoms_solvt = &top_solvt->atoms;
 621     if (0 == atoms_solvt->nr)
 622     {
 623         gmx_fatal(FARGS, "No solvent in %s, please check your input\n", filename);
 624     }
 625     sfree(filename);
 626     fprintf(stderr, "\n");
 627
 628     /* apply pbc for solvent configuration for whole molecules */
 629     rm_res_pbc(atoms_solvt, &x_solvt, box_solvt);
 630
 631     /* initialise distance arrays for solvent configuration */
 632     fprintf(stderr, "Initialising inter-atomic distances...\n");
 633     const std::vector<real> exclusionDistances(
 634             makeExclusionDistances(&top->atoms, aps, defaultDistance, scaleFactor));
 635     std::vector<real>       exclusionDistances_solvt(
 636             makeExclusionDistances(atoms_solvt, aps, defaultDistance, scaleFactor));
 637
 638     /* generate a new solvent configuration */
 639     fprintf(stderr, "Generating solvent configuration\n");
 640     t_pbc pbc;
 641     set_pbc(&pbc, ePBC, box);
 642     replicateSolventBox(atoms_solvt, &x_solvt, &v_solvt, &exclusionDistances_solvt,
 643                         box_solvt, box);
 644     if (ePBC != epbcNONE)
 645     {
 646         removeSolventBoxOverlap(atoms_solvt, &x_solvt, &v_solvt, &exclusionDistances_solvt, pbc);
 647     }
 648     if (top->atoms.nr > 0)
 649     {
 650         removeSoluteOverlap(atoms_solvt, &x_solvt, &v_solvt, &exclusionDistances_solvt, pbc,
 651                             *x, exclusionDistances, rshell);
 652     }
 653
 654     if (max_sol > 0 && atoms_solvt->nres > max_sol)
 655     {
 656         const int numberToRemove = atoms_solvt->nres - max_sol;
 657         removeExtraSolventMolecules(atoms_solvt, &x_solvt, &v_solvt, numberToRemove);
 658     }
 659
 660     /* Sort the solvent mixture, not the protein... */
 661     sort_molecule(&atoms_solvt, &x_solvt, &v_solvt);
 662
 663     // Merge the two configurations.
 664     x->insert(x->end(), x_solvt.begin(), x_solvt.end());
 665     if (!v->empty())
 666     {
 667         v->insert(v->end(), v_solvt.begin(), v_solvt.end());
 668     }
 669     {
 670         gmx::AtomsBuilder builder(&top->atoms, &top->symtab);
 671         builder.mergeAtoms(*atoms_solvt);
 672     }
 673     fprintf(stderr, "Generated solvent containing %d atoms in %d residues\n",
 674             atoms_solvt->nr, atoms_solvt->nres);
 675
 676     done_top(top_solvt);
 677     sfree(top_solvt);
 678 }
 679
 680 static void update_top(t_atoms *atoms, matrix box, int NFILE, t_filenm fnm[],
 681                        gmx_atomprop_t aps)
 682 {
 683     FILE       *fpin, *fpout;
 684     char        buf[STRLEN], buf2[STRLEN], *temp;
 685     const char *topinout;
 686     int         line;
 687     gmx_bool    bSystem, bMolecules, bSkip;
 688     int         i, nsol = 0;
 689     double      mtot;
 690     real        vol, mm;
 691
 692     for (i = 0; (i < atoms->nres); i++)
 693     {
 694         /* calculate number of SOLvent molecules */
 695         if ( (strcmp(*atoms->resinfo[i].name, "SOL") == 0) ||
 696              (strcmp(*atoms->resinfo[i].name, "WAT") == 0) ||
 697              (strcmp(*atoms->resinfo[i].name, "HOH") == 0) )
 698         {
 699             nsol++;
 700         }
 701     }
 702     mtot = 0;
 703     for (i = 0; (i < atoms->nr); i++)
 704     {
 705         gmx_atomprop_query(aps, epropMass,
 706                            *atoms->resinfo[atoms->atom[i].resind].name,
 707                            *atoms->atomname[i], &mm);
 708         mtot += mm;
 709     }
 710
 711     vol = det(box);
 712
 713     fprintf(stderr, "Volume                 :  %10g (nm^3)\n", vol);
 714     fprintf(stderr, "Density                :  %10g (g/l)\n",
 715             (mtot*1e24)/(AVOGADRO*vol));
 716     fprintf(stderr, "Number of SOL molecules:  %5d   \n\n", nsol);
 717
 718     /* open topology file and append sol molecules */
 719     topinout  = ftp2fn(efTOP, NFILE, fnm);
 720     if (ftp2bSet(efTOP, NFILE, fnm) )
 721     {
 722         char temporary_filename[STRLEN];
 723         strncpy(temporary_filename, "temp.topXXXXXX", STRLEN);
 724
 725         fprintf(stderr, "Processing topology\n");
 726         fpin    = gmx_ffopen(topinout, "r");
 727         fpout   = gmx_fopen_temporary(temporary_filename);
 728         line    = 0;
 729         bSystem = bMolecules = FALSE;
 730         while (fgets(buf, STRLEN, fpin))
 731         {
 732             bSkip = FALSE;
 733             line++;
 734             strcpy(buf2, buf);
 735             if ((temp = strchr(buf2, '\n')) != NULL)
 736             {
 737                 temp[0] = '\0';
 738             }
 739             ltrim(buf2);
 740             if (buf2[0] == '[')
 741             {
 742                 buf2[0] = ' ';
 743                 if ((temp = strchr(buf2, '\n')) != NULL)
 744                 {
 745                     temp[0] = '\0';
 746                 }
 747                 rtrim(buf2);
 748                 if (buf2[strlen(buf2)-1] == ']')
 749                 {
 750                     buf2[strlen(buf2)-1] = '\0';
 751                     ltrim(buf2);
 752                     rtrim(buf2);
 753                     bSystem    = (gmx_strcasecmp(buf2, "system") == 0);
 754                     bMolecules = (gmx_strcasecmp(buf2, "molecules") == 0);
 755                 }
 756             }
 757             else if (bSystem && nsol && (buf[0] != ';') )
 758             {
 759                 /* if sol present, append "in water" to system name */
 760                 rtrim(buf2);
 761                 if (buf2[0] && (!strstr(buf2, " water")) )
 762                 {
 763                     sprintf(buf, "%s in water\n", buf2);
 764                     bSystem = FALSE;
 765                 }
 766             }
 767             else if (bMolecules)
 768             {
 769                 /* check if this is a line with solvent molecules */
 770                 sscanf(buf, "%4095s", buf2);
 771                 if (strcmp(buf2, "SOL") == 0)
 772                 {
 773                     sscanf(buf, "%*4095s %20d", &i);
 774                     nsol -= i;
 775                     if (nsol < 0)
 776                     {
 777                         bSkip = TRUE;
 778                         nsol += i;
 779                     }
 780                 }
 781             }
 782             if (bSkip)
 783             {
 784                 if ((temp = strchr(buf, '\n')) != NULL)
 785                 {
 786                     temp[0] = '\0';
 787                 }
 788                 fprintf(stdout, "Removing line #%d '%s' from topology file (%s)\n",
 789                         line, buf, topinout);
 790             }
 791             else
 792             {
 793                 fprintf(fpout, "%s", buf);
 794             }
 795         }
 796         gmx_ffclose(fpin);
 797         if (nsol)
 798         {
 799             fprintf(stdout, "Adding line for %d solvent molecules to "
 800                     "topology file (%s)\n", nsol, topinout);
 801             fprintf(fpout, "%-15s %5d\n", "SOL", nsol);
 802         }
 803         gmx_ffclose(fpout);
 804         make_backup(topinout);
 805         gmx_file_rename(temporary_filename, topinout);
 806     }
 807 }
 808
 809 int gmx_solvate(int argc, char *argv[])
 810 {
 811     const char *desc[] = {
 812         "[THISMODULE] can do one of 2 things:[PAR]",
 813
 814         "1) Generate a box of solvent. Specify [TT]-cs[tt] and [TT]-box[tt].",
 815         "Or specify [TT]-cs[tt] and [TT]-cp[tt] with a structure file with",
 816         "a box, but without atoms.[PAR]",
 817
 818         "2) Solvate a solute configuration, e.g. a protein, in a bath of solvent ",
 819         "molecules. Specify [TT]-cp[tt] (solute) and [TT]-cs[tt] (solvent). ",
 820         "The box specified in the solute coordinate file ([TT]-cp[tt]) is used,",
 821         "unless [TT]-box[tt] is set.",
 822         "If you want the solute to be centered in the box,",
 823         "the program [gmx-editconf] has sophisticated options",
 824         "to change the box dimensions and center the solute.",
 825         "Solvent molecules are removed from the box where the ",
 826         "distance between any atom of the solute molecule(s) and any atom of ",
 827         "the solvent molecule is less than the sum of the scaled van der Waals",
 828         "radii of both atoms. A database ([TT]vdwradii.dat[tt]) of van der",
 829         "Waals radii is read by the program, and the resulting radii scaled",
 830         "by [TT]-scale[tt]. If radii are not found in the database, those"
 831         "atoms are assigned the (pre-scaled) distance [TT]-radius[tt].[PAR]",
 832
 833         "The default solvent is Simple Point Charge water (SPC), with coordinates ",
 834         "from [TT]$GMXLIB/spc216.gro[tt]. These coordinates can also be used",
 835         "for other 3-site water models, since a short equibilibration will remove",
 836         "the small differences between the models.",
 837         "Other solvents are also supported, as well as mixed solvents. The",
 838         "only restriction to solvent types is that a solvent molecule consists",
 839         "of exactly one residue. The residue information in the coordinate",
 840         "files is used, and should therefore be more or less consistent.",
 841         "In practice this means that two subsequent solvent molecules in the ",
 842         "solvent coordinate file should have different residue number.",
 843         "The box of solute is built by stacking the coordinates read from",
 844         "the coordinate file. This means that these coordinates should be ",
 845         "equlibrated in periodic boundary conditions to ensure a good",
 846         "alignment of molecules on the stacking interfaces.",
 847         "The [TT]-maxsol[tt] option simply adds only the first [TT]-maxsol[tt]",
 848         "solvent molecules and leaves out the rest that would have fitted",
 849         "into the box. This can create a void that can cause problems later.",
 850         "Choose your volume wisely.[PAR]",
 851
 852         "Setting [TT]-shell[tt] larger than zero will place a layer of water of",
 853         "the specified thickness (nm) around the solute. Hint: it is a good",
 854         "idea to put the protein in the center of a box first (using [gmx-editconf]).",
 855         "[PAR]",
 856
 857         "Finally, [THISMODULE] will optionally remove lines from your topology file in ",
 858         "which a number of solvent molecules is already added, and adds a ",
 859         "line with the total number of solvent molecules in your coordinate file."
 860     };
 861
 862     const char *bugs[] = {
 863         "Molecules must be whole in the initial configurations.",
 864     };
 865
 866     /* parameter data */
 867     gmx_bool       bProt, bBox;
 868     const char    *conf_prot, *confout;
 869     gmx_atomprop_t aps;
 870
 871     /* solute configuration data */
 872     t_topology *top;
 873     int         ePBC = -1;
 874     matrix      box;
 875
 876     t_filenm    fnm[] = {
 877         { efSTX, "-cp", "protein", ffOPTRD },
 878         { efSTX, "-cs", "spc216",  ffLIBRD},
 879         { efSTO, NULL,  NULL,      ffWRITE},
 880         { efTOP, NULL,  NULL,      ffOPTRW},
 881     };
 882 #define NFILE asize(fnm)
 883
 884     static real       defaultDistance = 0.105, r_shell = 0, scaleFactor = 0.57;
 885     static rvec       new_box         = {0.0, 0.0, 0.0};
 886     static gmx_bool   bReadV          = FALSE;
 887     static int        max_sol         = 0;
 888     gmx_output_env_t *oenv;
 889     t_pargs           pa[]              = {
 890         { "-box",    FALSE, etRVEC, {new_box},
 891           "Box size (in nm)" },
 892         { "-radius",   FALSE, etREAL, {&defaultDistance},
 893           "Default van der Waals distance"},
 894         { "-scale", FALSE, etREAL, {&scaleFactor},
 895           "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." },
 896         { "-shell",  FALSE, etREAL, {&r_shell},
 897           "Thickness of optional water layer around solute" },
 898         { "-maxsol", FALSE, etINT,  {&max_sol},
 899           "Maximum number of solvent molecules to add if they fit in the box. If zero (default) this is ignored" },
 900         { "-vel",    FALSE, etBOOL, {&bReadV},
 901           "Keep velocities from input solute and solvent" },
 902     };
 903
 904     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
 905                            asize(desc), desc, asize(bugs), bugs, &oenv))
 906     {
 907         return 0;
 908     }
 909
 910     const char *solventFileName = opt2fn("-cs", NFILE, fnm);
 911     bProt     = opt2bSet("-cp", NFILE, fnm);
 912     bBox      = opt2parg_bSet("-box", asize(pa), pa);
 913
 914     /* check input */
 915     if (!bProt && !bBox)
 916     {
 917         gmx_fatal(FARGS, "When no solute (-cp) is specified, "
 918                   "a box size (-box) must be specified");
 919     }
 920
 921     aps = gmx_atomprop_init();
 922
 923     std::vector<RVec> x;
 924     std::vector<RVec> v;
 925     snew(top, 1);
 926     if (bProt)
 927     {
 928         /* Generate a solute configuration */
 929         conf_prot = opt2fn("-cp", NFILE, fnm);
 930         readConformation(conf_prot, top, &x,
 931                          bReadV ? &v : NULL, &ePBC, box, "solute");
 932         if (bReadV && v.empty())
 933         {
 934             fprintf(stderr, "Note: no velocities found\n");
 935         }
 936         if (top->atoms.nr == 0)
 937         {
 938             fprintf(stderr, "Note: no atoms in %s\n", conf_prot);
 939             bProt = FALSE;
 940         }
 941     }
 942     if (bBox)
 943     {
 944         ePBC = epbcXYZ;
 945         clear_mat(box);
 946         box[XX][XX] = new_box[XX];
 947         box[YY][YY] = new_box[YY];
 948         box[ZZ][ZZ] = new_box[ZZ];
 949     }
 950     if (det(box) == 0)
 951     {
 952         gmx_fatal(FARGS, "Undefined solute box.\nCreate one with gmx editconf "
 953                   "or give explicit -box command line option");
 954     }
 955
 956     add_solv(solventFileName, top, &x, &v, ePBC, box,
 957              aps, defaultDistance, scaleFactor, r_shell, max_sol);
 958
 959     /* write new configuration 1 to file confout */
 960     confout = ftp2fn(efSTO, NFILE, fnm);
 961     fprintf(stderr, "Writing generated configuration to %s\n", confout);
 962     const char *outputTitle = (bProt ? *top->name : "Generated by gmx solvate");
 963     write_sto_conf(confout, outputTitle, &top->atoms, as_rvec_array(x.data()),
 964                    !v.empty() ? as_rvec_array(v.data()) : NULL, ePBC, box);
 965
 966     /* print size of generated configuration */
 967     fprintf(stderr, "\nOutput configuration contains %d atoms in %d residues\n",
 968             top->atoms.nr, top->atoms.nres);
 969     update_top(&top->atoms, box, NFILE, fnm, aps);
 970
 971     gmx_atomprop_destroy(aps);
 972     done_top(top);
 973     sfree(top);
 974     output_env_done(oenv);
 975     done_filenms(NFILE, fnm);
 976
 977     return 0;
 978 }