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