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