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, 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.
41 #include "gromacs/legacyheaders/typedefs.h"
51 econqCoord, /* Constrain coordinates (mass weighted) */
52 econqVeloc, /* Constrain velocities (mass weighted) */
53 econqDeriv, /* Constrain a derivative (mass weighted), *
54 * for instance velocity or acceleration, *
55 * constraint virial can not be calculated. */
56 econqDeriv_FlexCon, /* As econqDeriv, but only output flex. con. */
57 econqForce, /* Constrain forces (non mass-weighted) */
58 econqForceDispl /* Constrain forces (mass-weighted 1/0 for freeze) */
61 int n_flexible_constraints(struct gmx_constr *constr);
62 /* Returns the total number of flexible constraints in the system */
64 void too_many_constraint_warnings(int eConstrAlg, int warncount);
65 /* Generate a fatal error because of too many LINCS/SETTLE warnings */
67 gmx_shakedata_t shake_init();
68 /* Initializes and return the SHAKE data structure */
70 gmx_bool bshakef(FILE *log, /* Log file */
71 gmx_shakedata_t shaked, /* Total number of atoms */
72 real invmass[], /* Atomic masses */
73 int nblocks, /* The number of shake blocks */
74 int sblock[], /* The shake blocks */
75 t_idef *idef, /* The interaction def */
76 t_inputrec *ir, /* Input record */
77 rvec x_s[], /* Coords before update */
78 rvec prime[], /* Output coords */
79 t_nrnb *nrnb, /* Performance measure */
80 real *lagr, /* The Lagrange multipliers */
81 real lambda, /* FEP lambda */
82 real *dvdlambda, /* FEP force */
83 real invdt, /* 1/delta_t */
84 rvec *v, /* Also constrain v if v!=NULL */
85 gmx_bool bCalcVir, /* Calculate r x m delta_r */
86 tensor vir_r_m_dr, /* sum r x m delta_r */
87 gmx_bool bDumpOnError, /* Dump debugging stuff on error*/
88 int econq, /* which type of constrainint is occurring */
89 t_vetavars *vetavar); /* veta for pressure control */
90 /* Shake all the atoms blockwise. It is assumed that all the constraints
91 * in the idef->shakes field are sorted, to ascending block nr. The
92 * sblock array points into the idef->shakes.iatoms field, with block 0
94 * at sblock[0] and running to ( < ) sblock[1], block n running from
95 * sblock[n] to sblock[n+1]. Array sblock should be large enough.
96 * Return TRUE when OK, FALSE when shake-error
99 gmx_settledata_t settle_init(real mO, real mH, real invmO, real invmH,
101 /* Initializes and returns a structure with SETTLE parameters */
103 void csettle(gmx_settledata_t settled,
104 int nsettle, /* Number of settles */
105 t_iatom iatoms[], /* The settle iatom list */
106 const struct t_pbc *pbc, /* PBC data pointer, can be NULL */
107 real b4[], /* Old coordinates */
108 real after[], /* New coords, to be settled */
109 real invdt, /* 1/delta_t */
110 real *v, /* Also constrain v if v!=NULL */
111 int calcvir_atom_end, /* Calculate r x m delta_r up to this atom */
112 tensor vir_r_m_dr, /* sum r x m delta_r */
114 t_vetavars *vetavar /* variables for pressure control */
117 void settle_proj(gmx_settledata_t settled, int econq,
118 int nsettle, t_iatom iatoms[],
119 const struct t_pbc *pbc, /* PBC data pointer, can be NULL */
121 rvec *der, rvec *derp,
122 int CalcVirAtomEnd, tensor vir_r_m_dder,
123 t_vetavars *vetavar);
124 /* Analytical algorithm to subtract the components of derivatives
125 * of coordinates working on settle type constraint.
128 void cshake(atom_id iatom[], int ncon, int *nnit, int maxnit,
129 real dist2[], real xp[], real rij[], real m2[], real omega,
130 real invmass[], real tt[], real lagr[], int *nerror);
131 /* Regular iterative shake */
133 void crattle(atom_id iatom[], int ncon, int *nnit, int maxnit,
134 real dist2[], real vp[], real rij[], real m2[], real omega,
135 real invmass[], real tt[], real lagr[], int *nerror, real invdt, t_vetavars *vetavar);
137 gmx_bool constrain(FILE *log, gmx_bool bLog, gmx_bool bEner,
141 gmx_ekindata_t *ekind,
143 gmx_int64_t step, int delta_step,
146 rvec *x, rvec *xprime, rvec *min_proj,
147 gmx_bool bMolPBC, matrix box,
148 real lambda, real *dvdlambda,
149 rvec *v, tensor *vir,
150 t_nrnb *nrnb, int econq, gmx_bool bPscal, real veta, real vetanew);
152 * When econq=econqCoord constrains coordinates xprime using th
153 * directions in x, min_proj is not used.
155 * When econq=econqDeriv, calculates the components xprime in
156 * the constraint directions and subtracts these components from min_proj.
157 * So when min_proj=xprime, the constraint components are projected out.
159 * When econq=econqDeriv_FlexCon, the same is done as with econqDeriv,
160 * but only the components of the flexible constraints are stored.
162 * When bMolPBC=TRUE, assume that molecules might be broken: correct PBC.
164 * delta_step is used for determining the constraint reference lengths
165 * when lenA != lenB or will the pull code with a pulling rate.
166 * step + delta_step is the step at which the final configuration
167 * is meant to be; for update delta_step = 1.
169 * step_scaling can be used to update coordinates based on the time
170 * step multiplied by this factor. Thus, normally 1.0 is passed. The
171 * SD1 integrator uses 0.5 in one of its calls, to correct positions
172 * for half a step of changed velocities.
174 * If v!=NULL also constrain v by adding the constraint corrections / dt.
176 * If vir!=NULL calculate the constraint virial.
178 * if veta != NULL, constraints are done assuming isotropic pressure control
179 * (i.e. constraining rdot.r = (v + veta*r).r = 0 instead of v
181 * Init_constraints must have be called once, before calling constrain.
183 * Return TRUE if OK, FALSE in case of shake error
187 gmx_constr_t init_constraints(FILE *log,
188 gmx_mtop_t *mtop, t_inputrec *ir,
189 gmx_edsam_t ed, t_state *state,
191 /* Initialize constraints stuff */
193 void set_constraints(gmx_constr_t constr,
198 /* Set up all the local constraints for the node */
200 /* The at2con t_blocka struct returned by the routines below
201 * contains a list of constraints per atom.
202 * The F_CONSTRNC constraints in this structure number consecutively
203 * after the F_CONSTR constraints.
206 t_blocka make_at2con(int start, int natoms,
207 t_ilist *ilist, t_iparams *iparams,
208 gmx_bool bDynamics, int *nflexiblecons);
209 /* Returns a block struct to go from atoms to constraints */
211 const t_blocka *atom2constraints_moltype(gmx_constr_t constr);
212 /* Returns the an array of atom to constraints lists for the moltypes */
214 const int **atom2settle_moltype(gmx_constr_t constr);
215 /* Returns the an array of atom to settle for the moltypes */
217 #define constr_iatomptr(nconstr, iatom_constr, iatom_constrnc, con) ((con) < (nconstr) ? (iatom_constr)+(con)*3 : (iatom_constrnc)+(con-nconstr)*3)
218 /* Macro for getting the constraint iatoms for a constraint number con
219 * which comes from a list where F_CONSTR and F_CONSTRNC constraints
223 gmx_bool inter_charge_group_constraints(const gmx_mtop_t *mtop);
224 /* Returns if there are inter charge group constraints */
226 gmx_bool inter_charge_group_settles(const gmx_mtop_t *mtop);
227 /* Returns if there are inter charge group settles */
229 real *constr_rmsd_data(gmx_constr_t constr);
230 /* Return the data for determining constraint RMS relative deviations.
231 * Returns NULL when LINCS is not used.
234 real constr_rmsd(gmx_constr_t constr, gmx_bool bSD2);
235 /* Return the RMSD of the constraint, bSD2 selects the second SD step */
237 real *lincs_rmsd_data(gmx_lincsdata_t lincsd);
238 /* Return the data for determining constraint RMS relative deviations */
240 real lincs_rmsd(gmx_lincsdata_t lincsd, gmx_bool bSD2);
241 /* Return the RMSD of the constraint, bSD2 selects the second SD step */
243 gmx_lincsdata_t init_lincs(FILE *fplog, gmx_mtop_t *mtop,
244 int nflexcon_global, t_blocka *at2con,
245 gmx_bool bPLINCS, int nIter, int nProjOrder);
246 /* Initializes and returns the lincs data struct */
248 void set_lincs(t_idef *idef, t_mdatoms *md,
249 gmx_bool bDynamics, t_commrec *cr,
251 /* Initialize lincs stuff */
253 void set_lincs_matrix(gmx_lincsdata_t li, real *invmass, real lambda);
254 /* Sets the elements of the LINCS constraint coupling matrix */
256 real constr_r_max(FILE *fplog, gmx_mtop_t *mtop, t_inputrec *ir);
257 /* Returns an estimate of the maximum distance between atoms
258 * required for LINCS.
262 constrain_lincs(FILE *log, gmx_bool bLog, gmx_bool bEner,
265 gmx_lincsdata_t lincsd, t_mdatoms *md,
267 rvec *x, rvec *xprime, rvec *min_proj,
268 matrix box, struct t_pbc *pbc,
269 real lambda, real *dvdlambda,
271 gmx_bool bCalcVir, tensor vir_r_m_dr,
274 int maxwarn, int *warncount);
275 /* Returns if the constraining succeeded */
278 /* helper functions for andersen temperature control, because the
279 * gmx_constr construct is only defined in constr.c. Return the list
280 * of blocks (get_sblock) and the number of blocks (get_nblocks). */
282 int *get_sblock(struct gmx_constr *constr);
284 int get_nblocks(struct gmx_constr *constr);