2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
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"
61 rand_rot(int natoms, rvec x[], rvec v[], vec4 xrot[], vec4 vrot[], gmx::DefaultRandomEngine* rng, const rvec max_rot)
63 mat4 mt1, mt2, mr[DIM], mtemp1, mtemp2, mtemp3, mxtot, mvtot;
67 gmx::UniformRealDistribution<real> dist(-1.0, 1.0);
70 for (i = 0; (i < natoms); i++)
72 for (m = 0; (m < DIM); m++)
74 xcm[m] += x[i][m] / natoms; /* get center of mass of one molecule */
77 fprintf(stderr, "center of geometry: %f, %f, %f\n", xcm[0], xcm[1], xcm[2]);
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++)
83 phi = M_PI * max_rot[m] * dist(*rng) / 180;
84 gmx_mat4_init_rotation(m, phi, mr[m]);
86 gmx_mat4_init_translation(xcm[XX], xcm[YY], xcm[ZZ], mt2);
88 /* For gmx_mat4_mmul() we need to multiply in the opposite order
89 * compared to normal mathematical notation.
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);
97 for (i = 0; (i < natoms); i++)
99 gmx_mat4_transform_point(mxtot, x[i], xrot[i]);
100 gmx_mat4_transform_point(mvtot, v[i], vrot[i]);
104 int gmx_genconf(int argc, char* argv[])
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",
121 const char* bugs[] = { "The program should allow for random displacement of lattice points." };
124 rvec * x, *xx, *v; /* coordinates? */
128 matrix box, boxx; /* box length matrix */
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;
135 gmx_output_env_t* oenv;
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 */
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" }
156 if (!parse_common_args(
157 &argc, argv, 0, NFILE, fnm, asize(pa), pa, asize(desc), desc, asize(bugs), bugs, &oenv))
164 seed = static_cast<int>(gmx::makeRandomSeed());
166 gmx::DefaultRandomEngine rng(seed);
168 bTRX = ftp2bSet(efTRX, NFILE, fnm);
169 nx = gmx::roundToInt(nrbox[XX]);
170 ny = gmx::roundToInt(nrbox[YY]);
171 nz = gmx::roundToInt(nrbox[ZZ]);
173 if ((nx <= 0) || (ny <= 0) || (nz <= 0))
175 gmx_fatal(FARGS, "Number of boxes (-nbox) should be larger than zero");
178 vol = nx * ny * nz; /* calculate volume in grid points (= nr. molecules) */
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);
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? */
195 if (!read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &xx, boxx))
197 gmx_fatal(FARGS, "No atoms in trajectory %s", ftp2fn(efTRX, NFILE, fnm));
203 for (i = 0; i < natoms; i++)
205 copy_rvec(x[i], xx[i]);
210 for (k = 0; (k < nz); k++) /* loop over all gridpositions */
212 shift[ZZ] = k * (dist[ZZ] + box[ZZ][ZZ]);
214 for (j = 0; (j < ny); j++)
216 shift[YY] = j * (dist[YY] + box[YY][YY]) + k * box[ZZ][YY];
218 for (i = 0; (i < nx); i++)
220 shift[XX] = i * (dist[XX] + box[XX][XX]) + j * box[YY][XX] + k * box[ZZ][XX];
222 ndx = (i * ny * nz + j * nz + k) * natoms;
223 nrdx = (i * ny * nz + j * nz + k) * nres;
225 /* Random rotation on input coords */
228 rand_rot(natoms, xx, v, xrot, vrot, &rng, max_rot);
231 for (l = 0; (l < natoms); l++)
233 for (m = 0; (m < DIM); m++)
237 x[ndx + l][m] = xrot[l][m];
238 v[ndx + l][m] = vrot[l][m];
242 x[ndx + l][m] = xx[l][m];
243 v[ndx + l][m] = v[l][m];
246 if (pbcType == PbcType::Screw && i % 2 == 1)
248 /* Rotate around x axis */
249 for (m = YY; m <= ZZ; m++)
251 x[ndx + l][m] = box[YY][m] + box[ZZ][m] - x[ndx + l][m];
252 v[ndx + l][m] = -v[ndx + l][m];
255 for (m = 0; (m < DIM); m++)
257 x[ndx + l][m] += shift[m];
259 atoms.atom[ndx + l].resind = nrdx + atoms.atom[l].resind;
260 atoms.atomname[ndx + l] = atoms.atomname[l];
263 for (l = 0; (l < nres); l++)
265 atoms.resinfo[nrdx + l] = atoms.resinfo[l];
268 atoms.resinfo[nrdx + l].nr += nrdx;
273 if (!read_next_x(oenv, status, &t, xx, boxx) && ((i + 1) * (j + 1) * (k + 1) < vol))
275 gmx_fatal(FARGS, "Not enough frames in trajectory");
286 /* make box bigger */
287 for (m = 0; (m < DIM); m++)
289 box[m][m] += dist[m];
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)
296 /* With an even number of boxes in x we can forgot about the screw */
297 pbcType = PbcType::Xyz;
300 /*depending on how you look at it, this is either a nasty hack or the way it should work*/
303 for (i = 0; i < atoms.nres; i++)
305 atoms.resinfo[i].nr = i + 1;
309 write_sto_conf(opt2fn("-o", NFILE, fnm), *mtop.name, &atoms, x, v, pbcType, box);
317 output_env_done(oenv);