src/gromacs/legacyheaders/constr.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
  10  *
  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.
  15  *
  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.
  20  *
  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.
  25  *
  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.
  33  *
  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.
  36  */
  37
  38 #ifndef _constr_h
  39 #define _constr_h
  40
  41 #include "gromacs/legacyheaders/typedefs.h"
  42
  43 #ifdef __cplusplus
  44 extern "C" {
  45 #endif
  46
  47 struct t_pbc;
  48
  49 enum
  50 {
  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) */
  59 };
  60
  61 int n_flexible_constraints(struct gmx_constr *constr);
  62 /* Returns the total number of flexible constraints in the system */
  63
  64 void too_many_constraint_warnings(int eConstrAlg, int warncount);
  65 /* Generate a fatal error because of too many LINCS/SETTLE warnings */
  66
  67 gmx_shakedata_t shake_init();
  68 /* Initializes and return the SHAKE data structure */
  69
  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
  93  * starting
  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
  97  */
  98
  99 gmx_settledata_t settle_init(real mO, real mH, real invmO, real invmH,
 100                              real dOH, real dHH);
 101 /* Initializes and returns a structure with SETTLE parameters */
 102
 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            */
 113              int                *xerror,
 114              t_vetavars         *vetavar           /* variables for pressure control */
 115              );
 116
 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  */
 120                  rvec x[],
 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.
 126  */
 127
 128 void cshake(const atom_id iatom[], int ncon, int *nnit, int maxnit,
 129             const real dist2[], real xp[], const real rij[], const real m2[], real omega,
 130             const real invmass[], const real tt[], real lagr[], int *nerror);
 131 /* Regular iterative shake */
 132
 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);
 136
 137 gmx_bool constrain(FILE *log, gmx_bool bLog, gmx_bool bEner,
 138                    gmx_constr_t constr,
 139                    t_idef *idef,
 140                    t_inputrec *ir,
 141                    gmx_ekindata_t *ekind,
 142                    t_commrec *cr,
 143                    gmx_int64_t step, int delta_step,
 144                    real step_scaling,
 145                    t_mdatoms *md,
 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);
 151 /*
 152  * When econq=econqCoord constrains coordinates xprime using th
 153  * directions in x, min_proj is not used.
 154  *
 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.
 158  *
 159  * When econq=econqDeriv_FlexCon, the same is done as with econqDeriv,
 160  * but only the components of the flexible constraints are stored.
 161  *
 162  * When bMolPBC=TRUE, assume that molecules might be broken: correct PBC.
 163  *
 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.
 168  *
 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.
 173  *
 174  * If v!=NULL also constrain v by adding the constraint corrections / dt.
 175  *
 176  * If vir!=NULL calculate the constraint virial.
 177  *
 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
 180  *
 181  * Init_constraints must have be called once, before calling constrain.
 182  *
 183  * Return TRUE if OK, FALSE in case of shake error
 184  *
 185  */
 186
 187 gmx_constr_t init_constraints(FILE *log,
 188                               gmx_mtop_t *mtop, t_inputrec *ir,
 189                               gmx_edsam_t ed, t_state *state,
 190                               t_commrec *cr);
 191 /* Initialize constraints stuff */
 192
 193 void set_constraints(gmx_constr_t    constr,
 194                      gmx_localtop_t *top,
 195                      t_inputrec     *ir,
 196                      t_mdatoms      *md,
 197                      t_commrec      *cr);
 198 /* Set up all the local constraints for the node */
 199
 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.
 204  */
 205
 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 */
 210
 211 const t_blocka *atom2constraints_moltype(gmx_constr_t constr);
 212 /* Returns the an array of atom to constraints lists for the moltypes */
 213
 214 const int **atom2settle_moltype(gmx_constr_t constr);
 215 /* Returns the an array of atom to settle for the moltypes */
 216
 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
 220  * are concatenated.
 221  */
 222
 223 gmx_bool inter_charge_group_constraints(const gmx_mtop_t *mtop);
 224 /* Returns if there are inter charge group constraints */
 225
 226 gmx_bool inter_charge_group_settles(const gmx_mtop_t *mtop);
 227 /* Returns if there are inter charge group settles */
 228
 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.
 232  */
 233
 234 real constr_rmsd(gmx_constr_t constr, gmx_bool bSD2);
 235 /* Return the RMSD of the constraint, bSD2 selects the second SD step */
 236
 237 real *lincs_rmsd_data(gmx_lincsdata_t lincsd);
 238 /* Return the data for determining constraint RMS relative deviations */
 239
 240 real lincs_rmsd(gmx_lincsdata_t lincsd, gmx_bool bSD2);
 241 /* Return the RMSD of the constraint, bSD2 selects the second SD step */
 242
 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 */
 247
 248 void set_lincs(t_idef *idef, t_mdatoms *md,
 249                gmx_bool bDynamics, t_commrec *cr,
 250                gmx_lincsdata_t li);
 251 /* Initialize lincs stuff */
 252
 253 void set_lincs_matrix(gmx_lincsdata_t li, real *invmass, real lambda);
 254 /* Sets the elements of the LINCS constraint coupling matrix */
 255
 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.
 259  */
 260
 261 gmx_bool
 262 constrain_lincs(FILE *log, gmx_bool bLog, gmx_bool bEner,
 263                 t_inputrec *ir,
 264                 gmx_int64_t step,
 265                 gmx_lincsdata_t lincsd, t_mdatoms *md,
 266                 t_commrec *cr,
 267                 rvec *x, rvec *xprime, rvec *min_proj,
 268                 matrix box, struct t_pbc *pbc,
 269                 real lambda, real *dvdlambda,
 270                 real invdt, rvec *v,
 271                 gmx_bool bCalcVir, tensor vir_r_m_dr,
 272                 int econ,
 273                 t_nrnb *nrnb,
 274                 int maxwarn, int *warncount);
 275 /* Returns if the constraining succeeded */
 276
 277
 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).  */
 281
 282 int *get_sblock(struct gmx_constr *constr);
 283
 284 int get_nblocks(struct gmx_constr *constr);
 285
 286 #ifdef __cplusplus
 287 }
 288 #endif
 289 #endif