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 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * 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.
50 static char *RotStr = {"Enforced rotation:"};
53 static char s_vec[STRLEN];
56 static void string2dvec(char buf[], dvec nums)
58 if (sscanf(buf,"%lf%lf%lf",&nums[0],&nums[1],&nums[2]) != 3)
59 gmx_fatal(FARGS,"Expected three numbers at input line %s",buf);
63 extern char **read_rotparams(int *ninp_p,t_inpfile **inp_p,t_rot *rot,
71 char warn_buf[STRLEN];
78 /* read rotation parameters */
79 CTYPE("Output frequency for angle, torque and rotation potential energy for the whole group");
80 ITYPE("rot_nstrout", rot->nstrout, 100);
81 CTYPE("Output frequency for per-slab data (angles, torques and slab centers)");
82 ITYPE("rot_nstsout", rot->nstsout, 1000);
83 CTYPE("Number of rotation groups");
84 ITYPE("rot_ngroups", rot->ngrp,1);
88 gmx_fatal(FARGS,"rot_ngroups should be >= 1");
91 snew(rot->grp,rot->ngrp);
93 /* Read the rotation groups */
94 snew(grpbuf,rot->ngrp);
95 for(g=0; g<rot->ngrp; g++)
98 snew(grpbuf[g],STRLEN);
99 CTYPE("Rotation group name");
100 sprintf(buf,"rot_group%d",g);
101 STYPE(buf, grpbuf[g], "");
103 CTYPE("Rotation potential. Can be iso, iso-pf, pm, pm-pf, rm, rm-pf, rm2, rm2-pf, flex, flex-t, flex2, flex2-t");
104 sprintf(buf,"rot_type%d",g);
105 ETYPE(buf, rotg->eType, erotg_names);
107 CTYPE("Use mass-weighting of the rotation group positions");
108 sprintf(buf,"rot_massw%d",g);
109 ETYPE(buf, rotg->bMassW, yesno_names);
111 CTYPE("Rotation vector, will get normalized");
112 sprintf(buf,"rot_vec%d",g);
113 STYPE(buf, s_vec, "1.0 0.0 0.0");
114 string2dvec(s_vec,vec);
115 /* Normalize the rotation vector */
118 dsvmul(1.0/dnorm(vec),vec,vec);
122 sprintf(warn_buf, "rot_vec%d = 0", g);
123 warning_error(wi, warn_buf);
125 fprintf(stderr, "%s Group %d (%s) normalized rot. vector: %f %f %f\n",
126 RotStr, g, erotg_names[rotg->eType], vec[0], vec[1], vec[2]);
128 rotg->vec[m] = vec[m];
130 CTYPE("Pivot point for the potentials iso, pm, rm, and rm2 (nm)");
131 sprintf(buf,"rot_pivot%d",g);
132 STYPE(buf, s_vec, "0.0 0.0 0.0");
134 if ( (rotg->eType==erotgISO) || (rotg->eType==erotgPM) || (rotg->eType==erotgRM) || (rotg->eType==erotgRM2) )
135 string2dvec(s_vec,vec);
137 rotg->pivot[m] = vec[m];
139 CTYPE("Rotation rate (degree/ps) and force constant (kJ/(mol*nm^2))");
140 sprintf(buf,"rot_rate%d",g);
141 RTYPE(buf, rotg->rate, 0.0);
143 sprintf(buf,"rot_k%d",g);
144 RTYPE(buf, rotg->k, 0.0);
147 sprintf(warn_buf, "rot_k%d <= 0", g);
148 warning_note(wi, warn_buf);
151 CTYPE("Slab distance for flexible axis rotation (nm)");
152 sprintf(buf,"rot_slab_dist%d",g);
153 RTYPE(buf, rotg->slab_dist, 1.5);
154 if (rotg->slab_dist <= 0.0)
156 sprintf(warn_buf, "rot_slab_dist%d <= 0", g);
157 warning_error(wi, warn_buf);
160 CTYPE("Minimum value of Gaussian function for the force to be evaluated (for flex* potentials)");
161 sprintf(buf,"rot_min_gauss%d",g);
162 RTYPE(buf, rotg->min_gaussian, 1e-3);
163 if (rotg->min_gaussian <= 0.0)
165 sprintf(warn_buf, "rot_min_gauss%d <= 0", g);
166 warning_error(wi, warn_buf);
169 CTYPE("Value of additive constant epsilon' (nm^2) for rm2* and flex2* potentials");
170 sprintf(buf, "rot_eps%d",g);
171 RTYPE(buf, rotg->eps, 1e-4);
172 if ( (rotg->eps <= 0.0) && (rotg->eType==erotgRM2 || rotg->eType==erotgFLEX2) )
174 sprintf(warn_buf, "rot_eps%d <= 0", g);
175 warning_error(wi, warn_buf);
178 CTYPE("Fitting method to determine angle of rotation group (rmsd, norm, or potential)");
179 sprintf(buf,"rot_fit_method%d",g);
180 ETYPE(buf, rotg->eFittype, erotg_fitnames);
181 CTYPE("For fit type 'potential', nr. of angles around the reference for which the pot. is evaluated");
182 sprintf(buf,"rot_potfit_nsteps%d",g);
183 ITYPE(buf, rotg->PotAngle_nstep, 21);
184 if ( (rotg->eFittype==erotgFitPOT) && (rotg->PotAngle_nstep < 1) )
186 sprintf(warn_buf, "rot_potfit_nsteps%d < 1", g);
187 warning_error(wi, warn_buf);
189 CTYPE("For fit type 'potential', distance in degrees between two consecutive angles");
190 sprintf(buf,"rot_potfit_step%d",g);
191 RTYPE(buf, rotg->PotAngle_step, 0.25);
201 /* Check whether the box is unchanged */
202 static void check_box_unchanged(matrix f_box, matrix box, char fn[], warninp_t wi)
206 char warn_buf[STRLEN];
209 for (i=0; i<DIM; i++)
210 for (ii=0; ii<DIM; ii++)
211 if (f_box[i][ii] != box[i][ii])
215 sprintf(warn_buf, "%s Box size in reference file %s differs from actual box size!",
217 warning(wi, warn_buf);
218 pr_rvecs(stderr,0,"Your box is:",box ,3);
219 pr_rvecs(stderr,0,"Box in file:",f_box,3);
224 /* Extract the reference positions for the rotation group(s) */
225 extern void set_reference_positions(
226 t_rot *rot, gmx_mtop_t *mtop, rvec *x, matrix box,
227 const char *fn, gmx_bool bSet, warninp_t wi)
231 t_trnheader header; /* Header information of reference file */
232 char base[STRLEN],extension[STRLEN],reffile[STRLEN];
234 rvec f_box[3]; /* Box from reference file */
237 /* Base name and extension of the reference file: */
238 strncpy(base, fn, STRLEN - 1);
240 extpos = strrchr(base, '.');
241 strcpy(extension,extpos+1);
245 for (g=0; g<rot->ngrp; g++)
248 fprintf(stderr, "%s group %d has %d reference positions.\n",RotStr,g,rotg->nat);
249 snew(rotg->x_ref, rotg->nat);
251 /* Construct the name for the file containing the reference positions for this group: */
252 sprintf(reffile, "%s.%d.%s", base,g,extension);
254 /* If the base filename for the reference position files was explicitly set by
255 * the user, we issue a fatal error if the group file can not be found */
256 if (bSet && !gmx_fexist(reffile))
258 gmx_fatal(FARGS, "%s The file containing the reference positions was not found.\n"
259 "Expected the file '%s' for group %d.\n",
263 if (gmx_fexist(reffile))
265 fprintf(stderr, " Reading them from %s.\n", reffile);
266 read_trnheader(reffile, &header);
267 if (rotg->nat != header.natoms)
268 gmx_fatal(FARGS,"Number of atoms in file %s (%d) does not match the number of atoms in rotation group (%d)!\n",
269 reffile, header.natoms, rotg->nat);
270 read_trn(reffile, &header.step, &header.t, &header.lambda, f_box, &header.natoms, rotg->x_ref, NULL, NULL);
272 /* Check whether the box is unchanged and output a warning if not: */
273 check_box_unchanged(f_box,box,reffile,wi);
277 fprintf(stderr, " Saving them to %s.\n", reffile);
278 for(i=0; i<rotg->nat; i++)
281 copy_rvec(x[ii], rotg->x_ref[i]);
283 write_trn(reffile,g,0.0,0.0,box,rotg->nat,rotg->x_ref,NULL,NULL);
289 extern void make_rotation_groups(t_rot *rot,char **rotgnames,t_blocka *grps,char **gnames)
295 for (g=0; g<rot->ngrp; g++)
298 ig = search_string(rotgnames[g],grps->nr,gnames);
299 rotg->nat = grps->index[ig+1] - grps->index[ig];
303 fprintf(stderr,"Rotation group %d '%s' has %d atoms\n",g,rotgnames[g],rotg->nat);
304 snew(rotg->ind,rotg->nat);
305 for(i=0; i<rotg->nat; i++)
306 rotg->ind[i] = grps->a[grps->index[ig]+i];
309 gmx_fatal(FARGS,"Rotation group %d '%s' is empty",g,rotgnames[g]);