src/tools/gmx_genconf.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 #include "maths.h"
  43 #include "macros.h"
  44 #include "copyrite.h"
  45 #include "string2.h"
  46 #include "smalloc.h"
  47 #include "sysstuff.h"
  48 #include "confio.h"
  49 #include "statutil.h"
  50 #include "vec.h"
  51 #include "random.h"
  52 #include "3dview.h"
  53 #include "txtdump.h"
  54 #include "readinp.h"
  55 #include "names.h"
  56 #include "sortwater.h"
  57 #include "gmx_ana.h"
  58
  59 static void rand_rot(int natoms,rvec x[],rvec v[],vec4 xrot[],vec4 vrot[],
  60                      int *seed,rvec max_rot)
  61 {
  62   mat4 mt1,mt2,mr[DIM],mtemp1,mtemp2,mtemp3,mxtot,mvtot;
  63   rvec xcm;
  64   real phi;
  65   int  i,m;
  66
  67   clear_rvec(xcm);
  68   for(i=0; (i<natoms); i++)
  69     for(m=0; (m<DIM); m++) {
  70       xcm[m]+=x[i][m]/natoms;   /* get center of mass of one molecule  */
  71     }
  72   fprintf(stderr,"center of geometry: %f, %f, %f\n",xcm[0],xcm[1],xcm[2]);
  73
  74   translate(-xcm[XX],-xcm[YY],-xcm[ZZ],mt1);  /* move c.o.ma to origin */
  75   for(m=0; (m<DIM); m++) {
  76     phi=M_PI*max_rot[m]*(2*rando(seed) - 1)/180;
  77     rotate(m,phi,mr[m]);
  78   }
  79   translate(xcm[XX],xcm[YY],xcm[ZZ],mt2);
  80
  81   /* For mult_matrix we need to multiply in the opposite order
  82    * compared to normal mathematical notation.
  83    */
  84   mult_matrix(mtemp1,mt1,mr[XX]);
  85   mult_matrix(mtemp2,mr[YY],mr[ZZ]);
  86   mult_matrix(mtemp3,mtemp1,mtemp2);
  87   mult_matrix(mxtot,mtemp3,mt2);
  88   mult_matrix(mvtot,mr[XX],mtemp2);
  89
  90   for(i=0; (i<natoms); i++) {
  91     m4_op(mxtot,x[i],xrot[i]);
  92     m4_op(mvtot,v[i],vrot[i]);
  93   }
  94 }
  95
  96 static void move_x(int natoms,rvec x[],matrix box)
  97 {
  98   int  i,m;
  99   rvec xcm;
 100
 101   clear_rvec(xcm);
 102   for(i=0; (i<natoms); i++)
 103     for(m=0; (m<DIM); m++)
 104       xcm[m]+=x[i][m];
 105   for(m=0; (m<DIM); m++)
 106     xcm[m] = 0.5*box[m][m]-xcm[m]/natoms;
 107   for(i=0; (i<natoms); i++)
 108     for(m=0; (m<DIM); m++)
 109       x[i][m] += xcm[m];
 110 }
 111
 112 int gmx_genconf(int argc, char *argv[])
 113 {
 114   const char *desc[] = {
 115     "[TT]genconf[tt] multiplies a given coordinate file by simply stacking them",
 116     "on top of each other, like a small child playing with wooden blocks.",
 117     "The program makes a grid of [IT]user-defined[it]",
 118     "proportions ([TT]-nbox[tt]), ",
 119     "and interspaces the grid point with an extra space [TT]-dist[tt].[PAR]",
 120     "When option [TT]-rot[tt] is used the program does not check for overlap",
 121     "between molecules on grid points. It is recommended to make the box in",
 122     "the input file at least as big as the coordinates + ",
 123     "van der Waals radius.[PAR]",
 124     "If the optional trajectory file is given, conformations are not",
 125     "generated, but read from this file and translated appropriately to",
 126     "build the grid."
 127
 128   };
 129   const char *bugs[] = {
 130     "The program should allow for random displacement of lattice points." };
 131
 132   int     vol;
 133   t_atoms *atoms;       /* list with all atoms */
 134   char    title[STRLEN];
 135   rvec    *x,*xx,*v;        /* coordinates? */
 136   real    t;
 137   vec4    *xrot,*vrot;
 138   int     ePBC;
 139   matrix  box,boxx;          /* box length matrix */
 140   rvec    shift;
 141   int     natoms;       /* number of atoms in one molecule  */
 142   int     nres;         /* number of molecules? */
 143   int     i,j,k,l,m,ndx,nrdx,nx,ny,nz;
 144   t_trxstatus *status;
 145   gmx_bool    bTRX;
 146   output_env_t oenv;
 147
 148   t_filenm fnm[] = {
 149     { efSTX, "-f", "conf", ffREAD  },
 150     { efSTO, "-o", "out",  ffWRITE },
 151     { efTRX, "-trj",NULL,  ffOPTRD }
 152   };
 153 #define NFILE asize(fnm)
 154   static rvec nrbox    = {1,1,1};
 155   static int  seed     = 0;          /* seed for random number generator */
 156   static int  nmolat   = 3;
 157   static int  nblock   = 1;
 158   static gmx_bool bShuffle = FALSE;
 159   static gmx_bool bSort    = FALSE;
 160   static gmx_bool bRandom  = FALSE;      /* False: no random rotations */
 161   static gmx_bool bRenum   = TRUE;       /* renumber residues */
 162   static rvec dist     = {0,0,0};    /* space added between molecules ? */
 163   static rvec max_rot  = {180,180,180}; /* maximum rotation */
 164   t_pargs pa[] = {
 165     { "-nbox",   FALSE, etRVEC, {nrbox},   "Number of boxes" },
 166     { "-dist",   FALSE, etRVEC, {dist},    "Distance between boxes" },
 167     { "-seed",   FALSE, etINT,  {&seed},
 168       "Random generator seed, if 0 generated from the time" },
 169     { "-rot",    FALSE, etBOOL, {&bRandom}, "Randomly rotate conformations" },
 170     { "-shuffle",FALSE, etBOOL, {&bShuffle},"Random shuffling of molecules" },
 171     { "-sort",   FALSE, etBOOL, {&bSort},   "Sort molecules on X coord" },
 172     { "-block",  FALSE, etINT,  {&nblock},
 173       "Divide the box in blocks on this number of cpus" },
 174     { "-nmolat", FALSE, etINT,  {&nmolat},
 175       "Number of atoms per molecule, assumed to start from 0. "
 176       "If you set this wrong, it will screw up your system!" },
 177     { "-maxrot", FALSE, etRVEC, {max_rot}, "Maximum random rotation" },
 178     { "-renumber",FALSE,etBOOL, {&bRenum},  "Renumber residues" }
 179   };
 180
 181   CopyRight(stderr,argv[0]);
 182   parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
 183                     asize(desc),desc,asize(bugs),bugs,&oenv);
 184
 185   if (bRandom && (seed == 0))
 186     seed = make_seed();
 187
 188   bTRX = ftp2bSet(efTRX,NFILE,fnm);
 189   nx   = (int)(nrbox[XX]+0.5);
 190   ny   = (int)(nrbox[YY]+0.5);
 191   nz   = (int)(nrbox[ZZ]+0.5);
 192
 193   if ((nx <= 0) || (ny <= 0) || (nz <= 0))
 194     gmx_fatal(FARGS,"Number of boxes (-nbox) should be larger than zero");
 195   if ((nmolat <= 0) && bShuffle)
 196     gmx_fatal(FARGS,"Can not shuffle if the molecules only have %d atoms",
 197                 nmolat);
 198
 199   vol=nx*ny*nz;     /* calculate volume in grid points (= nr. molecules) */
 200
 201   get_stx_coordnum(opt2fn("-f",NFILE,fnm),&natoms);
 202   snew(atoms,1);
 203   /* make space for all the atoms */
 204   init_t_atoms(atoms,natoms*vol,FALSE);
 205   snew(x,natoms*vol);              /* get space for coordinates of all atoms */
 206   snew(xrot,natoms);               /* get space for rotation matrix? */
 207   snew(v,natoms*vol);              /* velocities. not really needed? */
 208   snew(vrot,natoms);
 209   /* set atoms->nr to the number in one box *
 210    * to avoid complaints in read_stx_conf   *
 211    */
 212   atoms->nr = natoms;
 213   read_stx_conf(opt2fn("-f",NFILE,fnm),title,atoms,x,v,&ePBC,box);
 214
 215   nres=atoms->nres;                /* nr of residues in one element? */
 216
 217   if (bTRX) {
 218     if (!read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&xx,boxx))
 219       gmx_fatal(FARGS,"No atoms in trajectory %s",ftp2fn(efTRX,NFILE,fnm));
 220   } else {
 221     snew(xx,natoms);
 222     for(i=0; i<natoms; i++) {
 223       copy_rvec(x[i],xx[i]);
 224     }
 225   }
 226
 227
 228   for(k=0; (k<nz); k++) {          /* loop over all gridpositions    */
 229     shift[ZZ]=k*(dist[ZZ]+box[ZZ][ZZ]);
 230
 231     for(j=0; (j<ny); j++) {
 232       shift[YY]=j*(dist[YY]+box[YY][YY])+k*box[ZZ][YY];
 233
 234       for(i=0; (i<nx); i++)  {
 235          shift[XX]=i*(dist[XX]+box[XX][XX])+j*box[YY][XX]+k*box[ZZ][XX];
 236
 237         ndx=(i*ny*nz+j*nz+k)*natoms;
 238         nrdx=(i*ny*nz+j*nz+k)*nres;
 239
 240         /* Random rotation on input coords */
 241         if (bRandom)
 242           rand_rot(natoms,xx,v,xrot,vrot,&seed,max_rot);
 243
 244         for(l=0; (l<natoms); l++) {
 245           for(m=0; (m<DIM); m++) {
 246             if (bRandom) {
 247               x[ndx+l][m] = xrot[l][m];
 248               v[ndx+l][m] = vrot[l][m];
 249             }
 250             else {
 251               x[ndx+l][m] = xx[l][m];
 252               v[ndx+l][m] = v[l][m];
 253             }
 254           }
 255           if (ePBC == epbcSCREW && i % 2 == 1) {
 256             /* Rotate around x axis */
 257             for(m=YY; m<=ZZ; m++) {
 258               x[ndx+l][m] = box[YY][m] + box[ZZ][m] - x[ndx+l][m];
 259               v[ndx+l][m] = -v[ndx+l][m];
 260             }
 261           }
 262           for(m=0; (m<DIM); m++) {
 263             x[ndx+l][m] += shift[m];
 264           }
 265           atoms->atom[ndx+l].resind = nrdx + atoms->atom[l].resind;
 266           atoms->atomname[ndx+l]=atoms->atomname[l];
 267         }
 268
 269         for(l=0; (l<nres); l++) {
 270           atoms->resinfo[nrdx+l] = atoms->resinfo[l];
 271           if (bRenum)
 272             atoms->resinfo[nrdx+l].nr += nrdx;
 273         }
 274         if (bTRX)
 275           if (!read_next_x(oenv,status,&t,natoms,xx,boxx) &&
 276               ((i+1)*(j+1)*(k+1) < vol))
 277             gmx_fatal(FARGS,"Not enough frames in trajectory");
 278       }
 279     }
 280   }
 281   if (bTRX)
 282     close_trj(status);
 283
 284   /* make box bigger */
 285   for(m=0; (m<DIM); m++)
 286     box[m][m] += dist[m];
 287   svmul(nx,box[XX],box[XX]);
 288   svmul(ny,box[YY],box[YY]);
 289   svmul(nz,box[ZZ],box[ZZ]);
 290   if (ePBC == epbcSCREW && nx % 2 == 0) {
 291     /* With an even number of boxes in x we can forgot about the screw */
 292     ePBC = epbcXYZ;
 293   }
 294
 295   /* move_x(natoms*vol,x,box); */          /* put atoms in box? */
 296
 297   atoms->nr*=vol;
 298   atoms->nres*=vol;
 299
 300   /*depending on how you look at it, this is either a nasty hack or the way it should work*/
 301   if (bRenum)
 302     for (i=0;i<atoms->nres;i++)
 303           atoms->resinfo[i].nr=i+1;
 304
 305
 306   if (bShuffle)
 307     randwater(0,atoms->nr/nmolat,nmolat,x,v,&seed);
 308   else if (bSort)
 309     sortwater(0,atoms->nr/nmolat,nmolat,x,v);
 310   else if (opt2parg_bSet("-block",asize(pa),pa))
 311     mkcompact(0,atoms->nr/nmolat,nmolat,x,v,nblock,box);
 312
 313   write_sto_conf(opt2fn("-o",NFILE,fnm),title,atoms,x,v,ePBC,box);
 314
 315   thanx(stderr);
 316
 317   return 0;
 318 }