src/gromacs/legacyheaders/constr.h

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  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.
  15
  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.
  20  *
  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.
  27  *
  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.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * Gromacs Runs On Most of All Computer Systems
  35  */
  36
  37 #ifndef _constr_h
  38 #define _constr_h
  39 #include "typedefs.h"
  40
  41 #ifdef __cplusplus
  42 extern "C" {
  43 #endif
  44
  45 enum
  46 {
  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) */
  55 };
  56
  57 int n_flexible_constraints(struct gmx_constr *constr);
  58 /* Returns the total number of flexible constraints in the system */
  59
  60 void too_many_constraint_warnings(int eConstrAlg,int warncount);
  61 /* Generate a fatal error because of too many LINCS/SETTLE warnings */
  62
  63 gmx_shakedata_t shake_init();
  64 /* Initializes and return the SHAKE data structure */
  65
  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
  91  * starting
  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
  95  */
  96
  97 gmx_settledata_t settle_init(real mO,real mH,real invmO,real invmH,
  98                                     real dOH,real dHH);
  99 /* Initializes and returns a structure with SETTLE parameters */
 100
 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            */
 110                     int *xerror,
 111                     t_vetavars *vetavar     /* variables for pressure control */
 112     );
 113
 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.
 121  */
 122
 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 */
 127
 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);
 131
 132 gmx_bool constrain(FILE *log,gmx_bool bLog,gmx_bool bEner,
 133                       gmx_constr_t constr,
 134                       t_idef *idef,
 135                       t_inputrec *ir,
 136                       gmx_ekindata_t *ekind,
 137                       t_commrec *cr,
 138                       gmx_large_int_t step,int delta_step,
 139                       t_mdatoms *md,
 140                       rvec *x,rvec *xprime,rvec *min_proj,matrix box,
 141                       real lambda,real *dvdlambda,
 142                       rvec *v,tensor *vir,
 143                       t_nrnb *nrnb,int econq, gmx_bool bPscal, real veta, real vetanew);
 144 /*
 145  * When econq=econqCoord constrains coordinates xprime using th
 146  * directions in x, min_proj is not used.
 147  *
 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.
 151  *
 152  * When econq=econqDeriv_FlexCon, the same is done as with econqDeriv,
 153  * but only the components of the flexible constraints are stored.
 154  *
 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.
 159  *
 160  * If v!=NULL also constrain v by adding the constraint corrections / dt.
 161  *
 162  * If vir!=NULL calculate the constraint virial.
 163  *
 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
 166  *
 167  * Init_constraints must have be called once, before calling constrain.
 168  *
 169  * Return TRUE if OK, FALSE in case of shake error
 170  *
 171  */
 172
 173 gmx_constr_t init_constraints(FILE *log,
 174                                      gmx_mtop_t *mtop,t_inputrec *ir,
 175                                      gmx_edsam_t ed,t_state *state,
 176                                      t_commrec *cr);
 177 /* Initialize constraints stuff */
 178
 179 void set_constraints(gmx_constr_t constr,
 180                                                         gmx_localtop_t *top,
 181                                                         t_inputrec *ir,
 182                                                         t_mdatoms *md,
 183                                                         t_commrec *cr);
 184 /* Set up all the local constraints for the node */
 185
 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.
 190  */
 191
 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 */
 196
 197 t_blocka *atom2constraints_moltype(gmx_constr_t constr);
 198 /* Returns the an arry of atom to constraints lists for the moltypes */
 199
 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
 203  * are concatenated.
 204  */
 205
 206 gmx_bool inter_charge_group_constraints(gmx_mtop_t *mtop);
 207 /* Returns if there are inter charge group constraints */
 208
 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.
 212  */
 213
 214 real constr_rmsd(gmx_constr_t constr,gmx_bool bSD2);
 215 /* Return the RMSD of the constraint, bSD2 selects the second SD step */
 216
 217 real *lincs_rmsd_data(gmx_lincsdata_t lincsd);
 218 /* Return the data for determining constraint RMS relative deviations */
 219
 220 real lincs_rmsd(gmx_lincsdata_t lincsd,gmx_bool bSD2);
 221 /* Return the RMSD of the constraint, bSD2 selects the second SD step */
 222
 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 */
 227
 228 void set_lincs(t_idef *idef,t_mdatoms *md,
 229                       gmx_bool bDynamics,t_commrec *cr,
 230                       gmx_lincsdata_t li);
 231 /* Initialize lincs stuff */
 232
 233 void set_lincs_matrix(gmx_lincsdata_t li,real *invmass,real lambda);
 234 /* Sets the elements of the LINCS constraint coupling matrix */
 235
 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.
 239  */
 240
 241 gmx_bool constrain_lincs(FILE *log,gmx_bool bLog,gmx_bool bEner,
 242                             t_inputrec *ir,
 243                             gmx_large_int_t step,
 244                             gmx_lincsdata_t lincsd,t_mdatoms *md,
 245                             t_commrec *cr,
 246                             rvec *x,rvec *xprime,rvec *min_proj,matrix box,
 247                             real lambda,real *dvdlambda,
 248                             real invdt,rvec *v,
 249                             gmx_bool bCalcVir,tensor rmdr,
 250                             int econ,
 251                             t_nrnb *nrnb,
 252                             int maxwarn,int *warncount);
 253 /* Returns if the constraining succeeded */
 254
 255
 256 /* helper functions for andersen temperature control, because the
 257  * gmx_constr construct is only defined in constr.c. Return the list
 258  * of blocks (get_sblock) and the number of blocks (get_nblocks).  */
 259
 260 int *get_sblock(struct gmx_constr *constr);
 261
 262 int get_nblocks(struct gmx_constr *constr);
 263
 264 #ifdef __cplusplus
 265 }
 266 #endif
 267 #endif