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,2013, 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.
41 #include "visibility.h"
43 #include "types/commrec.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) */
62 int n_flexible_constraints(struct gmx_constr
*constr
);
63 /* Returns the total number of flexible constraints in the system */
65 void too_many_constraint_warnings(int eConstrAlg
, int warncount
);
66 /* Generate a fatal error because of too many LINCS/SETTLE warnings */
68 gmx_shakedata_t
shake_init();
69 /* Initializes and return the SHAKE data structure */
71 gmx_bool
bshakef(FILE *log
, /* Log file */
72 gmx_shakedata_t shaked
, /* SHAKE data */
73 int natoms
, /* Total number of atoms */
74 real invmass
[], /* Atomic masses */
75 int nblocks
, /* The number of shake blocks */
76 int sblock
[], /* The shake blocks */
77 t_idef
*idef
, /* The interaction def */
78 t_inputrec
*ir
, /* Input record */
79 rvec x_s
[], /* Coords before update */
80 rvec prime
[], /* Output coords */
81 t_nrnb
*nrnb
, /* Performance measure */
82 real
*lagr
, /* The Lagrange multipliers */
83 real lambda
, /* FEP lambda */
84 real
*dvdlambda
, /* FEP force */
85 real invdt
, /* 1/delta_t */
86 rvec
*v
, /* Also constrain v if v!=NULL */
87 gmx_bool bCalcVir
, /* Calculate r x m delta_r */
88 tensor vir_r_m_dr
, /* sum r x m delta_r */
89 gmx_bool bDumpOnError
, /* Dump debugging stuff on error*/
90 int econq
, /* which type of constrainint is occurring */
91 t_vetavars
*vetavar
); /* veta for pressure control */
92 /* Shake all the atoms blockwise. It is assumed that all the constraints
93 * in the idef->shakes field are sorted, to ascending block nr. The
94 * sblock array points into the idef->shakes.iatoms field, with block 0
96 * at sblock[0] and running to ( < ) sblock[1], block n running from
97 * sblock[n] to sblock[n+1]. Array sblock should be large enough.
98 * Return TRUE when OK, FALSE when shake-error
101 gmx_settledata_t
settle_init(real mO
, real mH
, real invmO
, real invmH
,
103 /* Initializes and returns a structure with SETTLE parameters */
105 void csettle(gmx_settledata_t settled
,
106 int nsettle
, /* Number of settles */
107 t_iatom iatoms
[], /* The settle iatom list */
108 const t_pbc
*pbc
, /* PBC data pointer, can be NULL */
109 real b4
[], /* Old coordinates */
110 real after
[], /* New coords, to be settled */
111 real invdt
, /* 1/delta_t */
112 real
*v
, /* Also constrain v if v!=NULL */
113 int calcvir_atom_end
, /* Calculate r x m delta_r up to this atom */
114 tensor vir_r_m_dr
, /* sum r x m delta_r */
116 t_vetavars
*vetavar
/* variables for pressure control */
119 void settle_proj(FILE *fp
,
120 gmx_settledata_t settled
, int econq
,
121 int nsettle
, t_iatom iatoms
[],
122 const t_pbc
*pbc
, /* PBC data pointer, can be NULL */
124 rvec
*der
, rvec
*derp
,
125 int CalcVirAtomEnd
, tensor vir_r_m_dder
,
126 t_vetavars
*vetavar
);
127 /* Analytical algorithm to subtract the components of derivatives
128 * of coordinates working on settle type constraint.
131 void cshake(atom_id iatom
[], int ncon
, int *nnit
, int maxnit
,
132 real dist2
[], real xp
[], real rij
[], real m2
[], real omega
,
133 real invmass
[], real tt
[], real lagr
[], int *nerror
);
134 /* Regular iterative shake */
136 void crattle(atom_id iatom
[], int ncon
, int *nnit
, int maxnit
,
137 real dist2
[], real vp
[], real rij
[], real m2
[], real omega
,
138 real invmass
[], real tt
[], real lagr
[], int *nerror
, real invdt
, t_vetavars
*vetavar
);
140 gmx_bool
constrain(FILE *log
, gmx_bool bLog
, gmx_bool bEner
,
144 gmx_ekindata_t
*ekind
,
146 gmx_large_int_t step
, int delta_step
,
148 rvec
*x
, rvec
*xprime
, rvec
*min_proj
,
149 gmx_bool bMolPBC
, matrix box
,
150 real lambda
, real
*dvdlambda
,
151 rvec
*v
, tensor
*vir
,
152 t_nrnb
*nrnb
, int econq
, gmx_bool bPscal
, real veta
, real vetanew
);
154 * When econq=econqCoord constrains coordinates xprime using th
155 * directions in x, min_proj is not used.
157 * When econq=econqDeriv, calculates the components xprime in
158 * the constraint directions and subtracts these components from min_proj.
159 * So when min_proj=xprime, the constraint components are projected out.
161 * When econq=econqDeriv_FlexCon, the same is done as with econqDeriv,
162 * but only the components of the flexible constraints are stored.
164 * When bMolPBC=TRUE, assume that molecules might be broken: correct PBC.
166 * delta_step is used for determining the constraint reference lengths
167 * when lenA != lenB or will the pull code with a pulling rate.
168 * step + delta_step is the step at which the final configuration
169 * is meant to be; for update delta_step = 1.
171 * If v!=NULL also constrain v by adding the constraint corrections / dt.
173 * If vir!=NULL calculate the constraint virial.
175 * if veta != NULL, constraints are done assuming isotropic pressure control
176 * (i.e. constraining rdot.r = (v + veta*r).r = 0 instead of v
178 * Init_constraints must have be called once, before calling constrain.
180 * Return TRUE if OK, FALSE in case of shake error
185 gmx_constr_t
init_constraints(FILE *log
,
186 gmx_mtop_t
*mtop
, t_inputrec
*ir
,
187 gmx_edsam_t ed
, t_state
*state
,
189 /* Initialize constraints stuff */
192 void set_constraints(gmx_constr_t constr
,
197 /* Set up all the local constraints for the node */
199 /* The at2con t_blocka struct returned by the routines below
200 * contains a list of constraints per atom.
201 * The F_CONSTRNC constraints in this structure number consecutively
202 * after the F_CONSTR constraints.
205 t_blocka
make_at2con(int start
, int natoms
,
206 t_ilist
*ilist
, t_iparams
*iparams
,
207 gmx_bool bDynamics
, int *nflexiblecons
);
208 /* Returns a block struct to go from atoms to constraints */
210 const t_blocka
*atom2constraints_moltype(gmx_constr_t constr
);
211 /* Returns the an array of atom to constraints lists for the moltypes */
213 const int **atom2settle_moltype(gmx_constr_t constr
);
214 /* Returns the an array of atom to settle for the moltypes */
216 #define constr_iatomptr(nconstr, iatom_constr, iatom_constrnc, con) ((con) < (nconstr) ? (iatom_constr)+(con)*3 : (iatom_constrnc)+(con-nconstr)*3)
217 /* Macro for getting the constraint iatoms for a constraint number con
218 * which comes from a list where F_CONSTR and F_CONSTRNC constraints
222 gmx_bool
inter_charge_group_constraints(const gmx_mtop_t
*mtop
);
223 /* Returns if there are inter charge group constraints */
225 gmx_bool
inter_charge_group_settles(const gmx_mtop_t
*mtop
);
226 /* Returns if there are inter charge group settles */
228 real
*constr_rmsd_data(gmx_constr_t constr
);
229 /* Return the data for determining constraint RMS relative deviations.
230 * 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
,
264 gmx_large_int_t step
,
265 gmx_lincsdata_t lincsd
, t_mdatoms
*md
,
267 rvec
*x
, rvec
*xprime
, rvec
*min_proj
,
268 matrix box
, 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
);