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