src/gromacs/gmxpreprocess/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  * Copyright (c) 2013,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "genconf.h"
  38
  39 #include "gromacs/math/utilities.h"
  40 #include "macros.h"
  41 #include "smalloc.h"
  42 #include "sysstuff.h"
  43 #include "gromacs/fileio/confio.h"
  44 #include "gromacs/commandline/pargs.h"
  45 #include "vec.h"
  46 #include "gromacs/random/random.h"
  47 #include "gromacs/math/3dview.h"
  48 #include "txtdump.h"
  49 #include "readinp.h"
  50 #include "names.h"
  51 #include "sortwater.h"
  52 #include "gromacs/fileio/trxio.h"
  53
  54 static void rand_rot(int natoms, rvec x[], rvec v[], vec4 xrot[], vec4 vrot[],
  55                      gmx_rng_t rng, rvec max_rot)
  56 {
  57     mat4 mt1, mt2, mr[DIM], mtemp1, mtemp2, mtemp3, mxtot, mvtot;
  58     rvec xcm;
  59     real phi;
  60     int  i, m;
  61
  62     clear_rvec(xcm);
  63     for (i = 0; (i < natoms); i++)
  64     {
  65         for (m = 0; (m < DIM); m++)
  66         {
  67             xcm[m] += x[i][m]/natoms; /* get center of mass of one molecule  */
  68         }
  69     }
  70     fprintf(stderr, "center of geometry: %f, %f, %f\n", xcm[0], xcm[1], xcm[2]);
  71
  72     translate(-xcm[XX], -xcm[YY], -xcm[ZZ], mt1); /* move c.o.ma to origin */
  73     for (m = 0; (m < DIM); m++)
  74     {
  75         phi = M_PI*max_rot[m]*(2*gmx_rng_uniform_real(rng) - 1)/180;
  76         rotate(m, phi, mr[m]);
  77     }
  78     translate(xcm[XX], xcm[YY], xcm[ZZ], mt2);
  79
  80     /* For mult_matrix we need to multiply in the opposite order
  81      * compared to normal mathematical notation.
  82      */
  83     mult_matrix(mtemp1, mt1, mr[XX]);
  84     mult_matrix(mtemp2, mr[YY], mr[ZZ]);
  85     mult_matrix(mtemp3, mtemp1, mtemp2);
  86     mult_matrix(mxtot, mtemp3, mt2);
  87     mult_matrix(mvtot, mr[XX], mtemp2);
  88
  89     for (i = 0; (i < natoms); i++)
  90     {
  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     {
 104         for (m = 0; (m < DIM); m++)
 105         {
 106             xcm[m] += x[i][m];
 107         }
 108     }
 109     for (m = 0; (m < DIM); m++)
 110     {
 111         xcm[m] = 0.5*box[m][m]-xcm[m]/natoms;
 112     }
 113     for (i = 0; (i < natoms); i++)
 114     {
 115         for (m = 0; (m < DIM); m++)
 116         {
 117             x[i][m] += xcm[m];
 118         }
 119     }
 120 }
 121
 122 int gmx_genconf(int argc, char *argv[])
 123 {
 124     const char     *desc[] = {
 125         "[THISMODULE] multiplies a given coordinate file by simply stacking them",
 126         "on top of each other, like a small child playing with wooden blocks.",
 127         "The program makes a grid of [IT]user-defined[it]",
 128         "proportions ([TT]-nbox[tt]), ",
 129         "and interspaces the grid point with an extra space [TT]-dist[tt].[PAR]",
 130         "When option [TT]-rot[tt] is used the program does not check for overlap",
 131         "between molecules on grid points. It is recommended to make the box in",
 132         "the input file at least as big as the coordinates + ",
 133         "van der Waals radius.[PAR]",
 134         "If the optional trajectory file is given, conformations are not",
 135         "generated, but read from this file and translated appropriately to",
 136         "build the grid."
 137
 138     };
 139     const char     *bugs[] = {
 140         "The program should allow for random displacement of lattice points."
 141     };
 142
 143     int             vol;
 144     t_atoms        *atoms;      /* list with all atoms */
 145     char            title[STRLEN];
 146     rvec           *x, *xx, *v; /* coordinates? */
 147     real            t;
 148     vec4           *xrot, *vrot;
 149     int             ePBC;
 150     matrix          box, boxx; /* box length matrix */
 151     rvec            shift;
 152     int             natoms;    /* number of atoms in one molecule  */
 153     int             nres;      /* number of molecules? */
 154     int             i, j, k, l, m, ndx, nrdx, nx, ny, nz;
 155     t_trxstatus    *status;
 156     gmx_bool        bTRX;
 157     output_env_t    oenv;
 158     gmx_rng_t       rng;
 159
 160     t_filenm        fnm[] = {
 161         { efSTX, "-f", "conf", ffREAD  },
 162         { efSTO, "-o", "out",  ffWRITE },
 163         { efTRX, "-trj", NULL,  ffOPTRD }
 164     };
 165 #define NFILE asize(fnm)
 166     static rvec     nrbox    = {1, 1, 1};
 167     static int      seed     = 0;    /* seed for random number generator */
 168     static int      nmolat   = 3;
 169     static int      nblock   = 1;
 170     static gmx_bool bShuffle = FALSE;
 171     static gmx_bool bSort    = FALSE;
 172     static gmx_bool bRandom  = FALSE;           /* False: no random rotations */
 173     static gmx_bool bRenum   = TRUE;            /* renumber residues */
 174     static rvec     dist     = {0, 0, 0};       /* space added between molecules ? */
 175     static rvec     max_rot  = {180, 180, 180}; /* maximum rotation */
 176     t_pargs         pa[]     = {
 177         { "-nbox",   FALSE, etRVEC, {nrbox},   "Number of boxes" },
 178         { "-dist",   FALSE, etRVEC, {dist},    "Distance between boxes" },
 179         { "-seed",   FALSE, etINT,  {&seed},
 180           "Random generator seed, if 0 generated from the time" },
 181         { "-rot",    FALSE, etBOOL, {&bRandom}, "Randomly rotate conformations" },
 182         { "-shuffle", FALSE, etBOOL, {&bShuffle}, "Random shuffling of molecules" },
 183         { "-sort",   FALSE, etBOOL, {&bSort},   "Sort molecules on X coord" },
 184         { "-block",  FALSE, etINT,  {&nblock},
 185           "Divide the box in blocks on this number of cpus" },
 186         { "-nmolat", FALSE, etINT,  {&nmolat},
 187           "Number of atoms per molecule, assumed to start from 0. "
 188           "If you set this wrong, it will screw up your system!" },
 189         { "-maxrot", FALSE, etRVEC, {max_rot}, "Maximum random rotation" },
 190         { "-renumber", FALSE, etBOOL, {&bRenum},  "Renumber residues" }
 191     };
 192
 193     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
 194                            asize(desc), desc, asize(bugs), bugs, &oenv))
 195     {
 196         return 0;
 197     }
 198
 199     if (seed == 0)
 200     {
 201         rng = gmx_rng_init(gmx_rng_make_seed());
 202     }
 203     else
 204     {
 205         rng = gmx_rng_init(seed);
 206     }
 207
 208     bTRX = ftp2bSet(efTRX, NFILE, fnm);
 209     nx   = (int)(nrbox[XX]+0.5);
 210     ny   = (int)(nrbox[YY]+0.5);
 211     nz   = (int)(nrbox[ZZ]+0.5);
 212
 213     if ((nx <= 0) || (ny <= 0) || (nz <= 0))
 214     {
 215         gmx_fatal(FARGS, "Number of boxes (-nbox) should be larger than zero");
 216     }
 217     if ((nmolat <= 0) && bShuffle)
 218     {
 219         gmx_fatal(FARGS, "Can not shuffle if the molecules only have %d atoms",
 220                   nmolat);
 221     }
 222
 223     vol = nx*ny*nz; /* calculate volume in grid points (= nr. molecules) */
 224
 225     get_stx_coordnum(opt2fn("-f", NFILE, fnm), &natoms);
 226     snew(atoms, 1);
 227     /* make space for all the atoms */
 228     init_t_atoms(atoms, natoms*vol, FALSE);
 229     snew(x, natoms*vol);           /* get space for coordinates of all atoms */
 230     snew(xrot, natoms);            /* get space for rotation matrix? */
 231     snew(v, natoms*vol);           /* velocities. not really needed? */
 232     snew(vrot, natoms);
 233     /* set atoms->nr to the number in one box *
 234      * to avoid complaints in read_stx_conf   *
 235      */
 236     atoms->nr = natoms;
 237     read_stx_conf(opt2fn("-f", NFILE, fnm), title, atoms, x, v, &ePBC, box);
 238
 239     nres = atoms->nres;            /* nr of residues in one element? */
 240
 241     if (bTRX)
 242     {
 243         if (!read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &xx, boxx))
 244         {
 245             gmx_fatal(FARGS, "No atoms in trajectory %s", ftp2fn(efTRX, NFILE, fnm));
 246         }
 247     }
 248     else
 249     {
 250         snew(xx, natoms);
 251         for (i = 0; i < natoms; i++)
 252         {
 253             copy_rvec(x[i], xx[i]);
 254         }
 255     }
 256
 257
 258     for (k = 0; (k < nz); k++)     /* loop over all gridpositions    */
 259     {
 260         shift[ZZ] = k*(dist[ZZ]+box[ZZ][ZZ]);
 261
 262         for (j = 0; (j < ny); j++)
 263         {
 264             shift[YY] = j*(dist[YY]+box[YY][YY])+k*box[ZZ][YY];
 265
 266             for (i = 0; (i < nx); i++)
 267             {
 268                 shift[XX] = i*(dist[XX]+box[XX][XX])+j*box[YY][XX]+k*box[ZZ][XX];
 269
 270                 ndx  = (i*ny*nz+j*nz+k)*natoms;
 271                 nrdx = (i*ny*nz+j*nz+k)*nres;
 272
 273                 /* Random rotation on input coords */
 274                 if (bRandom)
 275                 {
 276                     rand_rot(natoms, xx, v, xrot, vrot, rng, max_rot);
 277                 }
 278
 279                 for (l = 0; (l < natoms); l++)
 280                 {
 281                     for (m = 0; (m < DIM); m++)
 282                     {
 283                         if (bRandom)
 284                         {
 285                             x[ndx+l][m] = xrot[l][m];
 286                             v[ndx+l][m] = vrot[l][m];
 287                         }
 288                         else
 289                         {
 290                             x[ndx+l][m] = xx[l][m];
 291                             v[ndx+l][m] = v[l][m];
 292                         }
 293                     }
 294                     if (ePBC == epbcSCREW && i % 2 == 1)
 295                     {
 296                         /* Rotate around x axis */
 297                         for (m = YY; m <= ZZ; m++)
 298                         {
 299                             x[ndx+l][m] = box[YY][m] + box[ZZ][m] - x[ndx+l][m];
 300                             v[ndx+l][m] = -v[ndx+l][m];
 301                         }
 302                     }
 303                     for (m = 0; (m < DIM); m++)
 304                     {
 305                         x[ndx+l][m] += shift[m];
 306                     }
 307                     atoms->atom[ndx+l].resind = nrdx + atoms->atom[l].resind;
 308                     atoms->atomname[ndx+l]    = atoms->atomname[l];
 309                 }
 310
 311                 for (l = 0; (l < nres); l++)
 312                 {
 313                     atoms->resinfo[nrdx+l] = atoms->resinfo[l];
 314                     if (bRenum)
 315                     {
 316                         atoms->resinfo[nrdx+l].nr += nrdx;
 317                     }
 318                 }
 319                 if (bTRX)
 320                 {
 321                     if (!read_next_x(oenv, status, &t, xx, boxx) &&
 322                         ((i+1)*(j+1)*(k+1) < vol))
 323                     {
 324                         gmx_fatal(FARGS, "Not enough frames in trajectory");
 325                     }
 326                 }
 327             }
 328         }
 329     }
 330     if (bTRX)
 331     {
 332         close_trj(status);
 333     }
 334
 335     /* make box bigger */
 336     for (m = 0; (m < DIM); m++)
 337     {
 338         box[m][m] += dist[m];
 339     }
 340     svmul(nx, box[XX], box[XX]);
 341     svmul(ny, box[YY], box[YY]);
 342     svmul(nz, box[ZZ], box[ZZ]);
 343     if (ePBC == epbcSCREW && nx % 2 == 0)
 344     {
 345         /* With an even number of boxes in x we can forgot about the screw */
 346         ePBC = epbcXYZ;
 347     }
 348
 349     /* move_x(natoms*vol,x,box); */          /* put atoms in box? */
 350
 351     atoms->nr   *= vol;
 352     atoms->nres *= vol;
 353
 354     /*depending on how you look at it, this is either a nasty hack or the way it should work*/
 355     if (bRenum)
 356     {
 357         for (i = 0; i < atoms->nres; i++)
 358         {
 359             atoms->resinfo[i].nr = i+1;
 360         }
 361     }
 362
 363
 364     if (bShuffle)
 365     {
 366         randwater(0, atoms->nr/nmolat, nmolat, x, v, rng);
 367     }
 368     else if (bSort)
 369     {
 370         sortwater(0, atoms->nr/nmolat, nmolat, x, v);
 371     }
 372     else if (opt2parg_bSet("-block", asize(pa), pa))
 373     {
 374         mkcompact(0, atoms->nr/nmolat, nmolat, x, v, nblock, box);
 375     }
 376     gmx_rng_destroy(rng);
 377
 378     write_sto_conf(opt2fn("-o", NFILE, fnm), title, atoms, x, v, ePBC, box);
 379
 380     return 0;
 381 }