src/gromacs/mdlib/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,2015,2016,2017,2018,2019, 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 /*! \libinternal \file
  38  * \brief Declares interface to constraint code.
  39  *
  40  * \author Berk Hess <hess@kth.se>
  41  * \author Mark Abraham <mark.j.abraham@gmail.com>
  42  * \ingroup module_mdlib
  43  * \inlibraryapi
  44  */
  45
  46 #ifndef GMX_MDLIB_CONSTR_H
  47 #define GMX_MDLIB_CONSTR_H
  48
  49 #include <cstdio>
  50
  51 #include "gromacs/math/vectypes.h"
  52 #include "gromacs/topology/idef.h"
  53 #include "gromacs/utility/arrayref.h"
  54 #include "gromacs/utility/basedefinitions.h"
  55 #include "gromacs/utility/classhelpers.h"
  56 #include "gromacs/utility/gmxassert.h"
  57 #include "gromacs/utility/real.h"
  58
  59 struct gmx_edsam;
  60 struct gmx_localtop_t;
  61 struct gmx_moltype_t;
  62 struct gmx_mtop_t;
  63 struct gmx_multisim_t;
  64 struct gmx_wallcycle;
  65 struct pull_t;
  66 struct t_blocka;
  67 struct t_commrec;
  68 struct t_ilist;
  69 struct t_inputrec;
  70 struct t_mdatoms;
  71 struct t_nrnb;
  72 struct t_pbc;
  73 class t_state;
  74
  75 namespace gmx
  76 {
  77 template<typename T>
  78 class ArrayRefWithPadding;
  79
  80 //! Describes supported flavours of constrained updates.
  81 enum class ConstraintVariable : int
  82 {
  83     Positions,     /* Constrain positions (mass weighted)             */
  84     Velocities,    /* Constrain velocities (mass weighted)            */
  85     Derivative,    /* Constrain a derivative (mass weighted),         *
  86                     * for instance velocity or acceleration,          *
  87                     * constraint virial can not be calculated.        */
  88     Deriv_FlexCon, /* As Derivative, but only output flex. con.       */
  89     Force,         /* Constrain forces (non mass-weighted)            */
  90     // TODO What does this do? Improve the comment.
  91     ForceDispl /* Like Force, but free particles will have mass
  92                 * 1 and frozen particles mass 0                   */
  93 };
  94
  95 /*! \libinternal
  96  * \brief Handles constraints */
  97 class Constraints
  98 {
  99 private:
 100     /*! \brief Constructor
 101      *
 102      * Private to enforce use of makeConstraints() factory
 103      * function. */
 104     Constraints(const gmx_mtop_t&     mtop,
 105                 const t_inputrec&     ir,
 106                 pull_t*               pull_work,
 107                 FILE*                 log,
 108                 const t_mdatoms&      md,
 109                 const t_commrec*      cr,
 110                 const gmx_multisim_t* ms,
 111                 t_nrnb*               nrnb,
 112                 gmx_wallcycle*        wcycle,
 113                 bool                  pbcHandlingRequired,
 114                 int                   numConstraints,
 115                 int                   numSettles);
 116
 117 public:
 118     /*! \brief This member type helps implement a factory
 119      * function, because its objects can access the private
 120      * constructor. */
 121     struct CreationHelper;
 122
 123     ~Constraints();
 124
 125     /*! \brief Returns the total number of flexible constraints in the system. */
 126     int numFlexibleConstraints() const;
 127
 128     /*! \brief Returns whether the system contains perturbed constraints */
 129     bool havePerturbedConstraints() const;
 130
 131     /*! \brief Set up all the local constraints for the domain.
 132      *
 133      * \todo Make this a callback that is called automatically
 134      * once a new domain has been made. */
 135     void setConstraints(const gmx_localtop_t& top, const t_mdatoms& md);
 136
 137     /*! \brief Applies constraints to coordinates.
 138      *
 139      * When econq=ConstraintVariable::Positions constrains
 140      * coordinates xprime using th directions in x, min_proj is
 141      * not used.
 142      *
 143      * When econq=ConstraintVariable::Derivative, calculates the
 144      * components xprime in the constraint directions and
 145      * subtracts these components from min_proj.  So when
 146      * min_proj=xprime, the constraint components are projected
 147      * out.
 148      *
 149      * When econq=ConstraintVariable::Deriv_FlexCon, the same is
 150      * done as with ConstraintVariable::Derivative, but only the
 151      * components of the flexible constraints are stored.
 152      *
 153      * delta_step is used for determining the constraint reference lengths
 154      * when lenA != lenB or will the pull code with a pulling rate.
 155      * step + delta_step is the step at which the final configuration
 156      * is meant to be; for update delta_step = 1.
 157      *
 158      * step_scaling can be used to update coordinates based on the time
 159      * step multiplied by this factor. Thus, normally 1.0 is passed. The
 160      * SD1 integrator uses 0.5 in one of its calls, to correct positions
 161      * for half a step of changed velocities.
 162      *
 163      * If v!=NULL also constrain v by adding the constraint corrections / dt.
 164      *
 165      * If vir!=NULL calculate the constraint virial.
 166      *
 167      * Return whether the application of constraints succeeded without error.
 168      */
 169     bool apply(bool               bLog,
 170                bool               bEner,
 171                int64_t            step,
 172                int                delta_step,
 173                real               step_scaling,
 174                rvec*              x,
 175                rvec*              xprime,
 176                rvec*              min_proj,
 177                const matrix       box,
 178                real               lambda,
 179                real*              dvdlambda,
 180                rvec*              v,
 181                tensor*            vir,
 182                ConstraintVariable econq);
 183     //! Links the essentialdynamics and constraint code.
 184     void saveEdsamPointer(gmx_edsam* ed);
 185     //! Getter for use by domain decomposition.
 186     ArrayRef<const t_blocka> atom2constraints_moltype() const;
 187     //! Getter for use by domain decomposition.
 188     ArrayRef<const std::vector<int>> atom2settle_moltype() const;
 189
 190     /*! \brief Return the data for reduction for determining
 191      * constraint RMS relative deviations, or an empty ArrayRef
 192      * when not supported for any active constraints. */
 193     ArrayRef<real> rmsdData() const;
 194     /*! \brief Return the RMSD of the constraints when available. */
 195     real rmsd() const;
 196
 197     /*! \brief Get the total number of constraints.
 198      *
 199      * \returns Total number of constraints, including SETTLE and LINCS/SHAKE constraints.
 200      */
 201     int numConstraintsTotal();
 202
 203 private:
 204     //! Implementation type.
 205     class Impl;
 206     //! Implementation object.
 207     PrivateImplPointer<Impl> impl_;
 208 };
 209
 210 /*! \brief Generate a fatal error because of too many LINCS/SETTLE warnings. */
 211 [[noreturn]] void too_many_constraint_warnings(int eConstrAlg, int warncount);
 212
 213 /*! \brief Returns whether constraint with parameter \p iparamsIndex is a flexible constraint */
 214 static inline bool isConstraintFlexible(const t_iparams* iparams, int iparamsIndex)
 215 {
 216     GMX_ASSERT(iparams != nullptr, "Need a valid iparams array");
 217
 218     return (iparams[iparamsIndex].constr.dA == 0 && iparams[iparamsIndex].constr.dB == 0);
 219 };
 220
 221 /* The at2con t_blocka struct returned by the routines below
 222  * contains a list of constraints per atom.
 223  * The F_CONSTRNC constraints in this structure number consecutively
 224  * after the F_CONSTR constraints.
 225  */
 226
 227 /*! \brief Tells make_at2con how to treat flexible constraints */
 228 enum class FlexibleConstraintTreatment
 229 {
 230     Include, //!< Include all flexible constraints
 231     Exclude  //!< Exclude all flexible constraints
 232 };
 233
 234 /*! \brief Returns the flexible constraint treatment depending on whether the integrator is dynamic */
 235 FlexibleConstraintTreatment flexibleConstraintTreatment(bool haveDynamicsIntegrator);
 236
 237 /*! \brief Returns a block struct to go from atoms to constraints
 238  *
 239  * The block struct will contain constraint indices with lower indices
 240  * directly matching the order in F_CONSTR and higher indices matching
 241  * the order in F_CONSTRNC offset by the number of constraints in F_CONSTR.
 242  *
 243  * \param[in]  moltype   The molecule data
 244  * \param[in]  iparams   Interaction parameters, can be null when
 245  * flexibleConstraintTreatment=Include \param[in]  flexibleConstraintTreatment  The flexible
 246  * constraint treatment, see enum above \returns a block struct with all constraints for each atom
 247  */
 248 t_blocka make_at2con(const gmx_moltype_t&           moltype,
 249                      gmx::ArrayRef<const t_iparams> iparams,
 250                      FlexibleConstraintTreatment    flexibleConstraintTreatment);
 251
 252 /*! \brief Returns a block struct to go from atoms to constraints
 253  *
 254  * The block struct will contain constraint indices with lower indices
 255  * directly matching the order in F_CONSTR and higher indices matching
 256  * the order in F_CONSTRNC offset by the number of constraints in F_CONSTR.
 257  *
 258  * \param[in]  numAtoms  The number of atoms to construct the list for
 259  * \param[in]  ilist     Interaction list, size F_NRE
 260  * \param[in]  iparams   Interaction parameters, can be null when
 261  * flexibleConstraintTreatment=Include \param[in]  flexibleConstraintTreatment  The flexible
 262  * constraint treatment, see enum above \returns a block struct with all constraints for each atom
 263  */
 264 t_blocka make_at2con(int                         numAtoms,
 265                      const t_ilist*              ilist,
 266                      const t_iparams*            iparams,
 267                      FlexibleConstraintTreatment flexibleConstraintTreatment);
 268
 269 /*! \brief Returns an array of atom to constraints lists for the moltypes */
 270 const t_blocka* atom2constraints_moltype(const Constraints* constr);
 271
 272 //! Return the number of flexible constraints in the \c ilist and \c iparams.
 273 int countFlexibleConstraints(const t_ilist* ilist, const t_iparams* iparams);
 274
 275 /*! \brief Returns the constraint iatoms for a constraint number con
 276  * which comes from a list where F_CONSTR and F_CONSTRNC constraints
 277  * are concatenated. */
 278 inline const int* constr_iatomptr(gmx::ArrayRef<const int> iatom_constr,
 279                                   gmx::ArrayRef<const int> iatom_constrnc,
 280                                   int                      con)
 281 {
 282     if (con * 3 < iatom_constr.ssize())
 283     {
 284         return iatom_constr.data() + con * 3;
 285     }
 286     else
 287     {
 288         return iatom_constrnc.data() + con * 3 - iatom_constr.size();
 289     }
 290 };
 291
 292 /*! \brief Returns whether there are inter charge group constraints */
 293 bool inter_charge_group_constraints(const gmx_mtop_t& mtop);
 294
 295 /*! \brief Returns whether there are inter charge group settles */
 296 bool inter_charge_group_settles(const gmx_mtop_t& mtop);
 297
 298 /*! \brief Constrain the initial coordinates and velocities */
 299 void do_constrain_first(FILE*                     log,
 300                         gmx::Constraints*         constr,
 301                         const t_inputrec*         inputrec,
 302                         const t_mdatoms*          md,
 303                         int                       natoms,
 304                         ArrayRefWithPadding<RVec> x,
 305                         ArrayRefWithPadding<RVec> v,
 306                         const matrix              box,
 307                         real                      lambda);
 308
 309 } // namespace gmx
 310
 311 #endif