2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
35 /*! \libinternal \file
36 * \brief Declares interface to SHAKE code.
38 * \author David van der Spoel <david.vanderspoel@icm.uu.se>
39 * \author Berk Hess <hess@kth.se>
40 * \author Mark Abraham <mark.j.abraham@gmail.com>
41 * \ingroup module_mdlib
45 #ifndef GMX_MDLIB_SHAKE_H
46 #define GMX_MDLIB_SHAKE_H
48 #include "gromacs/math/vec.h"
49 #include "gromacs/topology/block.h"
50 #include "gromacs/utility/real.h"
52 struct InteractionList;
53 class InteractionDefinitions;
63 enum class ConstraintVariable : int;
66 * \brief Working data for the SHAKE algorithm
70 //! Returns the number of SHAKE blocks */
71 int numShakeBlocks() const { return sblock.size() - 1; }
73 //! The reference constraint vectors
74 std::vector<RVec> rij;
75 //! The reduced mass of the two atoms in each constraint times 0.5
76 std::vector<real> half_of_reduced_mass;
77 //! Multiplicative tolerance on the difference in the square of the constrained distance
78 std::vector<real> distance_squared_tolerance;
79 //! The reference constraint distances squared
80 std::vector<real> constraint_distance_squared;
88 //! The SHAKE blocks, block i contains constraints sblock[i]/3 to sblock[i+1]/3 */
89 std::vector<int> sblock = { 0 };
90 /*! \brief Scaled Lagrange multiplier for each constraint.
92 * Value is -2 * eta from p. 336 of the paper, divided by the
93 * constraint distance. */
94 std::vector<real> scaled_lagrange_multiplier;
97 //! Make SHAKE blocks when not using DD.
98 void make_shake_sblock_serial(shakedata* shaked, InteractionDefinitions* idef, int numAtoms);
100 //! Make SHAKE blocks when using DD.
101 void make_shake_sblock_dd(shakedata* shaked, const InteractionList& ilcon);
103 /*! \brief Shake all the atoms blockwise. It is assumed that all the constraints
104 * in the idef->shakes field are sorted, to ascending block nr. The
105 * sblock array points into the idef->shakes.iatoms field, with block 0
107 * at sblock[0] and running to ( < ) sblock[1], block n running from
108 * sblock[n] to sblock[n+1]. Array sblock should be large enough.
109 * Return TRUE when OK, FALSE when shake-error
111 bool constrain_shake(FILE* log, /* Log file */
112 shakedata* shaked, /* Total number of atoms */
113 const real invmass[], /* Atomic masses */
114 const InteractionDefinitions& idef, /* The interaction def */
115 const t_inputrec& ir, /* Input record */
116 ArrayRef<const RVec> x_s, /* Coords before update */
117 ArrayRef<RVec> xprime, /* Output coords when constraining x */
118 ArrayRef<RVec> vprime, /* Output coords when constraining v */
119 const t_pbc* pbc, /* PBC information */
120 t_nrnb* nrnb, /* Performance measure */
121 real lambda, /* FEP lambda */
122 real* dvdlambda, /* FEP force */
123 real invdt, /* 1/delta_t */
124 ArrayRef<RVec> v, /* Also constrain v if not empty */
125 bool bCalcVir, /* Calculate r x m delta_r */
126 tensor vir_r_m_dr, /* sum r x m delta_r */
127 bool bDumpOnError, /* Dump debugging stuff on error*/
128 ConstraintVariable econq); /* which type of constraint is occurring */
130 /*! \brief Regular iterative shake */
131 void cshake(const int iatom[],
135 ArrayRef<const real> dist2,
138 ArrayRef<const RVec> rij,
139 ArrayRef<const real> half_of_reduced_mass,
141 const real invmass[],
142 ArrayRef<const real> distance_squared_tolerance,
143 ArrayRef<real> scaled_lagrange_multiplier,