src/gromacs/gmxana/gmx_genbox.cpp

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