Remove dependence of constraints on t_mdatoms
[alexxy/gromacs.git] / src / gromacs / mdlib / tests / constrtestdata.h
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
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.
8  *
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.
13  *
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.
18  *
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.
23  *
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.
31  *
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.
34  */
35 /*! \internal \file
36  * \brief SHAKE and LINCS tests header.
37  *
38  * Contains description and constructor for the test data accumulating object,
39  * declares CPU- and GPU-based functions used to apply SHAKE or LINCS on the
40  * test data.
41  *
42  * \author Artem Zhmurov <zhmurov@gmail.com>
43  * \ingroup module_mdlib
44  */
45
46 #ifndef GMX_MDLIB_TESTS_CONSTRTESTDATA_H
47 #define GMX_MDLIB_TESTS_CONSTRTESTDATA_H
48
49 #include <memory>
50 #include <vector>
51
52 #include "gromacs/gmxlib/nrnb.h"
53 #include "gromacs/gpu_utils/gpu_testutils.h"
54 #include "gromacs/math/paddedvector.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/math/vectypes.h"
57 #include "gromacs/mdlib/gmx_omp_nthreads.h"
58 #include "gromacs/mdlib/lincs.h"
59 #include "gromacs/mdlib/shake.h"
60 #include "gromacs/mdtypes/inputrec.h"
61 #include "gromacs/pbcutil/pbc.h"
62 #include "gromacs/topology/idef.h"
63 #include "gromacs/topology/ifunc.h"
64 #include "gromacs/topology/topology.h"
65
66 namespace gmx
67 {
68 namespace test
69 {
70
71 /* \brief
72  * Constraints test data structure.
73  *
74  * Structure to collect all the necessary data, including system coordinates and topology,
75  * constraints information, etc. The structure can be reset and reused.
76  */
77 class ConstraintsTestData
78 {
79 public:
80     //! Human-friendly name for a system
81     std::string title_;
82     //! Number of atoms
83     int numAtoms_;
84     //! Topology
85     gmx_mtop_t mtop_;
86     //! Masses
87     std::vector<real> masses_;
88     //! Inverse masses
89     std::vector<real> invmass_;
90     //! Input record (info that usually in .mdp file)
91     t_inputrec ir_;
92     //! Local topology
93     std::unique_ptr<InteractionDefinitions> idef_;
94     //! Computational time array (normally used to benchmark performance)
95     t_nrnb nrnb_;
96
97     //! Inverse timestep
98     real invdt_;
99     //! Number of flexible constraints
100     int nflexcon_ = 0;
101     //! Whether the virial should be computed
102     bool computeVirial_;
103     //! Scaled virial
104     tensor virialScaled_;
105     //! Scaled virial (reference values)
106     tensor virialScaledRef_;
107     //! If the free energy is computed
108     bool compute_dHdLambda_;
109     //! If there are atoms with perturbed mass
110     bool hasMassPerturbed_ = false;
111     //! Lambda value
112     real lambda_ = 0.0;
113     //! For free energy computation
114     real dHdLambda_;
115     //! For free energy computation (reference value)
116     real dHdLambdaRef_;
117
118     //! Coordinates before the timestep
119     PaddedVector<RVec> x_;
120     //! Coordinates after timestep, output for the constraints
121     PaddedVector<RVec> xPrime_;
122     //! Backup for coordinates (for reset)
123     PaddedVector<RVec> xPrime0_;
124     //! Intermediate set of coordinates (normally used for projection correction)
125     PaddedVector<RVec> xPrime2_;
126     //! Velocities
127     PaddedVector<RVec> v_;
128     //! Backup for velocities (for reset)
129     PaddedVector<RVec> v0_;
130
131     //! Constraints data (type1-i1-j1-type2-i2-j2-...)
132     std::vector<int> constraints_;
133     //! Target lengths for all constraint types
134     std::vector<real> constraintsR0_;
135
136     /*! \brief
137      * Constructor for the object with all parameters and variables needed by constraints algorithms.
138      *
139      * This constructor assembles stubs for all the data structures, required to initialize
140      * and apply LINCS and SHAKE constraints. The coordinates and velocities before constraining
141      * are saved to allow for reset. The constraints data are stored for testing after constraints
142      * were applied.
143      *
144      * \param[in]  title                Human-friendly name of the system.
145      * \param[in]  numAtoms             Number of atoms in the system.
146      * \param[in]  masses               Atom masses. Size of this vector should be equal to numAtoms.
147      * \param[in]  constraints          List of constraints, organized in triples of integers.
148      *                                  First integer is the index of type for a constraint, second
149      *                                  and third are the indices of constrained atoms. The types
150      *                                  of constraints should be sequential but not necessarily
151      *                                  start from zero (which is the way they normally are in
152      *                                  GROMACS).
153      * \param[in]  constraintsR0        Target values for bond lengths for bonds of each type. The
154      *                                  size of this vector should be equal to the total number of
155      *                                  unique types in constraints vector.
156      * \param[in]  computeVirial        Whether the virial should be computed.
157      * \param[in]  virialScaledRef      Reference values for scaled virial tensor.
158      * \param[in]  compute_dHdLambda    Whether free energy should be computed.
159      * \param[in]  dHdLambdaRef         Reference value for dHdLambda.
160      * \param[in]  initialTime          Initial time.
161      * \param[in]  timestep             Timestep.
162      * \param[in]  x                    Coordinates before integration step.
163      * \param[in]  xPrime               Coordinates after integration step, but before constraining.
164      * \param[in]  v                    Velocities before constraining.
165      * \param[in]  shakeTolerance       Target tolerance for SHAKE.
166      * \param[in]  shakeUseSOR          Use successive over-relaxation method for SHAKE iterations.
167      *                                  The general formula is:
168      *                                     x_n+1 = (1-omega)*x_n + omega*f(x_n),
169      *                                  where omega = 1 if SOR is off and may be < 1 if SOR is on.
170      * \param[in]  lincsNumIterations   Number of iterations used to compute the inverse matrix.
171      * \param[in]  lincsExpansionOrder  The order for algorithm that adjusts the direction of the
172      *                                  bond after constraints are applied.
173      * \param[in]  lincsWarnAngle       The threshold value for the change in bond angle. When
174      *                                  exceeded the program will issue a warning.
175      *
176      */
177     ConstraintsTestData(const std::string&       title,
178                         int                      numAtoms,
179                         std::vector<real>        masses,
180                         std::vector<int>         constraints,
181                         std::vector<real>        constraintsR0,
182                         bool                     computeVirial,
183                         tensor                   virialScaledRef,
184                         bool                     compute_dHdLambda,
185                         float                    dHdLambdaRef,
186                         real                     initialTime,
187                         real                     timestep,
188                         const std::vector<RVec>& x,
189                         const std::vector<RVec>& xPrime,
190                         const std::vector<RVec>& v,
191                         real                     shakeTolerance,
192                         gmx_bool                 shakeUseSOR,
193                         int                      lincsNumIterations,
194                         int                      lincsExpansionOrder,
195                         real                     lincsWarnAngle);
196
197     /*! \brief
198      * Reset the data structure so it can be reused.
199      *
200      * Set the coordinates and velocities back to their values before
201      * constraining. The scaled virial tensor and dHdLambda are zeroed.
202      *
203      */
204     void reset();
205 };
206
207 } // namespace test
208 } // namespace gmx
209
210 #endif // GMX_MDLIB_TESTS_CONSTRTESTDATA_H