include/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  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, 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.
  11  *
  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.
  16  *
  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.
  21  *
  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.
  26  *
  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.
  34  *
  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.
  37  */
  38
  39 #ifndef _constr_h
  40 #define _constr_h
  41 #include "visibility.h"
  42 #include "typedefs.h"
  43 #include "types/commrec.h"
  44
  45 #ifdef __cplusplus
  46 extern "C" {
  47 #endif
  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 GMX_LIBMD_EXPORT
  62 int n_flexible_constraints(struct gmx_constr *constr);
  63 /* Returns the total number of flexible constraints in the system */
  64
  65 void too_many_constraint_warnings(int eConstrAlg,int warncount);
  66 /* Generate a fatal error because of too many LINCS/SETTLE warnings */
  67
  68 gmx_shakedata_t shake_init();
  69 /* Initializes and return the SHAKE data structure */
  70
  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
  95  * starting
  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
  99  */
 100
 101 gmx_settledata_t settle_init(real mO,real mH,real invmO,real invmH,
 102                                     real dOH,real dHH);
 103 /* Initializes and returns a structure with SETTLE parameters */
 104
 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            */
 115              int *xerror,
 116              t_vetavars *vetavar     /* variables for pressure control */
 117     );
 118
 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  */
 123                  rvec x[],
 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.
 129  */
 130
 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 */
 135
 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);
 139
 140 gmx_bool constrain(FILE *log,gmx_bool bLog,gmx_bool bEner,
 141                    gmx_constr_t constr,
 142                    t_idef *idef,
 143                    t_inputrec *ir,
 144                    gmx_ekindata_t *ekind,
 145                    t_commrec *cr,
 146                    gmx_large_int_t step,int delta_step,
 147                    t_mdatoms *md,
 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);
 153 /*
 154  * When econq=econqCoord constrains coordinates xprime using th
 155  * directions in x, min_proj is not used.
 156  *
 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.
 160  *
 161  * When econq=econqDeriv_FlexCon, the same is done as with econqDeriv,
 162  * but only the components of the flexible constraints are stored.
 163  *
 164  * When bMolPBC=TRUE, assume that molecules might be broken: correct PBC.
 165  *
 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.
 170  *
 171  * If v!=NULL also constrain v by adding the constraint corrections / dt.
 172  *
 173  * If vir!=NULL calculate the constraint virial.
 174  *
 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
 177  *
 178  * Init_constraints must have be called once, before calling constrain.
 179  *
 180  * Return TRUE if OK, FALSE in case of shake error
 181  *
 182  */
 183
 184 GMX_LIBMD_EXPORT
 185 gmx_constr_t init_constraints(FILE *log,
 186                                      gmx_mtop_t *mtop,t_inputrec *ir,
 187                                      gmx_edsam_t ed,t_state *state,
 188                                      t_commrec *cr);
 189 /* Initialize constraints stuff */
 190
 191 GMX_LIBMD_EXPORT
 192 void set_constraints(gmx_constr_t constr,
 193                                                         gmx_localtop_t *top,
 194                                                         t_inputrec *ir,
 195                                                         t_mdatoms *md,
 196                                                         t_commrec *cr);
 197 /* Set up all the local constraints for the node */
 198
 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.
 203  */
 204
 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 */
 209
 210 const t_blocka *atom2constraints_moltype(gmx_constr_t constr);
 211 /* Returns the an array of atom to constraints lists for the moltypes */
 212
 213 const int **atom2settle_moltype(gmx_constr_t constr);
 214 /* Returns the an array of atom to settle for the moltypes */
 215
 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
 219  * are concatenated.
 220  */
 221
 222 gmx_bool inter_charge_group_constraints(const gmx_mtop_t *mtop);
 223 /* Returns if there are inter charge group constraints */
 224
 225 gmx_bool inter_charge_group_settles(const gmx_mtop_t *mtop);
 226 /* Returns if there are inter charge group settles */
 227
 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.
 231  */
 232
 233 GMX_LIBMD_EXPORT
 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_large_int_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,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