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