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