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.
39 #include "gromacs/fileio/readinp.h"
40 #include "gromacs/fileio/trrio.h"
41 #include "gromacs/fileio/warninp.h"
42 #include "gromacs/gmxpreprocess/readir.h"
43 #include "gromacs/math/vec.h"
44 #include "gromacs/math/vecdump.h"
45 #include "gromacs/mdtypes/inputrec.h"
46 #include "gromacs/mdtypes/md_enums.h"
47 #include "gromacs/topology/block.h"
48 #include "gromacs/utility/cstringutil.h"
49 #include "gromacs/utility/fatalerror.h"
50 #include "gromacs/utility/futil.h"
51 #include "gromacs/utility/smalloc.h"
53 static const char *RotStr = "Enforced rotation:";
56 static char s_vec[STRLEN];
59 static void string2dvec(char buf[], dvec nums)
61 if (sscanf(buf, "%lf%lf%lf", &nums[0], &nums[1], &nums[2]) != 3)
63 gmx_fatal(FARGS, "Expected three numbers at input line %s", buf);
68 extern char **read_rotparams(int *ninp_p, t_inpfile **inp_p, t_rot *rot,
76 char warn_buf[STRLEN];
83 /* read rotation parameters */
84 CTYPE("Output frequency for angle, torque and rotation potential energy for the whole group");
85 ITYPE("rot-nstrout", rot->nstrout, 100);
86 CTYPE("Output frequency for per-slab data (angles, torques and slab centers)");
87 ITYPE("rot-nstsout", rot->nstsout, 1000);
88 CTYPE("Number of rotation groups");
89 ITYPE("rot-ngroups", rot->ngrp, 1);
93 gmx_fatal(FARGS, "rot-ngroups should be >= 1");
96 snew(rot->grp, rot->ngrp);
98 /* Read the rotation groups */
99 snew(grpbuf, rot->ngrp);
100 for (g = 0; g < rot->ngrp; g++)
103 snew(grpbuf[g], STRLEN);
104 CTYPE("Rotation group name");
105 sprintf(buf, "rot-group%d", g);
106 STYPE(buf, grpbuf[g], "");
108 CTYPE("Rotation potential. Can be iso, iso-pf, pm, pm-pf, rm, rm-pf, rm2, rm2-pf, flex, flex-t, flex2, flex2-t");
109 sprintf(buf, "rot-type%d", g);
110 ETYPE(buf, rotg->eType, erotg_names);
112 CTYPE("Use mass-weighting of the rotation group positions");
113 sprintf(buf, "rot-massw%d", g);
114 ETYPE(buf, rotg->bMassW, yesno_names);
116 CTYPE("Rotation vector, will get normalized");
117 sprintf(buf, "rot-vec%d", g);
118 STYPE(buf, s_vec, "1.0 0.0 0.0");
119 string2dvec(s_vec, vec);
120 /* Normalize the rotation vector */
123 dsvmul(1.0/dnorm(vec), vec, vec);
127 sprintf(warn_buf, "rot-vec%d = 0", g);
128 warning_error(wi, warn_buf);
130 fprintf(stderr, "%s Group %d (%s) normalized rot. vector: %f %f %f\n",
131 RotStr, g, erotg_names[rotg->eType], vec[0], vec[1], vec[2]);
132 for (m = 0; m < DIM; m++)
134 rotg->vec[m] = vec[m];
137 CTYPE("Pivot point for the potentials iso, pm, rm, and rm2 (nm)");
138 sprintf(buf, "rot-pivot%d", g);
139 STYPE(buf, s_vec, "0.0 0.0 0.0");
141 if ( (rotg->eType == erotgISO) || (rotg->eType == erotgPM) || (rotg->eType == erotgRM) || (rotg->eType == erotgRM2) )
143 string2dvec(s_vec, vec);
145 for (m = 0; m < DIM; m++)
147 rotg->pivot[m] = vec[m];
150 CTYPE("Rotation rate (degree/ps) and force constant (kJ/(mol*nm^2))");
151 sprintf(buf, "rot-rate%d", g);
152 RTYPE(buf, rotg->rate, 0.0);
154 sprintf(buf, "rot-k%d", g);
155 RTYPE(buf, rotg->k, 0.0);
158 sprintf(warn_buf, "rot-k%d <= 0", g);
159 warning_note(wi, warn_buf);
162 CTYPE("Slab distance for flexible axis rotation (nm)");
163 sprintf(buf, "rot-slab-dist%d", g);
164 RTYPE(buf, rotg->slab_dist, 1.5);
165 if (rotg->slab_dist <= 0.0)
167 sprintf(warn_buf, "rot-slab-dist%d <= 0", g);
168 warning_error(wi, warn_buf);
171 CTYPE("Minimum value of Gaussian function for the force to be evaluated (for flex* potentials)");
172 sprintf(buf, "rot-min-gauss%d", g);
173 RTYPE(buf, rotg->min_gaussian, 1e-3);
174 if (rotg->min_gaussian <= 0.0)
176 sprintf(warn_buf, "rot-min-gauss%d <= 0", g);
177 warning_error(wi, warn_buf);
180 CTYPE("Value of additive constant epsilon' (nm^2) for rm2* and flex2* potentials");
181 sprintf(buf, "rot-eps%d", g);
182 RTYPE(buf, rotg->eps, 1e-4);
183 if ( (rotg->eps <= 0.0) && (rotg->eType == erotgRM2 || rotg->eType == erotgFLEX2) )
185 sprintf(warn_buf, "rot-eps%d <= 0", g);
186 warning_error(wi, warn_buf);
189 CTYPE("Fitting method to determine angle of rotation group (rmsd, norm, or potential)");
190 sprintf(buf, "rot-fit-method%d", g);
191 ETYPE(buf, rotg->eFittype, erotg_fitnames);
192 CTYPE("For fit type 'potential', nr. of angles around the reference for which the pot. is evaluated");
193 sprintf(buf, "rot-potfit-nsteps%d", g);
194 ITYPE(buf, rotg->PotAngle_nstep, 21);
195 if ( (rotg->eFittype == erotgFitPOT) && (rotg->PotAngle_nstep < 1) )
197 sprintf(warn_buf, "rot-potfit-nsteps%d < 1", g);
198 warning_error(wi, warn_buf);
200 CTYPE("For fit type 'potential', distance in degrees between two consecutive angles");
201 sprintf(buf, "rot-potfit-step%d", g);
202 RTYPE(buf, rotg->PotAngle_step, 0.25);
212 /* Check whether the box is unchanged */
213 static void check_box_unchanged(matrix f_box, matrix box, char fn[], warninp_t wi)
216 gmx_bool bSame = TRUE;
217 char warn_buf[STRLEN];
220 for (i = 0; i < DIM; i++)
222 for (ii = 0; ii < DIM; ii++)
224 if (f_box[i][ii] != box[i][ii])
232 sprintf(warn_buf, "%s Box size in reference file %s differs from actual box size!",
234 warning(wi, warn_buf);
235 pr_rvecs(stderr, 0, "Your box is:", box, 3);
236 pr_rvecs(stderr, 0, "Box in file:", f_box, 3);
241 /* Extract the reference positions for the rotation group(s) */
242 extern void set_reference_positions(
243 t_rot *rot, rvec *x, matrix box,
244 const char *fn, gmx_bool bSet, warninp_t wi)
248 gmx_trr_header_t header; /* Header information of reference file */
249 char base[STRLEN], extension[STRLEN], reffile[STRLEN];
251 rvec f_box[3]; /* Box from reference file */
254 /* Base name and extension of the reference file: */
255 strncpy(base, fn, STRLEN - 1);
256 base[STRLEN-1] = '\0';
257 extpos = strrchr(base, '.');
258 strcpy(extension, extpos+1);
262 for (g = 0; g < rot->ngrp; g++)
265 fprintf(stderr, "%s group %d has %d reference positions.\n", RotStr, g, rotg->nat);
266 snew(rotg->x_ref, rotg->nat);
268 /* Construct the name for the file containing the reference positions for this group: */
269 sprintf(reffile, "%s.%d.%s", base, g, extension);
271 /* If the base filename for the reference position files was explicitly set by
272 * the user, we issue a fatal error if the group file can not be found */
273 if (bSet && !gmx_fexist(reffile))
275 gmx_fatal(FARGS, "%s The file containing the reference positions was not found.\n"
276 "Expected the file '%s' for group %d.\n",
280 if (gmx_fexist(reffile))
282 fprintf(stderr, " Reading them from %s.\n", reffile);
283 gmx_trr_read_single_header(reffile, &header);
284 if (rotg->nat != header.natoms)
286 gmx_fatal(FARGS, "Number of atoms in file %s (%d) does not match the number of atoms in rotation group (%d)!\n",
287 reffile, header.natoms, rotg->nat);
289 gmx_trr_read_single_frame(reffile, &header.step, &header.t, &header.lambda, f_box, &header.natoms, rotg->x_ref, nullptr, nullptr);
291 /* Check whether the box is unchanged and output a warning if not: */
292 check_box_unchanged(f_box, box, reffile, wi);
296 fprintf(stderr, " Saving them to %s.\n", reffile);
297 for (i = 0; i < rotg->nat; i++)
300 copy_rvec(x[ii], rotg->x_ref[i]);
302 gmx_trr_write_single_frame(reffile, g, 0.0, 0.0, box, rotg->nat, rotg->x_ref, nullptr, nullptr);
308 extern void make_rotation_groups(t_rot *rot, char **rotgnames, t_blocka *grps, char **gnames)
314 for (g = 0; g < rot->ngrp; g++)
317 ig = search_string(rotgnames[g], grps->nr, gnames);
318 rotg->nat = grps->index[ig+1] - grps->index[ig];
322 fprintf(stderr, "Rotation group %d '%s' has %d atoms\n", g, rotgnames[g], rotg->nat);
323 snew(rotg->ind, rotg->nat);
324 for (i = 0; i < rotg->nat; i++)
326 rotg->ind[i] = grps->a[grps->index[ig]+i];
331 gmx_fatal(FARGS, "Rotation group %d '%s' is empty", g, rotgnames[g]);