src/gromacs/gmxpreprocess/genconf.cpp

   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,2015,2016, 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 "gmxpre.h"
  38
  39 #include "genconf.h"
  40
  41 #include "gromacs/commandline/pargs.h"
  42 #include "gromacs/fileio/confio.h"
  43 #include "gromacs/fileio/trxio.h"
  44 #include "gromacs/math/3dtransforms.h"
  45 #include "gromacs/math/utilities.h"
  46 #include "gromacs/math/vec.h"
  47 #include "gromacs/mdtypes/md_enums.h"
  48 #include "gromacs/pbcutil/pbc.h"
  49 #include "gromacs/random/threefry.h"
  50 #include "gromacs/random/uniformrealdistribution.h"
  51 #include "gromacs/topology/topology.h"
  52 #include "gromacs/utility/arraysize.h"
  53 #include "gromacs/utility/cstringutil.h"
  54 #include "gromacs/utility/fatalerror.h"
  55 #include "gromacs/utility/smalloc.h"
  56
  57 static void rand_rot(int natoms, rvec x[], rvec v[], vec4 xrot[], vec4 vrot[],
  58                      gmx::DefaultRandomEngine * rng, rvec max_rot)
  59 {
  60     mat4 mt1, mt2, mr[DIM], mtemp1, mtemp2, mtemp3, mxtot, mvtot;
  61     rvec xcm;
  62     real phi;
  63     int  i, m;
  64     gmx::UniformRealDistribution<real> dist(-1.0, 1.0);
  65
  66     clear_rvec(xcm);
  67     for (i = 0; (i < natoms); i++)
  68     {
  69         for (m = 0; (m < DIM); m++)
  70         {
  71             xcm[m] += x[i][m]/natoms; /* get center of mass of one molecule  */
  72         }
  73     }
  74     fprintf(stderr, "center of geometry: %f, %f, %f\n", xcm[0], xcm[1], xcm[2]);
  75
  76     /* move c.o.ma to origin */
  77     gmx_mat4_init_translation(-xcm[XX], -xcm[YY], -xcm[ZZ], mt1);
  78     for (m = 0; (m < DIM); m++)
  79     {
  80         phi = M_PI*max_rot[m]*dist(*rng)/180;
  81         gmx_mat4_init_rotation(m, phi, mr[m]);
  82     }
  83     gmx_mat4_init_translation(xcm[XX], xcm[YY], xcm[ZZ], mt2);
  84
  85     /* For gmx_mat4_mmul() we need to multiply in the opposite order
  86      * compared to normal mathematical notation.
  87      */
  88     gmx_mat4_mmul(mtemp1, mt1, mr[XX]);
  89     gmx_mat4_mmul(mtemp2, mr[YY], mr[ZZ]);
  90     gmx_mat4_mmul(mtemp3, mtemp1, mtemp2);
  91     gmx_mat4_mmul(mxtot, mtemp3, mt2);
  92     gmx_mat4_mmul(mvtot, mr[XX], mtemp2);
  93
  94     for (i = 0; (i < natoms); i++)
  95     {
  96         gmx_mat4_transform_point(mxtot, x[i], xrot[i]);
  97         gmx_mat4_transform_point(mvtot, v[i], vrot[i]);
  98     }
  99 }
 100
 101 int gmx_genconf(int argc, char *argv[])
 102 {
 103     const char       *desc[] = {
 104         "[THISMODULE] multiplies a given coordinate file by simply stacking them",
 105         "on top of each other, like a small child playing with wooden blocks.",
 106         "The program makes a grid of [IT]user-defined[it]",
 107         "proportions ([TT]-nbox[tt]), ",
 108         "and interspaces the grid point with an extra space [TT]-dist[tt].[PAR]",
 109         "When option [TT]-rot[tt] is used the program does not check for overlap",
 110         "between molecules on grid points. It is recommended to make the box in",
 111         "the input file at least as big as the coordinates + ",
 112         "van der Waals radius.[PAR]",
 113         "If the optional trajectory file is given, conformations are not",
 114         "generated, but read from this file and translated appropriately to",
 115         "build the grid."
 116
 117     };
 118     const char       *bugs[] = {
 119         "The program should allow for random displacement of lattice points."
 120     };
 121
 122     int               vol;
 123     t_atoms          *atoms;      /* list with all atoms */
 124     rvec             *x, *xx, *v; /* coordinates? */
 125     real              t;
 126     vec4             *xrot, *vrot;
 127     int               ePBC;
 128     matrix            box, boxx; /* box length matrix */
 129     rvec              shift;
 130     int               natoms;    /* number of atoms in one molecule  */
 131     int               nres;      /* number of molecules? */
 132     int               i, j, k, l, m, ndx, nrdx, nx, ny, nz;
 133     t_trxstatus      *status;
 134     gmx_bool          bTRX;
 135     gmx_output_env_t *oenv;
 136
 137     t_filenm          fnm[] = {
 138         { efSTX, "-f", "conf", ffREAD  },
 139         { efSTO, "-o", "out",  ffWRITE },
 140         { efTRX, "-trj", NULL,  ffOPTRD }
 141     };
 142 #define NFILE asize(fnm)
 143     static rvec       nrbox    = {1, 1, 1};
 144     static int        seed     = 0;               /* seed for random number generator */
 145     static gmx_bool   bRandom  = FALSE;           /* False: no random rotations */
 146     static gmx_bool   bRenum   = TRUE;            /* renumber residues */
 147     static rvec       dist     = {0, 0, 0};       /* space added between molecules ? */
 148     static rvec       max_rot  = {180, 180, 180}; /* maximum rotation */
 149     t_pargs           pa[]     = {
 150         { "-nbox",   FALSE, etRVEC, {nrbox},   "Number of boxes" },
 151         { "-dist",   FALSE, etRVEC, {dist},    "Distance between boxes" },
 152         { "-seed",   FALSE, etINT,  {&seed},
 153           "Random generator seed (0 means generate)" },
 154         { "-rot",    FALSE, etBOOL, {&bRandom}, "Randomly rotate conformations" },
 155         { "-maxrot", FALSE, etRVEC, {max_rot}, "Maximum random rotation" },
 156         { "-renumber", FALSE, etBOOL, {&bRenum},  "Renumber residues" }
 157     };
 158
 159     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
 160                            asize(desc), desc, asize(bugs), bugs, &oenv))
 161     {
 162         return 0;
 163     }
 164
 165     if (seed == 0)
 166     {
 167         seed = static_cast<int>(gmx::makeRandomSeed());
 168     }
 169     gmx::DefaultRandomEngine rng(seed);
 170
 171     bTRX = ftp2bSet(efTRX, NFILE, fnm);
 172     nx   = (int)(nrbox[XX]+0.5);
 173     ny   = (int)(nrbox[YY]+0.5);
 174     nz   = (int)(nrbox[ZZ]+0.5);
 175
 176     if ((nx <= 0) || (ny <= 0) || (nz <= 0))
 177     {
 178         gmx_fatal(FARGS, "Number of boxes (-nbox) should be larger than zero");
 179     }
 180
 181     vol = nx*ny*nz; /* calculate volume in grid points (= nr. molecules) */
 182
 183     t_topology *top;
 184     snew(top, 1);
 185     atoms = &top->atoms;
 186     read_tps_conf(opt2fn("-f", NFILE, fnm), top, &ePBC, &x, &v, box, FALSE);
 187     natoms = atoms->nr;
 188     nres   = atoms->nres;          /* nr of residues in one element? */
 189     /* make space for all the atoms */
 190     add_t_atoms(atoms, natoms*(vol-1), nres*(vol-1));
 191     srenew(x, natoms*vol);         /* get space for coordinates of all atoms */
 192     srenew(v, natoms*vol);         /* velocities. not really needed? */
 193     snew(xrot, natoms);            /* get space for rotation matrix? */
 194     snew(vrot, natoms);
 195
 196     if (bTRX)
 197     {
 198         if (!read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &xx, boxx))
 199         {
 200             gmx_fatal(FARGS, "No atoms in trajectory %s", ftp2fn(efTRX, NFILE, fnm));
 201         }
 202     }
 203     else
 204     {
 205         snew(xx, natoms);
 206         for (i = 0; i < natoms; i++)
 207         {
 208             copy_rvec(x[i], xx[i]);
 209         }
 210     }
 211
 212
 213     for (k = 0; (k < nz); k++)     /* loop over all gridpositions    */
 214     {
 215         shift[ZZ] = k*(dist[ZZ]+box[ZZ][ZZ]);
 216
 217         for (j = 0; (j < ny); j++)
 218         {
 219             shift[YY] = j*(dist[YY]+box[YY][YY])+k*box[ZZ][YY];
 220
 221             for (i = 0; (i < nx); i++)
 222             {
 223                 shift[XX] = i*(dist[XX]+box[XX][XX])+j*box[YY][XX]+k*box[ZZ][XX];
 224
 225                 ndx  = (i*ny*nz+j*nz+k)*natoms;
 226                 nrdx = (i*ny*nz+j*nz+k)*nres;
 227
 228                 /* Random rotation on input coords */
 229                 if (bRandom)
 230                 {
 231                     rand_rot(natoms, xx, v, xrot, vrot, &rng, max_rot);
 232                 }
 233
 234                 for (l = 0; (l < natoms); l++)
 235                 {
 236                     for (m = 0; (m < DIM); m++)
 237                     {
 238                         if (bRandom)
 239                         {
 240                             x[ndx+l][m] = xrot[l][m];
 241                             v[ndx+l][m] = vrot[l][m];
 242                         }
 243                         else
 244                         {
 245                             x[ndx+l][m] = xx[l][m];
 246                             v[ndx+l][m] = v[l][m];
 247                         }
 248                     }
 249                     if (ePBC == epbcSCREW && i % 2 == 1)
 250                     {
 251                         /* Rotate around x axis */
 252                         for (m = YY; m <= ZZ; m++)
 253                         {
 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                     {
 260                         x[ndx+l][m] += shift[m];
 261                     }
 262                     atoms->atom[ndx+l].resind = nrdx + atoms->atom[l].resind;
 263                     atoms->atomname[ndx+l]    = atoms->atomname[l];
 264                 }
 265
 266                 for (l = 0; (l < nres); l++)
 267                 {
 268                     atoms->resinfo[nrdx+l] = atoms->resinfo[l];
 269                     if (bRenum)
 270                     {
 271                         atoms->resinfo[nrdx+l].nr += nrdx;
 272                     }
 273                 }
 274                 if (bTRX)
 275                 {
 276                     if (!read_next_x(oenv, status, &t, xx, boxx) &&
 277                         ((i+1)*(j+1)*(k+1) < vol))
 278                     {
 279                         gmx_fatal(FARGS, "Not enough frames in trajectory");
 280                     }
 281                 }
 282             }
 283         }
 284     }
 285     if (bTRX)
 286     {
 287         close_trj(status);
 288     }
 289
 290     /* make box bigger */
 291     for (m = 0; (m < DIM); m++)
 292     {
 293         box[m][m] += dist[m];
 294     }
 295     svmul(nx, box[XX], box[XX]);
 296     svmul(ny, box[YY], box[YY]);
 297     svmul(nz, box[ZZ], box[ZZ]);
 298     if (ePBC == epbcSCREW && nx % 2 == 0)
 299     {
 300         /* With an even number of boxes in x we can forgot about the screw */
 301         ePBC = epbcXYZ;
 302     }
 303
 304     /*depending on how you look at it, this is either a nasty hack or the way it should work*/
 305     if (bRenum)
 306     {
 307         for (i = 0; i < atoms->nres; i++)
 308         {
 309             atoms->resinfo[i].nr = i+1;
 310         }
 311     }
 312
 313     write_sto_conf(opt2fn("-o", NFILE, fnm), *top->name, atoms, x, v, ePBC, box);
 314
 315     sfree(x);
 316     sfree(v);
 317     sfree(xrot);
 318     sfree(vrot);
 319     sfree(xx);
 320     done_top(top);
 321     sfree(top);
 322     done_filenms(NFILE, fnm);
 323     output_env_done(oenv);
 324
 325     return 0;
 326 }