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
39 #include "visibility.h"
41 #include "types/commrec.h"
49 econqCoord
, /* Constrain coordinates (mass weighted) */
50 econqVeloc
, /* Constrain velocities (mass weighted) */
51 econqDeriv
, /* Constrain a derivative (mass weighted), *
52 * for instance velocity or acceleration, *
53 * constraint virial can not be calculated. */
54 econqDeriv_FlexCon
, /* As econqDeriv, but only output flex. con. */
55 econqForce
, /* Constrain forces (non mass-weighted) */
56 econqForceDispl
/* Constrain forces (mass-weighted 1/0 for freeze) */
60 int n_flexible_constraints(struct gmx_constr
*constr
);
61 /* Returns the total number of flexible constraints in the system */
63 void too_many_constraint_warnings(int eConstrAlg
,int warncount
);
64 /* Generate a fatal error because of too many LINCS/SETTLE warnings */
66 gmx_shakedata_t
shake_init();
67 /* Initializes and return the SHAKE data structure */
69 gmx_bool
bshakef(FILE *log
, /* Log file */
70 gmx_shakedata_t shaked
, /* SHAKE data */
71 int natoms
, /* 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 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(FILE *fp
,
118 gmx_settledata_t settled
,int econq
,
119 int nsettle
, t_iatom iatoms
[],
120 const t_pbc
*pbc
, /* PBC data pointer, can be NULL */
122 rvec
*der
,rvec
*derp
,
123 int CalcVirAtomEnd
,tensor vir_r_m_dder
,
124 t_vetavars
*vetavar
);
125 /* Analytical algorithm to subtract the components of derivatives
126 * of coordinates working on settle type constraint.
129 void cshake(atom_id iatom
[],int ncon
,int *nnit
,int maxnit
,
130 real dist2
[],real xp
[],real rij
[],real m2
[],real omega
,
131 real invmass
[],real tt
[],real lagr
[],int *nerror
);
132 /* Regular iterative shake */
134 void crattle(atom_id iatom
[],int ncon
,int *nnit
,int maxnit
,
135 real dist2
[],real vp
[],real rij
[],real m2
[],real omega
,
136 real invmass
[],real tt
[],real lagr
[],int *nerror
,real invdt
,t_vetavars
*vetavar
);
138 gmx_bool
constrain(FILE *log
,gmx_bool bLog
,gmx_bool bEner
,
142 gmx_ekindata_t
*ekind
,
144 gmx_large_int_t step
,int delta_step
,
146 rvec
*x
,rvec
*xprime
,rvec
*min_proj
,
147 gmx_bool bMolPBC
,matrix box
,
148 real lambda
,real
*dvdlambda
,
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 * If v!=NULL also constrain v by adding the constraint corrections / dt.
171 * If vir!=NULL calculate the constraint virial.
173 * if veta != NULL, constraints are done assuming isotropic pressure control
174 * (i.e. constraining rdot.r = (v + veta*r).r = 0 instead of v
176 * Init_constraints must have be called once, before calling constrain.
178 * Return TRUE if OK, FALSE in case of shake error
183 gmx_constr_t
init_constraints(FILE *log
,
184 gmx_mtop_t
*mtop
,t_inputrec
*ir
,
185 gmx_edsam_t ed
,t_state
*state
,
187 /* Initialize constraints stuff */
190 void set_constraints(gmx_constr_t constr
,
195 /* Set up all the local constraints for the node */
197 /* The at2con t_blocka struct returned by the routines below
198 * contains a list of constraints per atom.
199 * The F_CONSTRNC constraints in this structure number consecutively
200 * after the F_CONSTR constraints.
203 t_blocka
make_at2con(int start
,int natoms
,
204 t_ilist
*ilist
,t_iparams
*iparams
,
205 gmx_bool bDynamics
,int *nflexiblecons
);
206 /* Returns a block struct to go from atoms to constraints */
208 const t_blocka
*atom2constraints_moltype(gmx_constr_t constr
);
209 /* Returns the an array of atom to constraints lists for the moltypes */
211 const int **atom2settle_moltype(gmx_constr_t constr
);
212 /* Returns the an array of atom to settle for the moltypes */
214 #define constr_iatomptr(nconstr,iatom_constr,iatom_constrnc,con) ((con) < (nconstr) ? (iatom_constr)+(con)*3 : (iatom_constrnc)+(con-nconstr)*3)
215 /* Macro for getting the constraint iatoms for a constraint number con
216 * which comes from a list where F_CONSTR and F_CONSTRNC constraints
220 gmx_bool
inter_charge_group_constraints(const gmx_mtop_t
*mtop
);
221 /* Returns if there are inter charge group constraints */
223 gmx_bool
inter_charge_group_settles(const gmx_mtop_t
*mtop
);
224 /* Returns if there are inter charge group settles */
226 real
*constr_rmsd_data(gmx_constr_t constr
);
227 /* Return the data for determining constraint RMS relative deviations.
228 * Returns NULL when LINCS is not used.
232 real
constr_rmsd(gmx_constr_t constr
,gmx_bool bSD2
);
233 /* Return the RMSD of the constraint, bSD2 selects the second SD step */
235 real
*lincs_rmsd_data(gmx_lincsdata_t lincsd
);
236 /* Return the data for determining constraint RMS relative deviations */
238 real
lincs_rmsd(gmx_lincsdata_t lincsd
,gmx_bool bSD2
);
239 /* Return the RMSD of the constraint, bSD2 selects the second SD step */
241 gmx_lincsdata_t
init_lincs(FILE *fplog
,gmx_mtop_t
*mtop
,
242 int nflexcon_global
,t_blocka
*at2con
,
243 gmx_bool bPLINCS
,int nIter
,int nProjOrder
);
244 /* Initializes and returns the lincs data struct */
246 void set_lincs(t_idef
*idef
,t_mdatoms
*md
,
247 gmx_bool bDynamics
,t_commrec
*cr
,
249 /* Initialize lincs stuff */
251 void set_lincs_matrix(gmx_lincsdata_t li
,real
*invmass
,real lambda
);
252 /* Sets the elements of the LINCS constraint coupling matrix */
254 real
constr_r_max(FILE *fplog
,gmx_mtop_t
*mtop
,t_inputrec
*ir
);
255 /* Returns an estimate of the maximum distance between atoms
256 * required for LINCS.
260 constrain_lincs(FILE *log
,gmx_bool bLog
,gmx_bool bEner
,
262 gmx_large_int_t step
,
263 gmx_lincsdata_t lincsd
,t_mdatoms
*md
,
265 rvec
*x
,rvec
*xprime
,rvec
*min_proj
,
266 matrix box
,t_pbc
*pbc
,
267 real lambda
,real
*dvdlambda
,
269 gmx_bool bCalcVir
,tensor vir_r_m_dr
,
272 int maxwarn
,int *warncount
);
273 /* Returns if the constraining succeeded */
276 /* helper functions for andersen temperature control, because the
277 * gmx_constr construct is only defined in constr.c. Return the list
278 * of blocks (get_sblock) and the number of blocks (get_nblocks). */
280 int *get_sblock(struct gmx_constr
*constr
);
282 int get_nblocks(struct gmx_constr
*constr
);