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, 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.
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.
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.
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.
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.
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.
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"
57 static void rand_rot(int natoms
, rvec x
[], rvec v
[], vec4 xrot
[], vec4 vrot
[],
58 gmx::DefaultRandomEngine
* rng
, rvec max_rot
)
60 mat4 mt1
, mt2
, mr
[DIM
], mtemp1
, mtemp2
, mtemp3
, mxtot
, mvtot
;
64 gmx::UniformRealDistribution
<real
> dist(-1.0, 1.0);
67 for (i
= 0; (i
< natoms
); i
++)
69 for (m
= 0; (m
< DIM
); m
++)
71 xcm
[m
] += x
[i
][m
]/natoms
; /* get center of mass of one molecule */
74 fprintf(stderr
, "center of geometry: %f, %f, %f\n", xcm
[0], xcm
[1], xcm
[2]);
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
++)
80 phi
= M_PI
*max_rot
[m
]*dist(*rng
)/180;
81 gmx_mat4_init_rotation(m
, phi
, mr
[m
]);
83 gmx_mat4_init_translation(xcm
[XX
], xcm
[YY
], xcm
[ZZ
], mt2
);
85 /* For gmx_mat4_mmul() we need to multiply in the opposite order
86 * compared to normal mathematical notation.
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
);
94 for (i
= 0; (i
< natoms
); i
++)
96 gmx_mat4_transform_point(mxtot
, x
[i
], xrot
[i
]);
97 gmx_mat4_transform_point(mvtot
, v
[i
], vrot
[i
]);
101 int gmx_genconf(int argc
, char *argv
[])
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",
118 const char *bugs
[] = {
119 "The program should allow for random displacement of lattice points."
123 t_atoms
*atoms
; /* list with all atoms */
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
;
138 { efSTX
, "-f", "conf", ffREAD
},
139 { efSTO
, "-o", "out", ffWRITE
},
140 { efTRX
, "-trj", nullptr, ffOPTRD
}
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 */
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" }
159 if (!parse_common_args(&argc
, argv
, 0, NFILE
, fnm
, asize(pa
), pa
,
160 asize(desc
), desc
, asize(bugs
), bugs
, &oenv
))
167 seed
= static_cast<int>(gmx::makeRandomSeed());
169 gmx::DefaultRandomEngine
rng(seed
);
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);
176 if ((nx
<= 0) || (ny
<= 0) || (nz
<= 0))
178 gmx_fatal(FARGS
, "Number of boxes (-nbox) should be larger than zero");
181 vol
= nx
*ny
*nz
; /* calculate volume in grid points (= nr. molecules) */
186 read_tps_conf(opt2fn("-f", NFILE
, fnm
), top
, &ePBC
, &x
, &v
, box
, FALSE
);
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? */
198 if (!read_first_x(oenv
, &status
, ftp2fn(efTRX
, NFILE
, fnm
), &t
, &xx
, boxx
))
200 gmx_fatal(FARGS
, "No atoms in trajectory %s", ftp2fn(efTRX
, NFILE
, fnm
));
206 for (i
= 0; i
< natoms
; i
++)
208 copy_rvec(x
[i
], xx
[i
]);
213 for (k
= 0; (k
< nz
); k
++) /* loop over all gridpositions */
215 shift
[ZZ
] = k
*(dist
[ZZ
]+box
[ZZ
][ZZ
]);
217 for (j
= 0; (j
< ny
); j
++)
219 shift
[YY
] = j
*(dist
[YY
]+box
[YY
][YY
])+k
*box
[ZZ
][YY
];
221 for (i
= 0; (i
< nx
); i
++)
223 shift
[XX
] = i
*(dist
[XX
]+box
[XX
][XX
])+j
*box
[YY
][XX
]+k
*box
[ZZ
][XX
];
225 ndx
= (i
*ny
*nz
+j
*nz
+k
)*natoms
;
226 nrdx
= (i
*ny
*nz
+j
*nz
+k
)*nres
;
228 /* Random rotation on input coords */
231 rand_rot(natoms
, xx
, v
, xrot
, vrot
, &rng
, max_rot
);
234 for (l
= 0; (l
< natoms
); l
++)
236 for (m
= 0; (m
< DIM
); m
++)
240 x
[ndx
+l
][m
] = xrot
[l
][m
];
241 v
[ndx
+l
][m
] = vrot
[l
][m
];
245 x
[ndx
+l
][m
] = xx
[l
][m
];
246 v
[ndx
+l
][m
] = v
[l
][m
];
249 if (ePBC
== epbcSCREW
&& i
% 2 == 1)
251 /* Rotate around x axis */
252 for (m
= YY
; m
<= ZZ
; m
++)
254 x
[ndx
+l
][m
] = box
[YY
][m
] + box
[ZZ
][m
] - x
[ndx
+l
][m
];
255 v
[ndx
+l
][m
] = -v
[ndx
+l
][m
];
258 for (m
= 0; (m
< DIM
); m
++)
260 x
[ndx
+l
][m
] += shift
[m
];
262 atoms
->atom
[ndx
+l
].resind
= nrdx
+ atoms
->atom
[l
].resind
;
263 atoms
->atomname
[ndx
+l
] = atoms
->atomname
[l
];
266 for (l
= 0; (l
< nres
); l
++)
268 atoms
->resinfo
[nrdx
+l
] = atoms
->resinfo
[l
];
271 atoms
->resinfo
[nrdx
+l
].nr
+= nrdx
;
276 if (!read_next_x(oenv
, status
, &t
, xx
, boxx
) &&
277 ((i
+1)*(j
+1)*(k
+1) < vol
))
279 gmx_fatal(FARGS
, "Not enough frames in trajectory");
290 /* make box bigger */
291 for (m
= 0; (m
< DIM
); m
++)
293 box
[m
][m
] += dist
[m
];
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)
300 /* With an even number of boxes in x we can forgot about the screw */
304 /*depending on how you look at it, this is either a nasty hack or the way it should work*/
307 for (i
= 0; i
< atoms
->nres
; i
++)
309 atoms
->resinfo
[i
].nr
= i
+1;
313 write_sto_conf(opt2fn("-o", NFILE
, fnm
), *top
->name
, atoms
, x
, v
, ePBC
, box
);
322 done_filenms(NFILE
, fnm
);
323 output_env_done(oenv
);