2 * This source code is part of
6 * GROningen MAchine for Chemical Simulations
8 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
9 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
10 * Copyright (c) 2001-2004, The GROMACS development team,
11 * check out http://www.gromacs.org for more information.
13 * This program is free software; you can redistribute it and/or
14 * modify it under the terms of the GNU General Public License
15 * as published by the Free Software Foundation; either version 2
16 * of the License, or (at your option) any later version.
18 * If you want to redistribute modifications, please consider that
19 * scientific software is very special. Version control is crucial -
20 * bugs must be traceable. We will be happy to consider code for
21 * inclusion in the official distribution, but derived work must not
22 * be called official GROMACS. Details are found in the README & COPYING
23 * files - if they are missing, get the official version at www.gromacs.org.
25 * To help us fund GROMACS development, we humbly ask that you cite
26 * the papers on the package - you can find them in the top README file.
28 * For more info, check our website at http://www.gromacs.org
31 * GROwing Monsters And Cloning Shrimps
41 #include "gromacs/fileio/futil.h"
42 #include "gromacs/fileio/trnio.h"
45 static char *RotStr = {"Enforced rotation:"};
48 static char s_vec[STRLEN];
51 static void string2dvec(char buf[], dvec nums)
53 if (sscanf(buf, "%lf%lf%lf", &nums[0], &nums[1], &nums[2]) != 3)
55 gmx_fatal(FARGS, "Expected three numbers at input line %s", buf);
60 extern char **read_rotparams(int *ninp_p, t_inpfile **inp_p, t_rot *rot,
68 char warn_buf[STRLEN];
75 /* read rotation parameters */
76 CTYPE("Output frequency for angle, torque and rotation potential energy for the whole group");
77 ITYPE("rot_nstrout", rot->nstrout, 100);
78 CTYPE("Output frequency for per-slab data (angles, torques and slab centers)");
79 ITYPE("rot_nstsout", rot->nstsout, 1000);
80 CTYPE("Number of rotation groups");
81 ITYPE("rot_ngroups", rot->ngrp, 1);
85 gmx_fatal(FARGS, "rot_ngroups should be >= 1");
88 snew(rot->grp, rot->ngrp);
90 /* Read the rotation groups */
91 snew(grpbuf, rot->ngrp);
92 for (g = 0; g < rot->ngrp; g++)
95 snew(grpbuf[g], STRLEN);
96 CTYPE("Rotation group name");
97 sprintf(buf, "rot_group%d", g);
98 STYPE(buf, grpbuf[g], "");
100 CTYPE("Rotation potential. Can be iso, iso-pf, pm, pm-pf, rm, rm-pf, rm2, rm2-pf, flex, flex-t, flex2, flex2-t");
101 sprintf(buf, "rot_type%d", g);
102 ETYPE(buf, rotg->eType, erotg_names);
104 CTYPE("Use mass-weighting of the rotation group positions");
105 sprintf(buf, "rot_massw%d", g);
106 ETYPE(buf, rotg->bMassW, yesno_names);
108 CTYPE("Rotation vector, will get normalized");
109 sprintf(buf, "rot_vec%d", g);
110 STYPE(buf, s_vec, "1.0 0.0 0.0");
111 string2dvec(s_vec, vec);
112 /* Normalize the rotation vector */
115 dsvmul(1.0/dnorm(vec), vec, vec);
119 sprintf(warn_buf, "rot_vec%d = 0", g);
120 warning_error(wi, warn_buf);
122 fprintf(stderr, "%s Group %d (%s) normalized rot. vector: %f %f %f\n",
123 RotStr, g, erotg_names[rotg->eType], vec[0], vec[1], vec[2]);
124 for (m = 0; m < DIM; m++)
126 rotg->vec[m] = vec[m];
129 CTYPE("Pivot point for the potentials iso, pm, rm, and rm2 (nm)");
130 sprintf(buf, "rot_pivot%d", g);
131 STYPE(buf, s_vec, "0.0 0.0 0.0");
133 if ( (rotg->eType == erotgISO) || (rotg->eType == erotgPM) || (rotg->eType == erotgRM) || (rotg->eType == erotgRM2) )
135 string2dvec(s_vec, vec);
137 for (m = 0; m < DIM; m++)
139 rotg->pivot[m] = vec[m];
142 CTYPE("Rotation rate (degree/ps) and force constant (kJ/(mol*nm^2))");
143 sprintf(buf, "rot_rate%d", g);
144 RTYPE(buf, rotg->rate, 0.0);
146 sprintf(buf, "rot_k%d", g);
147 RTYPE(buf, rotg->k, 0.0);
150 sprintf(warn_buf, "rot_k%d <= 0", g);
151 warning_note(wi, warn_buf);
154 CTYPE("Slab distance for flexible axis rotation (nm)");
155 sprintf(buf, "rot_slab_dist%d", g);
156 RTYPE(buf, rotg->slab_dist, 1.5);
157 if (rotg->slab_dist <= 0.0)
159 sprintf(warn_buf, "rot_slab_dist%d <= 0", g);
160 warning_error(wi, warn_buf);
163 CTYPE("Minimum value of Gaussian function for the force to be evaluated (for flex* potentials)");
164 sprintf(buf, "rot_min_gauss%d", g);
165 RTYPE(buf, rotg->min_gaussian, 1e-3);
166 if (rotg->min_gaussian <= 0.0)
168 sprintf(warn_buf, "rot_min_gauss%d <= 0", g);
169 warning_error(wi, warn_buf);
172 CTYPE("Value of additive constant epsilon' (nm^2) for rm2* and flex2* potentials");
173 sprintf(buf, "rot_eps%d", g);
174 RTYPE(buf, rotg->eps, 1e-4);
175 if ( (rotg->eps <= 0.0) && (rotg->eType == erotgRM2 || rotg->eType == erotgFLEX2) )
177 sprintf(warn_buf, "rot_eps%d <= 0", g);
178 warning_error(wi, warn_buf);
181 CTYPE("Fitting method to determine angle of rotation group (rmsd, norm, or potential)");
182 sprintf(buf, "rot_fit_method%d", g);
183 ETYPE(buf, rotg->eFittype, erotg_fitnames);
184 CTYPE("For fit type 'potential', nr. of angles around the reference for which the pot. is evaluated");
185 sprintf(buf, "rot_potfit_nsteps%d", g);
186 ITYPE(buf, rotg->PotAngle_nstep, 21);
187 if ( (rotg->eFittype == erotgFitPOT) && (rotg->PotAngle_nstep < 1) )
189 sprintf(warn_buf, "rot_potfit_nsteps%d < 1", g);
190 warning_error(wi, warn_buf);
192 CTYPE("For fit type 'potential', distance in degrees between two consecutive angles");
193 sprintf(buf, "rot_potfit_step%d", g);
194 RTYPE(buf, rotg->PotAngle_step, 0.25);
204 /* Check whether the box is unchanged */
205 static void check_box_unchanged(matrix f_box, matrix box, char fn[], warninp_t wi)
208 gmx_bool bSame = TRUE;
209 char warn_buf[STRLEN];
212 for (i = 0; i < DIM; i++)
214 for (ii = 0; ii < DIM; ii++)
216 if (f_box[i][ii] != box[i][ii])
224 sprintf(warn_buf, "%s Box size in reference file %s differs from actual box size!",
226 warning(wi, warn_buf);
227 pr_rvecs(stderr, 0, "Your box is:", box, 3);
228 pr_rvecs(stderr, 0, "Box in file:", f_box, 3);
233 /* Extract the reference positions for the rotation group(s) */
234 extern void set_reference_positions(
235 t_rot *rot, rvec *x, matrix box,
236 const char *fn, gmx_bool bSet, warninp_t wi)
240 t_trnheader header; /* Header information of reference file */
241 char base[STRLEN], extension[STRLEN], reffile[STRLEN];
243 rvec f_box[3]; /* Box from reference file */
246 /* Base name and extension of the reference file: */
247 strncpy(base, fn, STRLEN - 1);
248 base[STRLEN-1] = '\0';
249 extpos = strrchr(base, '.');
250 strcpy(extension, extpos+1);
254 for (g = 0; g < rot->ngrp; g++)
257 fprintf(stderr, "%s group %d has %d reference positions.\n", RotStr, g, rotg->nat);
258 snew(rotg->x_ref, rotg->nat);
260 /* Construct the name for the file containing the reference positions for this group: */
261 sprintf(reffile, "%s.%d.%s", base, g, extension);
263 /* If the base filename for the reference position files was explicitly set by
264 * the user, we issue a fatal error if the group file can not be found */
265 if (bSet && !gmx_fexist(reffile))
267 gmx_fatal(FARGS, "%s The file containing the reference positions was not found.\n"
268 "Expected the file '%s' for group %d.\n",
272 if (gmx_fexist(reffile))
274 fprintf(stderr, " Reading them from %s.\n", reffile);
275 read_trnheader(reffile, &header);
276 if (rotg->nat != header.natoms)
278 gmx_fatal(FARGS, "Number of atoms in file %s (%d) does not match the number of atoms in rotation group (%d)!\n",
279 reffile, header.natoms, rotg->nat);
281 read_trn(reffile, &header.step, &header.t, &header.lambda, f_box, &header.natoms, rotg->x_ref, NULL, NULL);
283 /* Check whether the box is unchanged and output a warning if not: */
284 check_box_unchanged(f_box, box, reffile, wi);
288 fprintf(stderr, " Saving them to %s.\n", reffile);
289 for (i = 0; i < rotg->nat; i++)
292 copy_rvec(x[ii], rotg->x_ref[i]);
294 write_trn(reffile, g, 0.0, 0.0, box, rotg->nat, rotg->x_ref, NULL, NULL);
300 extern void make_rotation_groups(t_rot *rot, char **rotgnames, t_blocka *grps, char **gnames)
306 for (g = 0; g < rot->ngrp; g++)
309 ig = search_string(rotgnames[g], grps->nr, gnames);
310 rotg->nat = grps->index[ig+1] - grps->index[ig];
314 fprintf(stderr, "Rotation group %d '%s' has %d atoms\n", g, rotgnames[g], rotg->nat);
315 snew(rotg->ind, rotg->nat);
316 for (i = 0; i < rotg->nat; i++)
318 rotg->ind[i] = grps->a[grps->index[ig]+i];
323 gmx_fatal(FARGS, "Rotation group %d '%s' is empty", g, rotgnames[g]);