src/tools/gmx_genbox.c

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