src/tools/gmx_genbox.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39
  40 #include "sysstuff.h"
  41 #include "typedefs.h"
  42 #include "smalloc.h"
  43 #include "string2.h"
  44 #include "physics.h"
  45 #include "confio.h"
  46 #include "copyrite.h"
  47 #include "txtdump.h"
  48 #include "math.h"
  49 #include "macros.h"
  50 #include "random.h"
  51 #include "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 "pdbio.h"
  61 #include "pbc.h"
  62 #include "gmx_ana.h"
  63
  64 #ifdef DEBUG
  65 void print_stat(rvec *x,int natoms,matrix box)
  66 {
  67   int i,m;
  68   rvec xmin,xmax;
  69   for(m=0;(m<DIM);m++) {
  70     xmin[m]=x[0][m];
  71     xmax[m]=x[0][m];
  72   }
  73   for(i=0;(i<natoms);i++) {
  74     for (m=0;m<DIM;m++) {
  75       xmin[m]=min(xmin[m],x[i][m]);
  76       xmax[m]=max(xmax[m],x[i][m]);
  77     }
  78   }
  79   for(m=0;(m<DIM);m++)
  80     fprintf(stderr,"DIM %d XMIN %8.3f XMAX %8.3f BOX %8.3f\n",
  81             m,xmin[m],xmax[m],box[m][m]);
  82 }
  83 #endif
  84
  85 static gmx_bool in_box(t_pbc *pbc,rvec x)
  86 {
  87   rvec box_center,dx;
  88   int  shift;
  89
  90   calc_box_center(ecenterTRIC,pbc->box,box_center);
  91
  92   shift = pbc_dx_aiuc(pbc,x,box_center,dx);
  93
  94   return (shift == CENTRAL);
  95 }
  96
  97 static void mk_vdw(t_atoms *a,real rvdw[],gmx_atomprop_t aps,
  98                    real r_distance)
  99 {
 100   int i;
 101
 102   /* initialise van der waals arrays of configuration */
 103   fprintf(stderr,"Initialising van der waals distances...\n");
 104   for(i=0; (i < a->nr); i++)
 105     if (!gmx_atomprop_query(aps,epropVDW,
 106                             *(a->resinfo[a->atom[i].resind].name),
 107                             *(a->atomname[i]),&(rvdw[i])))
 108       rvdw[i] = r_distance;
 109 }
 110
 111 typedef struct {
 112   char *name;
 113   int  natoms;
 114   int  nmol;
 115   int  i,i0;
 116   int  res0;
 117 } t_moltypes;
 118
 119 void sort_molecule(t_atoms **atoms_solvt,rvec *x,rvec *v,real *r)
 120 {
 121   int atnr,i,j,moltp=0,nrmoltypes,resi_o,resi_n,resnr;
 122   t_moltypes *moltypes;
 123   t_atoms *atoms,*newatoms;
 124   rvec *newx, *newv=NULL;
 125   real *newr;
 126
 127   fprintf(stderr,"Sorting configuration\n");
 128
 129   atoms = *atoms_solvt;
 130
 131   /* copy each residue from *atoms to a molecule in *molecule */
 132   moltypes=NULL;
 133   nrmoltypes=0;
 134   atnr=0;
 135   for (i=0; i<atoms->nr; i++) {
 136     if ( (i==0) || (atoms->atom[i].resind != atoms->atom[i-1].resind) ) {
 137       /* see if this was a molecule type we haven't had yet: */
 138       moltp=NOTSET;
 139       for (j=0; (j<nrmoltypes) && (moltp==NOTSET); j++)
 140         if (strcmp(*(atoms->resinfo[atoms->atom[i].resind].name),
 141                    moltypes[j].name)==0)
 142           moltp=j;
 143       if (moltp==NOTSET) {
 144         moltp=nrmoltypes;
 145         nrmoltypes++;
 146         srenew(moltypes,nrmoltypes);
 147         moltypes[moltp].name=*(atoms->resinfo[atoms->atom[i].resind].name);
 148         atnr = 0;
 149         while ((i+atnr<atoms->nr) &&
 150                (atoms->atom[i].resind == atoms->atom[i+atnr].resind))
 151           atnr++;
 152         moltypes[moltp].natoms=atnr;
 153         moltypes[moltp].nmol=0;
 154       }
 155       moltypes[moltp].nmol++;
 156     }
 157   }
 158
 159   fprintf(stderr,"Found %d%s molecule type%s:\n",
 160           nrmoltypes,nrmoltypes==1?"":" different",nrmoltypes==1?"":"s");
 161   for(j=0; j<nrmoltypes; j++) {
 162     if (j==0)
 163       moltypes[j].res0 = 0;
 164     else
 165       moltypes[j].res0 = moltypes[j-1].res0+moltypes[j-1].nmol;
 166     fprintf(stderr,"%7s (%4d atoms): %5d residues\n",
 167             moltypes[j].name,moltypes[j].natoms,moltypes[j].nmol);
 168   }
 169
 170   /* if we have only 1 moleculetype, we don't have to sort */
 171   if (nrmoltypes>1) {
 172     /* find out which molecules should go where: */
 173     moltypes[0].i = moltypes[0].i0 = 0;
 174     for(j=1; j<nrmoltypes; j++) {
 175       moltypes[j].i =
 176         moltypes[j].i0 =
 177         moltypes[j-1].i0+moltypes[j-1].natoms*moltypes[j-1].nmol;
 178     }
 179
 180     /* now put them there: */
 181     snew(newatoms,1);
 182     init_t_atoms(newatoms,atoms->nr,FALSE);
 183     newatoms->nres=atoms->nres;
 184     snew(newatoms->resinfo,atoms->nres);
 185     snew(newx,atoms->nr);
 186     if (v) snew(newv,atoms->nr);
 187     snew(newr,atoms->nr);
 188
 189     resi_n = 0;
 190     resnr = 1;
 191     j = 0;
 192     for(moltp=0; moltp<nrmoltypes; moltp++) {
 193       i = 0;
 194       while (i < atoms->nr) {
 195         resi_o = atoms->atom[i].resind;
 196         if (strcmp(*atoms->resinfo[resi_o].name,moltypes[moltp].name) == 0) {
 197           /* Copy the residue info */
 198           newatoms->resinfo[resi_n]    = atoms->resinfo[resi_o];
 199           newatoms->resinfo[resi_n].nr = resnr;
 200           /* Copy the atom info */
 201           do {
 202             newatoms->atom[j]        = atoms->atom[i];
 203             newatoms->atomname[j]    = atoms->atomname[i];
 204             newatoms->atom[j].resind = resi_n;
 205             copy_rvec(x[i],newx[j]);
 206             if (v != NULL) {
 207               copy_rvec(v[i],newv[j]);
 208             }
 209             newr[j] = r[i];
 210             i++;
 211             j++;
 212           } while (i < atoms->nr && atoms->atom[i].resind == resi_o);
 213           /* Increase the new residue counters */
 214           resi_n++;
 215           resnr++;
 216         } else {
 217           /* Skip this residue */
 218           do {
 219             i++;
 220           } while (i < atoms->nr && atoms->atom[i].resind == resi_o);
 221         }
 222       }
 223     }
 224
 225     /* put them back into the original arrays and throw away temporary arrays */
 226     sfree(atoms->atomname);
 227     sfree(atoms->resinfo);
 228     sfree(atoms->atom);
 229     sfree(atoms);
 230     *atoms_solvt = newatoms;
 231     for (i=0; i<(*atoms_solvt)->nr; i++) {
 232       copy_rvec(newx[i],x[i]);
 233       if (v) copy_rvec(newv[i],v[i]);
 234       r[i]=newr[i];
 235     }
 236     sfree(newx);
 237     if (v) sfree(newv);
 238     sfree(newr);
 239   }
 240   sfree(moltypes);
 241 }
 242
 243 void rm_res_pbc(t_atoms *atoms, rvec *x, matrix box)
 244 {
 245   int i,start,n,d,nat;
 246   rvec xcg;
 247
 248   start=0;
 249   nat=0;
 250   clear_rvec(xcg);
 251   for(n=0; n<atoms->nr; n++) {
 252     if (!is_hydrogen(*(atoms->atomname[n]))) {
 253       nat++;
 254       rvec_inc(xcg,x[n]);
 255     }
 256     if ( (n+1 == atoms->nr) ||
 257          (atoms->atom[n+1].resind != atoms->atom[n].resind) ) {
 258       /* if nat==0 we have only hydrogens in the solvent,
 259          we take last coordinate as cg */
 260       if (nat==0) {
 261         nat=1;
 262         copy_rvec(x[n],xcg);
 263       }
 264       svmul(1.0/nat,xcg,xcg);
 265       for(d=0; d<DIM; d++) {
 266         while (xcg[d]<0) {
 267           for(i=start; i<=n; i++)
 268             x[i][d]+=box[d][d];
 269           xcg[d]+=box[d][d];
 270         }
 271         while (xcg[d]>=box[d][d]) {
 272           for(i=start; i<=n; i++)
 273             x[i][d]-=box[d][d];
 274           xcg[d]-=box[d][d];
 275         }
 276       }
 277       start=n+1;
 278       nat=0;
 279       clear_rvec(xcg);
 280     }
 281   }
 282 }
 283
 284 static char *insert_mols(const char *mol_insrt,int nmol_insrt,int ntry,int seed,
 285                          t_atoms *atoms,rvec **x,real **r,int ePBC,matrix box,
 286                          gmx_atomprop_t aps,real r_distance,real rshell,
 287                          const output_env_t oenv)
 288 {
 289   t_pbc   pbc;
 290   static  char    *title_insrt;
 291   t_atoms atoms_insrt;
 292   rvec    *x_insrt,*x_n;
 293   real    *r_insrt;
 294   int     ePBC_insrt;
 295   matrix  box_insrt;
 296   int     i,mol,onr;
 297   real    alfa,beta,gamma;
 298   rvec    offset_x;
 299   int     trial;
 300
 301   set_pbc(&pbc,ePBC,box);
 302
 303   /* read number of atoms of insert molecules */
 304   get_stx_coordnum(mol_insrt,&atoms_insrt.nr);
 305   if (atoms_insrt.nr == 0)
 306     gmx_fatal(FARGS,"No molecule in %s, please check your input\n",mol_insrt);
 307   /* allocate memory for atom coordinates of insert molecules */
 308   snew(x_insrt,atoms_insrt.nr);
 309   snew(r_insrt,atoms_insrt.nr);
 310   snew(atoms_insrt.resinfo,atoms_insrt.nr);
 311   snew(atoms_insrt.atomname,atoms_insrt.nr);
 312   snew(atoms_insrt.atom,atoms_insrt.nr);
 313   atoms_insrt.pdbinfo = NULL;
 314   snew(x_n,atoms_insrt.nr);
 315   snew(title_insrt,STRLEN);
 316
 317   /* read residue number, residue names, atomnames, coordinates etc. */
 318   fprintf(stderr,"Reading molecule configuration \n");
 319   read_stx_conf(mol_insrt,title_insrt,&atoms_insrt,x_insrt,NULL,
 320                 &ePBC_insrt,box_insrt);
 321   fprintf(stderr,"%s\nContaining %d atoms in %d residue\n",
 322           title_insrt,atoms_insrt.nr,atoms_insrt.nres);
 323   srenew(atoms_insrt.resinfo,atoms_insrt.nres);
 324
 325   /* initialise van der waals arrays of insert molecules */
 326   mk_vdw(&atoms_insrt,r_insrt,aps,r_distance);
 327
 328   srenew(atoms->resinfo,(atoms->nres+nmol_insrt*atoms_insrt.nres));
 329   srenew(atoms->atomname,(atoms->nr+atoms_insrt.nr*nmol_insrt));
 330   srenew(atoms->atom,(atoms->nr+atoms_insrt.nr*nmol_insrt));
 331   srenew(*x,(atoms->nr+atoms_insrt.nr*nmol_insrt));
 332   srenew(*r,(atoms->nr+atoms_insrt.nr*nmol_insrt));
 333
 334   trial=mol=0;
 335   while ((mol < nmol_insrt) && (trial < ntry*nmol_insrt)) {
 336     fprintf(stderr,"\rTry %d",trial++);
 337     for (i=0;(i<atoms_insrt.nr);i++) {
 338       copy_rvec(x_insrt[i],x_n[i]);
 339     }
 340     alfa=2*M_PI*rando(&seed);
 341     beta=2*M_PI*rando(&seed);
 342     gamma=2*M_PI*rando(&seed);
 343     rotate_conf(atoms_insrt.nr,x_n,NULL,alfa,beta,gamma);
 344     offset_x[XX]=box[XX][XX]*rando(&seed);
 345     offset_x[YY]=box[YY][YY]*rando(&seed);
 346     offset_x[ZZ]=box[ZZ][ZZ]*rando(&seed);
 347     gen_box(0,atoms_insrt.nr,x_n,box_insrt,offset_x,TRUE);
 348     if (!in_box(&pbc,x_n[0]) || !in_box(&pbc,x_n[atoms_insrt.nr-1]))
 349       continue;
 350     onr=atoms->nr;
 351
 352     add_conf(atoms,x,NULL,r,FALSE,ePBC,box,TRUE,
 353              &atoms_insrt,x_n,NULL,r_insrt,FALSE,rshell,0,oenv);
 354
 355     if (atoms->nr==(atoms_insrt.nr+onr)) {
 356       mol++;
 357       fprintf(stderr," success (now %d atoms)!",atoms->nr);
 358     }
 359   }
 360   srenew(atoms->resinfo,  atoms->nres);
 361   srenew(atoms->atomname, atoms->nr);
 362   srenew(atoms->atom,     atoms->nr);
 363   srenew(*x,              atoms->nr);
 364   srenew(*r,              atoms->nr);
 365
 366   fprintf(stderr,"\n");
 367   /* print number of molecules added */
 368   fprintf(stderr,"Added %d molecules (out of %d requested) of %s\n",
 369           mol,nmol_insrt,*atoms_insrt.resinfo[0].name);
 370
 371   return title_insrt;
 372 }
 373
 374 static void add_solv(const char *fn,t_atoms *atoms,rvec **x,rvec **v,real **r,
 375                      int ePBC,matrix box,
 376                      gmx_atomprop_t aps,real r_distance,int *atoms_added,
 377                      int *residues_added,real rshell,int max_sol,
 378                      const output_env_t oenv)
 379 {
 380   int     i,nmol;
 381   ivec    n_box;
 382   char    filename[STRLEN];
 383   char    title_solvt[STRLEN];
 384   t_atoms *atoms_solvt;
 385   rvec    *x_solvt,*v_solvt=NULL;
 386   real    *r_solvt;
 387   int     ePBC_solvt;
 388   matrix  box_solvt;
 389   int     onr,onres;
 390   char    *lfn;
 391
 392   lfn = gmxlibfn(fn);
 393   strncpy(filename,lfn,STRLEN);
 394   sfree(lfn);
 395   snew(atoms_solvt,1);
 396   get_stx_coordnum(filename,&(atoms_solvt->nr));
 397   if (atoms_solvt->nr == 0)
 398     gmx_fatal(FARGS,"No solvent in %s, please check your input\n",filename);
 399   snew(x_solvt,atoms_solvt->nr);
 400   if (v) snew(v_solvt,atoms_solvt->nr);
 401   snew(r_solvt,atoms_solvt->nr);
 402   snew(atoms_solvt->resinfo,atoms_solvt->nr);
 403   snew(atoms_solvt->atomname,atoms_solvt->nr);
 404   snew(atoms_solvt->atom,atoms_solvt->nr);
 405   atoms_solvt->pdbinfo = NULL;
 406   fprintf(stderr,"Reading solvent configuration%s\n",
 407           v_solvt?" and velocities":"");
 408   read_stx_conf(filename,title_solvt,atoms_solvt,x_solvt,v_solvt,
 409                 &ePBC_solvt,box_solvt);
 410   fprintf(stderr,"\"%s\"\n",title_solvt);
 411   fprintf(stderr,"solvent configuration contains %d atoms in %d residues\n",
 412           atoms_solvt->nr,atoms_solvt->nres);
 413   fprintf(stderr,"\n");
 414
 415   /* apply pbc for solvent configuration for whole molecules */
 416   rm_res_pbc(atoms_solvt,x_solvt,box_solvt);
 417
 418   /* initialise van der waals arrays of solvent configuration */
 419   mk_vdw(atoms_solvt,r_solvt,aps,r_distance);
 420
 421   /* calculate the box multiplication factors n_box[0...DIM] */
 422   nmol=1;
 423   for (i=0; (i < DIM);i++) {
 424     n_box[i] = 1;
 425     while (n_box[i]*box_solvt[i][i] < box[i][i])
 426       n_box[i]++;
 427     nmol*=n_box[i];
 428   }
 429   fprintf(stderr,"Will generate new solvent configuration of %dx%dx%d boxes\n",
 430           n_box[XX],n_box[YY],n_box[ZZ]);
 431
 432   /* realloc atoms_solvt for the new solvent configuration */
 433   srenew(atoms_solvt->resinfo,atoms_solvt->nres*nmol);
 434   srenew(atoms_solvt->atomname,atoms_solvt->nr*nmol);
 435   srenew(atoms_solvt->atom,atoms_solvt->nr*nmol);
 436   srenew(x_solvt,atoms_solvt->nr*nmol);
 437   if (v_solvt) srenew(v_solvt,atoms_solvt->nr*nmol);
 438   srenew(r_solvt,atoms_solvt->nr*nmol);
 439
 440   /* generate a new solvent configuration */
 441   genconf(atoms_solvt,x_solvt,v_solvt,r_solvt,box_solvt,n_box);
 442
 443 #ifdef DEBUG
 444   print_stat(x_solvt,atoms_solvt->nr,box_solvt);
 445 #endif
 446
 447 #ifdef DEBUG
 448   print_stat(x_solvt,atoms_solvt->nr,box_solvt);
 449 #endif
 450   /* Sort the solvent mixture, not the protein... */
 451   sort_molecule(&atoms_solvt,x_solvt,v_solvt,r_solvt);
 452
 453   /* add the two configurations */
 454   onr=atoms->nr;
 455   onres=atoms->nres;
 456   add_conf(atoms,x,v,r,TRUE,ePBC,box,FALSE,
 457            atoms_solvt,x_solvt,v_solvt,r_solvt,TRUE,rshell,max_sol,oenv);
 458   *atoms_added=atoms->nr-onr;
 459   *residues_added=atoms->nres-onres;
 460
 461   sfree(x_solvt);
 462   sfree(r_solvt);
 463
 464   fprintf(stderr,"Generated solvent containing %d atoms in %d residues\n",
 465           *atoms_added,*residues_added);
 466 }
 467
 468 static char *read_prot(const char *confin,t_atoms *atoms,rvec **x,rvec **v,
 469                        real **r, int *ePBC,matrix box,gmx_atomprop_t aps,
 470                        real r_distance)
 471 {
 472   char *title;
 473   int  natoms;
 474
 475   snew(title,STRLEN);
 476   get_stx_coordnum(confin,&natoms);
 477
 478   /* allocate memory for atom coordinates of configuration 1 */
 479   snew(*x,natoms);
 480   if (v) snew(*v,natoms);
 481   snew(*r,natoms);
 482   init_t_atoms(atoms,natoms,FALSE);
 483
 484   /* read residue number, residue names, atomnames, coordinates etc. */
 485   fprintf(stderr,"Reading solute configuration%s\n",v?" and velocities":"");
 486   read_stx_conf(confin,title,atoms,*x,v?*v:NULL,ePBC,box);
 487   fprintf(stderr,"%s\nContaining %d atoms in %d residues\n",
 488           title,atoms->nr,atoms->nres);
 489
 490   /* initialise van der waals arrays of configuration 1 */
 491   mk_vdw(atoms,*r,aps,r_distance);
 492
 493   return title;
 494 }
 495
 496 static void update_top(t_atoms *atoms,matrix box,int NFILE,t_filenm fnm[],
 497                        gmx_atomprop_t aps)
 498 {
 499 #define TEMP_FILENM "temp.top"
 500   FILE   *fpin,*fpout;
 501   char   buf[STRLEN],buf2[STRLEN],*temp;
 502   const char *topinout;
 503   int    line;
 504   gmx_bool   bSystem,bMolecules,bSkip;
 505   int    i,nsol=0;
 506   double mtot;
 507   real   vol,mm;
 508
 509   for(i=0; (i<atoms->nres); i++) {
 510     /* calculate number of SOLvent molecules */
 511     if ( (strcmp(*atoms->resinfo[i].name,"SOL")==0) ||
 512          (strcmp(*atoms->resinfo[i].name,"WAT")==0) ||
 513          (strcmp(*atoms->resinfo[i].name,"HOH")==0) )
 514       nsol++;
 515   }
 516   mtot = 0;
 517   for(i=0; (i<atoms->nr); i++) {
 518     gmx_atomprop_query(aps,epropMass,
 519                        *atoms->resinfo[atoms->atom[i].resind].name,
 520                        *atoms->atomname[i],&mm);
 521     mtot += mm;
 522   }
 523
 524   vol=det(box);
 525
 526   fprintf(stderr,"Volume                 :  %10g (nm^3)\n",vol);
 527   fprintf(stderr,"Density                :  %10g (g/l)\n",
 528           (mtot*1e24)/(AVOGADRO*vol));
 529   fprintf(stderr,"Number of SOL molecules:  %5d   \n\n",nsol);
 530
 531   /* open topology file and append sol molecules */
 532   topinout  = ftp2fn(efTOP,NFILE,fnm);
 533   if ( ftp2bSet(efTOP,NFILE,fnm) ) {
 534     fprintf(stderr,"Processing topology\n");
 535     fpin = ffopen(topinout,"r");
 536     fpout= ffopen(TEMP_FILENM,"w");
 537     line=0;
 538     bSystem = bMolecules = FALSE;
 539     while (fgets(buf, STRLEN, fpin)) {
 540       bSkip=FALSE;
 541       line++;
 542       strcpy(buf2,buf);
 543       if ((temp=strchr(buf2,'\n')) != NULL)
 544         temp[0]='\0';
 545       ltrim(buf2);
 546       if (buf2[0]=='[') {
 547         buf2[0]=' ';
 548         if ((temp=strchr(buf2,'\n')) != NULL)
 549           temp[0]='\0';
 550         rtrim(buf2);
 551         if (buf2[strlen(buf2)-1]==']') {
 552           buf2[strlen(buf2)-1]='\0';
 553           ltrim(buf2);
 554           rtrim(buf2);
 555           bSystem=(gmx_strcasecmp(buf2,"system")==0);
 556           bMolecules=(gmx_strcasecmp(buf2,"molecules")==0);
 557         }
 558       } else if (bSystem && nsol && (buf[0]!=';') ) {
 559         /* if sol present, append "in water" to system name */
 560         rtrim(buf2);
 561         if (buf2[0] && (!strstr(buf2," water")) ) {
 562           sprintf(buf,"%s in water\n",buf2);
 563           bSystem = FALSE;
 564         }
 565       } else if (bMolecules) {
 566         /* check if this is a line with solvent molecules */
 567         sscanf(buf,"%s",buf2);
 568         if (strcmp(buf2,"SOL")==0) {
 569           sscanf(buf,"%*s %d",&i);
 570           nsol-=i;
 571           if (nsol<0) {
 572             bSkip=TRUE;
 573             nsol+=i;
 574           }
 575         }
 576       }
 577       if (bSkip) {
 578         if ((temp=strchr(buf,'\n')) != NULL)
 579           temp[0]='\0';
 580         fprintf(stdout,"Removing line #%d '%s' from topology file (%s)\n",
 581                 line,buf,topinout);
 582       } else
 583         fprintf(fpout,"%s",buf);
 584     }
 585     ffclose(fpin);
 586     if ( nsol ) {
 587       fprintf(stdout,"Adding line for %d solvent molecules to "
 588               "topology file (%s)\n",nsol,topinout);
 589       fprintf(fpout,"%-15s %5d\n","SOL",nsol);
 590     }
 591     ffclose(fpout);
 592     /* use ffopen to generate backup of topinout */
 593     fpout=ffopen(topinout,"w");
 594     ffclose(fpout);
 595     rename(TEMP_FILENM,topinout);
 596   }
 597 #undef TEMP_FILENM
 598 }
 599
 600 int gmx_genbox(int argc,char *argv[])
 601 {
 602   const char *desc[] = {
 603     "genbox can do one of 3 things:[PAR]",
 604
 605     "1) Generate a box of solvent. Specify [TT]-cs[tt] and [TT]-box[tt]. Or specify [TT]-cs[tt] and",
 606     "[TT]-cp[tt] with a structure file with a box, but without atoms.[PAR]",
 607
 608     "2) Solvate a solute configuration, e.g. a protein, in a bath of solvent ",
 609     "molecules. Specify [TT]-cp[tt] (solute) and [TT]-cs[tt] (solvent). ",
 610     "The box specified in the solute coordinate file ([TT]-cp[tt]) is used,",
 611     "unless [TT]-box[tt] is set.",
 612     "If you want the solute to be centered in the box,",
 613     "the program [TT]editconf[tt] has sophisticated options",
 614     "to change the box dimensions and center the solute.",
 615     "Solvent molecules are removed from the box where the ",
 616     "distance between any atom of the solute molecule(s) and any atom of ",
 617     "the solvent molecule is less than the sum of the van der Waals radii of ",
 618     "both atoms. A database ([TT]vdwradii.dat[tt]) of van der Waals radii is ",
 619     "read by the program, and atoms not in the database are ",
 620     "assigned a default distance [TT]-vdwd[tt].",
 621     "Note that this option will also influence the distances between",
 622     "solvent molecules if they contain atoms that are not in the database.",
 623     "[PAR]",
 624
 625     "3) Insert a number ([TT]-nmol[tt]) of extra molecules ([TT]-ci[tt]) ",
 626     "at random positions.",
 627     "The program iterates until [TT]nmol[tt] molecules",
 628     "have been inserted in the box. To test whether an insertion is ",
 629     "successful the same van der Waals criterium is used as for removal of ",
 630     "solvent molecules. When no appropriately ",
 631     "sized holes (holes that can hold an extra molecule) are available the ",
 632     "program tries for [TT]-nmol[tt] * [TT]-try[tt] times before giving up. ",
 633     "Increase [TT]-try[tt] if you have several small holes to fill.[PAR]",
 634
 635     "The default solvent is Simple Point Charge water (SPC), with coordinates ",
 636     "from [TT]$GMXLIB/spc216.gro[tt]. These coordinates can also be used",
 637     "for other 3-site water models, since a short equibilibration will remove",
 638     "the small differences between the models.",
 639     "Other solvents are also supported, as well as mixed solvents. The",
 640     "only restriction to solvent types is that a solvent molecule consists",
 641     "of exactly one residue. The residue information in the coordinate",
 642     "files is used, and should therefore be more or less consistent.",
 643     "In practice this means that two subsequent solvent molecules in the ",
 644     "solvent coordinate file should have different residue number.",
 645     "The box of solute is built by stacking the coordinates read from",
 646     "the coordinate file. This means that these coordinates should be ",
 647     "equlibrated in periodic boundary conditions to ensure a good",
 648     "alignment of molecules on the stacking interfaces.",
 649     "The [TT]-maxsol[tt] option simply adds only the first [TT]-maxsol[tt]",
 650     "solvent molecules and leaves out the rest would have fit into the box.",
 651     "[PAR]",
 652
 653     "The program can optionally rotate the solute molecule to align the",
 654     "longest molecule axis along a box edge. This way the amount of solvent",
 655     "molecules necessary is reduced.",
 656     "It should be kept in mind that this only works for",
 657     "short simulations, as e.g. an alpha-helical peptide in solution can ",
 658     "rotate over 90 degrees, within 500 ps. In general it is therefore ",
 659     "better to make a more or less cubic box.[PAR]",
 660
 661     "Setting [TT]-shell[tt] larger than zero will place a layer of water of",
 662     "the specified thickness (nm) around the solute. Hint: it is a good",
 663     "idea to put the protein in the center of a box first (using editconf).",
 664     "[PAR]",
 665
 666     "Finally, genbox will optionally remove lines from your topology file in ",
 667     "which a number of solvent molecules is already added, and adds a ",
 668     "line with the total number of solvent molecules in your coordinate file."
 669   };
 670
 671   const char *bugs[] = {
 672     "Molecules must be whole in the initial configurations.",
 673   };
 674
 675   /* parameter data */
 676   gmx_bool bSol,bProt,bBox;
 677   const char *conf_prot,*confout;
 678   int  bInsert;
 679   real *r;
 680   char *title_ins;
 681   gmx_atomprop_t aps;
 682
 683   /* protein configuration data */
 684   char    *title=NULL;
 685   t_atoms atoms;
 686   rvec    *x,*v=NULL;
 687   int     ePBC=-1;
 688   matrix  box;
 689   t_pbc   pbc;
 690
 691   /* other data types */
 692   int  atoms_added,residues_added;
 693
 694   t_filenm fnm[] = {
 695     { efSTX, "-cp", "protein", ffOPTRD },
 696     { efSTX, "-cs", "spc216",  ffLIBOPTRD},
 697     { efSTX, "-ci", "insert",  ffOPTRD},
 698     { efSTO, NULL,  NULL,      ffWRITE},
 699     { efTOP, NULL,  NULL,      ffOPTRW},
 700   };
 701 #define NFILE asize(fnm)
 702
 703   static int nmol_ins=0,nmol_try=10,seed=1997;
 704   static real r_distance=0.105,r_shell=0;
 705   static rvec new_box={0.0,0.0,0.0};
 706   static gmx_bool bReadV=FALSE;
 707   static int  max_sol = 0;
 708   output_env_t oenv;
 709   t_pargs pa[] = {
 710     { "-box",    FALSE, etRVEC, {new_box},
 711       "box size" },
 712     { "-nmol",   FALSE, etINT , {&nmol_ins},
 713       "no of extra molecules to insert" },
 714     { "-try",    FALSE, etINT , {&nmol_try},
 715       "try inserting -nmol*-try times" },
 716     { "-seed",   FALSE, etINT , {&seed},
 717       "random generator seed"},
 718     { "-vdwd",   FALSE, etREAL, {&r_distance},
 719       "default vdwaals distance"},
 720     { "-shell",  FALSE, etREAL, {&r_shell},
 721       "thickness of optional water layer around solute" },
 722     { "-maxsol", FALSE, etINT,  {&max_sol},
 723       "maximum number of solvent molecules to add if they fit in the box. If zero (default) this is ignored" },
 724     { "-vel",    FALSE, etBOOL, {&bReadV},
 725       "keep velocities from input solute and solvent" }
 726   };
 727
 728   CopyRight(stderr,argv[0]);
 729   parse_common_args(&argc,argv, PCA_BE_NICE,NFILE,fnm,asize(pa),pa,
 730                     asize(desc),desc,asize(bugs),bugs,&oenv);
 731
 732   bInsert   = opt2bSet("-ci",NFILE,fnm) && (nmol_ins > 0);
 733   bSol      = opt2bSet("-cs",NFILE,fnm);
 734   bProt     = opt2bSet("-cp",NFILE,fnm);
 735   bBox      = opt2parg_bSet("-box",asize(pa),pa);
 736
 737   /* check input */
 738   if (bInsert && nmol_ins<=0)
 739     gmx_fatal(FARGS,"When specifying inserted molecules (-ci), "
 740                 "-nmol must be larger than 0");
 741   if (!bProt && !bBox)
 742     gmx_fatal(FARGS,"When no solute (-cp) is specified, "
 743                 "a box size (-box) must be specified");
 744
 745   aps = gmx_atomprop_init();
 746
 747   if (bProt) {
 748     /*generate a solute configuration */
 749     conf_prot = opt2fn("-cp",NFILE,fnm);
 750     title = read_prot(conf_prot,&atoms,&x,bReadV?&v:NULL,&r,&ePBC,box,
 751                       aps,r_distance);
 752     if (bReadV && !v)
 753       fprintf(stderr,"Note: no velocities found\n");
 754     if (atoms.nr == 0) {
 755       fprintf(stderr,"Note: no atoms in %s\n",conf_prot);
 756       bProt = FALSE;
 757     }
 758   }
 759   if (!bProt) {
 760     atoms.nr=0;
 761     atoms.nres=0;
 762     atoms.resinfo=NULL;
 763     atoms.atomname=NULL;
 764     atoms.atom=NULL;
 765     atoms.pdbinfo=NULL;
 766     x=NULL;
 767     r=NULL;
 768   }
 769   if (bBox) {
 770     ePBC = epbcXYZ;
 771     clear_mat(box);
 772     box[XX][XX]=new_box[XX];
 773     box[YY][YY]=new_box[YY];
 774     box[ZZ][ZZ]=new_box[ZZ];
 775   }
 776   if (det(box) == 0)
 777     gmx_fatal(FARGS,"Undefined solute box.\nCreate one with editconf "
 778                 "or give explicit -box command line option");
 779
 780   /* add nmol_ins molecules of atoms_ins
 781      in random orientation at random place */
 782   if (bInsert)
 783     title_ins = insert_mols(opt2fn("-ci",NFILE,fnm),nmol_ins,nmol_try,seed,
 784                             &atoms,&x,&r,ePBC,box,aps,r_distance,r_shell,
 785                             oenv);
 786   else
 787     title_ins = strdup("Generated by genbox");
 788
 789   /* add solvent */
 790   if (bSol)
 791     add_solv(opt2fn("-cs",NFILE,fnm),&atoms,&x,v?&v:NULL,&r,ePBC,box,
 792              aps,r_distance,&atoms_added,&residues_added,r_shell,max_sol,
 793              oenv);
 794
 795   /* write new configuration 1 to file confout */
 796   confout = ftp2fn(efSTO,NFILE,fnm);
 797   fprintf(stderr,"Writing generated configuration to %s\n",confout);
 798   if (bProt) {
 799     write_sto_conf(confout,title,&atoms,x,v,ePBC,box);
 800     /* print box sizes and box type to stderr */
 801     fprintf(stderr,"%s\n",title);
 802   } else
 803     write_sto_conf(confout,title_ins,&atoms,x,v,ePBC,box);
 804
 805   /* print size of generated configuration */
 806   fprintf(stderr,"\nOutput configuration contains %d atoms in %d residues\n",
 807           atoms.nr,atoms.nres);
 808   update_top(&atoms,box,NFILE,fnm,aps);
 809
 810   gmx_atomprop_destroy(aps);
 811
 812   thanx(stderr);
 813
 814   return 0;
 815 }