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;
64 enum class ConstraintVariable : int;
67 * \brief Working data for the SHAKE algorithm
71 //! Returns the number of SHAKE blocks */
72 int numShakeBlocks() const { return sblock.size() - 1; }
74 //! The reference constraint vectors
75 std::vector<RVec> rij;
76 //! The reduced mass of the two atoms in each constraint times 0.5
77 std::vector<real> half_of_reduced_mass;
78 //! Multiplicative tolerance on the difference in the square of the constrained distance
79 std::vector<real> distance_squared_tolerance;
80 //! The reference constraint distances squared
81 std::vector<real> constraint_distance_squared;
89 //! The SHAKE blocks, block i contains constraints sblock[i]/3 to sblock[i+1]/3 */
90 std::vector<int> sblock = { 0 };
91 /*! \brief Scaled Lagrange multiplier for each constraint.
93 * Value is -2 * eta from p. 336 of the paper, divided by the
94 * constraint distance. */
95 std::vector<real> scaled_lagrange_multiplier;
98 //! Make SHAKE blocks when not using DD.
99 void make_shake_sblock_serial(shakedata* shaked, InteractionDefinitions* idef, const t_mdatoms& md);
101 //! Make SHAKE blocks when using DD.
102 void make_shake_sblock_dd(shakedata* shaked, const InteractionList& ilcon);
104 /*! \brief Shake all the atoms blockwise. It is assumed that all the constraints
105 * in the idef->shakes field are sorted, to ascending block nr. The
106 * sblock array points into the idef->shakes.iatoms field, with block 0
108 * at sblock[0] and running to ( < ) sblock[1], block n running from
109 * sblock[n] to sblock[n+1]. Array sblock should be large enough.
110 * Return TRUE when OK, FALSE when shake-error
112 bool constrain_shake(FILE* log, /* Log file */
113 shakedata* shaked, /* Total number of atoms */
114 const real invmass[], /* Atomic masses */
115 const InteractionDefinitions& idef, /* The interaction def */
116 const t_inputrec& ir, /* Input record */
117 ArrayRef<const RVec> x_s, /* Coords before update */
118 ArrayRef<RVec> xprime, /* Output coords when constraining x */
119 ArrayRef<RVec> vprime, /* Output coords when constraining v */
120 const t_pbc* pbc, /* PBC information */
121 t_nrnb* nrnb, /* Performance measure */
122 real lambda, /* FEP lambda */
123 real* dvdlambda, /* FEP force */
124 real invdt, /* 1/delta_t */
125 ArrayRef<RVec> v, /* Also constrain v if not empty */
126 bool bCalcVir, /* Calculate r x m delta_r */
127 tensor vir_r_m_dr, /* sum r x m delta_r */
128 bool bDumpOnError, /* Dump debugging stuff on error*/
129 ConstraintVariable econq); /* which type of constraint is occurring */
131 /*! \brief Regular iterative shake */
132 void cshake(const int iatom[],
136 ArrayRef<const real> dist2,
139 ArrayRef<const RVec> rij,
140 ArrayRef<const real> half_of_reduced_mass,
142 const real invmass[],
143 ArrayRef<const real> distance_squared_tolerance,
144 ArrayRef<real> scaled_lagrange_multiplier,