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 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "genconf.h"
  41
  42 #include "gromacs/commandline/pargs.h"
  43 #include "gromacs/fileio/confio.h"
  44 #include "gromacs/fileio/trxio.h"
  45 #include "gromacs/math/3dtransforms.h"
  46 #include "gromacs/math/units.h"
  47 #include "gromacs/math/utilities.h"
  48 #include "gromacs/math/vec.h"
  49 #include "gromacs/mdtypes/md_enums.h"
  50 #include "gromacs/pbcutil/pbc.h"
  51 #include "gromacs/random/threefry.h"
  52 #include "gromacs/random/uniformrealdistribution.h"
  53 #include "gromacs/topology/mtop_util.h"
  54 #include "gromacs/topology/topology.h"
  55 #include "gromacs/utility/arraysize.h"
  56 #include "gromacs/utility/cstringutil.h"
  57 #include "gromacs/utility/fatalerror.h"
  58 #include "gromacs/utility/smalloc.h"
  59
  60 static void
  61 rand_rot(int natoms, rvec x[], rvec v[], vec4 xrot[], vec4 vrot[], gmx::DefaultRandomEngine* rng, const rvec max_rot)
  62 {
  63     mat4                               mt1, mt2, mr[DIM], mtemp1, mtemp2, mtemp3, mxtot, mvtot;
  64     rvec                               xcm;
  65     real                               phi;
  66     int                                i, m;
  67     gmx::UniformRealDistribution<real> dist(-1.0, 1.0);
  68
  69     clear_rvec(xcm);
  70     for (i = 0; (i < natoms); i++)
  71     {
  72         for (m = 0; (m < DIM); m++)
  73         {
  74             xcm[m] += x[i][m] / natoms; /* get center of mass of one molecule  */
  75         }
  76     }
  77     fprintf(stderr, "center of geometry: %f, %f, %f\n", xcm[0], xcm[1], xcm[2]);
  78
  79     /* move c.o.ma to origin */
  80     gmx_mat4_init_translation(-xcm[XX], -xcm[YY], -xcm[ZZ], mt1);
  81     for (m = 0; (m < DIM); m++)
  82     {
  83         phi = M_PI * max_rot[m] * dist(*rng) / 180;
  84         gmx_mat4_init_rotation(m, phi, mr[m]);
  85     }
  86     gmx_mat4_init_translation(xcm[XX], xcm[YY], xcm[ZZ], mt2);
  87
  88     /* For gmx_mat4_mmul() we need to multiply in the opposite order
  89      * compared to normal mathematical notation.
  90      */
  91     gmx_mat4_mmul(mtemp1, mt1, mr[XX]);
  92     gmx_mat4_mmul(mtemp2, mr[YY], mr[ZZ]);
  93     gmx_mat4_mmul(mtemp3, mtemp1, mtemp2);
  94     gmx_mat4_mmul(mxtot, mtemp3, mt2);
  95     gmx_mat4_mmul(mvtot, mr[XX], mtemp2);
  96
  97     for (i = 0; (i < natoms); i++)
  98     {
  99         gmx_mat4_transform_point(mxtot, x[i], xrot[i]);
 100         gmx_mat4_transform_point(mvtot, v[i], vrot[i]);
 101     }
 102 }
 103
 104 int gmx_genconf(int argc, char* argv[])
 105 {
 106     const char* desc[] = {
 107         "[THISMODULE] multiplies a given coordinate file by simply stacking them",
 108         "on top of each other, like a small child playing with wooden blocks.",
 109         "The program makes a grid of [IT]user-defined[it]",
 110         "proportions ([TT]-nbox[tt]), ",
 111         "and interspaces the grid point with an extra space [TT]-dist[tt].[PAR]",
 112         "When option [TT]-rot[tt] is used the program does not check for overlap",
 113         "between molecules on grid points. It is recommended to make the box in",
 114         "the input file at least as big as the coordinates + ",
 115         "van der Waals radius.[PAR]",
 116         "If the optional trajectory file is given, conformations are not",
 117         "generated, but read from this file and translated appropriately to",
 118         "build the grid."
 119
 120     };
 121     const char* bugs[] = { "The program should allow for random displacement of lattice points." };
 122
 123     int               vol;
 124     rvec *            x, *xx, *v; /* coordinates? */
 125     real              t;
 126     vec4 *            xrot, *vrot;
 127     PbcType           pbcType;
 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     bool              bTRX;
 135     gmx_output_env_t* oenv;
 136
 137     t_filenm fnm[] = { { efSTX, "-f", "conf", ffREAD },
 138                        { efSTO, "-o", "out", ffWRITE },
 139                        { efTRX, "-trj", nullptr, ffOPTRD } };
 140 #define NFILE asize(fnm)
 141     rvec     nrbox   = { 1, 1, 1 };
 142     int      seed    = 0;                 /* seed for random number generator */
 143     gmx_bool bRandom = FALSE;             /* False: no random rotations */
 144     gmx_bool bRenum  = TRUE;              /* renumber residues */
 145     rvec     dist    = { 0, 0, 0 };       /* space added between molecules ? */
 146     rvec     max_rot = { 180, 180, 180 }; /* maximum rotation */
 147     t_pargs  pa[]    = {
 148         { "-nbox", FALSE, etRVEC, { nrbox }, "Number of boxes" },
 149         { "-dist", FALSE, etRVEC, { dist }, "Distance between boxes" },
 150         { "-seed", FALSE, etINT, { &seed }, "Random generator seed (0 means generate)" },
 151         { "-rot", FALSE, etBOOL, { &bRandom }, "Randomly rotate conformations" },
 152         { "-maxrot", FALSE, etRVEC, { max_rot }, "Maximum random rotation" },
 153         { "-renumber", FALSE, etBOOL, { &bRenum }, "Renumber residues" }
 154     };
 155
 156     if (!parse_common_args(
 157                 &argc, argv, 0, NFILE, fnm, asize(pa), pa, asize(desc), desc, asize(bugs), bugs, &oenv))
 158     {
 159         return 0;
 160     }
 161
 162     if (seed == 0)
 163     {
 164         seed = static_cast<int>(gmx::makeRandomSeed());
 165     }
 166     gmx::DefaultRandomEngine rng(seed);
 167
 168     bTRX = ftp2bSet(efTRX, NFILE, fnm);
 169     nx   = gmx::roundToInt(nrbox[XX]);
 170     ny   = gmx::roundToInt(nrbox[YY]);
 171     nz   = gmx::roundToInt(nrbox[ZZ]);
 172
 173     if ((nx <= 0) || (ny <= 0) || (nz <= 0))
 174     {
 175         gmx_fatal(FARGS, "Number of boxes (-nbox) should be larger than zero");
 176     }
 177
 178     vol = nx * ny * nz; /* calculate volume in grid points (= nr. molecules) */
 179
 180     gmx_mtop_t mtop;
 181     bool       haveTop = false;
 182     readConfAndTopology(opt2fn("-f", NFILE, fnm), &haveTop, &mtop, &pbcType, &x, &v, box);
 183     t_atoms atoms = gmx_mtop_global_atoms(mtop);
 184     natoms        = atoms.nr;
 185     nres          = atoms.nres; /* nr of residues in one element? */
 186     /* make space for all the atoms */
 187     add_t_atoms(&atoms, natoms * (vol - 1), nres * (vol - 1));
 188     srenew(x, natoms * vol); /* get space for coordinates of all atoms */
 189     srenew(v, natoms * vol); /* velocities. not really needed? */
 190     snew(xrot, natoms);      /* get space for rotation matrix? */
 191     snew(vrot, natoms);
 192
 193     if (bTRX)
 194     {
 195         if (!read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &xx, boxx))
 196         {
 197             gmx_fatal(FARGS, "No atoms in trajectory %s", ftp2fn(efTRX, NFILE, fnm));
 198         }
 199     }
 200     else
 201     {
 202         snew(xx, natoms);
 203         for (i = 0; i < natoms; i++)
 204         {
 205             copy_rvec(x[i], xx[i]);
 206         }
 207     }
 208
 209
 210     for (k = 0; (k < nz); k++) /* loop over all gridpositions    */
 211     {
 212         shift[ZZ] = k * (dist[ZZ] + box[ZZ][ZZ]);
 213
 214         for (j = 0; (j < ny); j++)
 215         {
 216             shift[YY] = j * (dist[YY] + box[YY][YY]) + k * box[ZZ][YY];
 217
 218             for (i = 0; (i < nx); i++)
 219             {
 220                 shift[XX] = i * (dist[XX] + box[XX][XX]) + j * box[YY][XX] + k * box[ZZ][XX];
 221
 222                 ndx  = (i * ny * nz + j * nz + k) * natoms;
 223                 nrdx = (i * ny * nz + j * nz + k) * nres;
 224
 225                 /* Random rotation on input coords */
 226                 if (bRandom)
 227                 {
 228                     rand_rot(natoms, xx, v, xrot, vrot, &rng, max_rot);
 229                 }
 230
 231                 for (l = 0; (l < natoms); l++)
 232                 {
 233                     for (m = 0; (m < DIM); m++)
 234                     {
 235                         if (bRandom)
 236                         {
 237                             x[ndx + l][m] = xrot[l][m];
 238                             v[ndx + l][m] = vrot[l][m];
 239                         }
 240                         else
 241                         {
 242                             x[ndx + l][m] = xx[l][m];
 243                             v[ndx + l][m] = v[l][m];
 244                         }
 245                     }
 246                     if (pbcType == PbcType::Screw && i % 2 == 1)
 247                     {
 248                         /* Rotate around x axis */
 249                         for (m = YY; m <= ZZ; m++)
 250                         {
 251                             x[ndx + l][m] = box[YY][m] + box[ZZ][m] - x[ndx + l][m];
 252                             v[ndx + l][m] = -v[ndx + l][m];
 253                         }
 254                     }
 255                     for (m = 0; (m < DIM); m++)
 256                     {
 257                         x[ndx + l][m] += shift[m];
 258                     }
 259                     atoms.atom[ndx + l].resind = nrdx + atoms.atom[l].resind;
 260                     atoms.atomname[ndx + l]    = atoms.atomname[l];
 261                 }
 262
 263                 for (l = 0; (l < nres); l++)
 264                 {
 265                     atoms.resinfo[nrdx + l] = atoms.resinfo[l];
 266                     if (bRenum)
 267                     {
 268                         atoms.resinfo[nrdx + l].nr += nrdx;
 269                     }
 270                 }
 271                 if (bTRX)
 272                 {
 273                     if (!read_next_x(oenv, status, &t, xx, boxx) && ((i + 1) * (j + 1) * (k + 1) < vol))
 274                     {
 275                         gmx_fatal(FARGS, "Not enough frames in trajectory");
 276                     }
 277                 }
 278             }
 279         }
 280     }
 281     if (bTRX)
 282     {
 283         close_trx(status);
 284     }
 285
 286     /* make box bigger */
 287     for (m = 0; (m < DIM); m++)
 288     {
 289         box[m][m] += dist[m];
 290     }
 291     svmul(nx, box[XX], box[XX]);
 292     svmul(ny, box[YY], box[YY]);
 293     svmul(nz, box[ZZ], box[ZZ]);
 294     if (pbcType == PbcType::Screw && nx % 2 == 0)
 295     {
 296         /* With an even number of boxes in x we can forgot about the screw */
 297         pbcType = PbcType::Xyz;
 298     }
 299
 300     /*depending on how you look at it, this is either a nasty hack or the way it should work*/
 301     if (bRenum)
 302     {
 303         for (i = 0; i < atoms.nres; i++)
 304         {
 305             atoms.resinfo[i].nr = i + 1;
 306         }
 307     }
 308
 309     write_sto_conf(opt2fn("-o", NFILE, fnm), *mtop.name, &atoms, x, v, pbcType, box);
 310
 311     sfree(x);
 312     sfree(v);
 313     sfree(xrot);
 314     sfree(vrot);
 315     sfree(xx);
 316     done_atom(&atoms);
 317     output_env_done(oenv);
 318
 319     return 0;
 320 }