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, 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 "gromacs/commandline/pargs.h"
  44 #include "gromacs/fileio/confio.h"
  45 #include "gromacs/fileio/pdbio.h"
  46 #include "gromacs/gmxlib/conformation-utilities.h"
  47 #include "gromacs/gmxpreprocess/addconf.h"
  48 #include "gromacs/gmxpreprocess/read-conformation.h"
  49 #include "gromacs/legacyheaders/macros.h"
  50 #include "gromacs/legacyheaders/names.h"
  51 #include "gromacs/legacyheaders/typedefs.h"
  52 #include "gromacs/math/vec.h"
  53 #include "gromacs/topology/atomprop.h"
  54 #include "gromacs/utility/cstringutil.h"
  55 #include "gromacs/utility/fatalerror.h"
  56 #include "gromacs/utility/futil.h"
  57 #include "gromacs/utility/smalloc.h"
  58
  59 #ifdef DEBUG
  60 static void print_stat(rvec *x, int natoms, matrix box)
  61 {
  62     int  i, m;
  63     rvec xmin, xmax;
  64     for (m = 0; (m < DIM); m++)
  65     {
  66         xmin[m] = x[0][m];
  67         xmax[m] = x[0][m];
  68     }
  69     for (i = 0; (i < natoms); i++)
  70     {
  71         for (m = 0; m < DIM; m++)
  72         {
  73             xmin[m] = min(xmin[m], x[i][m]);
  74             xmax[m] = max(xmax[m], x[i][m]);
  75         }
  76     }
  77     for (m = 0; (m < DIM); m++)
  78     {
  79         fprintf(stderr, "DIM %d XMIN %8.3f XMAX %8.3f BOX %8.3f\n",
  80                 m, xmin[m], xmax[m], box[m][m]);
  81     }
  82 }
  83 #endif
  84
  85 typedef struct {
  86     char *name;
  87     int   natoms;
  88     int   nmol;
  89     int   i, i0;
  90     int   res0;
  91 } t_moltypes;
  92
  93 static void sort_molecule(t_atoms **atoms_solvt, rvec *x, rvec *v, real *r)
  94 {
  95     int         atnr, i, j, moltp = 0, nrmoltypes, resi_o, resi_n, resnr;
  96     t_moltypes *moltypes;
  97     t_atoms    *atoms, *newatoms;
  98     rvec       *newx, *newv = NULL;
  99     real       *newr;
 100
 101     fprintf(stderr, "Sorting configuration\n");
 102
 103     atoms = *atoms_solvt;
 104
 105     /* copy each residue from *atoms to a molecule in *molecule */
 106     moltypes   = NULL;
 107     nrmoltypes = 0;
 108     for (i = 0; i < atoms->nr; i++)
 109     {
 110         if ( (i == 0) || (atoms->atom[i].resind != atoms->atom[i-1].resind) )
 111         {
 112             /* see if this was a molecule type we haven't had yet: */
 113             moltp = NOTSET;
 114             for (j = 0; (j < nrmoltypes) && (moltp == NOTSET); j++)
 115             {
 116                 /* cppcheck-suppress nullPointer
 117                  * moltypes is guaranteed to be allocated because otherwise
 118                  * nrmoltypes is 0. */
 119                 if (strcmp(*(atoms->resinfo[atoms->atom[i].resind].name), moltypes[j].name) == 0)
 120                 {
 121                     moltp = j;
 122                 }
 123             }
 124             if (moltp == NOTSET)
 125             {
 126                 moltp = nrmoltypes;
 127                 nrmoltypes++;
 128                 srenew(moltypes, nrmoltypes);
 129                 moltypes[moltp].name = *(atoms->resinfo[atoms->atom[i].resind].name);
 130                 atnr                 = 0;
 131                 while ((i+atnr < atoms->nr) &&
 132                        (atoms->atom[i].resind == atoms->atom[i+atnr].resind))
 133                 {
 134                     atnr++;
 135                 }
 136                 moltypes[moltp].natoms = atnr;
 137                 moltypes[moltp].nmol   = 0;
 138             }
 139             moltypes[moltp].nmol++;
 140         }
 141     }
 142
 143     fprintf(stderr, "Found %d%s molecule type%s:\n",
 144             nrmoltypes, nrmoltypes == 1 ? "" : " different", nrmoltypes == 1 ? "" : "s");
 145     for (j = 0; j < nrmoltypes; j++)
 146     {
 147         if (j == 0)
 148         {
 149             moltypes[j].res0 = 0;
 150         }
 151         else
 152         {
 153             moltypes[j].res0 = moltypes[j-1].res0+moltypes[j-1].nmol;
 154         }
 155         fprintf(stderr, "%7s (%4d atoms): %5d residues\n",
 156                 moltypes[j].name, moltypes[j].natoms, moltypes[j].nmol);
 157     }
 158
 159     /* if we have only 1 moleculetype, we don't have to sort */
 160     if (nrmoltypes > 1)
 161     {
 162         /* find out which molecules should go where: */
 163         moltypes[0].i = moltypes[0].i0 = 0;
 164         for (j = 1; j < nrmoltypes; j++)
 165         {
 166             moltypes[j].i      =
 167                 moltypes[j].i0 =
 168                     moltypes[j-1].i0+moltypes[j-1].natoms*moltypes[j-1].nmol;
 169         }
 170
 171         /* now put them there: */
 172         snew(newatoms, 1);
 173         init_t_atoms(newatoms, atoms->nr, FALSE);
 174         newatoms->nres = atoms->nres;
 175         snew(newatoms->resinfo, atoms->nres);
 176         snew(newx, atoms->nr);
 177         if (v)
 178         {
 179             snew(newv, atoms->nr);
 180         }
 181         snew(newr, atoms->nr);
 182
 183         resi_n = 0;
 184         resnr  = 1;
 185         j      = 0;
 186         for (moltp = 0; moltp < nrmoltypes; moltp++)
 187         {
 188             i = 0;
 189             while (i < atoms->nr)
 190             {
 191                 resi_o = atoms->atom[i].resind;
 192                 if (strcmp(*atoms->resinfo[resi_o].name, moltypes[moltp].name) == 0)
 193                 {
 194                     /* Copy the residue info */
 195                     newatoms->resinfo[resi_n]    = atoms->resinfo[resi_o];
 196                     newatoms->resinfo[resi_n].nr = resnr;
 197                     /* Copy the atom info */
 198                     do
 199                     {
 200                         newatoms->atom[j]        = atoms->atom[i];
 201                         newatoms->atomname[j]    = atoms->atomname[i];
 202                         newatoms->atom[j].resind = resi_n;
 203                         copy_rvec(x[i], newx[j]);
 204                         if (v != NULL)
 205                         {
 206                             copy_rvec(v[i], newv[j]);
 207                         }
 208                         newr[j] = r[i];
 209                         i++;
 210                         j++;
 211                     }
 212                     while (i < atoms->nr && atoms->atom[i].resind == resi_o);
 213                     /* Increase the new residue counters */
 214                     resi_n++;
 215                     resnr++;
 216                 }
 217                 else
 218                 {
 219                     /* Skip this residue */
 220                     do
 221                     {
 222                         i++;
 223                     }
 224                     while (i < atoms->nr && atoms->atom[i].resind == resi_o);
 225                 }
 226             }
 227         }
 228
 229         /* put them back into the original arrays and throw away temporary arrays */
 230         sfree(atoms->atomname);
 231         sfree(atoms->resinfo);
 232         sfree(atoms->atom);
 233         sfree(atoms);
 234         *atoms_solvt = newatoms;
 235         for (i = 0; i < (*atoms_solvt)->nr; i++)
 236         {
 237             copy_rvec(newx[i], x[i]);
 238             if (v)
 239             {
 240                 copy_rvec(newv[i], v[i]);
 241             }
 242             r[i] = newr[i];
 243         }
 244         sfree(newx);
 245         if (v)
 246         {
 247             sfree(newv);
 248         }
 249         sfree(newr);
 250     }
 251     sfree(moltypes);
 252 }
 253
 254 static void rm_res_pbc(t_atoms *atoms, rvec *x, matrix box)
 255 {
 256     int  i, start, n, d, nat;
 257     rvec xcg;
 258
 259     start = 0;
 260     nat   = 0;
 261     clear_rvec(xcg);
 262     for (n = 0; n < atoms->nr; n++)
 263     {
 264         if (!is_hydrogen(*(atoms->atomname[n])))
 265         {
 266             nat++;
 267             rvec_inc(xcg, x[n]);
 268         }
 269         if ( (n+1 == atoms->nr) ||
 270              (atoms->atom[n+1].resind != atoms->atom[n].resind) )
 271         {
 272             /* if nat==0 we have only hydrogens in the solvent,
 273                we take last coordinate as cg */
 274             if (nat == 0)
 275             {
 276                 nat = 1;
 277                 copy_rvec(x[n], xcg);
 278             }
 279             svmul(1.0/nat, xcg, xcg);
 280             for (d = 0; d < DIM; d++)
 281             {
 282                 while (xcg[d] < 0)
 283                 {
 284                     for (i = start; i <= n; i++)
 285                     {
 286                         x[i][d] += box[d][d];
 287                     }
 288                     xcg[d] += box[d][d];
 289                 }
 290                 while (xcg[d] >= box[d][d])
 291                 {
 292                     for (i = start; i <= n; i++)
 293                     {
 294                         x[i][d] -= box[d][d];
 295                     }
 296                     xcg[d] -= box[d][d];
 297                 }
 298             }
 299             start = n+1;
 300             nat   = 0;
 301             clear_rvec(xcg);
 302         }
 303     }
 304 }
 305
 306 /* Make a new configuration by adding boxes*/
 307 static void make_new_conformation(t_atoms *atoms, rvec *x, rvec *v, real *r, matrix box, ivec n_box)
 308 {
 309     int     i, ix, iy, iz, m, j, imol, offset;
 310     rvec    delta;
 311     int     nmol;
 312
 313     nmol = n_box[XX]*n_box[YY]*n_box[ZZ];
 314
 315     /*print message*/
 316     fprintf(stderr, "Generating configuration\n");
 317     imol = 0;
 318     for (ix = 0; (ix < n_box[XX]); ix++)
 319     {
 320         delta[XX] = ix*box[XX][XX];
 321         for (iy = 0; (iy < n_box[YY]); iy++)
 322         {
 323             delta[YY] = iy*box[YY][YY];
 324             for (iz = 0; (iz < n_box[ZZ]); iz++)
 325             {
 326                 delta[ZZ] = iz*box[ZZ][ZZ];
 327                 offset    = imol*atoms->nr;
 328                 for (i = 0; (i < atoms->nr); i++)
 329                 {
 330                     for (m = 0; (m < DIM); m++)
 331                     {
 332                         x[offset+i][m] = delta[m]+x[i][m];
 333                     }
 334                     if (v)
 335                     {
 336                         for (m = 0; (m < DIM); m++)
 337                         {
 338                             v[offset+i][m] = v[i][m];
 339                         }
 340                     }
 341                     r[offset+i] = r[i];
 342                 }
 343                 imol++;
 344             }
 345         }
 346     }
 347     for (i = 1; (i < nmol); i++)
 348     {
 349         int offs    = i*atoms->nr;
 350         int offsres = i*atoms->nres;
 351         for (j = 0; (j < atoms->nr); j++)
 352         {
 353             atoms->atomname[offs+j]                    = atoms->atomname[j];
 354             atoms->atom[offs+j].resind                 = atoms->atom[j].resind + offsres;
 355             atoms->resinfo[atoms->atom[offs+j].resind] =
 356                 atoms->resinfo[atoms->atom[j].resind];
 357             atoms->resinfo[atoms->atom[offs+j].resind].nr += offsres;
 358         }
 359     }
 360     atoms->nr   *= nmol;
 361     atoms->nres *= nmol;
 362     for (i = 0; i < DIM; i++)
 363     {
 364         for (j = 0; j < DIM; j++)
 365         {
 366             box[j][i] *= n_box[j];
 367         }
 368     }
 369 }
 370
 371 static void add_solv(const char *fn, t_atoms *atoms, rvec **x, rvec **v, real **exclusionDistances,
 372                      int ePBC, matrix box,
 373                      gmx_atomprop_t aps,
 374                      real defaultDistance, real scaleFactor, int *atoms_added,
 375                      int *residues_added, real rshell, int max_sol,
 376                      const output_env_t oenv)
 377 {
 378     int      i, nmol;
 379     ivec     n_box;
 380     char     filename[STRLEN];
 381     char     title_solvt[STRLEN];
 382     t_atoms *atoms_solvt;
 383     rvec    *x_solvt, *v_solvt = NULL;
 384     real    *exclusionDistances_solvt;
 385     int      ePBC_solvt;
 386     matrix   box_solvt;
 387     int      onr, onres;
 388     char    *lfn;
 389
 390     lfn = gmxlibfn(fn);
 391     strncpy(filename, lfn, STRLEN);
 392     sfree(lfn);
 393     {
 394         int natoms;
 395         get_stx_coordnum(filename, &natoms);
 396         if (0 == natoms)
 397         {
 398             gmx_fatal(FARGS, "No solvent in %s, please check your input\n", filename);
 399         }
 400         snew(atoms_solvt, 1);
 401         init_t_atoms(atoms_solvt, natoms, FALSE);
 402     }
 403     snew(x_solvt, atoms_solvt->nr);
 404     if (v)
 405     {
 406         snew(v_solvt, atoms_solvt->nr);
 407     }
 408     snew(exclusionDistances_solvt, atoms_solvt->nr);
 409     snew(atoms_solvt->resinfo, atoms_solvt->nr);
 410     snew(atoms_solvt->atomname, atoms_solvt->nr);
 411     snew(atoms_solvt->atom, atoms_solvt->nr);
 412     atoms_solvt->pdbinfo = NULL;
 413     fprintf(stderr, "Reading solvent configuration%s\n",
 414             v_solvt ? " and velocities" : "");
 415     read_stx_conf(filename, title_solvt, atoms_solvt, x_solvt, v_solvt,
 416                   &ePBC_solvt, box_solvt);
 417     fprintf(stderr, "\"%s\"\n", title_solvt);
 418     fprintf(stderr, "solvent configuration contains %d atoms in %d residues\n",
 419             atoms_solvt->nr, atoms_solvt->nres);
 420     fprintf(stderr, "\n");
 421
 422     /* apply pbc for solvent configuration for whole molecules */
 423     rm_res_pbc(atoms_solvt, x_solvt, box_solvt);
 424
 425     /* initialise distance arrays for solvent configuration */
 426     exclusionDistances_solvt = makeExclusionDistances(atoms_solvt, aps, defaultDistance, scaleFactor);
 427
 428     /* calculate the box multiplication factors n_box[0...DIM] */
 429     nmol = 1;
 430     for (i = 0; (i < DIM); i++)
 431     {
 432         n_box[i] = 1;
 433         while (n_box[i]*box_solvt[i][i] < box[i][i])
 434         {
 435             n_box[i]++;
 436         }
 437         nmol *= n_box[i];
 438     }
 439     fprintf(stderr, "Will generate new solvent configuration of %dx%dx%d boxes\n",
 440             n_box[XX], n_box[YY], n_box[ZZ]);
 441
 442     /* realloc atoms_solvt for the new solvent configuration */
 443     srenew(atoms_solvt->resinfo, atoms_solvt->nres*nmol);
 444     srenew(atoms_solvt->atomname, atoms_solvt->nr*nmol);
 445     srenew(atoms_solvt->atom, atoms_solvt->nr*nmol);
 446     srenew(x_solvt, atoms_solvt->nr*nmol);
 447     if (v_solvt)
 448     {
 449         srenew(v_solvt, atoms_solvt->nr*nmol);
 450     }
 451     srenew(exclusionDistances_solvt, atoms_solvt->nr*nmol);
 452
 453     /* generate a new solvent configuration */
 454     make_new_conformation(atoms_solvt, x_solvt, v_solvt, exclusionDistances_solvt, box_solvt, n_box);
 455
 456 #ifdef DEBUG
 457     print_stat(x_solvt, atoms_solvt->nr, box_solvt);
 458 #endif
 459
 460 #ifdef DEBUG
 461     print_stat(x_solvt, atoms_solvt->nr, box_solvt);
 462 #endif
 463     /* Sort the solvent mixture, not the protein... */
 464     sort_molecule(&atoms_solvt, x_solvt, v_solvt, exclusionDistances_solvt);
 465
 466     /* add the two configurations */
 467     onr   = atoms->nr;
 468     onres = atoms->nres;
 469     add_conf(atoms, x, v, exclusionDistances, TRUE, ePBC, box, FALSE,
 470              atoms_solvt, x_solvt, v_solvt, exclusionDistances_solvt, TRUE, rshell, max_sol, oenv);
 471     *atoms_added    = atoms->nr-onr;
 472     *residues_added = atoms->nres-onres;
 473
 474     sfree(x_solvt);
 475     sfree(exclusionDistances_solvt);
 476     done_atom(atoms_solvt);
 477     sfree(atoms_solvt);
 478
 479     fprintf(stderr, "Generated solvent containing %d atoms in %d residues\n",
 480             *atoms_added, *residues_added);
 481 }
 482
 483 static void update_top(t_atoms *atoms, matrix box, int NFILE, t_filenm fnm[],
 484                        gmx_atomprop_t aps)
 485 {
 486 #define TEMP_FILENM "temp.top"
 487     FILE       *fpin, *fpout;
 488     char        buf[STRLEN], buf2[STRLEN], *temp;
 489     const char *topinout;
 490     int         line;
 491     gmx_bool    bSystem, bMolecules, bSkip;
 492     int         i, nsol = 0;
 493     double      mtot;
 494     real        vol, mm;
 495
 496     for (i = 0; (i < atoms->nres); i++)
 497     {
 498         /* calculate number of SOLvent molecules */
 499         if ( (strcmp(*atoms->resinfo[i].name, "SOL") == 0) ||
 500              (strcmp(*atoms->resinfo[i].name, "WAT") == 0) ||
 501              (strcmp(*atoms->resinfo[i].name, "HOH") == 0) )
 502         {
 503             nsol++;
 504         }
 505     }
 506     mtot = 0;
 507     for (i = 0; (i < atoms->nr); i++)
 508     {
 509         gmx_atomprop_query(aps, epropMass,
 510                            *atoms->resinfo[atoms->atom[i].resind].name,
 511                            *atoms->atomname[i], &mm);
 512         mtot += mm;
 513     }
 514
 515     vol = det(box);
 516
 517     fprintf(stderr, "Volume                 :  %10g (nm^3)\n", vol);
 518     fprintf(stderr, "Density                :  %10g (g/l)\n",
 519             (mtot*1e24)/(AVOGADRO*vol));
 520     fprintf(stderr, "Number of SOL molecules:  %5d   \n\n", nsol);
 521
 522     /* open topology file and append sol molecules */
 523     topinout  = ftp2fn(efTOP, NFILE, fnm);
 524     if (ftp2bSet(efTOP, NFILE, fnm) )
 525     {
 526         fprintf(stderr, "Processing topology\n");
 527         fpin    = gmx_ffopen(topinout, "r");
 528         fpout   = gmx_ffopen(TEMP_FILENM, "w");
 529         line    = 0;
 530         bSystem = bMolecules = FALSE;
 531         while (fgets(buf, STRLEN, fpin))
 532         {
 533             bSkip = FALSE;
 534             line++;
 535             strcpy(buf2, buf);
 536             if ((temp = strchr(buf2, '\n')) != NULL)
 537             {
 538                 temp[0] = '\0';
 539             }
 540             ltrim(buf2);
 541             if (buf2[0] == '[')
 542             {
 543                 buf2[0] = ' ';
 544                 if ((temp = strchr(buf2, '\n')) != NULL)
 545                 {
 546                     temp[0] = '\0';
 547                 }
 548                 rtrim(buf2);
 549                 if (buf2[strlen(buf2)-1] == ']')
 550                 {
 551                     buf2[strlen(buf2)-1] = '\0';
 552                     ltrim(buf2);
 553                     rtrim(buf2);
 554                     bSystem    = (gmx_strcasecmp(buf2, "system") == 0);
 555                     bMolecules = (gmx_strcasecmp(buf2, "molecules") == 0);
 556                 }
 557             }
 558             else if (bSystem && nsol && (buf[0] != ';') )
 559             {
 560                 /* if sol present, append "in water" to system name */
 561                 rtrim(buf2);
 562                 if (buf2[0] && (!strstr(buf2, " water")) )
 563                 {
 564                     sprintf(buf, "%s in water\n", buf2);
 565                     bSystem = FALSE;
 566                 }
 567             }
 568             else if (bMolecules)
 569             {
 570                 /* check if this is a line with solvent molecules */
 571                 sscanf(buf, "%4095s", buf2);
 572                 if (strcmp(buf2, "SOL") == 0)
 573                 {
 574                     sscanf(buf, "%*4095s %20d", &i);
 575                     nsol -= i;
 576                     if (nsol < 0)
 577                     {
 578                         bSkip = TRUE;
 579                         nsol += i;
 580                     }
 581                 }
 582             }
 583             if (bSkip)
 584             {
 585                 if ((temp = strchr(buf, '\n')) != NULL)
 586                 {
 587                     temp[0] = '\0';
 588                 }
 589                 fprintf(stdout, "Removing line #%d '%s' from topology file (%s)\n",
 590                         line, buf, topinout);
 591             }
 592             else
 593             {
 594                 fprintf(fpout, "%s", buf);
 595             }
 596         }
 597         gmx_ffclose(fpin);
 598         if (nsol)
 599         {
 600             fprintf(stdout, "Adding line for %d solvent molecules to "
 601                     "topology file (%s)\n", nsol, topinout);
 602             fprintf(fpout, "%-15s %5d\n", "SOL", nsol);
 603         }
 604         gmx_ffclose(fpout);
 605         /* use gmx_ffopen to generate backup of topinout */
 606         fpout = gmx_ffopen(topinout, "w");
 607         gmx_ffclose(fpout);
 608         rename(TEMP_FILENM, topinout);
 609     }
 610 #undef TEMP_FILENM
 611 }
 612
 613 int gmx_solvate(int argc, char *argv[])
 614 {
 615     const char *desc[] = {
 616         "[THISMODULE] can do one of 2 things:[PAR]",
 617
 618         "1) Generate a box of solvent. Specify [TT]-cs[tt] and [TT]-box[tt].",
 619         "Or specify [TT]-cs[tt] and [TT]-cp[tt] with a structure file with",
 620         "a box, but without atoms.[PAR]",
 621
 622         "2) Solvate a solute configuration, e.g. a protein, in a bath of solvent ",
 623         "molecules. Specify [TT]-cp[tt] (solute) and [TT]-cs[tt] (solvent). ",
 624         "The box specified in the solute coordinate file ([TT]-cp[tt]) is used,",
 625         "unless [TT]-box[tt] is set.",
 626         "If you want the solute to be centered in the box,",
 627         "the program [gmx-editconf] has sophisticated options",
 628         "to change the box dimensions and center the solute.",
 629         "Solvent molecules are removed from the box where the ",
 630         "distance between any atom of the solute molecule(s) and any atom of ",
 631         "the solvent molecule is less than the sum of the scaled van der Waals",
 632         "radii of both atoms. A database ([TT]vdwradii.dat[tt]) of van der",
 633         "Waals radii is read by the program, and the resulting radii scaled",
 634         "by [TT]-scale[tt]. If radii are not found in the database, those"
 635         "atoms are assigned the (pre-scaled) distance [TT]-radius[tt].[PAR]",
 636
 637         "The default solvent is Simple Point Charge water (SPC), with coordinates ",
 638         "from [TT]$GMXLIB/spc216.gro[tt]. These coordinates can also be used",
 639         "for other 3-site water models, since a short equibilibration will remove",
 640         "the small differences between the models.",
 641         "Other solvents are also supported, as well as mixed solvents. The",
 642         "only restriction to solvent types is that a solvent molecule consists",
 643         "of exactly one residue. The residue information in the coordinate",
 644         "files is used, and should therefore be more or less consistent.",
 645         "In practice this means that two subsequent solvent molecules in the ",
 646         "solvent coordinate file should have different residue number.",
 647         "The box of solute is built by stacking the coordinates read from",
 648         "the coordinate file. This means that these coordinates should be ",
 649         "equlibrated in periodic boundary conditions to ensure a good",
 650         "alignment of molecules on the stacking interfaces.",
 651         "The [TT]-maxsol[tt] option simply adds only the first [TT]-maxsol[tt]",
 652         "solvent molecules and leaves out the rest that would have fitted",
 653         "into the box. This can create a void that can cause problems later.",
 654         "Choose your volume wisely.[PAR]",
 655
 656         "The program can optionally rotate the solute molecule to align the",
 657         "longest molecule axis along a box edge. This way the amount of solvent",
 658         "molecules necessary is reduced.",
 659         "It should be kept in mind that this only works for",
 660         "short simulations, as e.g. an alpha-helical peptide in solution can ",
 661         "rotate over 90 degrees, within 500 ps. In general it is therefore ",
 662         "better to make a more or less cubic box.[PAR]",
 663
 664         "Setting [TT]-shell[tt] larger than zero will place a layer of water of",
 665         "the specified thickness (nm) around the solute. Hint: it is a good",
 666         "idea to put the protein in the center of a box first (using [gmx-editconf]).",
 667         "[PAR]",
 668
 669         "Finally, [THISMODULE] will optionally remove lines from your topology file in ",
 670         "which a number of solvent molecules is already added, and adds a ",
 671         "line with the total number of solvent molecules in your coordinate file."
 672     };
 673
 674     const char *bugs[] = {
 675         "Molecules must be whole in the initial configurations.",
 676     };
 677
 678     /* parameter data */
 679     gmx_bool       bProt, bBox;
 680     const char    *conf_prot, *confout;
 681     real          *exclusionDistances = NULL;
 682     gmx_atomprop_t aps;
 683
 684     /* protein configuration data */
 685     char    *title = NULL;
 686     t_atoms *atoms;
 687     rvec    *x    = NULL, *v = NULL;
 688     int      ePBC = -1;
 689     matrix   box;
 690
 691     /* other data types */
 692     int      atoms_added, residues_added;
 693
 694     t_filenm fnm[] = {
 695         { efSTX, "-cp", "protein", ffOPTRD },
 696         { efSTX, "-cs", "spc216",  ffLIBRD},
 697         { efSTO, NULL,  NULL,      ffWRITE},
 698         { efTOP, NULL,  NULL,      ffOPTRW},
 699     };
 700 #define NFILE asize(fnm)
 701
 702     static real     defaultDistance = 0.105, r_shell = 0, scaleFactor = 0.57;
 703     static rvec     new_box         = {0.0, 0.0, 0.0};
 704     static gmx_bool bReadV          = FALSE;
 705     static int      max_sol         = 0;
 706     output_env_t    oenv;
 707     t_pargs         pa[]              = {
 708         { "-box",    FALSE, etRVEC, {new_box},
 709           "Box size (in nm)" },
 710         { "-radius",   FALSE, etREAL, {&defaultDistance},
 711           "Default van der Waals distance"},
 712         { "-scale", FALSE, etREAL, {&scaleFactor},
 713           "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." },
 714         { "-shell",  FALSE, etREAL, {&r_shell},
 715           "Thickness of optional water layer around solute" },
 716         { "-maxsol", FALSE, etINT,  {&max_sol},
 717           "Maximum number of solvent molecules to add if they fit in the box. If zero (default) this is ignored" },
 718         { "-vel",    FALSE, etBOOL, {&bReadV},
 719           "Keep velocities from input solute and solvent" },
 720     };
 721
 722     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
 723                            asize(desc), desc, asize(bugs), bugs, &oenv))
 724     {
 725         return 0;
 726     }
 727
 728     const char *solventFileName = opt2fn("-cs", NFILE, fnm);
 729     bProt     = opt2bSet("-cp", NFILE, fnm);
 730     bBox      = opt2parg_bSet("-box", asize(pa), pa);
 731
 732     /* check input */
 733     if (!bProt && !bBox)
 734     {
 735         gmx_fatal(FARGS, "When no solute (-cp) is specified, "
 736                   "a box size (-box) must be specified");
 737     }
 738
 739     aps = gmx_atomprop_init();
 740
 741     snew(atoms, 1);
 742     init_t_atoms(atoms, 0, FALSE);
 743     if (bProt)
 744     {
 745         /* Generate a solute configuration */
 746         conf_prot = opt2fn("-cp", NFILE, fnm);
 747         title     = readConformation(conf_prot, atoms, &x,
 748                                      bReadV ? &v : NULL, &ePBC, box);
 749         exclusionDistances = makeExclusionDistances(atoms, aps, defaultDistance, scaleFactor);
 750
 751         if (bReadV && !v)
 752         {
 753             fprintf(stderr, "Note: no velocities found\n");
 754         }
 755         if (atoms->nr == 0)
 756         {
 757             fprintf(stderr, "Note: no atoms in %s\n", conf_prot);
 758             bProt = FALSE;
 759         }
 760     }
 761     if (bBox)
 762     {
 763         ePBC = epbcXYZ;
 764         clear_mat(box);
 765         box[XX][XX] = new_box[XX];
 766         box[YY][YY] = new_box[YY];
 767         box[ZZ][ZZ] = new_box[ZZ];
 768     }
 769     if (det(box) == 0)
 770     {
 771         gmx_fatal(FARGS, "Undefined solute box.\nCreate one with gmx editconf "
 772                   "or give explicit -box command line option");
 773     }
 774
 775     add_solv(solventFileName, atoms, &x, v ? &v : NULL, &exclusionDistances, ePBC, box,
 776              aps, defaultDistance, scaleFactor, &atoms_added, &residues_added, r_shell, max_sol,
 777              oenv);
 778
 779     /* write new configuration 1 to file confout */
 780     confout = ftp2fn(efSTO, NFILE, fnm);
 781     fprintf(stderr, "Writing generated configuration to %s\n", confout);
 782     if (bProt)
 783     {
 784         write_sto_conf(confout, title, atoms, x, v, ePBC, box);
 785         /* print box sizes and box type to stderr */
 786         fprintf(stderr, "%s\n", title);
 787     }
 788     else
 789     {
 790         write_sto_conf(confout, "Generated by gmx solvate", atoms, x, v, ePBC, box);
 791     }
 792
 793     /* print size of generated configuration */
 794     fprintf(stderr, "\nOutput configuration contains %d atoms in %d residues\n",
 795             atoms->nr, atoms->nres);
 796     update_top(atoms, box, NFILE, fnm, aps);
 797
 798     gmx_atomprop_destroy(aps);
 799     sfree(exclusionDistances);
 800     sfree(x);
 801     sfree(v);
 802     done_atom(atoms);
 803     sfree(atoms);
 804
 805     return 0;
 806 }