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
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)
54 gmx_fatal(FARGS,"Expected three numbers at input line %s",buf);
58 extern char **read_rotparams(int *ninp_p,t_inpfile **inp_p,t_rot *rot,
66 char warn_buf[STRLEN];
73 /* read rotation parameters */
74 ITYPE("rot_nstrout", rot->nstrout, 100);
75 ITYPE("rot_nsttout", rot->nsttout, 1000);
76 CTYPE("Number of rotation groups");
77 ITYPE("rot_ngroups", rot->ngrp,1);
81 gmx_fatal(FARGS,"rot_ngroups should be >= 1");
84 snew(rot->grp,rot->ngrp);
86 /* Read the rotation groups */
87 snew(grpbuf,rot->ngrp);
88 for(g=0; g<rot->ngrp; g++)
91 snew(grpbuf[g],STRLEN);
92 CTYPE("Rotation group name");
93 sprintf(buf,"rot_group%d",g);
94 STYPE(buf, grpbuf[g], "");
96 CTYPE("Rotation potential. Can be iso, iso-pf, pm, pm-pf, rm, rm-pf, rm2, rm2-pf, flex, flex-t, flex2, flex2-t");
97 sprintf(buf,"rot_type%d",g);
98 ETYPE(buf, rotg->eType, erotg_names);
100 CTYPE("Use mass-weighting of the rotation group positions");
101 sprintf(buf,"rot_massw%d",g);
102 ETYPE(buf, rotg->bMassW, yesno_names);
104 CTYPE("Rotation vector, will get normalized");
105 sprintf(buf,"rot_vec%d",g);
106 STYPE(buf, s_vec, "1.0 0.0 0.0");
107 string2dvec(s_vec,vec);
108 /* Normalize the rotation vector */
111 dsvmul(1.0/dnorm(vec),vec,vec);
115 sprintf(warn_buf, "rot_vec%d = 0", g);
116 warning_error(wi, warn_buf);
118 fprintf(stderr, "%s Group %d (%s) normalized rot. vector: %f %f %f\n",
119 RotStr, g, erotg_names[rotg->eType], vec[0], vec[1], vec[2]);
121 rotg->vec[m] = vec[m];
123 CTYPE("Pivot point for iso, pm, rm, rm2 potential [nm]");
124 sprintf(buf,"rot_pivot%d",g);
125 STYPE(buf, s_vec, "0.0 0.0 0.0");
127 if ( (rotg->eType==erotgISO) || (rotg->eType==erotgPM) || (rotg->eType==erotgRM) || (rotg->eType==erotgRM2) )
128 string2dvec(s_vec,vec);
130 rotg->pivot[m] = vec[m];
132 CTYPE("Rotation rate [degree/ps] and force constant [kJ/(mol*nm^2)]");
133 sprintf(buf,"rot_rate%d",g);
134 RTYPE(buf, rotg->rate, 0.0);
136 sprintf(buf,"rot_k%d",g);
137 RTYPE(buf, rotg->k, 0.0);
140 sprintf(warn_buf, "rot_k%d <= 0", g);
141 warning_note(wi, warn_buf);
144 CTYPE("Slab distance for flexible rotation [nm] (flexible axis only)");
145 sprintf(buf,"rot_slab_dist%d",g);
146 RTYPE(buf, rotg->slab_dist, 1.5);
147 if (rotg->slab_dist <= 0.0)
149 sprintf(warn_buf, "rot_slab_dist%d <= 0", g);
150 warning_error(wi, warn_buf);
153 CTYPE("Minimum value of Gaussian for the force to be evaluated (for flex and flex2 potentials)");
154 sprintf(buf,"rot_min_gauss%d",g);
155 RTYPE(buf, rotg->min_gaussian, 1e-3);
156 if (rotg->min_gaussian <= 0.0)
158 sprintf(warn_buf, "rot_min_gauss%d <= 0", g);
159 warning_error(wi, warn_buf);
162 CTYPE("Value of additive constant epsilon' [nm^2] for rm2 and flex2 potentials");
163 sprintf(buf, "rot_eps%d",g);
164 RTYPE(buf, rotg->eps, 1e-4);
165 if ( (rotg->eps <= 0.0) && (rotg->eType==erotgRM2 || rotg->eType==erotgFLEX2) )
167 sprintf(warn_buf, "rot_eps%d <= 0", g);
168 warning_error(wi, warn_buf);
171 CTYPE("Fitting method to determine actual angle of rotation group (rmsd or norm) (flex and flex2 pot.)");
172 sprintf(buf,"rot_fit_method%d",g);
173 ETYPE(buf, rotg->eFittype, erotg_fitnames);
183 /* Check whether the box is unchanged */
184 static void check_box(matrix f_box, matrix box, char fn[], warninp_t wi)
188 char warn_buf[STRLEN];
191 for (i=0; i<DIM; i++)
192 for (ii=0; ii<DIM; ii++)
193 if (f_box[i][ii] != box[i][ii])
197 sprintf(warn_buf, "%s Box size in reference file %s differs from actual box size!",
199 warning(wi, warn_buf);
200 pr_rvecs(stderr,0,"Your box is:",box ,3);
201 pr_rvecs(stderr,0,"Box in file:",f_box,3);
206 /* Extract the reference positions for the rotation group(s) */
207 extern void set_reference_positions(
208 t_rot *rot, gmx_mtop_t *mtop, rvec *x, matrix box,
209 const char *fn, gmx_bool bSet, warninp_t wi)
213 t_trnheader header; /* Header information of reference file */
214 char base[STRLEN],extension[STRLEN],reffile[STRLEN];
216 rvec f_box[3]; /* Box from reference file */
219 /* Base name and extension of the reference file: */
220 strncpy(base, fn, STRLEN - 1);
221 extpos = strrchr(base, '.');
222 strcpy(extension,extpos+1);
226 for (g=0; g<rot->ngrp; g++)
229 fprintf(stderr, "%s group %d has %d reference positions.\n",RotStr,g,rotg->nat);
230 snew(rotg->x_ref, rotg->nat);
232 /* Construct the name for the file containing the reference positions for this group: */
233 sprintf(reffile, "%s.%d.%s", base,g,extension);
235 /* If the base filename for the reference position files was explicitly set by
236 * the user, we issue a fatal error if the group file can not be found */
237 if (bSet && !gmx_fexist(reffile))
239 gmx_fatal(FARGS, "%s The file containing the reference positions was not found.\n"
240 "Expected the file '%s' for group %d.\n",
244 if (gmx_fexist(reffile))
246 fprintf(stderr, " Reading them from %s.\n", reffile);
247 read_trnheader(reffile, &header);
248 if (rotg->nat != header.natoms)
249 gmx_fatal(FARGS,"Number of atoms in file %s (%d) does not match the number of atoms in rotation group (%d)!\n",
250 reffile, header.natoms, rotg->nat);
251 read_trn(reffile, &header.step, &header.t, &header.lambda, f_box, &header.natoms, rotg->x_ref, NULL, NULL);
253 /* Check whether the box is unchanged and output a warning if not: */
254 check_box(f_box,box,reffile,wi);
258 fprintf(stderr, " Saving them to %s.\n", reffile);
259 for(i=0; i<rotg->nat; i++)
262 copy_rvec(x[ii], rotg->x_ref[i]);
264 write_trn(reffile,g,0.0,0.0,box,rotg->nat,rotg->x_ref,NULL,NULL);
270 extern void make_rotation_groups(t_rot *rot,char **rotgnames,t_blocka *grps,char **gnames)
276 for (g=0; g<rot->ngrp; g++)
279 ig = search_string(rotgnames[g],grps->nr,gnames);
280 rotg->nat = grps->index[ig+1] - grps->index[ig];
284 fprintf(stderr,"Rotation group %d '%s' has %d atoms\n",g,rotgnames[g],rotg->nat);
285 snew(rotg->ind,rotg->nat);
286 for(i=0; i<rotg->nat; i++)
287 rotg->ind[i] = grps->a[grps->index[ig]+i];
290 gmx_fatal(FARGS,"Rotation group %d '%s' is empty",g,rotgnames[g]);