2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
37 /*! \libinternal \file
38 * \brief Declares interface to constraint code.
40 * \author Berk Hess <hess@kth.se>
41 * \author Mark Abraham <mark.j.abraham@gmail.com>
42 * \ingroup module_mdlib
46 #ifndef GMX_MDLIB_CONSTR_H
47 #define GMX_MDLIB_CONSTR_H
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"
60 struct gmx_localtop_t;
63 struct gmx_multisim_t;
78 class ArrayRefWithPadding;
80 //! Describes supported flavours of constrained updates.
81 enum class ConstraintVariable : int
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 */
96 * \brief Handles constraints */
100 /*! \brief Constructor
102 * Private to enforce use of makeConstraints() factory
104 Constraints(const gmx_mtop_t& mtop,
105 const t_inputrec& ir,
110 const gmx_multisim_t* ms,
112 gmx_wallcycle* wcycle,
113 bool pbcHandlingRequired,
118 /*! \brief This member type helps implement a factory
119 * function, because its objects can access the private
121 struct CreationHelper;
125 /*! \brief Returns the total number of flexible constraints in the system. */
126 int numFlexibleConstraints() const;
128 /*! \brief Returns whether the system contains perturbed constraints */
129 bool havePerturbedConstraints() const;
131 /*! \brief Set up all the local constraints for the domain.
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);
137 /*! \brief Applies constraints to coordinates.
139 * When econq=ConstraintVariable::Positions constrains
140 * coordinates xprime using th directions in x, min_proj is
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
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.
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.
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.
163 * If v!=NULL also constrain v by adding the constraint corrections / dt.
165 * If vir!=NULL calculate the constraint virial.
167 * Return whether the application of constraints succeeded without error.
169 bool apply(bool bLog,
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;
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. */
197 /*! \brief Get the total number of constraints.
199 * \returns Total number of constraints, including SETTLE and LINCS/SHAKE constraints.
201 int numConstraintsTotal();
204 //! Implementation type.
206 //! Implementation object.
207 PrivateImplPointer<Impl> impl_;
210 /*! \brief Generate a fatal error because of too many LINCS/SETTLE warnings. */
211 [[noreturn]] void too_many_constraint_warnings(int eConstrAlg, int warncount);
213 /*! \brief Returns whether constraint with parameter \p iparamsIndex is a flexible constraint */
214 static inline bool isConstraintFlexible(const t_iparams* iparams, int iparamsIndex)
216 GMX_ASSERT(iparams != nullptr, "Need a valid iparams array");
218 return (iparams[iparamsIndex].constr.dA == 0 && iparams[iparamsIndex].constr.dB == 0);
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.
227 /*! \brief Tells make_at2con how to treat flexible constraints */
228 enum class FlexibleConstraintTreatment
230 Include, //!< Include all flexible constraints
231 Exclude //!< Exclude all flexible constraints
234 /*! \brief Returns the flexible constraint treatment depending on whether the integrator is dynamic */
235 FlexibleConstraintTreatment flexibleConstraintTreatment(bool haveDynamicsIntegrator);
237 /*! \brief Returns a block struct to go from atoms to constraints
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.
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
248 t_blocka make_at2con(const gmx_moltype_t& moltype,
249 gmx::ArrayRef<const t_iparams> iparams,
250 FlexibleConstraintTreatment flexibleConstraintTreatment);
252 /*! \brief Returns a block struct to go from atoms to constraints
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.
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
264 t_blocka make_at2con(int numAtoms,
265 const t_ilist* ilist,
266 const t_iparams* iparams,
267 FlexibleConstraintTreatment flexibleConstraintTreatment);
269 /*! \brief Returns an array of atom to constraints lists for the moltypes */
270 const t_blocka* atom2constraints_moltype(const Constraints* constr);
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);
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,
282 if (con * 3 < iatom_constr.ssize())
284 return iatom_constr.data() + con * 3;
288 return iatom_constrnc.data() + con * 3 - iatom_constr.size();
292 /*! \brief Returns whether there are inter charge group constraints */
293 bool inter_charge_group_constraints(const gmx_mtop_t& mtop);
295 /*! \brief Returns whether there are inter charge group settles */
296 bool inter_charge_group_settles(const gmx_mtop_t& mtop);
298 /*! \brief Constrain the initial coordinates and velocities */
299 void do_constrain_first(FILE* log,
300 gmx::Constraints* constr,
301 const t_inputrec* inputrec,
304 ArrayRefWithPadding<RVec> x,
305 ArrayRefWithPadding<RVec> v,