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