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 "gromacs/fileio/confio.h"
  42 #include "txtdump.h"
  43 #include "readinp.h"
  44 #include "names.h"
  45 #include "sortwater.h"
  46 #include "gromacs/fileio/trxio.h"
  47
  48 #include "gromacs/commandline/pargs.h"
  49 #include "gromacs/math/3dtransforms.h"
  50 #include "gromacs/math/vec.h"
  51 #include "gromacs/random/random.h"
  52 #include "gromacs/utility/fatalerror.h"
  53 #include "gromacs/utility/smalloc.h"
  54
  55 static void rand_rot(int natoms, rvec x[], rvec v[], vec4 xrot[], vec4 vrot[],
  56                      gmx_rng_t rng, 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     {
  66         for (m = 0; (m < DIM); m++)
  67         {
  68             xcm[m] += x[i][m]/natoms; /* get center of mass of one molecule  */
  69         }
  70     }
  71     fprintf(stderr, "center of geometry: %f, %f, %f\n", xcm[0], xcm[1], xcm[2]);
  72
  73     /* move c.o.ma to origin */
  74     gmx_mat4_init_translation(-xcm[XX], -xcm[YY], -xcm[ZZ], mt1);
  75     for (m = 0; (m < DIM); m++)
  76     {
  77         phi = M_PI*max_rot[m]*(2*gmx_rng_uniform_real(rng) - 1)/180;
  78         gmx_mat4_init_rotation(m, phi, mr[m]);
  79     }
  80     gmx_mat4_init_translation(xcm[XX], xcm[YY], xcm[ZZ], mt2);
  81
  82     /* For gmx_mat4_mmul() we need to multiply in the opposite order
  83      * compared to normal mathematical notation.
  84      */
  85     gmx_mat4_mmul(mtemp1, mt1, mr[XX]);
  86     gmx_mat4_mmul(mtemp2, mr[YY], mr[ZZ]);
  87     gmx_mat4_mmul(mtemp3, mtemp1, mtemp2);
  88     gmx_mat4_mmul(mxtot, mtemp3, mt2);
  89     gmx_mat4_mmul(mvtot, mr[XX], mtemp2);
  90
  91     for (i = 0; (i < natoms); i++)
  92     {
  93         gmx_mat4_transform_point(mxtot, x[i], xrot[i]);
  94         gmx_mat4_transform_point(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         "[THISMODULE] 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     gmx_rng_t       rng;
 161
 162     t_filenm        fnm[] = {
 163         { efSTX, "-f", "conf", ffREAD  },
 164         { efSTO, "-o", "out",  ffWRITE },
 165         { efTRX, "-trj", NULL,  ffOPTRD }
 166     };
 167 #define NFILE asize(fnm)
 168     static rvec     nrbox    = {1, 1, 1};
 169     static int      seed     = 0;    /* seed for random number generator */
 170     static int      nmolat   = 3;
 171     static int      nblock   = 1;
 172     static gmx_bool bShuffle = FALSE;
 173     static gmx_bool bSort    = FALSE;
 174     static gmx_bool bRandom  = FALSE;           /* False: no random rotations */
 175     static gmx_bool bRenum   = TRUE;            /* renumber residues */
 176     static rvec     dist     = {0, 0, 0};       /* space added between molecules ? */
 177     static rvec     max_rot  = {180, 180, 180}; /* maximum rotation */
 178     t_pargs         pa[]     = {
 179         { "-nbox",   FALSE, etRVEC, {nrbox},   "Number of boxes" },
 180         { "-dist",   FALSE, etRVEC, {dist},    "Distance between boxes" },
 181         { "-seed",   FALSE, etINT,  {&seed},
 182           "Random generator seed, if 0 generated from the time" },
 183         { "-rot",    FALSE, etBOOL, {&bRandom}, "Randomly rotate conformations" },
 184         { "-shuffle", FALSE, etBOOL, {&bShuffle}, "Random shuffling of molecules" },
 185         { "-sort",   FALSE, etBOOL, {&bSort},   "Sort molecules on X coord" },
 186         { "-block",  FALSE, etINT,  {&nblock},
 187           "Divide the box in blocks on this number of cpus" },
 188         { "-nmolat", FALSE, etINT,  {&nmolat},
 189           "Number of atoms per molecule, assumed to start from 0. "
 190           "If you set this wrong, it will screw up your system!" },
 191         { "-maxrot", FALSE, etRVEC, {max_rot}, "Maximum random rotation" },
 192         { "-renumber", FALSE, etBOOL, {&bRenum},  "Renumber residues" }
 193     };
 194
 195     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
 196                            asize(desc), desc, asize(bugs), bugs, &oenv))
 197     {
 198         return 0;
 199     }
 200
 201     if (seed == 0)
 202     {
 203         rng = gmx_rng_init(gmx_rng_make_seed());
 204     }
 205     else
 206     {
 207         rng = gmx_rng_init(seed);
 208     }
 209
 210     bTRX = ftp2bSet(efTRX, NFILE, fnm);
 211     nx   = (int)(nrbox[XX]+0.5);
 212     ny   = (int)(nrbox[YY]+0.5);
 213     nz   = (int)(nrbox[ZZ]+0.5);
 214
 215     if ((nx <= 0) || (ny <= 0) || (nz <= 0))
 216     {
 217         gmx_fatal(FARGS, "Number of boxes (-nbox) should be larger than zero");
 218     }
 219     if ((nmolat <= 0) && bShuffle)
 220     {
 221         gmx_fatal(FARGS, "Can not shuffle if the molecules only have %d atoms",
 222                   nmolat);
 223     }
 224
 225     vol = nx*ny*nz; /* calculate volume in grid points (= nr. molecules) */
 226
 227     get_stx_coordnum(opt2fn("-f", NFILE, fnm), &natoms);
 228     snew(atoms, 1);
 229     /* make space for all the atoms */
 230     init_t_atoms(atoms, natoms*vol, FALSE);
 231     snew(x, natoms*vol);           /* get space for coordinates of all atoms */
 232     snew(xrot, natoms);            /* get space for rotation matrix? */
 233     snew(v, natoms*vol);           /* velocities. not really needed? */
 234     snew(vrot, natoms);
 235     /* set atoms->nr to the number in one box *
 236      * to avoid complaints in read_stx_conf   *
 237      */
 238     atoms->nr = natoms;
 239     read_stx_conf(opt2fn("-f", NFILE, fnm), title, atoms, x, v, &ePBC, box);
 240
 241     nres = atoms->nres;            /* nr of residues in one element? */
 242
 243     if (bTRX)
 244     {
 245         if (!read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &xx, boxx))
 246         {
 247             gmx_fatal(FARGS, "No atoms in trajectory %s", ftp2fn(efTRX, NFILE, fnm));
 248         }
 249     }
 250     else
 251     {
 252         snew(xx, natoms);
 253         for (i = 0; i < natoms; i++)
 254         {
 255             copy_rvec(x[i], xx[i]);
 256         }
 257     }
 258
 259
 260     for (k = 0; (k < nz); k++)     /* loop over all gridpositions    */
 261     {
 262         shift[ZZ] = k*(dist[ZZ]+box[ZZ][ZZ]);
 263
 264         for (j = 0; (j < ny); j++)
 265         {
 266             shift[YY] = j*(dist[YY]+box[YY][YY])+k*box[ZZ][YY];
 267
 268             for (i = 0; (i < nx); i++)
 269             {
 270                 shift[XX] = i*(dist[XX]+box[XX][XX])+j*box[YY][XX]+k*box[ZZ][XX];
 271
 272                 ndx  = (i*ny*nz+j*nz+k)*natoms;
 273                 nrdx = (i*ny*nz+j*nz+k)*nres;
 274
 275                 /* Random rotation on input coords */
 276                 if (bRandom)
 277                 {
 278                     rand_rot(natoms, xx, v, xrot, vrot, rng, max_rot);
 279                 }
 280
 281                 for (l = 0; (l < natoms); l++)
 282                 {
 283                     for (m = 0; (m < DIM); m++)
 284                     {
 285                         if (bRandom)
 286                         {
 287                             x[ndx+l][m] = xrot[l][m];
 288                             v[ndx+l][m] = vrot[l][m];
 289                         }
 290                         else
 291                         {
 292                             x[ndx+l][m] = xx[l][m];
 293                             v[ndx+l][m] = v[l][m];
 294                         }
 295                     }
 296                     if (ePBC == epbcSCREW && i % 2 == 1)
 297                     {
 298                         /* Rotate around x axis */
 299                         for (m = YY; m <= ZZ; m++)
 300                         {
 301                             x[ndx+l][m] = box[YY][m] + box[ZZ][m] - x[ndx+l][m];
 302                             v[ndx+l][m] = -v[ndx+l][m];
 303                         }
 304                     }
 305                     for (m = 0; (m < DIM); m++)
 306                     {
 307                         x[ndx+l][m] += shift[m];
 308                     }
 309                     atoms->atom[ndx+l].resind = nrdx + atoms->atom[l].resind;
 310                     atoms->atomname[ndx+l]    = atoms->atomname[l];
 311                 }
 312
 313                 for (l = 0; (l < nres); l++)
 314                 {
 315                     atoms->resinfo[nrdx+l] = atoms->resinfo[l];
 316                     if (bRenum)
 317                     {
 318                         atoms->resinfo[nrdx+l].nr += nrdx;
 319                     }
 320                 }
 321                 if (bTRX)
 322                 {
 323                     if (!read_next_x(oenv, status, &t, xx, boxx) &&
 324                         ((i+1)*(j+1)*(k+1) < vol))
 325                     {
 326                         gmx_fatal(FARGS, "Not enough frames in trajectory");
 327                     }
 328                 }
 329             }
 330         }
 331     }
 332     if (bTRX)
 333     {
 334         close_trj(status);
 335     }
 336
 337     /* make box bigger */
 338     for (m = 0; (m < DIM); m++)
 339     {
 340         box[m][m] += dist[m];
 341     }
 342     svmul(nx, box[XX], box[XX]);
 343     svmul(ny, box[YY], box[YY]);
 344     svmul(nz, box[ZZ], box[ZZ]);
 345     if (ePBC == epbcSCREW && nx % 2 == 0)
 346     {
 347         /* With an even number of boxes in x we can forgot about the screw */
 348         ePBC = epbcXYZ;
 349     }
 350
 351     /* move_x(natoms*vol,x,box); */          /* put atoms in box? */
 352
 353     atoms->nr   *= vol;
 354     atoms->nres *= vol;
 355
 356     /*depending on how you look at it, this is either a nasty hack or the way it should work*/
 357     if (bRenum)
 358     {
 359         for (i = 0; i < atoms->nres; i++)
 360         {
 361             atoms->resinfo[i].nr = i+1;
 362         }
 363     }
 364
 365
 366     if (bShuffle)
 367     {
 368         randwater(0, atoms->nr/nmolat, nmolat, x, v, rng);
 369     }
 370     else if (bSort)
 371     {
 372         sortwater(0, atoms->nr/nmolat, nmolat, x, v);
 373     }
 374     else if (opt2parg_bSet("-block", asize(pa), pa))
 375     {
 376         mkcompact(0, atoms->nr/nmolat, nmolat, x, v, nblock, box);
 377     }
 378     gmx_rng_destroy(rng);
 379
 380     write_sto_conf(opt2fn("-o", NFILE, fnm), title, atoms, x, v, ePBC, box);
 381
 382     return 0;
 383 }