src/gromacs/gmxana/gmx_genbox.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 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41
  42 #include "sysstuff.h"
  43 #include "typedefs.h"
  44 #include "smalloc.h"
  45 #include "string2.h"
  46 #include "physics.h"
  47 #include "gromacs/fileio/confio.h"
  48 #include "txtdump.h"
  49 #include <math.h>
  50 #include "macros.h"
  51 #include "gromacs/random/random.h"
  52 #include "gromacs/fileio/futil.h"
  53 #include "atomprop.h"
  54 #include "names.h"
  55 #include "vec.h"
  56 #include "gmx_fatal.h"
  57 #include "gromacs/commandline/pargs.h"
  58 #include "vec.h"
  59 #include "gbutil.h"
  60 #include "addconf.h"
  61 #include "gromacs/fileio/pdbio.h"
  62 #include "pbc.h"
  63 #include "gmx_ana.h"
  64 #include "xvgr.h"
  65
  66 #ifdef DEBUG
  67 void print_stat(rvec *x, int natoms, matrix box)
  68 {
  69     int  i, m;
  70     rvec xmin, xmax;
  71     for (m = 0; (m < DIM); m++)
  72     {
  73         xmin[m] = x[0][m];
  74         xmax[m] = x[0][m];
  75     }
  76     for (i = 0; (i < natoms); i++)
  77     {
  78         for (m = 0; m < DIM; m++)
  79         {
  80             xmin[m] = min(xmin[m], x[i][m]);
  81             xmax[m] = max(xmax[m], x[i][m]);
  82         }
  83     }
  84     for (m = 0; (m < DIM); m++)
  85     {
  86         fprintf(stderr, "DIM %d XMIN %8.3f XMAX %8.3f BOX %8.3f\n",
  87                 m, xmin[m], xmax[m], box[m][m]);
  88     }
  89 }
  90 #endif
  91
  92 static gmx_bool in_box(t_pbc *pbc, rvec x)
  93 {
  94     rvec box_center, dx;
  95     int  shift;
  96
  97     /* pbc_dx_aiuc only works correctly with the rectangular box center */
  98     calc_box_center(ecenterRECT, pbc->box, box_center);
  99
 100     shift = pbc_dx_aiuc(pbc, x, box_center, dx);
 101
 102     return (shift == CENTRAL);
 103 }
 104
 105 static void mk_vdw(t_atoms *a, real rvdw[], gmx_atomprop_t aps,
 106                    real r_distance, real r_scale)
 107 {
 108     int i;
 109
 110     /* initialise van der waals arrays of configuration */
 111     fprintf(stderr, "Initialising van der waals distances...\n");
 112     for (i = 0; (i < a->nr); i++)
 113     {
 114         if (!gmx_atomprop_query(aps, epropVDW,
 115                                 *(a->resinfo[a->atom[i].resind].name),
 116                                 *(a->atomname[i]), &(rvdw[i])))
 117         {
 118             rvdw[i] = r_distance;
 119         }
 120         else
 121         {
 122             rvdw[i] *= r_scale;
 123         }
 124     }
 125 }
 126
 127 typedef struct {
 128     char *name;
 129     int   natoms;
 130     int   nmol;
 131     int   i, i0;
 132     int   res0;
 133 } t_moltypes;
 134
 135 void sort_molecule(t_atoms **atoms_solvt, rvec *x, rvec *v, real *r)
 136 {
 137     int         atnr, i, j, moltp = 0, nrmoltypes, resi_o, resi_n, resnr;
 138     t_moltypes *moltypes;
 139     t_atoms    *atoms, *newatoms;
 140     rvec       *newx, *newv = NULL;
 141     real       *newr;
 142
 143     fprintf(stderr, "Sorting configuration\n");
 144
 145     atoms = *atoms_solvt;
 146
 147     /* copy each residue from *atoms to a molecule in *molecule */
 148     moltypes   = NULL;
 149     nrmoltypes = 0;
 150     for (i = 0; i < atoms->nr; i++)
 151     {
 152         if ( (i == 0) || (atoms->atom[i].resind != atoms->atom[i-1].resind) )
 153         {
 154             /* see if this was a molecule type we haven't had yet: */
 155             moltp = NOTSET;
 156             for (j = 0; (j < nrmoltypes) && (moltp == NOTSET); j++)
 157             {
 158                 if (strcmp(*(atoms->resinfo[atoms->atom[i].resind].name),
 159                            moltypes[j].name) == 0)
 160                 {
 161                     moltp = j;
 162                 }
 163             }
 164             if (moltp == NOTSET)
 165             {
 166                 moltp = nrmoltypes;
 167                 nrmoltypes++;
 168                 srenew(moltypes, nrmoltypes);
 169                 moltypes[moltp].name = *(atoms->resinfo[atoms->atom[i].resind].name);
 170                 atnr                 = 0;
 171                 while ((i+atnr < atoms->nr) &&
 172                        (atoms->atom[i].resind == atoms->atom[i+atnr].resind))
 173                 {
 174                     atnr++;
 175                 }
 176                 moltypes[moltp].natoms = atnr;
 177                 moltypes[moltp].nmol   = 0;
 178             }
 179             moltypes[moltp].nmol++;
 180         }
 181     }
 182
 183     fprintf(stderr, "Found %d%s molecule type%s:\n",
 184             nrmoltypes, nrmoltypes == 1 ? "" : " different", nrmoltypes == 1 ? "" : "s");
 185     for (j = 0; j < nrmoltypes; j++)
 186     {
 187         if (j == 0)
 188         {
 189             moltypes[j].res0 = 0;
 190         }
 191         else
 192         {
 193             moltypes[j].res0 = moltypes[j-1].res0+moltypes[j-1].nmol;
 194         }
 195         fprintf(stderr, "%7s (%4d atoms): %5d residues\n",
 196                 moltypes[j].name, moltypes[j].natoms, moltypes[j].nmol);
 197     }
 198
 199     /* if we have only 1 moleculetype, we don't have to sort */
 200     if (nrmoltypes > 1)
 201     {
 202         /* find out which molecules should go where: */
 203         moltypes[0].i = moltypes[0].i0 = 0;
 204         for (j = 1; j < nrmoltypes; j++)
 205         {
 206             moltypes[j].i      =
 207                 moltypes[j].i0 =
 208                     moltypes[j-1].i0+moltypes[j-1].natoms*moltypes[j-1].nmol;
 209         }
 210
 211         /* now put them there: */
 212         snew(newatoms, 1);
 213         init_t_atoms(newatoms, atoms->nr, FALSE);
 214         newatoms->nres = atoms->nres;
 215         snew(newatoms->resinfo, atoms->nres);
 216         snew(newx, atoms->nr);
 217         if (v)
 218         {
 219             snew(newv, atoms->nr);
 220         }
 221         snew(newr, atoms->nr);
 222
 223         resi_n = 0;
 224         resnr  = 1;
 225         j      = 0;
 226         for (moltp = 0; moltp < nrmoltypes; moltp++)
 227         {
 228             i = 0;
 229             while (i < atoms->nr)
 230             {
 231                 resi_o = atoms->atom[i].resind;
 232                 if (strcmp(*atoms->resinfo[resi_o].name, moltypes[moltp].name) == 0)
 233                 {
 234                     /* Copy the residue info */
 235                     newatoms->resinfo[resi_n]    = atoms->resinfo[resi_o];
 236                     newatoms->resinfo[resi_n].nr = resnr;
 237                     /* Copy the atom info */
 238                     do
 239                     {
 240                         newatoms->atom[j]        = atoms->atom[i];
 241                         newatoms->atomname[j]    = atoms->atomname[i];
 242                         newatoms->atom[j].resind = resi_n;
 243                         copy_rvec(x[i], newx[j]);
 244                         if (v != NULL)
 245                         {
 246                             copy_rvec(v[i], newv[j]);
 247                         }
 248                         newr[j] = r[i];
 249                         i++;
 250                         j++;
 251                     }
 252                     while (i < atoms->nr && atoms->atom[i].resind == resi_o);
 253                     /* Increase the new residue counters */
 254                     resi_n++;
 255                     resnr++;
 256                 }
 257                 else
 258                 {
 259                     /* Skip this residue */
 260                     do
 261                     {
 262                         i++;
 263                     }
 264                     while (i < atoms->nr && atoms->atom[i].resind == resi_o);
 265                 }
 266             }
 267         }
 268
 269         /* put them back into the original arrays and throw away temporary arrays */
 270         sfree(atoms->atomname);
 271         sfree(atoms->resinfo);
 272         sfree(atoms->atom);
 273         sfree(atoms);
 274         *atoms_solvt = newatoms;
 275         for (i = 0; i < (*atoms_solvt)->nr; i++)
 276         {
 277             copy_rvec(newx[i], x[i]);
 278             if (v)
 279             {
 280                 copy_rvec(newv[i], v[i]);
 281             }
 282             r[i] = newr[i];
 283         }
 284         sfree(newx);
 285         if (v)
 286         {
 287             sfree(newv);
 288         }
 289         sfree(newr);
 290     }
 291     sfree(moltypes);
 292 }
 293
 294 void rm_res_pbc(t_atoms *atoms, rvec *x, matrix box)
 295 {
 296     int  i, start, n, d, nat;
 297     rvec xcg;
 298
 299     start = 0;
 300     nat   = 0;
 301     clear_rvec(xcg);
 302     for (n = 0; n < atoms->nr; n++)
 303     {
 304         if (!is_hydrogen(*(atoms->atomname[n])))
 305         {
 306             nat++;
 307             rvec_inc(xcg, x[n]);
 308         }
 309         if ( (n+1 == atoms->nr) ||
 310              (atoms->atom[n+1].resind != atoms->atom[n].resind) )
 311         {
 312             /* if nat==0 we have only hydrogens in the solvent,
 313                we take last coordinate as cg */
 314             if (nat == 0)
 315             {
 316                 nat = 1;
 317                 copy_rvec(x[n], xcg);
 318             }
 319             svmul(1.0/nat, xcg, xcg);
 320             for (d = 0; d < DIM; d++)
 321             {
 322                 while (xcg[d] < 0)
 323                 {
 324                     for (i = start; i <= n; i++)
 325                     {
 326                         x[i][d] += box[d][d];
 327                     }
 328                     xcg[d] += box[d][d];
 329                 }
 330                 while (xcg[d] >= box[d][d])
 331                 {
 332                     for (i = start; i <= n; i++)
 333                     {
 334                         x[i][d] -= box[d][d];
 335                     }
 336                     xcg[d] -= box[d][d];
 337                 }
 338             }
 339             start = n+1;
 340             nat   = 0;
 341             clear_rvec(xcg);
 342         }
 343     }
 344 }
 345
 346 /* This is a (maybe) slow workaround to avoid the neighbor searching in addconf.c, which
 347  * leaks memory (May 2012). The function could be deleted as soon as the momory leaks
 348  * in addconf.c are fixed.
 349  * However, when inserting a small molecule in a system containing not too many atoms,
 350  * allPairsDistOk is probably even faster than addconf.c
 351  */
 352 static gmx_bool
 353 allPairsDistOk(t_atoms *atoms, rvec *x, real *r,
 354                int ePBC, matrix box,
 355                t_atoms *atoms_insrt, rvec *x_n, real *r_insrt)
 356 {
 357     int   i, j;
 358     rvec  dx;
 359     real  n2, r2;
 360     t_pbc pbc;
 361
 362     set_pbc(&pbc, ePBC, box);
 363     for (i = 0; i < atoms->nr; i++)
 364     {
 365         for (j = 0; j < atoms_insrt->nr; j++)
 366         {
 367             pbc_dx(&pbc, x[i], x_n[j], dx);
 368             n2 = norm2(dx);
 369             r2 = sqr(r[i]+r_insrt[j]);
 370             if (n2 < r2)
 371             {
 372                 return FALSE;
 373             }
 374         }
 375     }
 376     return TRUE;
 377 }
 378
 379 /* enum for random rotations of inserted solutes */
 380 enum {
 381     en_rot, en_rotXYZ, en_rotZ, en_rotNone, en_NR
 382 };
 383
 384 static char *insert_mols(const char *mol_insrt, int nmol_insrt, int ntry, int seed,
 385                          t_atoms *atoms, rvec **x, real **r, int ePBC, matrix box,
 386                          gmx_atomprop_t aps,
 387                          real r_distance, real r_scale, real rshell,
 388                          const output_env_t oenv,
 389                          const char* posfn, const rvec deltaR, int enum_rot,
 390                          gmx_bool bCheckAllPairDist)
 391 {
 392     t_pbc            pbc;
 393     static  char    *title_insrt;
 394     t_atoms          atoms_insrt;
 395     rvec            *x_insrt, *x_n;
 396     real            *r_insrt;
 397     int              ePBC_insrt;
 398     matrix           box_insrt;
 399     int              i, mol, onr, ncol;
 400     real             alfa = 0., beta = 0., gamma = 0.;
 401     rvec             offset_x;
 402     int              trial;
 403     double         **rpos;
 404     gmx_rng_t        rng;
 405
 406     rng = gmx_rng_init(seed);
 407     set_pbc(&pbc, ePBC, box);
 408
 409     /* read number of atoms of insert molecules */
 410     get_stx_coordnum(mol_insrt, &atoms_insrt.nr);
 411     if (atoms_insrt.nr == 0)
 412     {
 413         gmx_fatal(FARGS, "No molecule in %s, please check your input\n", mol_insrt);
 414     }
 415     /* allocate memory for atom coordinates of insert molecules */
 416     snew(x_insrt, atoms_insrt.nr);
 417     snew(r_insrt, atoms_insrt.nr);
 418     snew(atoms_insrt.resinfo, atoms_insrt.nr);
 419     snew(atoms_insrt.atomname, atoms_insrt.nr);
 420     snew(atoms_insrt.atom, atoms_insrt.nr);
 421     atoms_insrt.pdbinfo = NULL;
 422     snew(x_n, atoms_insrt.nr);
 423     snew(title_insrt, STRLEN);
 424
 425     /* read residue number, residue names, atomnames, coordinates etc. */
 426     fprintf(stderr, "Reading molecule configuration \n");
 427     read_stx_conf(mol_insrt, title_insrt, &atoms_insrt, x_insrt, NULL,
 428                   &ePBC_insrt, box_insrt);
 429     fprintf(stderr, "%s\nContaining %d atoms in %d residue\n",
 430             title_insrt, atoms_insrt.nr, atoms_insrt.nres);
 431     srenew(atoms_insrt.resinfo, atoms_insrt.nres);
 432
 433     /* initialise van der waals arrays of insert molecules */
 434     mk_vdw(&atoms_insrt, r_insrt, aps, r_distance, r_scale);
 435
 436     /* With -ip, take nmol_insrt from file posfn */
 437     if (posfn != NULL)
 438     {
 439         nmol_insrt = read_xvg(posfn, &rpos, &ncol);
 440         if (ncol != 3)
 441         {
 442             gmx_fatal(FARGS, "Expected 3 columns (x/y/z coordinates) in file %s\n", ncol, posfn);
 443         }
 444         fprintf(stderr, "Read %d positions from file %s\n\n", nmol_insrt, posfn);
 445     }
 446
 447     srenew(atoms->resinfo, (atoms->nres+nmol_insrt*atoms_insrt.nres));
 448     srenew(atoms->atomname, (atoms->nr+atoms_insrt.nr*nmol_insrt));
 449     srenew(atoms->atom, (atoms->nr+atoms_insrt.nr*nmol_insrt));
 450     srenew(*x, (atoms->nr+atoms_insrt.nr*nmol_insrt));
 451     srenew(*r, (atoms->nr+atoms_insrt.nr*nmol_insrt));
 452
 453     trial = mol = 0;
 454     while ((mol < nmol_insrt) && (trial < ntry*nmol_insrt))
 455     {
 456         fprintf(stderr, "\rTry %d", trial++);
 457         for (i = 0; (i < atoms_insrt.nr); i++)
 458         {
 459             copy_rvec(x_insrt[i], x_n[i]);
 460         }
 461         switch (enum_rot)
 462         {
 463             case en_rotXYZ:
 464                 alfa  = 2*M_PI * gmx_rng_uniform_real(rng);
 465                 beta  = 2*M_PI * gmx_rng_uniform_real(rng);
 466                 gamma = 2*M_PI * gmx_rng_uniform_real(rng);
 467                 break;
 468             case en_rotZ:
 469                 alfa  = beta = 0.;
 470                 gamma = 2*M_PI * gmx_rng_uniform_real(rng);
 471                 break;
 472             case en_rotNone:
 473                 alfa = beta = gamma = 0.;
 474                 break;
 475         }
 476         if (enum_rot == en_rotXYZ || (enum_rot == en_rotZ))
 477         {
 478             rotate_conf(atoms_insrt.nr, x_n, NULL, alfa, beta, gamma);
 479         }
 480         if (posfn == NULL)
 481         {
 482             /* insert at random positions */
 483             offset_x[XX] = box[XX][XX] * gmx_rng_uniform_real(rng);
 484             offset_x[YY] = box[YY][YY] * gmx_rng_uniform_real(rng);
 485             offset_x[ZZ] = box[ZZ][ZZ] * gmx_rng_uniform_real(rng);
 486             gen_box(0, atoms_insrt.nr, x_n, box_insrt, offset_x, TRUE);
 487             if (!in_box(&pbc, x_n[0]) || !in_box(&pbc, x_n[atoms_insrt.nr-1]))
 488             {
 489                 continue;
 490             }
 491         }
 492         else
 493         {
 494             /* Insert at positions taken from option -ip file */
 495             offset_x[XX] = rpos[XX][mol] + deltaR[XX]*(2 * gmx_rng_uniform_real(rng)-1);
 496             offset_x[YY] = rpos[YY][mol] + deltaR[YY]*(2 * gmx_rng_uniform_real(rng)-1);
 497             offset_x[ZZ] = rpos[ZZ][mol] + deltaR[ZZ]*(2 * gmx_rng_uniform_real(rng)-1);
 498             for (i = 0; i < atoms_insrt.nr; i++)
 499             {
 500                 rvec_inc(x_n[i], offset_x);
 501             }
 502         }
 503
 504         onr = atoms->nr;
 505
 506         /* This is a (maybe) slow workaround to avoid too many calls of add_conf, which
 507          * leaks memory (status May 2012). If the momory leaks in add_conf() are fixed,
 508          * this check could be removed. Note, however, that allPairsDistOk is probably
 509          * even faster than add_conf() when inserting a small molecule into a moderately
 510          * small system.
 511          */
 512         if (bCheckAllPairDist && !allPairsDistOk(atoms, *x, *r, ePBC, box, &atoms_insrt, x_n, r_insrt))
 513         {
 514             continue;
 515         }
 516
 517         add_conf(atoms, x, NULL, r, FALSE, ePBC, box, TRUE,
 518                  &atoms_insrt, x_n, NULL, r_insrt, FALSE, rshell, 0, oenv);
 519
 520         if (atoms->nr == (atoms_insrt.nr+onr))
 521         {
 522             mol++;
 523             fprintf(stderr, " success (now %d atoms)!\n", atoms->nr);
 524         }
 525     }
 526     gmx_rng_destroy(rng);
 527     srenew(atoms->resinfo,  atoms->nres);
 528     srenew(atoms->atomname, atoms->nr);
 529     srenew(atoms->atom,     atoms->nr);
 530     srenew(*x,              atoms->nr);
 531     srenew(*r,              atoms->nr);
 532
 533     fprintf(stderr, "\n");
 534     /* print number of molecules added */
 535     fprintf(stderr, "Added %d molecules (out of %d requested) of %s\n",
 536             mol, nmol_insrt, *atoms_insrt.resinfo[0].name);
 537
 538     return title_insrt;
 539 }
 540
 541 static void add_solv(const char *fn, t_atoms *atoms, rvec **x, rvec **v, real **r,
 542                      int ePBC, matrix box,
 543                      gmx_atomprop_t aps,
 544                      real r_distance, real r_scale, int *atoms_added,
 545                      int *residues_added, real rshell, int max_sol,
 546                      const output_env_t oenv)
 547 {
 548     int      i, nmol;
 549     ivec     n_box;
 550     char     filename[STRLEN];
 551     char     title_solvt[STRLEN];
 552     t_atoms *atoms_solvt;
 553     rvec    *x_solvt, *v_solvt = NULL;
 554     real    *r_solvt;
 555     int      ePBC_solvt;
 556     matrix   box_solvt;
 557     int      onr, onres;
 558     char    *lfn;
 559
 560     lfn = gmxlibfn(fn);
 561     strncpy(filename, lfn, STRLEN);
 562     sfree(lfn);
 563     snew(atoms_solvt, 1);
 564     get_stx_coordnum(filename, &(atoms_solvt->nr));
 565     if (atoms_solvt->nr == 0)
 566     {
 567         gmx_fatal(FARGS, "No solvent in %s, please check your input\n", filename);
 568     }
 569     snew(x_solvt, atoms_solvt->nr);
 570     if (v)
 571     {
 572         snew(v_solvt, atoms_solvt->nr);
 573     }
 574     snew(r_solvt, atoms_solvt->nr);
 575     snew(atoms_solvt->resinfo, atoms_solvt->nr);
 576     snew(atoms_solvt->atomname, atoms_solvt->nr);
 577     snew(atoms_solvt->atom, atoms_solvt->nr);
 578     atoms_solvt->pdbinfo = NULL;
 579     fprintf(stderr, "Reading solvent configuration%s\n",
 580             v_solvt ? " and velocities" : "");
 581     read_stx_conf(filename, title_solvt, atoms_solvt, x_solvt, v_solvt,
 582                   &ePBC_solvt, box_solvt);
 583     fprintf(stderr, "\"%s\"\n", title_solvt);
 584     fprintf(stderr, "solvent configuration contains %d atoms in %d residues\n",
 585             atoms_solvt->nr, atoms_solvt->nres);
 586     fprintf(stderr, "\n");
 587
 588     /* apply pbc for solvent configuration for whole molecules */
 589     rm_res_pbc(atoms_solvt, x_solvt, box_solvt);
 590
 591     /* initialise van der waals arrays of solvent configuration */
 592     mk_vdw(atoms_solvt, r_solvt, aps, r_distance, r_scale);
 593
 594     /* calculate the box multiplication factors n_box[0...DIM] */
 595     nmol = 1;
 596     for (i = 0; (i < DIM); i++)
 597     {
 598         n_box[i] = 1;
 599         while (n_box[i]*box_solvt[i][i] < box[i][i])
 600         {
 601             n_box[i]++;
 602         }
 603         nmol *= n_box[i];
 604     }
 605     fprintf(stderr, "Will generate new solvent configuration of %dx%dx%d boxes\n",
 606             n_box[XX], n_box[YY], n_box[ZZ]);
 607
 608     /* realloc atoms_solvt for the new solvent configuration */
 609     srenew(atoms_solvt->resinfo, atoms_solvt->nres*nmol);
 610     srenew(atoms_solvt->atomname, atoms_solvt->nr*nmol);
 611     srenew(atoms_solvt->atom, atoms_solvt->nr*nmol);
 612     srenew(x_solvt, atoms_solvt->nr*nmol);
 613     if (v_solvt)
 614     {
 615         srenew(v_solvt, atoms_solvt->nr*nmol);
 616     }
 617     srenew(r_solvt, atoms_solvt->nr*nmol);
 618
 619     /* generate a new solvent configuration */
 620     genconf(atoms_solvt, x_solvt, v_solvt, r_solvt, box_solvt, n_box);
 621
 622 #ifdef DEBUG
 623     print_stat(x_solvt, atoms_solvt->nr, box_solvt);
 624 #endif
 625
 626 #ifdef DEBUG
 627     print_stat(x_solvt, atoms_solvt->nr, box_solvt);
 628 #endif
 629     /* Sort the solvent mixture, not the protein... */
 630     sort_molecule(&atoms_solvt, x_solvt, v_solvt, r_solvt);
 631
 632     /* add the two configurations */
 633     onr   = atoms->nr;
 634     onres = atoms->nres;
 635     add_conf(atoms, x, v, r, TRUE, ePBC, box, FALSE,
 636              atoms_solvt, x_solvt, v_solvt, r_solvt, TRUE, rshell, max_sol, oenv);
 637     *atoms_added    = atoms->nr-onr;
 638     *residues_added = atoms->nres-onres;
 639
 640     sfree(x_solvt);
 641     sfree(r_solvt);
 642
 643     fprintf(stderr, "Generated solvent containing %d atoms in %d residues\n",
 644             *atoms_added, *residues_added);
 645 }
 646
 647 static char *read_prot(const char *confin, t_atoms *atoms, rvec **x, rvec **v,
 648                        real **r, int *ePBC, matrix box, gmx_atomprop_t aps,
 649                        real r_distance, real r_scale)
 650 {
 651     char *title;
 652     int   natoms;
 653
 654     snew(title, STRLEN);
 655     get_stx_coordnum(confin, &natoms);
 656
 657     /* allocate memory for atom coordinates of configuration 1 */
 658     snew(*x, natoms);
 659     if (v)
 660     {
 661         snew(*v, natoms);
 662     }
 663     snew(*r, natoms);
 664     init_t_atoms(atoms, natoms, FALSE);
 665
 666     /* read residue number, residue names, atomnames, coordinates etc. */
 667     fprintf(stderr, "Reading solute configuration%s\n", v ? " and velocities" : "");
 668     read_stx_conf(confin, title, atoms, *x, v ? *v : NULL, ePBC, box);
 669     fprintf(stderr, "%s\nContaining %d atoms in %d residues\n",
 670             title, atoms->nr, atoms->nres);
 671
 672     /* initialise van der waals arrays of configuration 1 */
 673     mk_vdw(atoms, *r, aps, r_distance, r_scale);
 674
 675     return title;
 676 }
 677
 678 static void update_top(t_atoms *atoms, matrix box, int NFILE, t_filenm fnm[],
 679                        gmx_atomprop_t aps)
 680 {
 681 #define TEMP_FILENM "temp.top"
 682     FILE       *fpin, *fpout;
 683     char        buf[STRLEN], buf2[STRLEN], *temp;
 684     const char *topinout;
 685     int         line;
 686     gmx_bool    bSystem, bMolecules, bSkip;
 687     int         i, nsol = 0;
 688     double      mtot;
 689     real        vol, mm;
 690
 691     for (i = 0; (i < atoms->nres); i++)
 692     {
 693         /* calculate number of SOLvent molecules */
 694         if ( (strcmp(*atoms->resinfo[i].name, "SOL") == 0) ||
 695              (strcmp(*atoms->resinfo[i].name, "WAT") == 0) ||
 696              (strcmp(*atoms->resinfo[i].name, "HOH") == 0) )
 697         {
 698             nsol++;
 699         }
 700     }
 701     mtot = 0;
 702     for (i = 0; (i < atoms->nr); i++)
 703     {
 704         gmx_atomprop_query(aps, epropMass,
 705                            *atoms->resinfo[atoms->atom[i].resind].name,
 706                            *atoms->atomname[i], &mm);
 707         mtot += mm;
 708     }
 709
 710     vol = det(box);
 711
 712     fprintf(stderr, "Volume                 :  %10g (nm^3)\n", vol);
 713     fprintf(stderr, "Density                :  %10g (g/l)\n",
 714             (mtot*1e24)/(AVOGADRO*vol));
 715     fprintf(stderr, "Number of SOL molecules:  %5d   \n\n", nsol);
 716
 717     /* open topology file and append sol molecules */
 718     topinout  = ftp2fn(efTOP, NFILE, fnm);
 719     if (ftp2bSet(efTOP, NFILE, fnm) )
 720     {
 721         fprintf(stderr, "Processing topology\n");
 722         fpin    = ffopen(topinout, "r");
 723         fpout   = ffopen(TEMP_FILENM, "w");
 724         line    = 0;
 725         bSystem = bMolecules = FALSE;
 726         while (fgets(buf, STRLEN, fpin))
 727         {
 728             bSkip = FALSE;
 729             line++;
 730             strcpy(buf2, buf);
 731             if ((temp = strchr(buf2, '\n')) != NULL)
 732             {
 733                 temp[0] = '\0';
 734             }
 735             ltrim(buf2);
 736             if (buf2[0] == '[')
 737             {
 738                 buf2[0] = ' ';
 739                 if ((temp = strchr(buf2, '\n')) != NULL)
 740                 {
 741                     temp[0] = '\0';
 742                 }
 743                 rtrim(buf2);
 744                 if (buf2[strlen(buf2)-1] == ']')
 745                 {
 746                     buf2[strlen(buf2)-1] = '\0';
 747                     ltrim(buf2);
 748                     rtrim(buf2);
 749                     bSystem    = (gmx_strcasecmp(buf2, "system") == 0);
 750                     bMolecules = (gmx_strcasecmp(buf2, "molecules") == 0);
 751                 }
 752             }
 753             else if (bSystem && nsol && (buf[0] != ';') )
 754             {
 755                 /* if sol present, append "in water" to system name */
 756                 rtrim(buf2);
 757                 if (buf2[0] && (!strstr(buf2, " water")) )
 758                 {
 759                     sprintf(buf, "%s in water\n", buf2);
 760                     bSystem = FALSE;
 761                 }
 762             }
 763             else if (bMolecules)
 764             {
 765                 /* check if this is a line with solvent molecules */
 766                 sscanf(buf, "%4095s", buf2);
 767                 if (strcmp(buf2, "SOL") == 0)
 768                 {
 769                     sscanf(buf, "%*4095s %20d", &i);
 770                     nsol -= i;
 771                     if (nsol < 0)
 772                     {
 773                         bSkip = TRUE;
 774                         nsol += i;
 775                     }
 776                 }
 777             }
 778             if (bSkip)
 779             {
 780                 if ((temp = strchr(buf, '\n')) != NULL)
 781                 {
 782                     temp[0] = '\0';
 783                 }
 784                 fprintf(stdout, "Removing line #%d '%s' from topology file (%s)\n",
 785                         line, buf, topinout);
 786             }
 787             else
 788             {
 789                 fprintf(fpout, "%s", buf);
 790             }
 791         }
 792         ffclose(fpin);
 793         if (nsol)
 794         {
 795             fprintf(stdout, "Adding line for %d solvent molecules to "
 796                     "topology file (%s)\n", nsol, topinout);
 797             fprintf(fpout, "%-15s %5d\n", "SOL", nsol);
 798         }
 799         ffclose(fpout);
 800         /* use ffopen to generate backup of topinout */
 801         fpout = ffopen(topinout, "w");
 802         ffclose(fpout);
 803         rename(TEMP_FILENM, topinout);
 804     }
 805 #undef TEMP_FILENM
 806 }
 807
 808 int gmx_genbox(int argc, char *argv[])
 809 {
 810     const char *desc[] = {
 811         "[THISMODULE] can do one of 4 things:[PAR]",
 812
 813         "1) Generate a box of solvent. Specify [TT]-cs[tt] and [TT]-box[tt]. Or specify [TT]-cs[tt] and",
 814         "[TT]-cp[tt] with a structure file with a box, but without atoms.[PAR]",
 815
 816         "2) Solvate a solute configuration, e.g. a protein, in a bath of solvent ",
 817         "molecules. Specify [TT]-cp[tt] (solute) and [TT]-cs[tt] (solvent). ",
 818         "The box specified in the solute coordinate file ([TT]-cp[tt]) is used,",
 819         "unless [TT]-box[tt] is set.",
 820         "If you want the solute to be centered in the box,",
 821         "the program [TT]editconf[tt] has sophisticated options",
 822         "to change the box dimensions and center the solute.",
 823         "Solvent molecules are removed from the box where the ",
 824         "distance between any atom of the solute molecule(s) and any atom of ",
 825         "the solvent molecule is less than the sum of the van der Waals radii of ",
 826         "both atoms. A database ([TT]vdwradii.dat[tt]) of van der Waals radii is ",
 827         "read by the program, and atoms not in the database are ",
 828         "assigned a default distance [TT]-vdwd[tt].",
 829         "Note that this option will also influence the distances between",
 830         "solvent molecules if they contain atoms that are not in the database.",
 831         "[PAR]",
 832
 833         "3) Insert a number ([TT]-nmol[tt]) of extra molecules ([TT]-ci[tt]) ",
 834         "at random positions.",
 835         "The program iterates until [TT]nmol[tt] molecules",
 836         "have been inserted in the box. To test whether an insertion is ",
 837         "successful the same van der Waals criterium is used as for removal of ",
 838         "solvent molecules. When no appropriately-sized ",
 839         "holes (holes that can hold an extra molecule) are available, the ",
 840         "program tries for [TT]-nmol[tt] * [TT]-try[tt] times before giving up. ",
 841         "Increase [TT]-try[tt] if you have several small holes to fill.",
 842         "Option [TT]-rot[tt] defines if the molecules are randomly oriented.",
 843         "[PAR]",
 844
 845         "4) Insert a number of molecules ([TT]-ci[tt]) at positions defined in",
 846         "positions.dat ([TT]-ip[tt]). positions.dat should have 3 columns (x/y/z),",
 847         "that give the displacements compared to the input molecule position ([TT]-ci[tt]).",
 848         "Hence, if positions.dat should contain the absolut positions, the molecule ",
 849         "must be centered to 0/0/0 before using genbox (use, e.g., editconf -center).",
 850         "Comments in positions.dat starting with # are ignored. Option [TT]-dr[tt]",
 851         "defines the maximally allowed displacements during insertial trials.",
 852         "[TT]-try[tt] and [TT]-rot[tt] work as in mode (3) (see above)",
 853         "[PAR]",
 854
 855         "If you need to do more than one of the above operations, it can be",
 856         "best to call [THISMODULE] separately for each operation, so that",
 857         "you are sure of the order in which the operations occur.[PAR]",
 858
 859         "The default solvent is Simple Point Charge water (SPC), with coordinates ",
 860         "from [TT]$GMXLIB/spc216.gro[tt]. These coordinates can also be used",
 861         "for other 3-site water models, since a short equibilibration will remove",
 862         "the small differences between the models.",
 863         "Other solvents are also supported, as well as mixed solvents. The",
 864         "only restriction to solvent types is that a solvent molecule consists",
 865         "of exactly one residue. The residue information in the coordinate",
 866         "files is used, and should therefore be more or less consistent.",
 867         "In practice this means that two subsequent solvent molecules in the ",
 868         "solvent coordinate file should have different residue number.",
 869         "The box of solute is built by stacking the coordinates read from",
 870         "the coordinate file. This means that these coordinates should be ",
 871         "equlibrated in periodic boundary conditions to ensure a good",
 872         "alignment of molecules on the stacking interfaces.",
 873         "The [TT]-maxsol[tt] option simply adds only the first [TT]-maxsol[tt]",
 874         "solvent molecules and leaves out the rest that would have fitted",
 875         "into the box. This can create a void that can cause problems later.",
 876         "Choose your volume wisely.[PAR]",
 877
 878         "The program can optionally rotate the solute molecule to align the",
 879         "longest molecule axis along a box edge. This way the amount of solvent",
 880         "molecules necessary is reduced.",
 881         "It should be kept in mind that this only works for",
 882         "short simulations, as e.g. an alpha-helical peptide in solution can ",
 883         "rotate over 90 degrees, within 500 ps. In general it is therefore ",
 884         "better to make a more or less cubic box.[PAR]",
 885
 886         "Setting [TT]-shell[tt] larger than zero will place a layer of water of",
 887         "the specified thickness (nm) around the solute. Hint: it is a good",
 888         "idea to put the protein in the center of a box first (using [gmx-editconf]).",
 889         "[PAR]",
 890
 891         "Finally, [THISMODULE] will optionally remove lines from your topology file in ",
 892         "which a number of solvent molecules is already added, and adds a ",
 893         "line with the total number of solvent molecules in your coordinate file."
 894     };
 895
 896     const char *bugs[] = {
 897         "Molecules must be whole in the initial configurations.",
 898         "Many repeated neighbor searchings with -ci blows up the allocated memory. "
 899         "Option -allpair avoids this using all-to-all distance checks (slow for large systems)"
 900     };
 901
 902     /* parameter data */
 903     gmx_bool       bSol, bProt, bBox;
 904     const char    *conf_prot, *confout;
 905     int            bInsert;
 906     real          *r;
 907     char          *title_ins;
 908     gmx_atomprop_t aps;
 909
 910     /* protein configuration data */
 911     char    *title = NULL;
 912     t_atoms  atoms;
 913     rvec    *x, *v = NULL;
 914     int      ePBC = -1;
 915     matrix   box;
 916
 917     /* other data types */
 918     int      atoms_added, residues_added;
 919
 920     t_filenm fnm[] = {
 921         { efSTX, "-cp", "protein", ffOPTRD },
 922         { efSTX, "-cs", "spc216",  ffLIBOPTRD},
 923         { efSTX, "-ci", "insert",  ffOPTRD},
 924         { efDAT, "-ip", "positions",  ffOPTRD},
 925         { efSTO, NULL,  NULL,      ffWRITE},
 926         { efTOP, NULL,  NULL,      ffOPTRW},
 927     };
 928 #define NFILE asize(fnm)
 929
 930     static int      nmol_ins   = 0, nmol_try = 10, seed = 1997, enum_rot;
 931     static real     r_distance = 0.105, r_shell = 0, r_scale = 0.57;
 932     static rvec     new_box    = {0.0, 0.0, 0.0}, deltaR = {0.0, 0.0, 0.0};
 933     static gmx_bool bReadV     = FALSE, bCheckAllPairDist = FALSE;
 934     static int      max_sol    = 0;
 935     output_env_t    oenv;
 936     const char     *enum_rot_string[] = {NULL, "xyz", "z", "none", NULL};
 937     t_pargs         pa[]              = {
 938         { "-box",    FALSE, etRVEC, {new_box},
 939           "Box size" },
 940         { "-nmol",   FALSE, etINT, {&nmol_ins},
 941           "Number of extra molecules to insert" },
 942         { "-try",    FALSE, etINT, {&nmol_try},
 943           "Try inserting [TT]-nmol[tt] times [TT]-try[tt] times" },
 944         { "-seed",   FALSE, etINT, {&seed},
 945           "Random generator seed"},
 946         { "-vdwd",   FALSE, etREAL, {&r_distance},
 947           "Default van der Waals distance"},
 948         { "-vdwscale", FALSE, etREAL, {&r_scale},
 949           "HIDDENScale 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." },
 950         { "-shell",  FALSE, etREAL, {&r_shell},
 951           "Thickness of optional water layer around solute" },
 952         { "-maxsol", FALSE, etINT,  {&max_sol},
 953           "Maximum number of solvent molecules to add if they fit in the box. If zero (default) this is ignored" },
 954         { "-vel",    FALSE, etBOOL, {&bReadV},
 955           "Keep velocities from input solute and solvent" },
 956         { "-dr",    FALSE, etRVEC, {deltaR},
 957           "Allowed displacement in x/y/z from positions in [TT]-ip[tt] file" },
 958         { "-rot", FALSE,  etENUM, {enum_rot_string},
 959           "rotate inserted molecules randomly" },
 960         { "-allpair",    FALSE, etBOOL, {&bCheckAllPairDist},
 961           "Avoid momory leaks during neighbor searching with option -ci. May be slow for large systems." },
 962     };
 963
 964     if (!parse_common_args(&argc, argv, PCA_BE_NICE, NFILE, fnm, asize(pa), pa,
 965                            asize(desc), desc, asize(bugs), bugs, &oenv))
 966     {
 967         return 0;
 968     }
 969
 970     bInsert   = opt2bSet("-ci", NFILE, fnm) && ((nmol_ins > 0) || opt2bSet("-ip", NFILE, fnm));
 971     bSol      = opt2bSet("-cs", NFILE, fnm);
 972     bProt     = opt2bSet("-cp", NFILE, fnm);
 973     bBox      = opt2parg_bSet("-box", asize(pa), pa);
 974     enum_rot  = nenum(enum_rot_string);
 975
 976     /* check input */
 977     if (bInsert && (nmol_ins <= 0 && !opt2bSet("-ip", NFILE, fnm)))
 978     {
 979         gmx_fatal(FARGS, "When specifying inserted molecules (-ci), "
 980                   "-nmol must be larger than 0 or positions must be given with -ip");
 981     }
 982     if (!bProt && !bBox)
 983     {
 984         gmx_fatal(FARGS, "When no solute (-cp) is specified, "
 985                   "a box size (-box) must be specified");
 986     }
 987
 988     aps = gmx_atomprop_init();
 989
 990     if (bProt)
 991     {
 992         /*generate a solute configuration */
 993         conf_prot = opt2fn("-cp", NFILE, fnm);
 994         title     = read_prot(conf_prot, &atoms, &x, bReadV ? &v : NULL, &r, &ePBC, box,
 995                               aps, r_distance, r_scale);
 996         if (bReadV && !v)
 997         {
 998             fprintf(stderr, "Note: no velocities found\n");
 999         }
1000         if (atoms.nr == 0)
1001         {
1002             fprintf(stderr, "Note: no atoms in %s\n", conf_prot);
1003             bProt = FALSE;
1004         }
1005     }
1006     if (!bProt)
1007     {
1008         atoms.nr       = 0;
1009         atoms.nres     = 0;
1010         atoms.resinfo  = NULL;
1011         atoms.atomname = NULL;
1012         atoms.atom     = NULL;
1013         atoms.pdbinfo  = NULL;
1014         x              = NULL;
1015         r              = NULL;
1016     }
1017     if (bBox)
1018     {
1019         ePBC = epbcXYZ;
1020         clear_mat(box);
1021         box[XX][XX] = new_box[XX];
1022         box[YY][YY] = new_box[YY];
1023         box[ZZ][ZZ] = new_box[ZZ];
1024     }
1025     if (det(box) == 0)
1026     {
1027         gmx_fatal(FARGS, "Undefined solute box.\nCreate one with editconf "
1028                   "or give explicit -box command line option");
1029     }
1030
1031     /* add nmol_ins molecules of atoms_ins
1032        in random orientation at random place */
1033     if (bInsert)
1034     {
1035         title_ins = insert_mols(opt2fn("-ci", NFILE, fnm), nmol_ins, nmol_try, seed,
1036                                 &atoms, &x, &r, ePBC, box, aps,
1037                                 r_distance, r_scale, r_shell,
1038                                 oenv, opt2fn_null("-ip", NFILE, fnm), deltaR, enum_rot,
1039                                 bCheckAllPairDist);
1040     }
1041     else
1042     {
1043         title_ins = strdup("Generated by genbox");
1044     }
1045
1046     /* add solvent */
1047     if (bSol)
1048     {
1049         add_solv(opt2fn("-cs", NFILE, fnm), &atoms, &x, v ? &v : NULL, &r, ePBC, box,
1050                  aps, r_distance, r_scale, &atoms_added, &residues_added, r_shell, max_sol,
1051                  oenv);
1052     }
1053
1054     /* write new configuration 1 to file confout */
1055     confout = ftp2fn(efSTO, NFILE, fnm);
1056     fprintf(stderr, "Writing generated configuration to %s\n", confout);
1057     if (bProt)
1058     {
1059         write_sto_conf(confout, title, &atoms, x, v, ePBC, box);
1060         /* print box sizes and box type to stderr */
1061         fprintf(stderr, "%s\n", title);
1062     }
1063     else
1064     {
1065         write_sto_conf(confout, title_ins, &atoms, x, v, ePBC, box);
1066     }
1067
1068     sfree(title_ins);
1069
1070     /* print size of generated configuration */
1071     fprintf(stderr, "\nOutput configuration contains %d atoms in %d residues\n",
1072             atoms.nr, atoms.nres);
1073     update_top(&atoms, box, NFILE, fnm, aps);
1074
1075     gmx_atomprop_destroy(aps);
1076
1077     return 0;
1078 }