/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2018, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include <cmath>
#include <algorithm>
+#include <unordered_map>
#include <vector>
#include <gtest/gtest.h>
#include "testutils/refdata.h"
#include "testutils/testasserts.h"
-#include "constrdata.h"
-
namespace gmx
{
namespace test
{
-/*! \brief Test fixture for SHAKE and LINCS
- *
- * The fixture uses three simple and one real-life system for testing.
- * Simple systems include:
- * 1.1. Two atoms, connected with one constrain.
- * 1.2. Three atoms, connected consequently with two constraints
- * 1.3. Four atoms, connected by two independent constraints.
- * 1.4. Three atoms, connected by three constraints in a triangle.
- * 1.5. Four atoms, connected by three consequentive constraints.
- * The real-life system is a Cys-Val-Trp peptide with constraints
- * on:
- * 2.1. Bonds, containing hydrogens (SHAKE and LINCS)
- * 2.2. All bonds (SHAKE and LINCS).
- * 2.3. Angles with hydrogens and all bonds (SHAKE).
- * 2.4. All angles and all bonds (disabled).
+/*! \brief
+ * Constraints test data structure.
*
- * For all systems, the final length of the constraints is tested ageinst the
- * reference values. For systems 2.1 and 2.2, the exact final coordinates of
- * all atoms are also tested against the reference values.
+ * Structure to collect all the necessary data, including system coordinates and topology,
+ * constraints information, etc. The structure can be reset and reused.
*/
-class ConstraintsTest : public ::testing::Test
+struct ConstraintsTestData
{
public:
+ std::string title_; //!< Human-friendly name for a system
+ int nAtom_; //!< Number of atoms
+ gmx_mtop_t mtop_; //!< Topology
+ std::vector<real> masses_; //!< Masses
+ std::vector<real> invmass_; //!< Inverse masses
+ t_commrec cr_; //!< Communication record
+ t_inputrec ir_; //!< Input record (info that usually in .mdp file)
+ t_idef idef_; //!< Local topology
+ t_mdatoms md_; //!< MD atoms
+ gmx_multisim_t ms_; //!< Multisim data
+ t_nrnb nrnb_; //!< Computational time array (normally used in benchmarks)
+
+ real invdt_; //!< Inverse timestep
+ int nflexcon_ = 0; //!< Number of flexible constraints
+ bool computeVirial_; //!< Whether the virial should be computed
+ tensor virialScaled_; //!< Scaled virial
+ tensor virialScaledRef_; //!< Scaled virial (reference values)
+ bool compute_dHdLambda_; //!< If the free energy is computed
+ real dHdLambda_; //!< For free energy computation
+ real dHdLambdaRef_; //!< For free energy computation (reference value)
+
+ std::vector<RVec> x_; //!< Coordinates before the timestep
+ std::vector<RVec> xPrime_; //!< Coordinates after timestep, output for the constraints
+ std::vector<RVec> xPrime0_; //!< Backup for coordinates (for reset)
+ std::vector<RVec> xPrime2_; //!< Intermediate set of coordinates used by LINCS and
+ //!< SHAKE for different purposes
+ std::vector<RVec> v_; //!< Velocities
+ std::vector<RVec> v0_; //!< Backup for velocities (for reset)
- /*! \brief
- * Cleaning up the memory.
- */
- void TearDown() override
- {
- sfree(idef.il[F_CONSTR].iatoms);
- sfree(idef.iparams);
- sfree(x);
- sfree(xprime);
- sfree(xprime2);
- sfree(v);
- gmx_fio_close(log);
- }
+ // Fields to store constraints data for testing
+ std::vector<int> constraints_; //!< Constraints data (type1-i1-j1-type2-i2-j2-...)
+ std::vector<real> constraintsR0_; //!< Target lengths for all constraint types
/*! \brief
- * Method used to initialize atoms, constraints and their parameters.
+ * Constructor for the object with all parameters and variables needed by constraints algorithms.
*
- * This method constructs stubs for all the data structures, required to initialize
- * and apply LINCS and SHAKE constraints. The constraints data is also saved as a field
- * to check if the target distances are achieved after the constraints are applied.
+ * This constructor assembles stubs for all the data structures, required to initialize
+ * and apply LINCS and SHAKE constraints. The coordinates and velocities before constraining
+ * are saved to allow for reset. The constraints data are stored for testing after constraints
+ * were applied.
*
- * \param[in] natom Number of atoms in the system.
- * \param[in] masses Atom masses. Size of this vector should be equal to natom.
- * \param[in] constraints List of constraints, organized in tripples of integers.
- * First integer is the index of type for a constrain, second
- * and third are the indices of constrained atoms. The types
- * of constraints should be sequential but not nesseseraly
- * start from zero (which is the way they normally are in
- * GROMACS).
- * \param[in] constraintsR0 Target values for bond lengths for bonds of each type. The
- * size of this vector should be equal to the total number of
- * unique types in constraints vector.
- * \param[in] firstType The index of first interaction type used for constraints
- * (i.e. the lowest number that appear as type in constraints
- * vector). This is not fail-safe, since technically the
- * numbering of the types may not be sequantial, but (i) this
- * is the way things are in GROMACS, (ii) everything is fine
- * in the examples provided for the test.
- * \param[in] epbc Type of PBC (epbcXYZ / epbcNONE are used in this test,
- * see src/gromacs/pbcutil/pbc.h for details).
- * \param[in] box The PBC box.
- * \param[in] init_t Initial time.
- * \param[in] delta_t Timestep.
- * \param[in] coordinates Initial (unconstrained) coordinates.
+ * \param[in] title Human-friendly name of the system.
+ * \param[in] nAtom Number of atoms in the system.
+ * \param[in] masses Atom masses. Size of this vector should be equal to nAtom.
+ * \param[in] constraints List of constraints, organized in triples of integers.
+ * First integer is the index of type for a constraint, second
+ * and third are the indices of constrained atoms. The types
+ * of constraints should be sequential but not necessarily
+ * start from zero (which is the way they normally are in
+ * GROMACS).
+ * \param[in] constraintsR0 Target values for bond lengths for bonds of each type. The
+ * size of this vector should be equal to the total number of
+ * unique types in constraints vector.
+ * \param[in] computeVirial Whether the virial should be computed.
+ * \param[in] virialScaledRef Reference values for scaled virial tensor.
+ * \param[in] compute_dHdLambda Whether free energy should be computed.
+ * \param[in] dHdLambdaRef Reference value for dHdLambda.
+ * \param[in] initialTime Initial time.
+ * \param[in] timestep Timestep.
+ * \param[in] x Coordinates before integration step.
+ * \param[in] xPrime Coordinates after integration step, but before constraining.
+ * \param[in] v Velocities before constraining.
+ * \param[in] shakeTolerance Target tolerance for SHAKE.
+ * \param[in] shakeUseSOR Use successive over-relaxation method for SHAKE iterations.
+ * The general formula is:
+ * x_n+1 = (1-omega)*x_n + omega*f(x_n),
+ * where omega = 1 if SOR is off and may be < 1 if SOR is on.
+ * \param[in] nLincsIter Number of iterations used to compute the inverse matrix.
+ * \param[in] nProjOrder The order for algorithm that adjusts the direction of the
+ * bond after constraints are applied.
+ * \param[in] lincsWarnAngle The threshold value for the change in bond angle. When
+ * exceeded the program will issue a warning.
*
*/
- void initSystem(int natom, std::vector<real> masses,
- std::vector<int> constraints, std::vector<real> constraintsR0, int firstType,
- int epbc, const matrix box,
- real init_t, real delta_t,
- std::vector<RVec> coordinates)
-
+ ConstraintsTestData(const std::string &title,
+ int nAtom, std::vector<real> masses,
+ std::vector<int> constraints, std::vector<real> constraintsR0,
+ bool computeVirial, tensor virialScaledRef,
+ bool compute_dHdLambda, float dHdLambdaRef,
+ real initialTime, real timestep,
+ const std::vector<RVec> &x, const std::vector<RVec> &xPrime, const std::vector<RVec> &v,
+ real shakeTolerance, gmx_bool shakeUseSOR,
+ int nLincsIter, int nProjOrder, real lincsWarnAngle)
{
+ title_ = title; // Human-friendly name of the system
+ nAtom_ = nAtom; // Number of atoms
- this->n = natom; // Number of atoms
-
- invmass.resize(n); // Vector of inverce masses
+ // Masses of atoms
+ masses_ = masses;
+ invmass_.resize(nAtom); // Vector of inverse masses
- for (int i = 0; i < n; i++)
+ for (int i = 0; i < nAtom; i++)
{
- invmass.at(i) = 1.0/masses.at(i);
+ invmass_[i] = 1.0/masses.at(i);
}
// Saving constraints to check if they are satisfied after algorithm was applied
- this->constraints = constraints; // Constraints indexes (in type-i-j format)
- this->constraintsR0 = constraintsR0; // Euilibrium distances for each type of constraint
- this->firstType = firstType; // The first type index used for constraints
+ constraints_ = constraints; // Constraints indexes (in type-i-j format)
+ constraintsR0_ = constraintsR0; // Equilibrium distances for each type of constraint
- // PBC initialization
- this->epbc = epbc;
- for (int i = 0; i < DIM; i++)
+ invdt_ = 1.0/timestep; // Inverse timestep
+
+ // Communication record
+ cr_.nnodes = 1;
+ cr_.dd = nullptr;
+
+ // Multisim data
+ ms_.sim = 0;
+ ms_.nsim = 1;
+
+ // Input record - data that usually comes from configuration file (.mdp)
+ ir_.efep = 0;
+ ir_.init_t = initialTime;
+ ir_.delta_t = timestep;
+ ir_.eI = 0;
+
+ // MD atoms data
+ md_.nMassPerturbed = 0;
+ md_.lambda = 0.0;
+ md_.invmass = invmass_.data();
+ md_.nr = nAtom;
+ md_.homenr = nAtom;
+
+ // Virial evaluation
+ computeVirial_ = computeVirial;
+ if (computeVirial)
{
- for (int j = 0; j < DIM; j++)
+ for (int i = 0; i < DIM; i++)
{
- this->box[i][j] = box[i][j]; // Periodic box
+ for (int j = 0; j < DIM; j++)
+ {
+ virialScaled_[i][j] = 0;
+ virialScaledRef_[i][j] = virialScaledRef[i][j];
+ }
}
}
- set_pbc(&pbc, epbc, box);
-
- this->invdt = 1.0/delta_t; // Inverse timestep
-
- // Communication record
- cr.nnodes = 1;
- cr.dd = nullptr;
- // Input record - data that usualy comes from configurtion file (.mdp)
- ir.efep = 0;
- ir.init_t = init_t;
- ir.delta_t = delta_t;
- ir.eI = 0;
- // MD atoms data
- md.nMassPerturbed = 0;
- md.lambda = 0.0;
- md.invmass = invmass.data();
- md.nr = n;
- md.homenr = n;
+ // Free energy evaluation
+ compute_dHdLambda_ = compute_dHdLambda;
+ dHdLambda_ = 0;
+ if (compute_dHdLambda_)
+ {
+ ir_.efep = efepYES;
+ dHdLambdaRef_ = dHdLambdaRef;
+ }
+ else
+ {
+ ir_.efep = efepNO;
+ dHdLambdaRef_ = 0;
+ }
- // Constraints and their parameters in old data format (local topology)
+ // Constraints and their parameters (local topology)
for (int i = 0; i < F_NRE; i++)
{
- idef.il[i].nr = 0;
+ idef_.il[i].nr = 0;
}
- idef.il[F_CONSTR].nr = constraints.size();
+ idef_.il[F_CONSTR].nr = constraints.size();
- snew(idef.il[F_CONSTR].iatoms, constraints.size());
+ snew(idef_.il[F_CONSTR].iatoms, constraints.size());
int maxType = 0;
for (unsigned i = 0; i < constraints.size(); i++)
{
maxType = constraints.at(i);
}
}
- idef.il[F_CONSTR].iatoms[i] = constraints.at(i);
+ idef_.il[F_CONSTR].iatoms[i] = constraints.at(i);
}
- snew(idef.iparams, maxType + 1);
+ snew(idef_.iparams, maxType + 1);
for (unsigned i = 0; i < constraints.size()/3; i++)
{
- idef.iparams[constraints.at(3*i)].constr.dA = constraintsR0.at(constraints.at(3*i) - firstType);
- idef.iparams[constraints.at(3*i)].constr.dB = constraintsR0.at(constraints.at(3*i) - firstType);
+ idef_.iparams[constraints.at(3*i)].constr.dA = constraintsR0.at(constraints.at(3*i));
+ idef_.iparams[constraints.at(3*i)].constr.dB = constraintsR0.at(constraints.at(3*i));
}
-
- // Constraints and their parameters in new data format (global topology)
- InteractionList ilist;
- ilist.iatoms.resize(constraints.size());
+ // Constraints and their parameters (global topology)
+ InteractionList interactionList;
+ interactionList.iatoms.resize(constraints.size());
for (unsigned i = 0; i < constraints.size(); i++)
{
- ilist.iatoms.at(i) = constraints.at(i);
+ interactionList.iatoms.at(i) = constraints.at(i);
}
- InteractionList ilistEmpty;
- ilistEmpty.iatoms.resize(0);
-
- gmx_moltype_t moltype;
- moltype.atoms.nr = n;
- moltype.ilist.at(F_CONSTR) = ilist;
- moltype.ilist.at(F_CONSTRNC) = ilistEmpty;
- mtop.moltype.push_back(moltype);
-
- gmx_molblock_t molblock;
- molblock.type = 0;
- molblock.nmol = 1;
- mtop.molblock.push_back(molblock);
-
- mtop.natoms = n;
- mtop.ffparams.iparams.resize(maxType + 1);
+ InteractionList interactionListEmpty;
+ interactionListEmpty.iatoms.resize(0);
+
+ gmx_moltype_t molType;
+ molType.atoms.nr = nAtom;
+ molType.ilist.at(F_CONSTR) = interactionList;
+ molType.ilist.at(F_CONSTRNC) = interactionListEmpty;
+ mtop_.moltype.push_back(molType);
+
+ gmx_molblock_t molBlock;
+ molBlock.type = 0;
+ molBlock.nmol = 1;
+ mtop_.molblock.push_back(molBlock);
+
+ mtop_.natoms = nAtom;
+ mtop_.ffparams.iparams.resize(maxType + 1);
for (int i = 0; i <= maxType; i++)
{
- mtop.ffparams.iparams.at(i) = idef.iparams[i];
+ mtop_.ffparams.iparams.at(i) = idef_.iparams[i];
}
- mtop.bIntermolecularInteractions = false;
+ mtop_.bIntermolecularInteractions = false;
+
+ // Coordinates and velocities
+ x_ = x;
+ xPrime_ = xPrime;
+ xPrime0_ = xPrime;
+ xPrime2_ = xPrime;
- // Log file may be here. Redirect it somewhere usefull?
- log = gmx_fio_open("constraintstest.log", "w");
+ v_ = v;
+ v0_ = v;
- // Set the coordinates in convinient format
- snew(x, n);
- snew(xprime, n);
- snew(xprime2, n);
+ // SHAKE-specific parameters
+ ir_.shake_tol = shakeTolerance;
+ ir_.bShakeSOR = shakeUseSOR;
- for (int i = 0; i < n; i++)
+ // LINCS-specific parameters
+ ir_.nLincsIter = nLincsIter;
+ ir_.nProjOrder = nProjOrder;
+ ir_.LincsWarnAngle = lincsWarnAngle;
+ }
+
+ /*! \brief
+ * Reset the data structure so it can be reused.
+ *
+ * Set the coordinates and velocities back to their values before
+ * constraining. The scaled virial tensor and dHdLambda are zeroed.
+ *
+ */
+ void reset()
+ {
+ xPrime_ = xPrime0_;
+ xPrime2_ = xPrime0_;
+ v_ = v0_;
+
+ if (computeVirial_)
{
- for (int j = 0; j < DIM; j++)
+ for (int i = 0; i < DIM; i++)
{
- x[i][j] = coordinates.at(i)[j];
- xprime[i][j] = coordinates.at(i)[j];
+ for (int j = 0; j < DIM; j++)
+ {
+ virialScaled_[i][j] = 0;
+ }
}
}
+ dHdLambda_ = 0;
+ }
+
+ /*! \brief
+ * Cleaning up the memory.
+ */
+ ~ConstraintsTestData()
+ {
+ sfree(idef_.il[F_CONSTR].iatoms);
+ sfree(idef_.iparams);
+ }
- // Velocities (not used)
- snew(v, n);
+};
+/*! \brief The two-dimensional parameter space for test.
+ *
+ * The test will run for all possible combinations of accessible
+ * values of the:
+ * 1. PBC setup ("PBCNONE" or "PBCXYZ")
+ * 2. The algorithm ("SHAKE" or "LINCS").
+ */
+typedef std::tuple<std::string, std::string> ConstraintsTestParameters;
+/*! \brief Test fixture for constraints.
+ *
+ * The fixture uses following test systems:
+ * 1. Two atoms, connected with one constraint (e.g. NH).
+ * 2. Three atoms, connected consequently with two constraints (e.g. CH2).
+ * 3. Three atoms, constrained to the fourth atom (e.g. CH3).
+ * 4. Four atoms, connected by two independent constraints.
+ * 5. Three atoms, connected by three constraints in a triangle
+ * (e.g. H2O with constrained H-O-H angle).
+ * 6. Four atoms, connected by three consequential constraints.
+ *
+ * For all systems, the final lengths of the constraints are tested against the
+ * reference values, the direction of each constraint is checked.
+ * Test also verifies that the center of mass has not been
+ * shifted by the constraints and that its velocity has not changed.
+ * For some systems, the value for scaled virial tensor is checked against
+ * pre-computed data.
+ */
+class ConstraintsTest : public ::testing::TestWithParam<ConstraintsTestParameters>
+{
+ public:
+ //! PBC setups
+ std::unordered_map <std::string, t_pbc> pbcs_;
+ //! Algorithms (SHAKE and LINCS)
+ std::unordered_map <std::string, void(*)(ConstraintsTestData *testData, t_pbc pbc)> algorithms_;
+
+ /*! \brief Test setup function.
+ *
+ * Setting up the pbcs and algorithms. Note, that corresponding string keywords
+ * have to be explicitly added at the end of this file when the tests are called.
+ *
+ */
+ void SetUp() override
+ {
+
+ //
+ // PBC initialization
+ //
+ t_pbc pbc;
+
+ // Infinitely small box
+ matrix boxNone = { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} };
+ set_pbc(&pbc, epbcNONE, boxNone);
+ pbcs_["PBCNone"] = pbc;
+
+ // Rectangular box
+ matrix boxXyz = { {10.0, 0.0, 0.0}, {0.0, 20.0, 0.0}, {0.0, 0.0, 15.0} };
+ set_pbc(&pbc, epbcXYZ, boxXyz);
+ pbcs_["PBCXYZ"] = pbc;
+
+ //
+ // Algorithms
+ //
+ // SHAKE
+ algorithms_["SHAKE"] = applyShake;
+ // LINCS
+ algorithms_["LINCS"] = applyLincs;
+
+ }
+
+ /*! \brief
+ * Initialize and apply SHAKE constraints.
+ *
+ * \param[in] testData Test data structure.
+ * \param[in] pbc Periodic boundary data (not used in SHAKE).
+ */
+ static void applyShake(ConstraintsTestData *testData, t_pbc gmx_unused pbc)
+ {
+ shakedata* shaked = shake_init();
+ make_shake_sblock_serial(shaked, &testData->idef_, testData->md_);
+ bool success = constrain_shake(
+ nullptr,
+ shaked,
+ testData->invmass_.data(),
+ testData->idef_,
+ testData->ir_,
+ as_rvec_array(testData->x_.data()),
+ as_rvec_array(testData->xPrime_.data()),
+ as_rvec_array(testData->xPrime2_.data()),
+ &testData->nrnb_,
+ testData->md_.lambda,
+ &testData->dHdLambda_,
+ testData->invdt_,
+ as_rvec_array(testData->v_.data()),
+ testData->computeVirial_,
+ testData->virialScaled_,
+ false,
+ gmx::ConstraintVariable::Positions);
+ EXPECT_TRUE(success) << "Test failed with a false return value in SHAKE.";
+ done_shake(shaked);
}
/*! \brief
- * This method initializes and applies LINCS constraints.
+ * Initialize and apply LINCS constraints.
*
- * \param[in] nLincsIter Number of iterations used to compute the inverse matrix.
- * \param[in] nProjOrder The order for algorithm that adjusts the direction of the
- * bond after constraints are applied.
- * \param[in] LincsWarnAngle The value for the change in bond angle after which the
- * program will issue a warning.
+ * \param[in] testData Test data structure.
+ * \param[in] pbc Periodic boundary data.
*/
- void applyConstraintsLincs(int nLincsIter, int nProjOrder, real LincsWarnAngle)
+ static void applyLincs(ConstraintsTestData *testData, t_pbc pbc)
{
Lincs *lincsd;
int maxwarn = 100;
int warncount_lincs = 0;
- ir.nLincsIter = nLincsIter;
- ir.nProjOrder = nProjOrder;
- ir.LincsWarnAngle = LincsWarnAngle;
gmx_omp_nthreads_set(emntLINCS, 1);
// Make blocka structure for faster LINCS setup
std::vector<t_blocka> at2con_mt;
- at2con_mt.reserve(mtop.moltype.size());
- for (const gmx_moltype_t &moltype : mtop.moltype)
+ at2con_mt.reserve(testData->mtop_.moltype.size());
+ for (const gmx_moltype_t &moltype : testData->mtop_.moltype)
{
// This function is in constr.cpp
at2con_mt.push_back(make_at2con(moltype,
- mtop.ffparams.iparams,
- flexibleConstraintTreatment(EI_DYNAMICS(ir.eI))));
+ testData->mtop_.ffparams.iparams,
+ flexibleConstraintTreatment(EI_DYNAMICS(testData->ir_.eI))));
}
- // Initialie LINCS
- lincsd = init_lincs(gmx_fio_getfp(log), mtop,
- nflexcon, at2con_mt,
+ // Initialize LINCS
+ lincsd = init_lincs(nullptr, testData->mtop_,
+ testData->nflexcon_, at2con_mt,
false,
- ir.nLincsIter, ir.nProjOrder);
- set_lincs(idef, md, EI_DYNAMICS(ir.eI), &cr, lincsd);
+ testData->ir_.nLincsIter, testData->ir_.nProjOrder);
+ set_lincs(testData->idef_, testData->md_, EI_DYNAMICS(testData->ir_.eI), &testData->cr_, lincsd);
// Evaluate constraints
- bool bOK = constrain_lincs(false, ir, 0, lincsd, md,
- &cr,
- nullptr,
- x, xprime, xprime2,
- box, &pbc, md.lambda, dvdlambda,
- invdt, v,
- computeVirial, vir_r_m_dr,
- gmx::ConstraintVariable::Positions, &nrnb,
- maxwarn, &warncount_lincs);
- EXPECT_TRUE(bOK) << "Something went wrong: LINCS returned false value.";
+ bool sucess = constrain_lincs(false, testData->ir_, 0, lincsd, testData->md_,
+ &testData->cr_,
+ &testData->ms_,
+ as_rvec_array(testData->x_.data()),
+ as_rvec_array(testData->xPrime_.data()),
+ as_rvec_array(testData->xPrime2_.data()),
+ pbc.box, &pbc,
+ testData->md_.lambda, &testData->dHdLambda_,
+ testData->invdt_,
+ as_rvec_array(testData->v_.data()),
+ testData->computeVirial_, testData->virialScaled_,
+ gmx::ConstraintVariable::Positions,
+ &testData->nrnb_,
+ maxwarn, &warncount_lincs);
+ EXPECT_TRUE(sucess) << "Test failed with a false return value in LINCS.";
EXPECT_EQ(warncount_lincs, 0) << "There were warnings in LINCS.";
- for (unsigned int i = 0; i < mtop.moltype.size(); i++)
+ for (unsigned int i = 0; i < testData->mtop_.moltype.size(); i++)
{
sfree(at2con_mt.at(i).index);
sfree(at2con_mt.at(i).a);
}
/*! \brief
- * Initialize and apply SHAKE constraints.
+ * The test on the final length of constrained bonds.
*
- * \param[in] shake_tol Target tolerance for SHAKE.
- * \param[in] bShakeSOR Use successive over-relaxation method for SHAKE iterative process.
- * The general formula is:
- * x_n+1 = (1-omega)*x_n + omega*f(x_n),
- * where omega = 1 if SOR is off and may be < 1 if SOR is on.
+ * Goes through all the constraints and checks if the final length of all the constraints is equal
+ * to the target length with provided tolerance.
*
+ * \param[in] tolerance Allowed tolerance in final lengths.
+ * \param[in] testData Test data structure.
+ * \param[in] pbc Periodic boundary data.
*/
- void applyConstraintsShake(real shake_tol, gmx_bool bShakeSOR)
+ void checkConstrainsLength(FloatingPointTolerance tolerance, const ConstraintsTestData &testData, t_pbc pbc)
{
- ir.shake_tol = shake_tol;
- ir.bShakeSOR = bShakeSOR;
-
- shakedata* shaked = shake_init();
- make_shake_sblock_serial(shaked, &idef, md);
- bool bOK = constrain_shake(
- gmx_fio_getfp(log),
- shaked,
- invmass.data(),
- idef,
- ir,
- x,
- xprime,
- xprime2,
- &nrnb,
- md.lambda,
- dvdlambda,
- invdt,
- v,
- computeVirial,
- vir_r_m_dr,
- false,
- gmx::ConstraintVariable::Positions);
- EXPECT_TRUE(bOK) << "Something went wrong: SHAKE returned false value.";
- done_shake(shaked); // Not yet implemented. We will have memory leaks without this.
+ // Test if all the constraints are satisfied
+ for (unsigned c = 0; c < testData.constraints_.size()/3; c++)
+ {
+ real r0 = testData.constraintsR0_.at(testData.constraints_.at(3*c));
+ int i = testData.constraints_.at(3*c + 1);
+ int j = testData.constraints_.at(3*c + 2);
+ RVec xij0, xij1;
+ real d0, d1;
+ if (pbc.ePBC == epbcXYZ)
+ {
+ pbc_dx_aiuc(&pbc, testData.x_[i], testData.x_[j], xij0);
+ pbc_dx_aiuc(&pbc, testData.xPrime_[i], testData.xPrime_[j], xij1);
+ }
+ else
+ {
+ rvec_sub(testData.x_[i], testData.x_[j], xij0);
+ rvec_sub(testData.xPrime_[i], testData.xPrime_[j], xij1);
+ }
+ d0 = norm(xij0);
+ d1 = norm(xij1);
+ EXPECT_REAL_EQ_TOL(r0, d1, tolerance) << gmx::formatString(
+ "rij = %f, which is not equal to r0 = %f for constraint #%u, between atoms %d and %d"
+ " (before constraining rij was %f).", d1, r0, c, i, j, d0);
+ }
}
/*! \brief
* The test on the final length of constrained bonds.
*
- * Goes through all the constraints and checks if the final length of all the constraints is equal
- * to the target lenght with provided tolerance.
+ * Goes through all the constraints and checks if the direction of constraint has not changed
+ * by the algorithm (i.e. the constraints algorithm arrived to the solution that is closest
+ * to the initial system conformation).
*
- * \param[in] tolerance Allowed tolerance in final lengths.
+ * \param[in] testData Test data structure.
+ * \param[in] pbc Periodic boundary data.
*/
-
- void checkConstrained(FloatingPointTolerance tolerance)
+ void checkConstrainsDirection(const ConstraintsTestData &testData, t_pbc pbc)
{
- // Test if all the constraints are satisfied
- for (unsigned c = 0; c < constraints.size()/3; c++)
+ for (unsigned c = 0; c < testData.constraints_.size()/3; c++)
{
- real r0 = constraintsR0.at(constraints.at(3*c) - firstType);
- int i = constraints.at(3*c + 1);
- int j = constraints.at(3*c + 2);
- rvec v0, v1;
- real d0, d1;
- if (this->epbc == epbcXYZ)
+ int i = testData.constraints_.at(3*c + 1);
+ int j = testData.constraints_.at(3*c + 2);
+ RVec xij0, xij1;
+ if (pbc.ePBC == epbcXYZ)
{
- pbc_dx_aiuc(&pbc, x[i], x[j], v0);
- pbc_dx_aiuc(&pbc, xprime[i], xprime[j], v1);
+ pbc_dx_aiuc(&pbc, testData.x_[i], testData.x_[j], xij0);
+ pbc_dx_aiuc(&pbc, testData.xPrime_[i], testData.xPrime_[j], xij1);
}
else
{
- rvec_sub(x[i], x[j], v0);
- rvec_sub(xprime[i], xprime[j], v1);
+ rvec_sub(testData.x_[i], testData.x_[j], xij0);
+ rvec_sub(testData.xPrime_[i], testData.xPrime_[j], xij1);
}
- d0 = norm(v0);
- d1 = norm(v1);
- EXPECT_REAL_EQ_TOL(r0, d1, tolerance) << "rij = " << d1 << ", which is not equal to r0 = " << r0
- << " for constraint #" << c << ", between atoms " << i << " and " << j
- << " (before lincs rij was " << d0 << ").";
+
+ real dot = xij0.dot(xij1);
+
+ EXPECT_GE(dot, 0.0) << gmx::formatString(
+ "The constraint %u changed direction. Constraining algorithm might have returned the wrong root "
+ "of the constraints equation.", c);
+
+ }
+ }
+
+ /*! \brief
+ * The test on the coordinates of the center of the mass (COM) of the system.
+ *
+ * Checks if the center of mass has not been shifted by the constraints. Note,
+ * that this test does not take into account the periodic boundary conditions.
+ * Hence it will not work should the constraints decide to move atoms across
+ * PBC borders.
+ *
+ * \param[in] tolerance Allowed tolerance in COM coordinates.
+ * \param[in] testData Test data structure.
+ */
+ void checkCOMCoordinates(FloatingPointTolerance tolerance, const ConstraintsTestData &testData)
+ {
+
+ RVec comPrime0({0.0, 0.0, 0.0});
+ RVec comPrime({0.0, 0.0, 0.0});
+ for (int i = 0; i < testData.nAtom_; i++)
+ {
+ comPrime0 += testData.masses_[i]*testData.xPrime0_[i];
+ comPrime += testData.masses_[i]*testData.xPrime_[i];
+ }
+
+ comPrime0 /= testData.nAtom_;
+ comPrime /= testData.nAtom_;
+
+ EXPECT_REAL_EQ_TOL(comPrime[XX], comPrime0[XX], tolerance)
+ << "Center of mass was shifted by constraints in x-direction.";
+ EXPECT_REAL_EQ_TOL(comPrime[YY], comPrime0[YY], tolerance)
+ << "Center of mass was shifted by constraints in y-direction.";
+ EXPECT_REAL_EQ_TOL(comPrime[ZZ], comPrime0[ZZ], tolerance)
+ << "Center of mass was shifted by constraints in z-direction.";
+
+ }
+
+ /*! \brief
+ * The test on the velocity of the center of the mass (COM) of the system.
+ *
+ * Checks if the velocity of the center of mass has not changed.
+ *
+ * \param[in] tolerance Allowed tolerance in COM velocity components.
+ * \param[in] testData Test data structure.
+ */
+ void checkCOMVelocity(FloatingPointTolerance tolerance, const ConstraintsTestData &testData)
+ {
+
+ RVec comV0({0.0, 0.0, 0.0});
+ RVec comV({0.0, 0.0, 0.0});
+ for (int i = 0; i < testData.nAtom_; i++)
+ {
+ comV0 += testData.masses_[i]*testData.v0_[i];
+ comV += testData.masses_[i]*testData.v_[i];
}
+ comV0 /= testData.nAtom_;
+ comV /= testData.nAtom_;
+
+ EXPECT_REAL_EQ_TOL(comV[XX], comV0[XX], tolerance)
+ << "Velocity of the center of mass in x-direction has been changed by constraints.";
+ EXPECT_REAL_EQ_TOL(comV[YY], comV0[YY], tolerance)
+ << "Velocity of the center of mass in y-direction has been changed by constraints.";
+ EXPECT_REAL_EQ_TOL(comV[ZZ], comV0[ZZ], tolerance)
+ << "Velocity of the center of mass in z-direction has been changed by constraints.";
}
/*! \brief
- * The test on the final coordinates.
+ * The test on the final coordinates (not used).
*
* Goes through all atoms and checks if the final positions correspond to the
* provided reference set of coordinates.
*
- * \param[in] finalCoordinates The reference set of coordinates.
- * \param[in] finalCoordinatesTolerance Tolerance for the coordinaes test.
+ * \param[in] xPrimeRef The reference set of coordinates.
+ * \param[in] tolerance Tolerance for the coordinates test.
+ * \param[in] testData Test data structure.
*/
- void checkFinalCoordinates(std::vector<RVec> finalCoordinates, FloatingPointTolerance finalCoordinatesTolerance)
+ void checkFinalCoordinates(std::vector<RVec> xPrimeRef, FloatingPointTolerance tolerance,
+ const ConstraintsTestData &testData)
{
- for (int i = 0; i < n; i++)
+ for (int i = 0; i < testData.nAtom_; i++)
{
for (int d = 0; d < DIM; d++)
{
- EXPECT_REAL_EQ_TOL(finalCoordinates.at(i)[d], xprime[i][d], finalCoordinatesTolerance) <<
- "Coordinates after constrains were applied differ from these in the reference set for atom #" << i << ".";
+ EXPECT_REAL_EQ_TOL(xPrimeRef.at(i)[d], testData.xPrime_[i][d], tolerance) << gmx::formatString(
+ "Coordinates after constrains were applied differ from these in the reference set for atom #%d.", i);
}
}
}
+ /*! \brief
+ * The test of virial tensor.
+ *
+ * Checks if the values in the scaled virial tensor are equal to pre-computed values.
+ *
+ * \param[in] tolerance Tolerance for the tensor values.
+ * \param[in] testData Test data structure.
+ */
+ void checkVirialTensor(FloatingPointTolerance tolerance, const ConstraintsTestData &testData)
+ {
+ for (int i = 0; i < DIM; i++)
+ {
+ for (int j = 0; j < DIM; j++)
+ {
+ EXPECT_REAL_EQ_TOL(testData.virialScaledRef_[i][j], testData.virialScaled_[i][j],
+ tolerance) << gmx::formatString(
+ "Values in virial tensor at [%d][%d] are not within the tolerance from reference value.", i, j);
+ }
+ }
+ }
- private:
- int n; // Number of atoms
- gmx_mtop_t mtop; // Topology
- std::vector<real> invmass; // Inverse masses
- int epbc; // PBC used (epbcNONE/epbcXYZ)
- matrix box; // PBC box
- t_pbc pbc; // PBC object
- t_commrec cr; // Communication record
- t_inputrec ir; // Input record (info that usually in .mdp file)
- t_idef idef; // Local topology
- t_mdatoms md; // MD atoms
- t_nrnb nrnb;
-
- real invdt = 1.0/0.001; // Inverse timestep
- int nflexcon = 0; // Number of flexible constraints
- real *dvdlambda = nullptr; // Free energy computation stub (not tested)
- bool computeVirial = false; // If the virials should be conputed (not tested)
- tensor vir_r_m_dr; // Virials stuff
- rvec *x; // Coordinates before constraints are applied
- rvec *xprime; // Coordinates after constraints are applied
- rvec *xprime2; // Intermediate set of coordinates used by LINCS and
- // SHAKE for different purposes
- rvec *v; // Velocities (can also be constrained, this is not tested)
- t_fileio *log; // log file (now set to "constraintstest.log")
+ /*! \brief
+ * The test for FEP (not used).
+ *
+ * Checks if the value of dH/dLambda is equal to the reference value.
+ * \todo Add tests for dHdLambda values.
+ *
+ * \param[in] dHdLambdaRef Reference value.
+ * \param[in] tolerance Tolerance.
+ * \param[in] testData Test data structure.
+ */
+ void checkFEP(const real dHdLambdaRef, FloatingPointTolerance tolerance,
+ const ConstraintsTestData &testData)
+ {
+ EXPECT_REAL_EQ_TOL(dHdLambdaRef, testData.dHdLambda_, tolerance) <<
+ "Computed value for dV/dLambda is not equal to the reference value. ";
+ }
- // Fields to store constraints data for testing
- int firstType; // First interaction type used for constraints
- std::vector<int> constraints; // Constraints data (type1-i1-j1-type2-i2-j2-...)
- std::vector<real> constraintsR0; // Target lengths for al the constrain
-};
-/*
- *
- * Tests that the interatomic distance along constraints correspond to the reference distances.
- * Performed on basic systems of up to four constaints.
- *
- */
-/*
- * Simple bond test for SHAKE
- */
-TEST_F(ConstraintsTest, ShakeOneBond)
-{
+};
+
+TEST_P(ConstraintsTest, SingleConstraint){
+ std::string title = "one constraint (e.g. OH)";
+ int nAtom = 2;
+
+ std::vector<real> masses = {1.0, 12.0};
+ std::vector<int> constraints = {0, 0, 1};
+ std::vector<real> constraintsR0 = {0.1};
+
+ real oneTenthOverSqrtTwo = 0.1_real / std::sqrt(2.0_real);
+
+ std::vector<RVec> x = {{ 0.0, oneTenthOverSqrtTwo, 0.0 },
+ { oneTenthOverSqrtTwo, 0.0, 0.0 }};
+ std::vector<RVec> xPrime = {{ 0.01, 0.08, 0.01 },
+ { 0.06, 0.01, -0.01 }};
+ std::vector<RVec> v = {{ 1.0, 2.0, 3.0 },
+ { 3.0, 2.0, 1.0 }};
+
+ tensor virialScaledRef = {{-5.58e-04, 5.58e-04, 0.00e+00 },
+ { 5.58e-04, -5.58e-04, 0.00e+00 },
+ { 0.00e+00, 0.00e+00, 0.00e+00 } };
+
+ real shakeTolerance = 0.0001;
+ gmx_bool shakeUseSOR = false;
+
+ int lincsNIter = 1;
+ int lincsNProjOrder = 4;
+ real lincsWarnAngle = 30.0;
+
+ std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>
+ (title, nAtom, masses,
+ constraints, constraintsR0,
+ true, virialScaledRef,
+ false, 0,
+ real(0.0), real(0.001),
+ x, xPrime, v,
+ shakeTolerance, shakeUseSOR,
+ lincsNIter, lincsNProjOrder, lincsWarnAngle);
+ std::string pbcName;
+ std::string algorithmName;
+ std::tie(pbcName, algorithmName) = GetParam();
+ t_pbc pbc = pbcs_.at(pbcName);
+
+ // Apply constraints
+ algorithms_.at(algorithmName)(testData.get(), pbc);
+
+ checkConstrainsLength(absoluteTolerance(0.0002), *testData, pbc);
+ checkConstrainsDirection(*testData, pbc);
+ checkCOMCoordinates(absoluteTolerance(0.0001), *testData);
+ checkCOMVelocity(absoluteTolerance(0.0001), *testData);
+
+ checkVirialTensor(absoluteTolerance(0.0001), *testData);
- FloatingPointTolerance tolerance = relativeToleranceAsFloatingPoint(0.1, 0.0001);
- initSystem(2, c_oneBondMasses,
- c_oneBondConstraints, c_oneBondConstraintsR0, c_oneBondConstraintsFirstType,
- epbcNONE, c_infinitesimalBox,
- real(0.0), real(0.001), c_oneBondCoordinates);
- applyConstraintsShake(0.0001, false);
- checkConstrained(tolerance);
}
-/*
- * Simple bond test for LINCS
- */
-TEST_F(ConstraintsTest, LincsOneBond)
-{
+TEST_P(ConstraintsTest, TwoDisjointConstraints){
- FloatingPointTolerance tolerance = relativeToleranceAsFloatingPoint(0.1, 0.0000002);
- initSystem(2, c_oneBondMasses,
- c_oneBondConstraints, c_oneBondConstraintsR0, c_oneBondConstraintsFirstType,
- epbcNONE, c_infinitesimalBox,
- real(0.0), real(0.001), c_oneBondCoordinates);
- applyConstraintsLincs(1, 4, real(30.0));
- checkConstrained(tolerance);
-}
+ std::string title = "two disjoint constraints";
+ int nAtom = 4;
+ std::vector<real> masses = {0.5, 1.0/3.0, 0.25, 1.0};
+ std::vector<int> constraints = {0, 0, 1, 1, 2, 3};
+ std::vector<real> constraintsR0 = {2.0, 1.0};
-/*
- * Two disjoint bonds test for SHAKE
- */
-TEST_F(ConstraintsTest, ShakeTwoDisjointBonds)
-{
- FloatingPointTolerance tolerance = relativeToleranceAsFloatingPoint(0.1, 0.000001);
- initSystem(4, c_twoDJBondsMasses,
- c_twoDJBondsConstraints, c_twoDJBondsConstraintsR0, c_twoDJBondsConstraintsFirstType,
- epbcNONE, c_infinitesimalBox,
- real(0.0), real(0.001), c_twoDJBondsCoordinates);
+ std::vector<RVec> x = {{ 2.50, -3.10, 15.70 },
+ { 0.51, -3.02, 15.55 },
+ { -0.50, -3.00, 15.20 },
+ { -1.51, -2.95, 15.05 }};
- applyConstraintsShake(0.000001, true);
- checkConstrained(tolerance);
-}
+ std::vector<RVec> xPrime = {{ 2.50, -3.10, 15.70 },
+ { 0.51, -3.02, 15.55 },
+ { -0.50, -3.00, 15.20 },
+ { -1.51, -2.95, 15.05 }};
-/*
- * Two disjoint bonds test for LINCS
- */
-TEST_F(ConstraintsTest, LincsTwoDisjointBonds)
-{
+ std::vector<RVec> v = {{ 0.0, 1.0, 0.0 },
+ { 1.0, 0.0, 0.0 },
+ { 0.0, 0.0, 1.0 },
+ { 0.0, 0.0, 0.0 }};
- FloatingPointTolerance tolerance = relativeToleranceAsFloatingPoint(0.1, 0.0000002);
- initSystem(4, c_twoDJBondsMasses,
- c_twoDJBondsConstraints, c_twoDJBondsConstraintsR0, c_twoDJBondsConstraintsFirstType,
- epbcNONE, c_infinitesimalBox,
- real(0.0), real(0.001), c_twoDJBondsCoordinates);
+ tensor virialScaledRef = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
- applyConstraintsLincs(1, 4, real(30.0));
- checkConstrained(tolerance);
-}
+ real shakeTolerance = 0.0001;
+ gmx_bool shakeUseSOR = false;
-/*
- * Two consequentive constraints test for SHAKE.
- */
-TEST_F(ConstraintsTest, ShakeTwoBonds)
-{
+ int lincsNIter = 1;
+ int lincsNProjOrder = 4;
+ real lincsWarnAngle = 30.0;
- FloatingPointTolerance tolerance = relativeToleranceAsFloatingPoint(0.1, 0.0001);
- initSystem(3, c_twoBondsMasses,
- c_twoBondsConstraints, c_twoBondsConstraintsR0, c_twoBondsConstraintsFirstType,
- epbcNONE, c_infinitesimalBox,
- real(0.0), real(0.001), c_twoBondsCoordinates);
+ std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>
+ (title, nAtom, masses,
+ constraints, constraintsR0,
+ false, virialScaledRef,
+ false, 0,
+ real(0.0), real(0.001),
+ x, xPrime, v,
+ shakeTolerance, shakeUseSOR,
+ lincsNIter, lincsNProjOrder, lincsWarnAngle);
- applyConstraintsShake(0.0001, true);
- checkConstrained(tolerance);
-}
+ std::string pbcName;
+ std::string algorithmName;
+ std::tie(pbcName, algorithmName) = GetParam();
+ t_pbc pbc = pbcs_.at(pbcName);
-/*
- * Two consequentive constraints test for LINCS.
- */
-TEST_F(ConstraintsTest, LincsTwoBonds)
-{
+ // Apply constraints
+ algorithms_.at(algorithmName)(testData.get(), pbc);
- FloatingPointTolerance tolerance = relativeToleranceAsFloatingPoint(0.1, 0.000001);
- initSystem(3, c_twoBondsMasses,
- c_twoBondsConstraints, c_twoBondsConstraintsR0, c_twoBondsConstraintsFirstType,
- epbcNONE, c_infinitesimalBox,
- real(0.0), real(0.001), c_twoBondsCoordinates);
+ checkConstrainsLength(absoluteTolerance(0.0002), *testData, pbc);
+ checkConstrainsDirection(*testData, pbc);
+ checkCOMCoordinates(absoluteTolerance(0.0001), *testData);
+ checkCOMVelocity(absoluteTolerance(0.0001), *testData);
- applyConstraintsLincs(1, 4, real(30.0));
- checkConstrained(tolerance);
}
-/*
- * Triangle of constraints test for SHAKE
- */
-TEST_F(ConstraintsTest, ShakeTriangleOfBonds)
-{
+TEST_P(ConstraintsTest, ThreeSequentialConstraints){
- FloatingPointTolerance tolerance = relativeToleranceAsFloatingPoint(0.1, 0.0001);
- initSystem(3, c_triangleMasses,
- c_triangleConstraints, c_triangleConstraintsR0, c_triangleConstraintsFirstType,
- epbcNONE, c_infinitesimalBox,
- real(0.0), real(0.001), c_triangleCoordinates);
- applyConstraintsShake(0.0001, true);
- checkConstrained(tolerance);
-}
+ std::string title = "three atoms, connected longitudinally (e.g. CH2)";
+ int nAtom = 3;
+ std::vector<real> masses = {1.0, 12.0, 16.0 };
+ std::vector<int> constraints = {0, 0, 1, 1, 1, 2};
+ std::vector<real> constraintsR0 = {0.1, 0.2};
-/*
- * Triangle of constraints test for LINCS
- */
-TEST_F(ConstraintsTest, LincsTriangleOfBonds)
-{
+ real oneTenthOverSqrtTwo = 0.1_real / std::sqrt(2.0_real);
+ real twoTenthsOverSqrtThree = 0.2_real / std::sqrt(3.0_real);
- FloatingPointTolerance tolerance = relativeToleranceAsFloatingPoint(0.1, 0.0001);
- initSystem(3, c_triangleMasses,
- c_triangleConstraints, c_triangleConstraintsR0, c_triangleConstraintsFirstType,
- epbcNONE, c_infinitesimalBox,
- real(0.0), real(0.001), c_triangleCoordinates);
- applyConstraintsLincs(4, 4, real(30.0));
- checkConstrained(tolerance);
-}
+ std::vector<RVec> x = {{ oneTenthOverSqrtTwo, oneTenthOverSqrtTwo, 0.0 },
+ { 0.0, 0.0, 0.0 },
+ { twoTenthsOverSqrtThree, twoTenthsOverSqrtThree, twoTenthsOverSqrtThree }};
-/*
- * Three consecutive constraints test for SHAKE.
- */
-TEST_F(ConstraintsTest, ShakeThreeConsequativeConstraints)
-{
+ std::vector<RVec> xPrime = {{ 0.08, 0.07, 0.01 },
+ { -0.02, 0.01, -0.02 },
+ { 0.10, 0.12, 0.11 }};
- FloatingPointTolerance tolerance = relativeToleranceAsFloatingPoint(0.1, 0.001);
- initSystem(4, c_threeBondsMasses,
- c_threeBondsConstraints, c_threeBondsConstraintsR0, c_threeBondsConstraintsFirstType,
- epbcNONE, c_infinitesimalBox,
- real(0.0), real(0.001), c_threeBondsCoordinates);
- applyConstraintsShake(0.0001, true);
- checkConstrained(tolerance);
-}
+ std::vector<RVec> v = {{ 1.0, 0.0, 0.0 },
+ { 0.0, 1.0, 0.0 },
+ { 0.0, 0.0, 1.0 }};
-/*
- * Three consecutive constraints test for LINCS.
- */
-TEST_F(ConstraintsTest, LincsThreeConsequativeConstraints)
-{
+ tensor virialScaledRef = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
- FloatingPointTolerance tolerance = relativeToleranceAsFloatingPoint(0.1, 0.001);
- initSystem(4, c_threeBondsMasses,
- c_threeBondsConstraints, c_threeBondsConstraintsR0, c_threeBondsConstraintsFirstType,
- epbcNONE, c_infinitesimalBox,
- real(0.0), real(0.001), c_threeBondsCoordinates);
- applyConstraintsLincs(4, 4, real(30.0));
- checkConstrained(tolerance);
-}
+ real shakeTolerance = 0.0001;
+ gmx_bool shakeUseSOR = false;
+ int lincsNIter = 1;
+ int lincsNProjOrder = 4;
+ real lincsWarnAngle = 30.0;
-/*
- *
- * Tests of both final interatomic distances and final coordinates.
- * Peformed on the peptide Cys-Val-Trp, with constraints along HBonds,
- * AllBonds (for both SHAKE and LINCS), HAngles and AllAlngles (for
- * SHAKE only).
- */
+ std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>
+ (title, nAtom, masses,
+ constraints, constraintsR0,
+ false, virialScaledRef,
+ false, 0,
+ real(0.0), real(0.001),
+ x, xPrime, v,
+ shakeTolerance, shakeUseSOR,
+ lincsNIter, lincsNProjOrder, lincsWarnAngle);
-/*
- * All bonds that involve hydrogen atom (HBonds) are constrained by SHAKE. Both final length
- * of each constraint and final coordinates for each atom are tested.
- */
-TEST_F(ConstraintsTest, ShakeCysValTrpHBonds)
-{
+ std::string pbcName;
+ std::string algorithmName;
+ std::tie(pbcName, algorithmName) = GetParam();
+ t_pbc pbc = pbcs_.at(pbcName);
+
+ // Apply constraints
+ algorithms_.at(algorithmName)(testData.get(), pbc);
- FloatingPointTolerance tolerance = relativeToleranceAsFloatingPoint(0.1, 0.0001);
- FloatingPointTolerance coordinatesTolerance = absoluteTolerance(0.0001);
+ checkConstrainsLength(absoluteTolerance(0.0002), *testData, pbc);
+ checkConstrainsDirection(*testData, pbc);
+ checkCOMCoordinates(absoluteTolerance(0.0001), *testData);
+ checkCOMVelocity(absoluteTolerance(0.0001), *testData);
- initSystem(54, c_cvwMasses,
- c_cvwHBondsConstraints, c_cvwHBondsConstraintsR0, c_cvwHBondsConstraintsFirstType,
- epbcXYZ, c_cvwBox,
- real(0.0), real(0.001), c_cvwInitialCoordinates);
- applyConstraintsShake(0.0001, true);
- checkConstrained(tolerance);
- checkFinalCoordinates(c_cvwFinalCoordinatesHBonds, coordinatesTolerance);
}
-/*
- * All bonds that involve hydrogen atom (HBonds) are constrained by LINCS. Both final length
- * of each constraint and final coordinates for each atom are tested.
- */
-TEST_F(ConstraintsTest, LincsCysValTrpHBonds)
-{
+TEST_P(ConstraintsTest, ThreeConstraintsWithCentralAtom){
- FloatingPointTolerance tolerance = relativeToleranceAsFloatingPoint(0.1, 0.00001);
- FloatingPointTolerance coordinatesTolerance = absoluteTolerance(0.0001);
+ std::string title = "three atoms, connected to the central atom (e.g. CH3)";
+ int nAtom = 4;
+ std::vector<real> masses = {12.0, 1.0, 1.0, 1.0 };
+ std::vector<int> constraints = {0, 0, 1, 0, 0, 2, 0, 0, 3};
+ std::vector<real> constraintsR0 = {0.1};
- initSystem(54, c_cvwMasses,
- c_cvwHBondsConstraints, c_cvwHBondsConstraintsR0, c_cvwHBondsConstraintsFirstType,
- epbcXYZ, c_cvwBox,
- real(0.0), real(0.001), c_cvwInitialCoordinates);
- applyConstraintsLincs(1, 4, real(30.0));
- checkConstrained(tolerance);
- checkFinalCoordinates(c_cvwFinalCoordinatesHBonds, coordinatesTolerance);
-}
-/*
- * All bonds are constrained by SHAKE. Both final length of each constraint and final coordinates
- * for each atom are tested.
- */
-TEST_F(ConstraintsTest, ShakeCysValTrpAllBonds)
-{
+ std::vector<RVec> x = {{ 0.00, 0.00, 0.00 },
+ { 0.10, 0.00, 0.00 },
+ { 0.00, -0.10, 0.00 },
+ { 0.00, 0.00, 0.10 }};
+
+ std::vector<RVec> xPrime = {{ 0.004, 0.009, -0.010 },
+ { 0.110, -0.006, 0.003 },
+ {-0.007, -0.102, -0.007 },
+ {-0.005, 0.011, 0.102 }};
+
+ std::vector<RVec> v = {{ 1.0, 0.0, 0.0 },
+ { 1.0, 0.0, 0.0 },
+ { 1.0, 0.0, 0.0 },
+ { 1.0, 0.0, 0.0 }};
+
+ tensor virialScaledRef = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
+
+ real shakeTolerance = 0.0001;
+ gmx_bool shakeUseSOR = false;
- FloatingPointTolerance tolerance = relativeToleranceAsFloatingPoint(0.1, 0.0001);
- FloatingPointTolerance coordinatesTolerance = absoluteTolerance(0.001);
+ int lincsNIter = 1;
+ int lincsNProjOrder = 4;
+ real lincsWarnAngle = 30.0;
+
+ std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>
+ (title, nAtom, masses,
+ constraints, constraintsR0,
+ false, virialScaledRef,
+ false, 0,
+ real(0.0), real(0.001),
+ x, xPrime, v,
+ shakeTolerance, shakeUseSOR,
+ lincsNIter, lincsNProjOrder, lincsWarnAngle);
+
+ std::string pbcName;
+ std::string algorithmName;
+ std::tie(pbcName, algorithmName) = GetParam();
+ t_pbc pbc = pbcs_.at(pbcName);
+
+ // Apply constraints
+ algorithms_.at(algorithmName)(testData.get(), pbc);
+
+ checkConstrainsLength(absoluteTolerance(0.0002), *testData, pbc);
+ checkConstrainsDirection(*testData, pbc);
+ checkCOMCoordinates(absoluteTolerance(0.0001), *testData);
+ checkCOMVelocity(absoluteTolerance(0.0001), *testData);
- initSystem(54, c_cvwMasses,
- c_cvwAllBondsConstraints, c_cvwAllBondsConstraintsR0, c_cvwAllBondsConstraintsFirstType,
- epbcXYZ, c_cvwBox,
- real(0.0), real(0.001), c_cvwInitialCoordinates);
- applyConstraintsShake(0.0001, true);
- checkConstrained(tolerance);
- checkFinalCoordinates(c_cvwFinalCoordinatesAllBonds, coordinatesTolerance);
}
-/*
- * All bonds are constrained by LINCS. Both final length of each constraint and final coordinates
- * for each atom are tested.
- */
-TEST_F(ConstraintsTest, LincsCysValTrpAllBonds)
-{
+TEST_P(ConstraintsTest, FourSequentialConstraints){
- FloatingPointTolerance tolerance = relativeToleranceAsFloatingPoint(0.1, 0.0001);
- FloatingPointTolerance coordinatesTolerance = absoluteTolerance(0.001);
+ std::string title = "four atoms, connected longitudinally";
+ int nAtom = 4;
+ std::vector<real> masses = {0.5, 1.0/3.0, 0.25, 1.0};
+ std::vector<int> constraints = {0, 0, 1, 1, 1, 2, 2, 2, 3};
+ std::vector<real> constraintsR0 = {2.0, 1.0, 1.0};
+
+
+ std::vector<RVec> x = {{ 2.50, -3.10, 15.70 },
+ { 0.51, -3.02, 15.55 },
+ { -0.50, -3.00, 15.20 },
+ { -1.51, -2.95, 15.05 }};
+
+ std::vector<RVec> xPrime = {{ 2.50, -3.10, 15.70 },
+ { 0.51, -3.02, 15.55 },
+ { -0.50, -3.00, 15.20 },
+ { -1.51, -2.95, 15.05 }};
+
+ std::vector<RVec> v = {{ 0.0, 0.0, 2.0 },
+ { 0.0, 0.0, 3.0 },
+ { 0.0, 0.0, -4.0 },
+ { 0.0, 0.0, -1.0 }};
+
+ tensor virialScaledRef = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
+
+ real shakeTolerance = 0.0001;
+ gmx_bool shakeUseSOR = false;
+
+ int lincsNIter = 4;
+ int lincsNProjOrder = 8;
+ real lincsWarnAngle = 30.0;
+
+ std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>
+ (title, nAtom, masses,
+ constraints, constraintsR0,
+ false, virialScaledRef,
+ false, 0,
+ real(0.0), real(0.001),
+ x, xPrime, v,
+ shakeTolerance, shakeUseSOR,
+ lincsNIter, lincsNProjOrder, lincsWarnAngle);
+
+ std::string pbcName;
+ std::string algorithmName;
+ std::tie(pbcName, algorithmName) = GetParam();
+ t_pbc pbc = pbcs_.at(pbcName);
+
+ // Apply constraints
+ algorithms_.at(algorithmName)(testData.get(), pbc);
+
+ checkConstrainsLength(absoluteTolerance(0.0002), *testData, pbc);
+ checkConstrainsDirection(*testData, pbc);
+ checkCOMCoordinates(absoluteTolerance(0.0001), *testData);
+ checkCOMVelocity(absoluteTolerance(0.0001), *testData);
- initSystem(54, c_cvwMasses,
- c_cvwAllBondsConstraints, c_cvwAllBondsConstraintsR0, c_cvwAllBondsConstraintsFirstType,
- epbcXYZ, c_cvwBox,
- real(0.0), real(0.001), c_cvwInitialCoordinates);
- applyConstraintsLincs(4, 4, real(30.0));
- checkConstrained(tolerance);
- checkFinalCoordinates(c_cvwFinalCoordinatesAllBonds, coordinatesTolerance);
}
-/*
- * H angles and all bonds are constrained by SHAKE. Only final lengths of constraints are tested.
- */
-TEST_F(ConstraintsTest, ShakeCysValTrpHAngles)
-{
+TEST_P(ConstraintsTest, TriangleOfConstraints){
+
+ std::string title = "basic triangle (tree atoms, connected to each other)";
+ int nAtom = 3;
+ std::vector<real> masses = {1.0, 1.0, 1.0};
+ std::vector<int> constraints = {0, 0, 1, 2, 0, 2, 1, 1, 2};
+ std::vector<real> constraintsR0 = {0.1, 0.1, 0.1};
- FloatingPointTolerance tolerance = relativeToleranceAsFloatingPoint(0.1, 0.0001);
+ real oneTenthOverSqrtTwo = 0.1_real / std::sqrt(2.0_real);
+
+ std::vector<RVec> x = {{ oneTenthOverSqrtTwo, 0.0, 0.0 },
+ { 0.0, oneTenthOverSqrtTwo, 0.0 },
+ { 0.0, 0.0, oneTenthOverSqrtTwo }};
+
+ std::vector<RVec> xPrime = {{ 0.09, -0.02, 0.01 },
+ { -0.02, 0.10, -0.02 },
+ { 0.03, -0.01, 0.07 }};
+
+ std::vector<RVec> v = {{ 1.0, 1.0, 1.0 },
+ { -2.0, -2.0, -2.0 },
+ { 1.0, 1.0, 1.0 }};
+
+ tensor virialScaledRef = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
+
+ real shakeTolerance = 0.0001;
+ gmx_bool shakeUseSOR = false;
+
+ int lincsNIter = 1;
+ int lincsNProjOrder = 4;
+ real lincsWarnAngle = 30.0;
+
+ std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>
+ (title, nAtom, masses,
+ constraints, constraintsR0,
+ false, virialScaledRef,
+ false, 0,
+ real(0.0), real(0.001),
+ x, xPrime, v,
+ shakeTolerance, shakeUseSOR,
+ lincsNIter, lincsNProjOrder, lincsWarnAngle);
+
+ std::string pbcName;
+ std::string algorithmName;
+ std::tie(pbcName, algorithmName) = GetParam();
+ t_pbc pbc = pbcs_.at(pbcName);
+
+ // Apply constraints
+ algorithms_.at(algorithmName)(testData.get(), pbc);
+
+ checkConstrainsLength(absoluteTolerance(0.0002), *testData, pbc);
+ checkConstrainsDirection(*testData, pbc);
+ checkCOMCoordinates(absoluteTolerance(0.0001), *testData);
+ checkCOMVelocity(absoluteTolerance(0.0001), *testData);
- initSystem(54, c_cvwMasses,
- c_cvwHAnglesConstraints, c_cvwHAnglesConstraintsR0, c_cvwHAnglesConstraintsFirstType,
- epbcXYZ, c_cvwBox,
- real(0.0), real(0.001), c_cvwInitialCoordinates);
- applyConstraintsShake(0.0001, true);
- checkConstrained(tolerance);
}
-/*
- * All angles and all bonds are constrained by SHAKE. Only final lengths of constraints are tested.
- */
-/*TEST_F(ConstraintsTest, ShakeCysValTrpAllAngles)
- {
-
- FloatingPointTolerance tolerance = relativeToleranceAsFloatingPoint(0.1, 0.001);
-
- initSystem(54, c_cvwMasses,
- c_cvwAllAnglesConstraints, c_cvwAllAnglesConstraintsR0, c_cvwAllAnglesConstraintsFirstType,
- epbcXYZ, c_cvwBox,
- real(0.0), real(0.001), c_cvwInitialCoordinates);
- applyConstraintsShake(0.000001, true);
- checkConstrained(tolerance);
- }
- */
+
+INSTANTIATE_TEST_CASE_P(WithParameters, ConstraintsTest,
+ ::testing::Combine(::testing::Values("PBCNone", "PBCXYZ"),
+ ::testing::Values("SHAKE", "LINCS")));
+
} // namespace test
} // namespace gmx
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 2018, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-/*! \internal \file
- * \brief Data used by LINCS tests.
- *
- * \author Artem Zhmurov <zhmurov@gmail.com>
- * \ingroup module_mdlib
- */
-
-#include "gromacs/math/vectypes.h"
-
-namespace gmx
-{
-
-namespace test
-{
-
-//! A box with zero sides to check if PBC are actually disabled.
-static const matrix c_infinitesimalBox = { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} };
-
-/*
- * System of two connected atoms.
- */
-//! System of two connected atoms: one constraint (type, i, j).
-static const std::vector<int> c_oneBondConstraints( {0, 0, 1} );
-//! System of two connected atoms: equilibrium distance.
-static const std::vector<real> c_oneBondConstraintsR0( {0.1} );
-//! System of two connected atoms: starting type id for constraints.
-static const int c_oneBondConstraintsFirstType = 0;
-//! System of two connected atoms: coordinates.
-static const std::vector<RVec> c_oneBondCoordinates(
- {{
- { 0.10, -0.13, 0.05 },
- { -0.13, 0.05, -0.10 }
- }} );
-//! System of two connected atoms: masses.
-static const std::vector<real> c_oneBondMasses ( {1.0, 12.0} );
-
-/*
- * Two disjoint bonds (taken from the SHAKE test).
- */
-//! Two disjoint bonds: two constraint (type1, i1, j1, type2, i2, j2).
-static const std::vector<int> c_twoDJBondsConstraints( {0, 0, 1, 1, 2, 3} );
-//! Two disjoint bonds: equilibrium distances.
-static const std::vector<real> c_twoDJBondsConstraintsR0( {2.0, 1.0} );
-//! Two disjoint bonds: starting type id for constraints.
-static const int c_twoDJBondsConstraintsFirstType = 0;
-//! Two disjoint bonds: coordinates.
-static const std::vector<RVec> c_twoDJBondsCoordinates(
- {{
- { 2.50, -3.10, 15.70 },
- { 0.51, -3.02, 15.55 },
- { -0.50, -3.00, 15.20 },
- { -1.51, -2.95, 15.05 }
- }} );
-//! Two disjoint bonds: masses.
-static const std::vector<real> c_twoDJBondsMasses ( {0.5, 1.0/3.0, 0.25, 1.0} );
-
-
-/*
- * Three atoms, connected longitudaly.
- */
-//! Three atoms, connected longitudaly: two constraints (type1, i1, j1, type2, i2, j2).
-static const std::vector<int> c_twoBondsConstraints( {0, 0, 1, 1, 1, 2} );
-//! Three atoms, connected longitudaly: two distances.
-static const std::vector<real> c_twoBondsConstraintsR0( {0.1, 0.2} );
-//! Three atoms, connected longitudaly: starting type id for constraints.
-static const int c_twoBondsConstraintsFirstType = 0;
-//! Three atoms, connected longitudaly: coordinates.
-static const std::vector<RVec> c_twoBondsCoordinates(
- {{
- { 0.08, 0.05, 0.01 },
- { -0.02, 0.03, -0.02 },
- { 0.14, 0.02, 0.01 }
- }} );
-//! Three atoms, connected longitudaly: masses.
-static const std::vector<real> c_twoBondsMasses ( {1.0, 12.0, 16.0 } );
-
-/*
- * Four atoms, connected longitudaly (taken from SHAKE test).
- */
-//! Four atoms, connected longitudaly: two constraints (type1, i1, j1, type2, i2, j2).
-static const std::vector<int> c_threeBondsConstraints( {0, 0, 1, 1, 1, 2, 2, 2, 3} );
-//! Four atoms, connected longitudaly: two distances.
-static const std::vector<real> c_threeBondsConstraintsR0( {2.0, 1.0, 1.0} );
-//! Four atoms, connected longitudaly: starting type id for constraints.
-static const int c_threeBondsConstraintsFirstType = 0;
-//! Four atoms, connected longitudaly: coordinates.
-static const std::vector<RVec> c_threeBondsCoordinates(
- {{
- { 2.50, -3.10, 15.70 },
- { 0.51, -3.02, 15.55 },
- { -0.50, -3.00, 15.20 },
- { -1.51, -2.95, 15.05 }
- }} );
-//! Four atoms, connected longitudaly: masses.
-static const std::vector<real> c_threeBondsMasses ( {0.5, 1.0/3.0, 0.25, 1.0} );
-
-/*
- * Basic triangle (tree atoms, connected with each other).
- */
-//! Basic triangle: three constraints (type1, i1, j1, type2, i2, j2, type3, i3, j3).
-static const std::vector<int> c_triangleConstraints( {0, 0, 1, 2, 0, 2, 1, 1, 2} );
-//! Basic triangle: euilibrium distances.
-static const std::vector<real> c_triangleConstraintsR0( {0.1, 0.1, 0.1} );
-//! Basic triangle: starting type id for constraints.
-static const int c_triangleConstraintsFirstType = 0;
-//! Basic triangle: coordinates.
-static const std::vector<RVec> c_triangleCoordinates(
- {{
- { 0.09, -0.02, 0.01 },
- { -0.02, 0.10, -0.02 },
- { 0.03, -0.01, 0.07 }
- }} );
-//! Basic triangle: masses.
-static const std::vector<real> c_triangleMasses ( {1.0, 1.0, 1.0} );
-
-
-
-/*
- * Real-life system: Cys-Val-Trp peptide.
- */
-//! CVW peptide: periodic box.
-static const matrix c_cvwBox = {{real(2.570950), 0, 0}, {0, real(2.570950), 0}, {0, 0, real(2.570950)}};
-
-
-/*
- * Constraints only on covalent bonds with hydrogens.
- */
-//! CVW peptide: constraints on bonds with hydrogens (type1, i1, j1, type2, i2, j2,...).
-static const std::vector<int> c_cvwHBondsConstraints(
- {465, 0, 1, 465, 0, 2, 465, 0, 3, 466, 4, 5, 467, 6, 7, 467, 6, 8, 468, 9, 10, 469, 13, 14,
- 466, 15, 16, 467, 17, 18, 467, 19, 20, 467, 19, 21, 467, 19, 22, 467, 23, 24, 467, 23, 25, 467, 23, 26,
- 469, 29, 30, 466, 31, 32, 467, 33, 34, 467, 33, 35, 466, 37, 38, 470, 39, 40, 466, 43, 44, 466, 45, 46,
- 466, 47, 48, 466, 49, 50} );
-//! CVW peptide: equilibrium distances (one for each type of HBond).
-static const std::vector<real> c_cvwHBondsConstraintsR0(
- {0.104000, 0.108000, 0.111100, 0.132500, 0.099700, 0.097600} );
-//! CVW peptide: first type id used for constraints when HBonds are constrained.
-static const int c_cvwHBondsConstraintsFirstType = 465;
-
-
-/*
- * Constraints on all covalent bonds.
- */
-//! CVW peptide: constraints on all bonds (type1, i1, j1, type2, i2, j2,...).
-static const std::vector<int> c_cvwAllBondsConstraints(
- {447, 0, 1, 447, 0, 2, 447, 0, 3, 448, 0, 4, 449, 4, 5, 450, 4, 6, 451, 4, 11, 452, 6, 7,
- 452, 6, 8, 453, 6, 9, 454, 9, 10, 455, 11, 12, 456, 11, 13, 457, 13, 14, 458, 13, 15, 449, 15, 16,
- 459, 15, 17, 451, 15, 27, 452, 17, 18, 450, 17, 19, 450, 17, 23, 452, 19, 20, 452, 19, 21, 452, 19, 22,
- 452, 23, 24, 452, 23, 25, 452, 23, 26, 455, 27, 28, 456, 27, 29, 457, 29, 30, 458, 29, 31, 449, 31, 32,
- 450, 31, 33, 460, 31, 51, 452, 33, 34, 452, 33, 35, 461, 33, 36, 462, 36, 37, 463, 36, 42, 449, 37, 38,
- 464, 37, 39, 465, 39, 40, 466, 39, 41, 467, 41, 42, 468, 41, 47, 468, 42, 43, 449, 43, 44, 466, 43, 45,
- 449, 45, 46, 466, 45, 49, 449, 47, 48, 466, 47, 49, 449, 49, 50, 469, 51, 52, 469, 51, 53} );
-//! CVW peptide: equilibrium distances (one for each type of bond).
-static const std::vector<real> c_cvwAllBondsConstraintsR0(
- {0.104000, 0.148000, 0.108000, 0.153800, 0.149000, 0.111100, 0.181800, 0.132500, 0.123000, 0.134500,
- 0.099700, 0.143000, 0.150000, 0.152200, 0.151000, 0.136500, 0.144000, 0.137000, 0.097600, 0.137500,
- 0.140000, 0.136800, 0.126000} );
-//! CVW peptide: first type id used for constraints when all bonds are constrained.
-static const int c_cvwAllBondsConstraintsFirstType = 447;
-
-
-/*
- * Constraints on all covalent bonds and all angles with hydrogens.
- */
-//! CVW peptide: constraints on angles with hydrogens and all bonds (type1, i1, j1, type2, i2, j2,...).
-static const std::vector<int> c_cvwHAnglesConstraints(
- {444, 1, 2, 444, 1, 3, 444, 2, 3, 445, 7, 8, 446, 20, 21, 446, 20, 22, 446, 21, 22, 446, 24, 25,
- 446, 24, 26, 446, 25, 26, 445, 34, 35, 447, 0, 1, 447, 0, 2, 447, 0, 3, 448, 0, 4, 449, 4, 5,
- 450, 4, 6, 451, 4, 11, 452, 6, 7, 452, 6, 8, 453, 6, 9, 454, 9, 10, 455, 11, 12, 456, 11, 13,
- 457, 13, 14, 458, 13, 15, 449, 15, 16, 459, 15, 17, 451, 15, 27, 452, 17, 18, 450, 17, 19, 450, 17, 23,
- 452, 19, 20, 452, 19, 21, 452, 19, 22, 452, 23, 24, 452, 23, 25, 452, 23, 26, 455, 27, 28, 456, 27, 29,
- 457, 29, 30, 458, 29, 31, 449, 31, 32, 450, 31, 33, 460, 31, 51, 452, 33, 34, 452, 33, 35, 461, 33, 36,
- 462, 36, 37, 463, 36, 42, 449, 37, 38, 464, 37, 39, 465, 39, 40, 466, 39, 41, 467, 41, 42, 468, 41, 47,
- 468, 42, 43, 449, 43, 44, 466, 43, 45, 449, 45, 46, 466, 45, 49, 449, 47, 48, 466, 47, 49, 449, 49, 50,
- 469, 51, 52, 469, 51, 53} );
-//! CVW peptide: equilibrium distances (one for each type of constraint).
-static const std::vector<real> c_cvwHAnglesConstraintsR0(
- {0.169861, 0.180896, 0.180218, 0.104000, 0.148000, 0.108000, 0.153800, 0.149000, 0.111100, 0.181800,
- 0.132500, 0.123000, 0.134500, 0.099700, 0.143000, 0.150000, 0.152200, 0.151000, 0.136500, 0.144000,
- 0.137000, 0.097600, 0.137500, 0.140000, 0.136800, 0.126000} );
-//! CVW peptide: first type id used.
-static const int c_cvwHAnglesConstraintsFirstType = 444;
-
-// Constrain all bonds and angles.
-//! CVW peptide: constraints on all angles and all bonds (type1, i1, j1, type2, i2, j2,...).
-static const std::vector<int> c_cvwAllAnglesConstraints(
- {402, 1, 2, 403, 52, 53, 404, 31, 53, 404, 31, 52, 405, 47, 50, 405, 45, 50, 406, 45, 47, 405, 48, 49,
- 407, 41, 49, 408, 41, 48, 405, 46, 49, 406, 43, 49, 405, 43, 46, 405, 44, 45, 407, 42, 45, 408, 42, 44,
- 409, 41, 43, 410, 36, 43, 411, 36, 41, 409, 42, 47, 412, 39, 47, 413, 39, 42, 414, 40, 41, 415, 37, 41,
- 416, 37, 40, 417, 38, 39, 418, 36, 39, 419, 36, 38, 420, 37, 42, 421, 33, 42, 422, 33, 37, 423, 35, 36,
- 423, 34, 36, 424, 34, 35, 425, 31, 36, 426, 31, 35, 426, 31, 34, 427, 33, 51, 428, 32, 51, 429, 32, 33,
- 430, 29, 51, 431, 29, 33, 432, 29, 32, 433, 30, 31, 434, 27, 31, 435, 27, 30, 436, 28, 29, 437, 15, 29,
- 438, 15, 28, 439, 25, 26, 439, 24, 26, 439, 24, 25, 426, 17, 26, 426, 17, 25, 426, 17, 24, 439, 21, 22,
- 439, 20, 22, 439, 20, 21, 426, 17, 22, 426, 17, 21, 426, 17, 20, 440, 19, 23, 441, 18, 23, 441, 18, 19,
- 442, 15, 23, 442, 15, 19, 443, 15, 18, 444, 17, 27, 445, 16, 27, 446, 16, 17, 447, 13, 27, 448, 13, 17,
- 432, 13, 16, 433, 14, 15, 434, 11, 15, 435, 11, 14, 436, 12, 13, 437, 4, 13, 438, 4, 12, 449, 6, 10,
- 450, 8, 9, 450, 7, 9, 424, 7, 8, 451, 4, 9, 426, 4, 8, 426, 4, 7, 452, 6, 11, 445, 5, 11,
- 429, 5, 6, 453, 0, 11, 454, 0, 6, 455, 0, 5, 456, 3, 4, 456, 2, 4, 402, 2, 3, 456, 1, 4,
- 402, 1, 3, 457, 0, 1, 457, 0, 2, 457, 0, 3, 458, 0, 4, 459, 4, 5, 460, 4, 6, 461, 4, 11,
- 462, 6, 7, 462, 6, 8, 463, 6, 9, 464, 9, 10, 465, 11, 12, 466, 11, 13, 467, 13, 14, 468, 13, 15,
- 459, 15, 16, 469, 15, 17, 461, 15, 27, 462, 17, 18, 460, 17, 19, 460, 17, 23, 462, 19, 20, 462, 19, 21,
- 462, 19, 22, 462, 23, 24, 462, 23, 25, 462, 23, 26, 465, 27, 28, 466, 27, 29, 467, 29, 30, 468, 29, 31,
- 459, 31, 32, 460, 31, 33, 470, 31, 51, 462, 33, 34, 462, 33, 35, 471, 33, 36, 472, 36, 37, 473, 36, 42,
- 459, 37, 38, 474, 37, 39, 475, 39, 40, 476, 39, 41, 477, 41, 42, 478, 41, 47, 478, 42, 43, 459, 43, 44,
- 476, 43, 45, 459, 45, 46, 476, 45, 49, 459, 47, 48, 476, 47, 49, 459, 49, 50, 479, 51, 52, 479, 51, 53} );
-//! CVW peptide: equilibrium distances (one for each type of constraint).
-static const std::vector<real> c_cvwAllAnglesConstraintsR0(
- {0.169861, 0.222503, 0.238845, 0.213120, 0.238157, 0.235121, 0.214562, 0.242100, 0.255127, 0.228896,
- 0.249204, 0.223650, 0.210257, 0.222075, 0.209794, 0.217730, 0.224038, 0.217273, 0.226106, 0.260490,
- 0.259998, 0.215277, 0.180896, 0.255631, 0.218499, 0.247561, 0.214016, 0.217310, 0.237362, 0.248280,
- 0.204103, 0.208169, 0.240360, 0.206488, 0.225825, 0.241196, 0.237083, 0.180218, 0.257975, 0.218499,
- 0.246566, 0.215168, 0.241897, 0.211207, 0.213951, 0.234753, 0.245062, 0.234108, 0.245088, 0.279474,
- 0.244987, 0.243289, 0.247242, 0.207800, 0.207355, 0.104000, 0.148000, 0.108000, 0.153800, 0.149000,
- 0.111100, 0.181800, 0.132500, 0.123000, 0.134500, 0.099700, 0.143000, 0.150000, 0.152200, 0.151000,
- 0.136500, 0.144000, 0.137000, 0.097600, 0.137500, 0.140000, 0.136800, 0.126000} );
-//! CVW peptide: first type id used.
-static const int c_cvwAllAnglesConstraintsFirstType = 402;
-
-
-//! Coordinates of all 54 atoms of the Cys-Val-Trp system before constraints were applied
-static const std::vector<RVec> c_cvwInitialCoordinates(
- {{
- { 1.746000, 1.404000, 1.052000 },
- { 1.806000, 1.450000, 0.979000 },
- { 1.672000, 1.473000, 1.086000 },
- { 1.814000, 1.368000, 1.123000 },
- { 1.681000, 1.283000, 0.992000 },
- { 1.612000, 1.309000, 0.913000 },
- { 1.801000, 1.204000, 0.919000 },
- { 1.772000, 1.100000, 0.896000 },
- { 1.814000, 1.266000, 0.827000 },
- { 1.974000, 1.193000, 0.998000 },
- { 2.058000, 1.204000, 0.896000 },
- { 1.594000, 1.204000, 1.095000 },
- { 1.568000, 1.083000, 1.091000 },
- { 1.533000, 1.270000, 1.201000 },
- { 1.540000, 1.372000, 1.213000 },
- { 1.458000, 1.206000, 1.310000 },
- { 1.457000, 1.094000, 1.307000 },
- { 1.532000, 1.230000, 1.433000 },
- { 1.527000, 1.340000, 1.449000 },
- { 1.478000, 1.148000, 1.544000 },
- { 1.544000, 1.157000, 1.633000 },
- { 1.378000, 1.176000, 1.581000 },
- { 1.461000, 1.043000, 1.518000 },
- { 1.685000, 1.181000, 1.411000 },
- { 1.747000, 1.243000, 1.345000 },
- { 1.730000, 1.190000, 1.515000 },
- { 1.688000, 1.075000, 1.377000 },
- { 1.312000, 1.253000, 1.324000 },
- { 1.268000, 1.367000, 1.318000 },
- { 1.224000, 1.149000, 1.344000 },
- { 1.264000, 1.058000, 1.349000 },
- { 1.095000, 1.170000, 1.405000 },
- { 1.103000, 1.263000, 1.454000 },
- { 0.975000, 1.166000, 1.297000 },
- { 1.025000, 1.113000, 1.214000 },
- { 0.887000, 1.110000, 1.333000 },
- { 0.924000, 1.306000, 1.265000 },
- { 0.952000, 1.386000, 1.165000 },
- { 1.017000, 1.361000, 1.079000 },
- { 0.882000, 1.507000, 1.174000 },
- { 0.911000, 1.589000, 1.129000 },
- { 0.819000, 1.515000, 1.296000 },
- { 0.840000, 1.389000, 1.355000 },
- { 0.777000, 1.364000, 1.480000 },
- { 0.785000, 1.264000, 1.522000 },
- { 0.698000, 1.458000, 1.538000 },
- { 0.640000, 1.433000, 1.626000 },
- { 0.746000, 1.617000, 1.350000 },
- { 0.719000, 1.711000, 1.307000 },
- { 0.682000, 1.582000, 1.472000 },
- { 0.631000, 1.668000, 1.515000 },
- { 1.070000, 1.060000, 1.508000 },
- { 1.150000, 0.960000, 1.515000 },
- { 0.966000, 1.066000, 1.569000 }
- }});
-
-//! Coordinates of 54 atoms of the Cys-Val-Trp system after H-Bonds were constrained
-static const std::vector<RVec> c_cvwFinalCoordinatesHBonds(
- {{
- { 1.745950, 1.404137, 1.052042 },
- { 1.805376, 1.449522, 0.979759 },
- { 1.673804, 1.471318, 1.085171 },
- { 1.813515, 1.368257, 1.122494 },
- { 1.680997, 1.283001, 0.991996 },
- { 1.612038, 1.308986, 0.913044 },
- { 1.801019, 1.204075, 0.918974 },
- { 1.771832, 1.099398, 0.895867 },
- { 1.813938, 1.265703, 0.827441 },
- { 1.974002, 1.193000, 0.997998 },
- { 2.057943, 1.203992, 0.896070 },
- { 1.594000, 1.204000, 1.095000 },
- { 1.568000, 1.083000, 1.091000 },
- { 1.533015, 1.270216, 1.201025 },
- { 1.539794, 1.369004, 1.212648 },
- { 1.457997, 1.205687, 1.309992 },
- { 1.457033, 1.097730, 1.307100 },
- { 1.531999, 1.230013, 1.433002 },
- { 1.527007, 1.339845, 1.448977 },
- { 1.478082, 1.148102, 1.544007 },
- { 1.543998, 1.157000, 1.632997 },
- { 1.377262, 1.176207, 1.581273 },
- { 1.460769, 1.041573, 1.517647 },
- { 1.685022, 1.180937, 1.411234 },
- { 1.747673, 1.243673, 1.344284 },
- { 1.729067, 1.189813, 1.512844 },
- { 1.687993, 1.075258, 1.377083 },
- { 1.312000, 1.253000, 1.324000 },
- { 1.268000, 1.367000, 1.318000 },
- { 1.223995, 1.149011, 1.343999 },
- { 1.264064, 1.057854, 1.349008 },
- { 1.094985, 1.169824, 1.404907 },
- { 1.103180, 1.265097, 1.455105 },
- { 0.975024, 1.166056, 1.297020 },
- { 1.025283, 1.112700, 1.213531 },
- { 0.886430, 1.109637, 1.333233 },
- { 0.924000, 1.306000, 1.265000 },
- { 0.952121, 1.385953, 1.164840 },
- { 1.015558, 1.361555, 1.080908 },
- { 0.882007, 1.507018, 1.173990 },
- { 0.910909, 1.588743, 1.129141 },
- { 0.819000, 1.515000, 1.296000 },
- { 0.840000, 1.389000, 1.355000 },
- { 0.777004, 1.363946, 1.480023 },
- { 0.784949, 1.264642, 1.521730 },
- { 0.697987, 1.457994, 1.538020 },
- { 0.640158, 1.433068, 1.625761 },
- { 0.746023, 1.616921, 1.350036 },
- { 0.718729, 1.711945, 1.306568 },
- { 0.681970, 1.582051, 1.472026 },
- { 0.631363, 1.667388, 1.514694 },
- { 1.070000, 1.060000, 1.508000 },
- { 1.150000, 0.960000, 1.515000 },
- { 0.966000, 1.066000, 1.569000 }
- }});
-
-//! Coordinates of 54 atoms of the Cys-Val-Trp system after all bonds were constrained
-static const std::vector<RVec> c_cvwFinalCoordinatesAllBonds(
- {{
- { 1.744192, 1.400876, 1.050408 },
- { 1.804620, 1.448942, 0.980679 },
- { 1.674799, 1.470390, 1.084714 },
- { 1.812763, 1.368655, 1.121708 },
- { 1.683948, 1.276865, 0.996242 },
- { 1.616170, 1.307429, 0.917774 },
- { 1.803010, 1.207954, 0.927155 },
- { 1.773417, 1.105083, 0.897124 },
- { 1.813404, 1.263159, 0.831216 },
- { 1.970550, 1.193223, 0.996389 },
- { 2.057327, 1.203912, 0.896817 },
- { 1.595071, 1.207624, 1.094104 },
- { 1.568932, 1.087338, 1.091143 },
- { 1.533602, 1.264980, 1.199359 },
- { 1.539437, 1.363790, 1.212034 },
- { 1.455744, 1.209606, 1.306025 },
- { 1.457067, 1.101535, 1.307202 },
- { 1.536585, 1.230178, 1.430802 },
- { 1.527027, 1.339398, 1.448912 },
- { 1.477277, 1.146856, 1.545698 },
- { 1.544449, 1.157061, 1.633605 },
- { 1.376383, 1.176453, 1.581598 },
- { 1.460628, 1.040705, 1.517432 },
- { 1.681448, 1.182094, 1.411715 },
- { 1.746716, 1.242717, 1.345302 },
- { 1.728696, 1.189739, 1.511987 },
- { 1.687955, 1.076581, 1.377507 },
- { 1.313454, 1.250569, 1.323999 },
- { 1.268736, 1.365093, 1.318100 },
- { 1.224082, 1.152044, 1.344591 },
- { 1.262989, 1.060300, 1.348874 },
- { 1.093660, 1.168520, 1.401086 },
- { 1.102950, 1.262423, 1.453696 },
- { 0.977797, 1.168815, 1.299821 },
- { 1.024283, 1.113760, 1.215191 },
- { 0.888605, 1.111021, 1.332343 },
- { 0.921772, 1.305848, 1.269492 },
- { 0.951364, 1.387941, 1.164422 },
- { 1.015236, 1.361678, 1.081333 },
- { 0.882251, 1.505804, 1.174960 },
- { 0.910507, 1.587606, 1.129765 },
- { 0.819059, 1.513487, 1.296914 },
- { 0.840064, 1.387880, 1.355340 },
- { 0.778667, 1.365803, 1.475722 },
- { 0.784710, 1.267623, 1.520478 },
- { 0.698397, 1.459946, 1.535866 },
- { 0.641444, 1.433622, 1.623810 },
- { 0.745363, 1.614637, 1.352384 },
- { 0.719425, 1.709521, 1.307677 },
- { 0.683959, 1.580100, 1.470510 },
- { 0.632793, 1.664977, 1.513489 },
- { 1.072860, 1.060117, 1.505945 },
- { 1.149574, 0.960533, 1.514963 },
- { 0.964441, 1.066090, 1.569914 }
- }});
-
-//! Masses of atoms in the Cys-Val-Trp system
-static const std::vector<real> c_cvwMasses ( {
- 14.007000,
- 1.008000,
- 1.008000,
- 1.008000,
- 12.011000,
- 1.008000,
- 12.011000,
- 1.008000,
- 1.008000,
- 32.060001,
- 1.008000,
- 12.011000,
- 15.999000,
- 14.007000,
- 1.008000,
- 12.011000,
- 1.008000,
- 12.011000,
- 1.008000,
- 12.011000,
- 1.008000,
- 1.008000,
- 1.008000,
- 12.011000,
- 1.008000,
- 1.008000,
- 1.008000,
- 12.011000,
- 15.999000,
- 14.007000,
- 1.008000,
- 12.011000,
- 1.008000,
- 12.011000,
- 1.008000,
- 1.008000,
- 12.011000,
- 12.011000,
- 1.008000,
- 14.007000,
- 1.008000,
- 12.011000,
- 12.011000,
- 12.011000,
- 1.008000,
- 12.011000,
- 1.008000,
- 12.011000,
- 1.008000,
- 12.011000,
- 1.008000,
- 12.011000,
- 15.999400,
- 15.999400
- } );
-
-
-} // namespace test
-
-} // namespace gmx
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff