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,2018, 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
,
74 char warn_buf
[STRLEN
];
78 /* read rotation parameters */
79 printStringNoNewline(ninp_p
, inp_p
, "Output frequency for angle, torque and rotation potential energy for the whole group");
80 rot
->nstrout
= get_eint(ninp_p
, inp_p
, "rot-nstrout", 100, wi
);
81 printStringNoNewline(ninp_p
, inp_p
, "Output frequency for per-slab data (angles, torques and slab centers)");
82 rot
->nstsout
= get_eint(ninp_p
, inp_p
, "rot-nstsout", 1000, wi
);
83 printStringNoNewline(ninp_p
, inp_p
, "Number of rotation groups");
84 rot
->ngrp
= get_eint(ninp_p
, inp_p
, "rot-ngroups", 1, wi
);
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 printStringNoNewline(ninp_p
, inp_p
, "Rotation group name");
100 sprintf(buf
, "rot-group%d", g
);
101 setStringEntry(ninp_p
, inp_p
, buf
, grpbuf
[g
], "");
103 printStringNoNewline(ninp_p
, inp_p
, "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 rotg
->eType
= get_eenum(ninp_p
, inp_p
, buf
, erotg_names
);
107 printStringNoNewline(ninp_p
, inp_p
, "Use mass-weighting of the rotation group positions");
108 sprintf(buf
, "rot-massw%d", g
);
109 rotg
->bMassW
= get_eenum(ninp_p
, inp_p
, buf
, yesno_names
);
111 printStringNoNewline(ninp_p
, inp_p
, "Rotation vector, will get normalized");
112 sprintf(buf
, "rot-vec%d", g
);
113 setStringEntry(ninp_p
, inp_p
, 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]);
127 for (m
= 0; m
< DIM
; m
++)
129 rotg
->vec
[m
] = vec
[m
];
132 printStringNoNewline(ninp_p
, inp_p
, "Pivot point for the potentials iso, pm, rm, and rm2 (nm)");
133 sprintf(buf
, "rot-pivot%d", g
);
134 setStringEntry(ninp_p
, inp_p
, buf
, s_vec
, "0.0 0.0 0.0");
136 if ( (rotg
->eType
== erotgISO
) || (rotg
->eType
== erotgPM
) || (rotg
->eType
== erotgRM
) || (rotg
->eType
== erotgRM2
) )
138 string2dvec(s_vec
, vec
);
140 for (m
= 0; m
< DIM
; m
++)
142 rotg
->pivot
[m
] = vec
[m
];
145 printStringNoNewline(ninp_p
, inp_p
, "Rotation rate (degree/ps) and force constant (kJ/(mol*nm^2))");
146 sprintf(buf
, "rot-rate%d", g
);
147 rotg
->rate
= get_ereal(ninp_p
, inp_p
, buf
, 0.0, wi
);
149 sprintf(buf
, "rot-k%d", g
);
150 rotg
->k
= get_ereal(ninp_p
, inp_p
, buf
, 0.0, wi
);
153 sprintf(warn_buf
, "rot-k%d <= 0", g
);
154 warning_note(wi
, warn_buf
);
157 printStringNoNewline(ninp_p
, inp_p
, "Slab distance for flexible axis rotation (nm)");
158 sprintf(buf
, "rot-slab-dist%d", g
);
159 rotg
->slab_dist
= get_ereal(ninp_p
, inp_p
, buf
, 1.5, wi
);
160 if (rotg
->slab_dist
<= 0.0)
162 sprintf(warn_buf
, "rot-slab-dist%d <= 0", g
);
163 warning_error(wi
, warn_buf
);
166 printStringNoNewline(ninp_p
, inp_p
, "Minimum value of Gaussian function for the force to be evaluated (for flex* potentials)");
167 sprintf(buf
, "rot-min-gauss%d", g
);
168 rotg
->min_gaussian
= get_ereal(ninp_p
, inp_p
, buf
, 1e-3, wi
);
169 if (rotg
->min_gaussian
<= 0.0)
171 sprintf(warn_buf
, "rot-min-gauss%d <= 0", g
);
172 warning_error(wi
, warn_buf
);
175 printStringNoNewline(ninp_p
, inp_p
, "Value of additive constant epsilon' (nm^2) for rm2* and flex2* potentials");
176 sprintf(buf
, "rot-eps%d", g
);
177 rotg
->eps
= get_ereal(ninp_p
, inp_p
, buf
, 1e-4, wi
);
178 if ( (rotg
->eps
<= 0.0) && (rotg
->eType
== erotgRM2
|| rotg
->eType
== erotgFLEX2
) )
180 sprintf(warn_buf
, "rot-eps%d <= 0", g
);
181 warning_error(wi
, warn_buf
);
184 printStringNoNewline(ninp_p
, inp_p
, "Fitting method to determine angle of rotation group (rmsd, norm, or potential)");
185 sprintf(buf
, "rot-fit-method%d", g
);
186 rotg
->eFittype
= get_eenum(ninp_p
, inp_p
, buf
, erotg_fitnames
);
187 printStringNoNewline(ninp_p
, inp_p
, "For fit type 'potential', nr. of angles around the reference for which the pot. is evaluated");
188 sprintf(buf
, "rot-potfit-nsteps%d", g
);
189 rotg
->PotAngle_nstep
= get_eint(ninp_p
, inp_p
, buf
, 21, wi
);
190 if ( (rotg
->eFittype
== erotgFitPOT
) && (rotg
->PotAngle_nstep
< 1) )
192 sprintf(warn_buf
, "rot-potfit-nsteps%d < 1", g
);
193 warning_error(wi
, warn_buf
);
195 printStringNoNewline(ninp_p
, inp_p
, "For fit type 'potential', distance in degrees between two consecutive angles");
196 sprintf(buf
, "rot-potfit-step%d", g
);
197 rotg
->PotAngle_step
= get_ereal(ninp_p
, inp_p
, buf
, 0.25, wi
);
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 gmx_trr_header_t 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 gmx_trr_read_single_header(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 gmx_trr_read_single_frame(reffile
, &header
.step
, &header
.t
, &header
.lambda
, f_box
, &header
.natoms
, rotg
->x_ref
, nullptr, nullptr);
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 gmx_trr_write_single_frame(reffile
, g
, 0.0, 0.0, box
, rotg
->nat
, rotg
->x_ref
, nullptr, nullptr);
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
]);