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