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