1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * Gromacs Runs On Most of All Computer Systems
47 econqCoord, /* Constrain coordinates (mass weighted) */
48 econqVeloc, /* Constrain velocities (mass weighted) */
49 econqDeriv, /* Constrain a derivative (mass weighted), *
50 * for instance velocity or acceleration, *
51 * constraint virial can not be calculated. */
52 econqDeriv_FlexCon, /* As econqDeriv, but only output flex. con. */
53 econqForce, /* Constrain forces (non mass-weighted) */
54 econqForceDispl /* Constrain forces (mass-weighted 1/0 for freeze) */
57 int n_flexible_constraints(struct gmx_constr *constr);
58 /* Returns the total number of flexible constraints in the system */
60 void too_many_constraint_warnings(int eConstrAlg,int warncount);
61 /* Generate a fatal error because of too many LINCS/SETTLE warnings */
63 gmx_shakedata_t shake_init();
64 /* Initializes and return the SHAKE data structure */
66 gmx_bool bshakef(FILE *log, /* Log file */
67 gmx_shakedata_t shaked, /* SHAKE data */
68 int natoms, /* Total number of atoms */
69 real invmass[], /* Atomic masses */
70 int nblocks, /* The number of shake blocks */
71 int sblock[], /* The shake blocks */
72 t_idef *idef, /* The interaction def */
73 t_inputrec *ir, /* Input record */
74 matrix box, /* The box */
75 rvec x_s[], /* Coords before update */
76 rvec prime[], /* Output coords */
77 t_nrnb *nrnb, /* Performance measure */
78 real *lagr, /* The Lagrange multipliers */
79 real lambda, /* FEP lambda */
80 real *dvdlambda, /* FEP force */
81 real invdt, /* 1/delta_t */
82 rvec *v, /* Also constrain v if v!=NULL */
83 gmx_bool bCalcVir, /* Calculate r x m delta_r */
84 tensor rmdr, /* sum r x m delta_r */
85 gmx_bool bDumpOnError, /* Dump debugging stuff on error*/
86 int econq, /* which type of constrainint is occurring */
87 t_vetavars *vetavar); /* veta for pressure control */
88 /* Shake all the atoms blockwise. It is assumed that all the constraints
89 * in the idef->shakes field are sorted, to ascending block nr. The
90 * sblock array points into the idef->shakes.iatoms field, with block 0
92 * at sblock[0] and running to ( < ) sblock[1], block n running from
93 * sblock[n] to sblock[n+1]. Array sblock should be large enough.
94 * Return TRUE when OK, FALSE when shake-error
97 gmx_settledata_t settle_init(real mO,real mH,real invmO,real invmH,
99 /* Initializes and returns a structure with SETTLE parameters */
101 void csettle(gmx_settledata_t settled,
102 int nsettle, /* Number of settles */
103 t_iatom iatoms[], /* The settle iatom list */
104 real b4[], /* Old coordinates */
105 real after[], /* New coords, to be settled */
106 real invdt, /* 1/delta_t */
107 real *v, /* Also constrain v if v!=NULL */
108 gmx_bool bCalcVir, /* Calculate r x m delta_r */
109 tensor rmdr, /* sum r x m delta_r */
111 t_vetavars *vetavar /* variables for pressure control */
114 void settle_proj(FILE *fp,
115 gmx_settledata_t settled,int econq,
116 int nsettle, t_iatom iatoms[],rvec x[],
117 rvec *der,rvec *derp,
118 gmx_bool bCalcVir,tensor rmdder, t_vetavars *vetavar);
119 /* Analytical algorithm to subtract the components of derivatives
120 * of coordinates working on settle type constraint.
123 void cshake(atom_id iatom[],int ncon,int *nnit,int maxnit,
124 real dist2[],real xp[],real rij[],real m2[],real omega,
125 real invmass[],real tt[],real lagr[],int *nerror);
126 /* Regular iterative shake */
128 void crattle(atom_id iatom[],int ncon,int *nnit,int maxnit,
129 real dist2[],real vp[],real rij[],real m2[],real omega,
130 real invmass[],real tt[],real lagr[],int *nerror,real invdt,t_vetavars *vetavar);
132 gmx_bool constrain(FILE *log,gmx_bool bLog,gmx_bool bEner,
136 gmx_ekindata_t *ekind,
138 gmx_large_int_t step,int delta_step,
140 rvec *x,rvec *xprime,rvec *min_proj,matrix box,
141 real lambda,real *dvdlambda,
143 t_nrnb *nrnb,int econq, gmx_bool bPscal, real veta, real vetanew);
145 * When econq=econqCoord constrains coordinates xprime using th
146 * directions in x, min_proj is not used.
148 * When econq=econqDeriv, calculates the components xprime in
149 * the constraint directions and subtracts these components from min_proj.
150 * So when min_proj=xprime, the constraint components are projected out.
152 * When econq=econqDeriv_FlexCon, the same is done as with econqDeriv,
153 * but only the components of the flexible constraints are stored.
155 * delta_step is used for determining the constraint reference lengths
156 * when lenA != lenB or will the pull code with a pulling rate.
157 * step + delta_step is the step at which the final configuration
158 * is meant to be; for update delta_step = 1.
160 * If v!=NULL also constrain v by adding the constraint corrections / dt.
162 * If vir!=NULL calculate the constraint virial.
164 * if veta != NULL, constraints are done assuming isotropic pressure control
165 * (i.e. constraining rdot.r = (v + veta*r).r = 0 instead of v
167 * Init_constraints must have be called once, before calling constrain.
169 * Return TRUE if OK, FALSE in case of shake error
173 gmx_constr_t init_constraints(FILE *log,
174 gmx_mtop_t *mtop,t_inputrec *ir,
175 gmx_edsam_t ed,t_state *state,
177 /* Initialize constraints stuff */
179 void set_constraints(gmx_constr_t constr,
184 /* Set up all the local constraints for the node */
186 /* The at2con t_blocka struct returned by the routines below
187 * contains a list of constraints per atom.
188 * The F_CONSTRNC constraints in this structure number consecutively
189 * after the F_CONSTR constraints.
192 t_blocka make_at2con(int start,int natoms,
193 t_ilist *ilist,t_iparams *iparams,
194 gmx_bool bDynamics,int *nflexiblecons);
195 /* Returns a block struct to go from atoms to constraints */
197 t_blocka *atom2constraints_moltype(gmx_constr_t constr);
198 /* Returns the an arry of atom to constraints lists for the moltypes */
200 #define constr_iatomptr(nconstr,iatom_constr,iatom_constrnc,con) ((con) < (nconstr) ? (iatom_constr)+(con)*3 : (iatom_constrnc)+(con-nconstr)*3)
201 /* Macro for getting the constraint iatoms for a constraint number con
202 * which comes from a list where F_CONSTR and F_CONSTRNC constraints
206 gmx_bool inter_charge_group_constraints(gmx_mtop_t *mtop);
207 /* Returns if there are inter charge group constraints */
209 real *constr_rmsd_data(gmx_constr_t constr);
210 /* Return the data for determining constraint RMS relative deviations.
211 * Returns NULL when LINCS is not used.
214 real constr_rmsd(gmx_constr_t constr,gmx_bool bSD2);
215 /* Return the RMSD of the constraint, bSD2 selects the second SD step */
217 real *lincs_rmsd_data(gmx_lincsdata_t lincsd);
218 /* Return the data for determining constraint RMS relative deviations */
220 real lincs_rmsd(gmx_lincsdata_t lincsd,gmx_bool bSD2);
221 /* Return the RMSD of the constraint, bSD2 selects the second SD step */
223 gmx_lincsdata_t init_lincs(FILE *fplog,gmx_mtop_t *mtop,
224 int nflexcon_global,t_blocka *at2con,
225 gmx_bool bPLINCS,int nIter,int nProjOrder);
226 /* Initializes and returns the lincs data struct */
228 void set_lincs(t_idef *idef,t_mdatoms *md,
229 gmx_bool bDynamics,t_commrec *cr,
231 /* Initialize lincs stuff */
233 void set_lincs_matrix(gmx_lincsdata_t li,real *invmass,real lambda);
234 /* Sets the elements of the LINCS constraint coupling matrix */
236 real constr_r_max(FILE *fplog,gmx_mtop_t *mtop,t_inputrec *ir);
237 /* Returns an estimate of the maximum distance between atoms
238 * required for LINCS.
241 gmx_bool constrain_lincs(FILE *log,gmx_bool bLog,gmx_bool bEner,
243 gmx_large_int_t step,
244 gmx_lincsdata_t lincsd,t_mdatoms *md,
246 rvec *x,rvec *xprime,rvec *min_proj,matrix box,
247 real lambda,real *dvdlambda,
249 gmx_bool bCalcVir,tensor rmdr,
252 int maxwarn,int *warncount);
253 /* Returns if the constraining succeeded */