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     "genconf 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 off 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,status=-1;
 141   bool    bTRX;
 142   output_env_t oenv;
 143
 144   t_filenm fnm[] = {
 145     { efSTX, "-f", "conf", ffREAD  },
 146     { efSTO, "-o", "out",  ffWRITE },
 147     { efTRX, "-trj",NULL,  ffOPTRD }
 148   };
 149 #define NFILE asize(fnm)
 150   static rvec nrbox    = {1,1,1};
 151   static int  seed     = 0;          /* seed for random number generator */
 152   static int  nmolat   = 3;
 153   static int  nblock   = 1;
 154   static bool bShuffle = FALSE;
 155   static bool bSort    = FALSE;
 156   static bool bRandom  = FALSE;      /* False: no random rotations */
 157   static bool bRenum   = TRUE;       /* renumber residues */
 158   static rvec dist     = {0,0,0};    /* space added between molecules ? */
 159   static rvec max_rot  = {180,180,180}; /* maximum rotation */
 160   t_pargs pa[] = {
 161     { "-nbox",   FALSE, etRVEC, {nrbox},   "Number of boxes" },
 162     { "-dist",   FALSE, etRVEC, {dist},    "Distance between boxes" },
 163     { "-seed",   FALSE, etINT,  {&seed},
 164       "Random generator seed, if 0 generated from the time" },
 165     { "-rot",    FALSE, etBOOL, {&bRandom}, "Randomly rotate conformations" },
 166     { "-shuffle",FALSE, etBOOL, {&bShuffle},"Random shuffling of molecules" },
 167     { "-sort",   FALSE, etBOOL, {&bSort},   "Sort molecules on X coord" },
 168     { "-block",  FALSE, etINT,  {&nblock},
 169       "Divide the box in blocks on this number of cpus" },
 170     { "-nmolat", FALSE, etINT,  {&nmolat},
 171       "Number of atoms per molecule, assumed to start from 0. "
 172       "If you set this wrong, it will screw up your system!" },
 173     { "-maxrot", FALSE, etRVEC, {max_rot}, "Maximum random rotation" },
 174     { "-renumber",FALSE,etBOOL, {&bRenum},  "Renumber residues" }
 175   };
 176
 177   CopyRight(stderr,argv[0]);
 178   parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
 179                     asize(desc),desc,asize(bugs),bugs,&oenv);
 180
 181   if (bRandom && (seed == 0))
 182     seed = make_seed();
 183
 184   bTRX = ftp2bSet(efTRX,NFILE,fnm);
 185   nx   = (int)(nrbox[XX]+0.5);
 186   ny   = (int)(nrbox[YY]+0.5);
 187   nz   = (int)(nrbox[ZZ]+0.5);
 188
 189   if ((nx <= 0) || (ny <= 0) || (nz <= 0))
 190     gmx_fatal(FARGS,"Number of boxes (-nbox) should be larger than zero");
 191   if ((nmolat <= 0) && bShuffle)
 192     gmx_fatal(FARGS,"Can not shuffle if the molecules only have %d atoms",
 193                 nmolat);
 194
 195   vol=nx*ny*nz;     /* calculate volume in grid points (= nr. molecules) */
 196
 197   get_stx_coordnum(opt2fn("-f",NFILE,fnm),&natoms);
 198   snew(atoms,1);
 199   /* make space for all the atoms */
 200   init_t_atoms(atoms,natoms*vol,FALSE);
 201   snew(x,natoms*vol);              /* get space for coordinates of all atoms */
 202   snew(xrot,natoms);               /* get space for rotation matrix? */
 203   snew(v,natoms*vol);              /* velocities. not really needed? */
 204   snew(vrot,natoms);
 205   /* set atoms->nr to the number in one box *
 206    * to avoid complaints in read_stx_conf   *
 207    */
 208   atoms->nr = natoms;
 209   read_stx_conf(opt2fn("-f",NFILE,fnm),title,atoms,x,v,&ePBC,box);
 210
 211   nres=atoms->nres;                /* nr of residues in one element? */
 212
 213   if (bTRX) {
 214     if (!read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&xx,boxx))
 215       gmx_fatal(FARGS,"No atoms in trajectory %s",ftp2fn(efTRX,NFILE,fnm));
 216   } else {
 217     snew(xx,natoms);
 218     for(i=0; i<natoms; i++) {
 219       copy_rvec(x[i],xx[i]);
 220     }
 221   }
 222
 223
 224   for(k=0; (k<nz); k++) {          /* loop over all gridpositions    */
 225     shift[ZZ]=k*(dist[ZZ]+box[ZZ][ZZ]);
 226
 227     for(j=0; (j<ny); j++) {
 228       shift[YY]=j*(dist[YY]+box[YY][YY])+k*box[ZZ][YY];
 229
 230       for(i=0; (i<nx); i++)  {
 231          shift[XX]=i*(dist[XX]+box[XX][XX])+j*box[YY][XX]+k*box[ZZ][XX];
 232
 233         ndx=(i*ny*nz+j*nz+k)*natoms;
 234         nrdx=(i*ny*nz+j*nz+k)*nres;
 235
 236         /* Random rotation on input coords */
 237         if (bRandom)
 238           rand_rot(natoms,xx,v,xrot,vrot,&seed,max_rot);
 239
 240         for(l=0; (l<natoms); l++) {
 241           for(m=0; (m<DIM); m++) {
 242             if (bRandom) {
 243               x[ndx+l][m] = xrot[l][m];
 244               v[ndx+l][m] = vrot[l][m];
 245             }
 246             else {
 247               x[ndx+l][m] = xx[l][m];
 248               v[ndx+l][m] = v[l][m];
 249             }
 250           }
 251           if (ePBC == epbcSCREW && i % 2 == 1) {
 252             /* Rotate around x axis */
 253             for(m=YY; m<=ZZ; m++) {
 254               x[ndx+l][m] = box[YY][m] + box[ZZ][m] - x[ndx+l][m];
 255               v[ndx+l][m] = -v[ndx+l][m];
 256             }
 257           }
 258           for(m=0; (m<DIM); m++) {
 259             x[ndx+l][m] += shift[m];
 260           }
 261           atoms->atom[ndx+l].resind = nrdx + atoms->atom[l].resind;
 262           atoms->atomname[ndx+l]=atoms->atomname[l];
 263         }
 264
 265         for(l=0; (l<nres); l++) {
 266           atoms->resinfo[nrdx+l] = atoms->resinfo[l];
 267           if (bRenum)
 268             atoms->resinfo[nrdx+l].nr += nrdx;
 269         }
 270         if (bTRX)
 271           if (!read_next_x(oenv,status,&t,natoms,xx,boxx) &&
 272               ((i+1)*(j+1)*(k+1) < vol))
 273             gmx_fatal(FARGS,"Not enough frames in trajectory");
 274       }
 275     }
 276   }
 277   if (bTRX)
 278     close_trj(status);
 279
 280   /* make box bigger */
 281   for(m=0; (m<DIM); m++)
 282     box[m][m] += dist[m];
 283   svmul(nx,box[XX],box[XX]);
 284   svmul(ny,box[YY],box[YY]);
 285   svmul(nz,box[ZZ],box[ZZ]);
 286   if (ePBC == epbcSCREW && nx % 2 == 0) {
 287     /* With an even number of boxes in x we can forgot about the screw */
 288     ePBC = epbcXYZ;
 289   }
 290
 291   /* move_x(natoms*vol,x,box); */          /* put atoms in box? */
 292
 293   atoms->nr*=vol;
 294   atoms->nres*=vol;
 295
 296   /*depending on how you look at it, this is either a nasty hack or the way it should work*/
 297   if (bRenum)
 298     for (i=0;i<atoms->nres;i++)
 299           atoms->resinfo[i].nr=i+1;
 300
 301
 302   if (bShuffle)
 303     randwater(0,atoms->nr/nmolat,nmolat,x,v,&seed);
 304   else if (bSort)
 305     sortwater(0,atoms->nr/nmolat,nmolat,x,v);
 306   else if (opt2parg_bSet("-block",asize(pa),pa))
 307     mkcompact(0,atoms->nr/nmolat,nmolat,x,v,nblock,box);
 308
 309   write_sto_conf(opt2fn("-o",NFILE,fnm),title,atoms,x,v,ePBC,box);
 310
 311   thanx(stderr);
 312
 313   return 0;
 314 }