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