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