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.
36 * \brief SHAKE and LINCS tests.
38 * \todo Better tests for virial are needed.
39 * \todo Tests for bigger systems to test threads synchronization,
40 * reduction, etc. on the GPU.
41 * \todo Tests for algorithms for derivatives.
42 * \todo Free-energy perturbation tests
44 * \author Artem Zhmurov <zhmurov@gmail.com>
45 * \ingroup module_mdlib
50 #include "constrtestdata.h"
52 #include "gromacs/utility/smalloc.h"
53 #include "gromacs/utility/stringutil.h"
60 ConstraintsTestData::ConstraintsTestData(const std::string
& title
,
62 std::vector
<real
> masses
,
63 std::vector
<int> constraints
,
64 std::vector
<real
> constraintsR0
,
66 tensor virialScaledRef
,
67 bool compute_dHdLambda
,
71 const std::vector
<RVec
>& x
,
72 const std::vector
<RVec
>& xPrime
,
73 const std::vector
<RVec
>& v
,
76 int lincsNumIterations
,
77 int lincsExpansionOrder
,
80 title_
= title
; // Human-friendly name of the system
81 numAtoms_
= numAtoms
; // Number of atoms
85 invmass_
.resize(numAtoms
); // Vector of inverse masses
87 for (int i
= 0; i
< numAtoms
; i
++)
89 invmass_
[i
] = 1.0 / masses
.at(i
);
92 // Saving constraints to check if they are satisfied after algorithm was applied
93 constraints_
= constraints
; // Constraints indices (in type-i-j format)
94 constraintsR0_
= constraintsR0
; // Equilibrium distances for each type of constraint
96 invdt_
= 1.0 / timestep
; // Inverse timestep
98 // Input record - data that usually comes from configuration file (.mdp)
100 ir_
.init_t
= initialTime
;
101 ir_
.delta_t
= timestep
;
105 md_
.nMassPerturbed
= 0;
107 md_
.invmass
= invmass_
.data();
109 md_
.homenr
= numAtoms
;
112 computeVirial_
= computeVirial
;
115 for (int i
= 0; i
< DIM
; i
++)
117 for (int j
= 0; j
< DIM
; j
++)
119 virialScaled_
[i
][j
] = 0;
120 virialScaledRef_
[i
][j
] = virialScaledRef
[i
][j
];
126 // Free energy evaluation
127 compute_dHdLambda_
= compute_dHdLambda
;
129 if (compute_dHdLambda_
)
132 dHdLambdaRef_
= dHdLambdaRef
;
141 for (index i
= 0; i
< ssize(constraints
); i
+= 3)
143 if (maxType
< constraints
.at(i
))
145 maxType
= constraints
.at(i
);
148 auto& iparams
= mtop_
.ffparams
.iparams
;
149 iparams
.resize(maxType
+ 1);
150 for (index i
= 0; i
< ssize(constraints
) / 3; i
++)
152 iparams
[constraints
.at(3 * i
)].constr
.dA
= constraintsR0
.at(constraints
.at(3 * i
));
153 iparams
[constraints
.at(3 * i
)].constr
.dB
= constraintsR0
.at(constraints
.at(3 * i
));
155 idef_
= std::make_unique
<InteractionDefinitions
>(mtop_
.ffparams
);
156 for (index i
= 0; i
< ssize(constraints
); i
++)
158 idef_
->il
[F_CONSTR
].iatoms
.push_back(constraints
.at(i
));
161 // Constraints and their parameters (global topology)
162 InteractionList interactionList
;
163 interactionList
.iatoms
.resize(constraints
.size());
164 std::copy(constraints
.begin(), constraints
.end(), interactionList
.iatoms
.begin());
165 InteractionList interactionListEmpty
;
166 interactionListEmpty
.iatoms
.resize(0);
168 gmx_moltype_t molType
;
169 molType
.atoms
.nr
= numAtoms
;
170 molType
.ilist
.at(F_CONSTR
) = interactionList
;
171 molType
.ilist
.at(F_CONSTRNC
) = interactionListEmpty
;
172 mtop_
.moltype
.push_back(molType
);
174 gmx_molblock_t molBlock
;
177 mtop_
.molblock
.push_back(molBlock
);
179 mtop_
.natoms
= numAtoms
;
180 mtop_
.bIntermolecularInteractions
= false;
182 // Coordinates and velocities
183 x_
.resizeWithPadding(numAtoms
);
184 xPrime_
.resizeWithPadding(numAtoms
);
185 xPrime0_
.resizeWithPadding(numAtoms
);
186 xPrime2_
.resizeWithPadding(numAtoms
);
188 v_
.resizeWithPadding(numAtoms
);
189 v0_
.resizeWithPadding(numAtoms
);
191 std::copy(x
.begin(), x
.end(), x_
.begin());
192 std::copy(xPrime
.begin(), xPrime
.end(), xPrime_
.begin());
193 std::copy(xPrime
.begin(), xPrime
.end(), xPrime0_
.begin());
194 std::copy(xPrime
.begin(), xPrime
.end(), xPrime2_
.begin());
196 std::copy(v
.begin(), v
.end(), v_
.begin());
197 std::copy(v
.begin(), v
.end(), v0_
.begin());
199 // SHAKE-specific parameters
200 ir_
.shake_tol
= shakeTolerance
;
201 ir_
.bShakeSOR
= shakeUseSOR
;
203 // LINCS-specific parameters
204 ir_
.nLincsIter
= lincsNumIterations
;
205 ir_
.nProjOrder
= lincsExpansionOrder
;
206 ir_
.LincsWarnAngle
= lincsWarnAngle
;
210 * Reset the data structure so it can be reused.
212 * Set the coordinates and velocities back to their values before
213 * constraining. The scaled virial tensor and dHdLambda are zeroed.
216 void ConstraintsTestData::reset()
224 for (int i
= 0; i
< DIM
; i
++)
226 for (int j
= 0; j
< DIM
; j
++)
228 virialScaled_
[i
][j
] = 0;