/* TODO Specialize this routine into init-time and loop-time versions?
e.g. bReadEkin is only true when restoring from checkpoint */
-void compute_globals(FILE *fplog, gmx_global_stat *gstat, t_commrec *cr, t_inputrec *ir,
+void compute_globals(gmx_global_stat *gstat, t_commrec *cr, const t_inputrec *ir,
t_forcerec *fr, gmx_ekindata_t *ekind,
- rvec *x, rvec *v, matrix box, real vdwLambda, t_mdatoms *mdatoms,
+ const rvec *x, const rvec *v, const matrix box,
+ real vdwLambda, const t_mdatoms *mdatoms,
t_nrnb *nrnb, t_vcm *vcm, gmx_wallcycle_t wcycle,
gmx_enerdata_t *enerd, tensor force_vir, tensor shake_vir, tensor total_vir,
tensor pres, rvec mu_tot, gmx::Constraints *constr,
gmx::SimulationSignaller *signalCoordinator,
- matrix lastbox, int *totalNumberOfBondedInteractions,
- gmx_bool *bSumEkinhOld, int flags)
+ const matrix lastbox, int *totalNumberOfBondedInteractions,
+ gmx_bool *bSumEkinhOld, const int flags)
{
gmx_bool bEner, bPres, bTemp;
gmx_bool bStopCM, bGStat,
/* Calculate center of mass velocity if necessary, also parallellized */
if (bStopCM)
{
- calc_vcm_grp(0, mdatoms->homenr, mdatoms,
- x, v, vcm);
+ calc_vcm_grp(*mdatoms, x, v, vcm);
}
if (bTemp || bStopCM || bPres || bEner || bConstrain || bCheckNumberOfBondedInteractions)
}
}
- /* Do center of mass motion removal */
- if (bStopCM)
- {
- check_cm_grp(fplog, vcm, ir, 1);
- /* At initialization, do not pass x with acceleration-correction mode
- * to avoid (incorrect) correction of the initial coordinates.
- */
- rvec *xPtr = nullptr;
- if (vcm->mode == ecmANGULAR || (vcm->mode == ecmLINEAR_ACCELERATION_CORRECTION && !(flags & CGLO_INITIALIZATION)))
- {
- xPtr = x;
- }
- do_stopcm_grp(*mdatoms,
- xPtr, v, *vcm);
- inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
- }
-
if (bEner)
{
/* Calculate the amplitude of the cosine velocity profile */
* passed to compute_globals in md.c and global_stat.
*/
-/* we are initializing and not yet in the actual MD loop */
-#define CGLO_INITIALIZATION (1u<<1u)
/* we are computing the kinetic energy from average velocities */
#define CGLO_EKINAVEVEL (1u<<2u)
/* we are removing the center of mass momenta */
int multisim_min(const gmx_multisim_t *ms, int nmin, int n);
/* Set an appropriate value for n across the whole multi-simulation */
-void compute_globals(FILE *fplog, gmx_global_stat *gstat, t_commrec *cr, t_inputrec *ir,
+
+/* Compute global variables during integration
+ *
+ * Coordinates x are needed for kinetic energy calculation with cosine accelation
+ * and for COM removal with rotational and acceleration correction modes.
+ * Velocities v are needed for kinetic energy calculation and for COM removal.
+ */
+void compute_globals(gmx_global_stat *gstat, t_commrec *cr, const t_inputrec *ir,
t_forcerec *fr, gmx_ekindata_t *ekind,
- rvec *x, rvec *v, matrix box, real vdwLambda, t_mdatoms *mdatoms,
+ const rvec *x, const rvec *v, const matrix box,
+ real vdwLambda, const t_mdatoms *mdatoms,
t_nrnb *nrnb, t_vcm *vcm, gmx_wallcycle_t wcycle,
gmx_enerdata_t *enerd, tensor force_vir, tensor shake_vir, tensor total_vir,
tensor pres, rvec mu_tot, gmx::Constraints *constr,
gmx::SimulationSignaller *signalCoordinator,
- matrix lastbox, int *totalNumberOfBondedInteractions,
+ const matrix lastbox, int *totalNumberOfBondedInteractions,
gmx_bool *bSumEkinhOld, int flags);
-/* Compute global variables during integration */
#endif
# the research papers on the package. Check out http://www.gromacs.org.
file(GLOB MDLIB_TEST_SOURCES
- constr_impl.cpp
+ constrtestrunners.cpp
+ constrtestdata.cpp
)
if (GMX_USE_CUDA)
file(GLOB MDLIB_TEST_CUDA_SOURCES
- constr_impl.cu
+ constrtestrunners.cu
settle_runners.cu
)
endif()
* \author Artem Zhmurov <zhmurov@gmail.com>
* \ingroup module_mdlib
*/
-
#include "gmxpre.h"
-#include "gromacs/mdlib/constr.h"
-
#include "config.h"
#include <assert.h>
-#include <cmath>
-
-#include <algorithm>
#include <unordered_map>
#include <vector>
#include <gtest/gtest.h>
-#include "gromacs/fileio/gmxfio.h"
-#include "gromacs/gmxlib/nrnb.h"
-#include "gromacs/math/paddedvector.h"
#include "gromacs/math/vec.h"
#include "gromacs/math/vectypes.h"
-#include "gromacs/mdlib/lincs.h"
-#include "gromacs/mdlib/shake.h"
-#include "gromacs/mdrunutility/multisim.h"
-#include "gromacs/mdtypes/commrec.h"
-#include "gromacs/mdtypes/inputrec.h"
-#include "gromacs/mdtypes/mdatom.h"
#include "gromacs/pbcutil/pbc.h"
-#include "gromacs/topology/block.h"
-#include "gromacs/topology/idef.h"
-#include "gromacs/topology/ifunc.h"
-#include "gromacs/topology/topology.h"
-#include "gromacs/utility/smalloc.h"
#include "gromacs/utility/stringutil.h"
-#include "gromacs/utility/unique_cptr.h"
-#include "testutils/refdata.h"
#include "testutils/testasserts.h"
-#include "constr_impl.h"
+#include "constrtestdata.h"
+#include "constrtestrunners.h"
namespace gmx
{
namespace test
{
-
-ConstraintsTestData::ConstraintsTestData(const std::string &title,
- int numAtoms, 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 lincsNumIterations, int lincsExpansionOrder, real lincsWarnAngle)
-{
- // This is to trick Gerrit
- {
- {
- title_ = title; // Human-friendly name of the system
- numAtoms_ = numAtoms; // Number of atoms
-
- // Masses of atoms
- masses_ = masses;
- invmass_.resize(numAtoms); // Vector of inverse masses
-
- for (int i = 0; i < numAtoms; i++)
- {
- invmass_[i] = 1.0/masses.at(i);
- }
-
- // Saving constraints to check if they are satisfied after algorithm was applied
- constraints_ = constraints; // Constraints indices (in type-i-j format)
- constraintsR0_ = constraintsR0; // Equilibrium distances for each type of constraint
-
- 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 = numAtoms;
- md_.homenr = numAtoms;
-
- // Virial evaluation
- computeVirial_ = computeVirial;
- if (computeVirial)
- {
- for (int i = 0; i < DIM; i++)
- {
- for (int j = 0; j < DIM; j++)
- {
- virialScaled_[i][j] = 0;
- virialScaledRef_[i][j] = virialScaledRef[i][j];
- }
- }
- }
-
-
- // 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 (local topology)
- for (int i = 0; i < F_NRE; i++)
- {
- idef_.il[i].nr = 0;
- }
- idef_.il[F_CONSTR].nr = constraints.size();
-
- snew(idef_.il[F_CONSTR].iatoms, constraints.size());
- int maxType = 0;
- for (unsigned i = 0; i < constraints.size(); i++)
- {
- if (i % 3 == 0)
- {
- if (maxType < constraints.at(i))
- {
- maxType = constraints.at(i);
- }
- }
- idef_.il[F_CONSTR].iatoms[i] = constraints.at(i);
- }
- 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));
- idef_.iparams[constraints.at(3*i)].constr.dB = constraintsR0.at(constraints.at(3*i));
- }
-
- // Constraints and their parameters (global topology)
- InteractionList interactionList;
- interactionList.iatoms.resize(constraints.size());
- for (unsigned i = 0; i < constraints.size(); i++)
- {
- interactionList.iatoms.at(i) = constraints.at(i);
- }
- InteractionList interactionListEmpty;
- interactionListEmpty.iatoms.resize(0);
-
- gmx_moltype_t molType;
- molType.atoms.nr = numAtoms;
- 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 = numAtoms;
- mtop_.ffparams.iparams.resize(maxType + 1);
- for (int i = 0; i <= maxType; i++)
- {
- mtop_.ffparams.iparams.at(i) = idef_.iparams[i];
- }
- mtop_.bIntermolecularInteractions = false;
-
- // Coordinates and velocities
- x_.resizeWithPadding(numAtoms);
- xPrime_.resizeWithPadding(numAtoms);
- xPrime0_.resizeWithPadding(numAtoms);
- xPrime2_.resizeWithPadding(numAtoms);
-
- v_.resizeWithPadding(numAtoms);
- v0_.resizeWithPadding(numAtoms);
-
- std::copy(x.begin(), x.end(), x_.begin());
- std::copy(xPrime.begin(), xPrime.end(), xPrime_.begin());
- std::copy(xPrime.begin(), xPrime.end(), xPrime0_.begin());
- std::copy(xPrime.begin(), xPrime.end(), xPrime2_.begin());
-
- std::copy(v.begin(), v.end(), v_.begin());
- std::copy(v.begin(), v.end(), v0_.begin());
-
- // SHAKE-specific parameters
- ir_.shake_tol = shakeTolerance;
- ir_.bShakeSOR = shakeUseSOR;
-
- // LINCS-specific parameters
- ir_.nLincsIter = lincsNumIterations;
- ir_.nProjOrder = lincsExpansionOrder;
- 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 ConstraintsTestData::reset()
-{
- xPrime_ = xPrime0_;
- xPrime2_ = xPrime0_;
- v_ = v0_;
-
- if (computeVirial_)
- {
- for (int i = 0; i < DIM; i++)
- {
- for (int j = 0; j < DIM; j++)
- {
- virialScaled_[i][j] = 0;
- }
- }
- }
- dHdLambda_ = 0;
-}
-
-/*! \brief
- * Cleaning up the memory.
- */
-ConstraintsTestData::~ConstraintsTestData()
-{
- sfree(idef_.il[F_CONSTR].iatoms);
- sfree(idef_.iparams);
-}
-
namespace
{
*/
typedef std::tuple<std::string, std::string> ConstraintsTestParameters;
-//! Names of all availible algorithms
-std::vector<std::string> algorithmsNames;
+//! Names of all availible runners
+std::vector<std::string> runnersNames;
-//! Method that fills and returns algorithmNames to the test macros.
-std::vector<std::string> getAlgorithmsNames()
+//! Method that fills and returns runnersNames to the test macros.
+std::vector<std::string> getRunnersNames()
{
- algorithmsNames.emplace_back("SHAKE");
- algorithmsNames.emplace_back("LINCS");
+ runnersNames.emplace_back("SHAKE");
+ runnersNames.emplace_back("LINCS");
if (GMX_GPU == GMX_GPU_CUDA && canComputeOnGpu())
{
- algorithmsNames.emplace_back("LINCS_CUDA");
+ runnersNames.emplace_back("LINCS_CUDA");
}
- return algorithmsNames;
+ return runnersNames;
}
/*! \brief Test fixture for constraints.
lincsNIter, lincslincsExpansionOrder, lincsWarnAngle);
std::string pbcName;
- std::string algorithmName;
- std::tie(pbcName, algorithmName) = GetParam();
- t_pbc pbc = pbcs_.at(pbcName);
+ std::string runnerName;
+ std::tie(pbcName, runnerName) = GetParam();
+ t_pbc pbc = pbcs_.at(pbcName);
// Apply constraints
- algorithms_.at(algorithmName)(testData.get(), pbc);
+ algorithms_.at(runnerName)(testData.get(), pbc);
checkConstrainsLength(absoluteTolerance(0.0002), *testData, pbc);
checkConstrainsDirection(*testData, pbc);
INSTANTIATE_TEST_CASE_P(WithParameters, ConstraintsTest,
::testing::Combine(::testing::Values("PBCNone", "PBCXYZ"),
- ::testing::ValuesIn(getAlgorithmsNames())));
+ ::testing::ValuesIn(getRunnersNames())));
} // namespace
} // namespace test
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * 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.
+ *
+ * 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 SHAKE and LINCS tests.
+ *
+ * \todo Better tests for virial are needed.
+ * \todo Tests for bigger systems to test threads synchronization,
+ * reduction, etc. on the GPU.
+ * \todo Tests for algorithms for derivatives.
+ * \todo Free-energy perturbation tests
+ *
+ * \author Artem Zhmurov <zhmurov@gmail.com>
+ * \ingroup module_mdlib
+ */
+
+#include "gmxpre.h"
+
+#include "constrtestdata.h"
+
+#include "gromacs/utility/smalloc.h"
+#include "gromacs/utility/stringutil.h"
+
+namespace gmx
+{
+namespace test
+{
+
+ConstraintsTestData::ConstraintsTestData(const std::string &title,
+ int numAtoms, 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 lincsNumIterations, int lincsExpansionOrder, real lincsWarnAngle)
+{
+ title_ = title; // Human-friendly name of the system
+ numAtoms_ = numAtoms; // Number of atoms
+
+ // Masses of atoms
+ masses_ = masses;
+ invmass_.resize(numAtoms); // Vector of inverse masses
+
+ for (int i = 0; i < numAtoms; i++)
+ {
+ invmass_[i] = 1.0/masses.at(i);
+ }
+
+ // Saving constraints to check if they are satisfied after algorithm was applied
+ constraints_ = constraints; // Constraints indices (in type-i-j format)
+ constraintsR0_ = constraintsR0; // Equilibrium distances for each type of constraint
+
+ 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 = numAtoms;
+ md_.homenr = numAtoms;
+
+ // Virial evaluation
+ computeVirial_ = computeVirial;
+ if (computeVirial)
+ {
+ for (int i = 0; i < DIM; i++)
+ {
+ for (int j = 0; j < DIM; j++)
+ {
+ virialScaled_[i][j] = 0;
+ virialScaledRef_[i][j] = virialScaledRef[i][j];
+ }
+ }
+ }
+
+
+ // 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 (local topology)
+ for (int i = 0; i < F_NRE; i++)
+ {
+ idef_.il[i].nr = 0;
+ }
+ idef_.il[F_CONSTR].nr = constraints.size();
+
+ snew(idef_.il[F_CONSTR].iatoms, constraints.size());
+ int maxType = 0;
+ for (unsigned i = 0; i < constraints.size(); i++)
+ {
+ if (i % 3 == 0)
+ {
+ if (maxType < constraints.at(i))
+ {
+ maxType = constraints.at(i);
+ }
+ }
+ idef_.il[F_CONSTR].iatoms[i] = constraints.at(i);
+ }
+ 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));
+ idef_.iparams[constraints.at(3*i)].constr.dB = constraintsR0.at(constraints.at(3*i));
+ }
+
+ // Constraints and their parameters (global topology)
+ InteractionList interactionList;
+ interactionList.iatoms.resize(constraints.size());
+ for (unsigned i = 0; i < constraints.size(); i++)
+ {
+ interactionList.iatoms.at(i) = constraints.at(i);
+ }
+ InteractionList interactionListEmpty;
+ interactionListEmpty.iatoms.resize(0);
+
+ gmx_moltype_t molType;
+ molType.atoms.nr = numAtoms;
+ 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 = numAtoms;
+ mtop_.ffparams.iparams.resize(maxType + 1);
+ for (int i = 0; i <= maxType; i++)
+ {
+ mtop_.ffparams.iparams.at(i) = idef_.iparams[i];
+ }
+ mtop_.bIntermolecularInteractions = false;
+
+ // Coordinates and velocities
+ x_.resizeWithPadding(numAtoms);
+ xPrime_.resizeWithPadding(numAtoms);
+ xPrime0_.resizeWithPadding(numAtoms);
+ xPrime2_.resizeWithPadding(numAtoms);
+
+ v_.resizeWithPadding(numAtoms);
+ v0_.resizeWithPadding(numAtoms);
+
+ std::copy(x.begin(), x.end(), x_.begin());
+ std::copy(xPrime.begin(), xPrime.end(), xPrime_.begin());
+ std::copy(xPrime.begin(), xPrime.end(), xPrime0_.begin());
+ std::copy(xPrime.begin(), xPrime.end(), xPrime2_.begin());
+
+ std::copy(v.begin(), v.end(), v_.begin());
+ std::copy(v.begin(), v.end(), v0_.begin());
+
+ // SHAKE-specific parameters
+ ir_.shake_tol = shakeTolerance;
+ ir_.bShakeSOR = shakeUseSOR;
+
+ // LINCS-specific parameters
+ ir_.nLincsIter = lincsNumIterations;
+ ir_.nProjOrder = lincsExpansionOrder;
+ 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 ConstraintsTestData::reset()
+{
+ xPrime_ = xPrime0_;
+ xPrime2_ = xPrime0_;
+ v_ = v0_;
+
+ if (computeVirial_)
+ {
+ for (int i = 0; i < DIM; i++)
+ {
+ for (int j = 0; j < DIM; j++)
+ {
+ virialScaled_[i][j] = 0;
+ }
+ }
+ }
+ dHdLambda_ = 0;
+}
+
+/*! \brief
+ * Cleaning up the memory.
+ */
+ConstraintsTestData::~ConstraintsTestData()
+{
+ sfree(idef_.il[F_CONSTR].iatoms);
+ sfree(idef_.iparams);
+}
+
+} // namespace test
+} // namespace gmx
* \ingroup module_mdlib
*/
-#ifndef GMX_MDLIB_TESTS_CONSTR_IMPL_H
-#define GMX_MDLIB_TESTS_CONSTR_IMPL_H
+#ifndef GMX_MDLIB_TESTS_CONSTRTESTDATA_H
+#define GMX_MDLIB_TESTS_CONSTRTESTDATA_H
-#include "gmxpre.h"
-
-#include <assert.h>
-
-#include <cmath>
-
-#include <algorithm>
-#include <unordered_map>
#include <vector>
-#include "gromacs/fileio/gmxfio.h"
#include "gromacs/gmxlib/nrnb.h"
-#include "gromacs/gmxlib/nonbonded/nonbonded.h"
#include "gromacs/gpu_utils/gpu_testutils.h"
#include "gromacs/math/paddedvector.h"
#include "gromacs/math/vec.h"
#include "gromacs/math/vectypes.h"
-#include "gromacs/mdlib/constr.h"
#include "gromacs/mdlib/gmx_omp_nthreads.h"
#include "gromacs/mdlib/lincs.h"
#include "gromacs/mdlib/shake.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/mdatom.h"
#include "gromacs/pbcutil/pbc.h"
-#include "gromacs/topology/block.h"
#include "gromacs/topology/idef.h"
#include "gromacs/topology/ifunc.h"
#include "gromacs/topology/topology.h"
-#include "gromacs/utility/smalloc.h"
-#include "gromacs/utility/stringutil.h"
-#include "gromacs/utility/unique_cptr.h"
namespace gmx
{
~ConstraintsTestData();
};
-/*! \brief Apply SHAKE constraints to the test data.
- */
-void applyShake(ConstraintsTestData *testData, t_pbc pbc);
-/*! \brief Apply LINCS constraints to the test data.
- */
-void applyLincs(ConstraintsTestData *testData, t_pbc pbc);
-/*! \brief Apply CUDA version of LINCS constraints to the test data.
- *
- * All the data is copied to the GPU device, then LINCS is applied and
- * the resulting coordinates are copied back.
- */
-void applyLincsCuda(ConstraintsTestData *testData, t_pbc pbc);
-
} // namespace test
} // namespace gmx
-#endif // GMX_MDLIB_TESTS_CONSTR_IMPL_H
+#endif // GMX_MDLIB_TESTS_CONSTRTESTDATA_H
#include "gmxpre.h"
-#include "constr_impl.h"
+#include "constrtestrunners.h"
#include "config.h"
#include <cmath>
#include <algorithm>
-#include <unordered_map>
#include <vector>
#include <gtest/gtest.h>
*/
#include "gmxpre.h"
-#include "constr_impl.h"
+#include "constrtestrunners.h"
#include <assert.h>
#include <cmath>
#include <algorithm>
-#include <unordered_map>
#include <vector>
#include "gromacs/gpu_utils/devicebuffer.cuh"
#include "gromacs/gpu_utils/gpu_utils.h"
-#include "gromacs/math/vec.h"
-#include "gromacs/math/vectypes.h"
-#include "gromacs/mdlib/constr.h"
#include "gromacs/mdlib/lincs_cuda.cuh"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/utility/unique_cptr.h"
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * 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.
+ *
+ * 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 SHAKE and LINCS tests header.
+ *
+ * Contains description and constructor for the test data accumulating object,
+ * declares CPU- and GPU-based functions used to apply SHAKE or LINCS on the
+ * test data.
+ *
+ * \author Artem Zhmurov <zhmurov@gmail.com>
+ * \ingroup module_mdlib
+ */
+
+#ifndef GMX_MDLIB_TESTS_CONSTRTESTRUNNERS_H
+#define GMX_MDLIB_TESTS_CONSTRTESTRUNNERS_H
+
+#include "constrtestdata.h"
+
+struct t_pbc;
+
+namespace gmx
+{
+namespace test
+{
+
+/*! \brief Apply SHAKE constraints to the test data.
+ */
+void applyShake(ConstraintsTestData *testData, t_pbc pbc);
+/*! \brief Apply LINCS constraints to the test data.
+ */
+void applyLincs(ConstraintsTestData *testData, t_pbc pbc);
+/*! \brief Apply CUDA version of LINCS constraints to the test data.
+ *
+ * All the data is copied to the GPU device, then LINCS is applied and
+ * the resulting coordinates are copied back.
+ */
+void applyLincsCuda(ConstraintsTestData *testData, t_pbc pbc);
+
+} // namespace test
+} // namespace gmx
+
+#endif // GMX_MDLIB_TESTS_CONSTRTESTRUNNERS_H
}
void calc_ke_part(
- rvec *x, rvec *v, matrix box,
+ const rvec *x, const rvec *v, const matrix box,
const t_grpopts *opts, const t_mdatoms *md,
gmx_ekindata_t *ekind, t_nrnb *nrnb, gmx_bool bEkinAveVel)
{
/* Return TRUE if OK, FALSE in case of Shake Error */
void calc_ke_part(
- rvec *x, rvec *v, matrix box,
+ const rvec *x, const rvec *v, const matrix box,
const t_grpopts *opts, const t_mdatoms *md,
gmx_ekindata_t *ekind, t_nrnb *nrnb, gmx_bool bEkinAveVel);
/*
#include "gromacs/utility/gmxomp.h"
#include "gromacs/utility/smalloc.h"
-t_vcm::t_vcm(const SimulationGroups &groups, const t_inputrec &ir)
+t_vcm::t_vcm(const SimulationGroups &groups, const t_inputrec &ir) :
+ integratorConservesMomentum(!EI_RANDOM(ir.eI))
{
mode = (ir.nstcomm > 0) ? ir.comm_mode : ecmNO;
ndim = ndof_com(&ir);
}
/* Center of mass code for groups */
-void calc_vcm_grp(int start, int homenr, t_mdatoms *md,
- rvec x[], rvec v[], t_vcm *vcm)
+void calc_vcm_grp(const t_mdatoms &md,
+ const rvec x[], const rvec v[], t_vcm *vcm)
{
int nthreads = gmx_omp_nthreads_get(emntDefault);
if (vcm->mode != ecmNO)
}
#pragma omp for schedule(static)
- for (int i = start; i < start+homenr; i++)
+ for (int i = 0; i < md.homenr; i++)
{
int g = 0;
- real m0 = md->massT[i];
- if (md->cVCM)
+ real m0 = md.massT[i];
+ if (md.cVCM)
{
- g = md->cVCM[i];
+ g = md.cVCM[i];
}
t_vcm_thread *vcm_t = &vcm->thread_vcm[t*vcm->stride + g];
/* Calculate linear momentum */
}
}
-void do_stopcm_grp(const t_mdatoms &mdatoms,
- rvec x[], rvec v[], const t_vcm &vcm)
+static void do_stopcm_grp(const t_mdatoms &mdatoms,
+ rvec x[], rvec v[], const t_vcm &vcm)
{
if (vcm.mode != ecmNO)
{
}
}
-void check_cm_grp(FILE *fp, t_vcm *vcm, t_inputrec *ir, real Temp_Max)
+/* Processes VCM after reduction over ranks and prints warning with high VMC and fp != nullptr */
+static void process_and_check_cm_grp(FILE *fp, t_vcm *vcm, real Temp_Max)
{
int m, g;
real ekcm, ekrot, tm, tm_1, Temp_cm;
if (vcm->mode == ecmANGULAR)
{
ekrot = 0.5*iprod(vcm->group_j[g], vcm->group_w[g]);
- if ((ekrot > 1) && fp && !EI_RANDOM(ir->eI))
+ // TODO: Change absolute energy comparison to relative
+ if ((ekrot > 1) && fp && vcm->integratorConservesMomentum)
{
/* if we have an integrator that may not conserve momenta, skip */
tm = vcm->group_mass[g];
}
}
}
+
+void process_and_stopcm_grp(FILE *fplog,
+ t_vcm *vcm,
+ const t_mdatoms &mdatoms,
+ rvec x[], rvec v[])
+{
+ if (vcm->mode != ecmNO)
+ {
+ // TODO: Replace fixed temperature of 1 by a system value
+ process_and_check_cm_grp(fplog, vcm, 1);
+
+ do_stopcm_grp(mdatoms, x, v, *vcm);
+ }
+}
ivec *nFreeze; /* Tells whether dimensions are frozen per freeze group */
std::vector<t_vcm_thread> thread_vcm; /* Temporary data per thread and group */
+ // Tell whether the integrator conserves momentum
+ bool integratorConservesMomentum;
+
t_vcm(const SimulationGroups &groups, const t_inputrec &ir);
~t_vcm();
};
/* Do a per group center of mass things */
-void calc_vcm_grp(int start, int homenr, t_mdatoms *md,
- rvec x[], rvec v[], t_vcm *vcm);
+void calc_vcm_grp(const t_mdatoms &md,
+ const rvec x[], const rvec v[], t_vcm *vcm);
/* Set the COM velocity to zero and potentially correct the COM position.
*
- * With linear modes nullptr can be passed for x.
+ * Processes the kinetic energy reduced over MPI before removing COM motion.
+ * With mode linear, nullptr can be passed for x.
+ * With acceleration correction nullptr should be passed for x at initialization
+ * and a pointer to the coordinates at normal MD steps.
+ * When fplog != nullptr, a warning is printed to fplog with high COM velocity.
*/
-void do_stopcm_grp(const t_mdatoms &mdatoms,
- rvec x[], rvec v[], const t_vcm &vcm);
-
-void check_cm_grp(FILE *fp, t_vcm *vcm, t_inputrec *ir, real Temp_Max);
+void process_and_stopcm_grp(FILE *fplog, t_vcm *vcm, const t_mdatoms &mdatoms,
+ rvec x[], rvec v[]);
#endif
restore_ekinstate_from_state(cr, ekind, &state_global->ekinstate);
}
- unsigned int cglo_flags = (CGLO_INITIALIZATION | CGLO_TEMPERATURE | CGLO_GSTAT
+ unsigned int cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
| (EI_VV(ir->eI) ? CGLO_PRESSURE : 0)
| (EI_VV(ir->eI) ? CGLO_CONSTRAINT : 0)
| (hasReadEkinState ? CGLO_READEKIN : 0));
cglo_flags_iteration |= CGLO_STOPCM;
cglo_flags_iteration &= ~CGLO_TEMPERATURE;
}
- compute_globals(fplog, gstat, cr, ir, fr, ekind,
+ compute_globals(gstat, cr, ir, fr, ekind,
state->x.rvec_array(), state->v.rvec_array(), state->box, state->lambda[efptVDW],
mdatoms, nrnb, &vcm,
nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
constr, &nullSignaller, state->box,
&totalNumberOfBondedInteractions, &bSumEkinhOld, cglo_flags_iteration
| (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0));
+ if (cglo_flags_iteration & CGLO_STOPCM)
+ {
+ /* At initialization, do not pass x with acceleration-correction mode
+ * to avoid (incorrect) correction of the initial coordinates.
+ */
+ rvec *xPtr = nullptr;
+ if (vcm.mode != ecmLINEAR_ACCELERATION_CORRECTION)
+ {
+ xPtr = state->x.rvec_array();
+ }
+ process_and_stopcm_grp(fplog, &vcm, *mdatoms, xPtr, state->v.rvec_array());
+ inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
+ }
}
checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
top_global, &top, state->x.rvec_array(), state->box,
kinetic energy calculation. This minimized excess variables, but
perhaps loses some logic?*/
- compute_globals(fplog, gstat, cr, ir, fr, ekind,
+ compute_globals(gstat, cr, ir, fr, ekind,
state->x.rvec_array(), state->v.rvec_array(), state->box, state->lambda[efptVDW],
mdatoms, nrnb, &vcm,
nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
/* We need the kinetic energy at minus the half step for determining
* the full step kinetic energy and possibly for T-coupling.*/
/* This may not be quite working correctly yet . . . . */
- compute_globals(fplog, gstat, cr, ir, fr, ekind,
+ compute_globals(gstat, cr, ir, fr, ekind,
state->x.rvec_array(), state->v.rvec_array(), state->box, state->lambda[efptVDW],
mdatoms, nrnb, &vcm,
wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
if (bGStat || do_per_step(step-1, nstglobalcomm))
{
wallcycle_stop(wcycle, ewcUPDATE);
- compute_globals(fplog, gstat, cr, ir, fr, ekind,
+ compute_globals(gstat, cr, ir, fr, ekind,
state->x.rvec_array(), state->v.rvec_array(), state->box, state->lambda[efptVDW],
mdatoms, nrnb, &vcm,
wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
top_global, &top, state->x.rvec_array(), state->box,
&shouldCheckNumberOfBondedInteractions);
+ if (bStopCM)
+ {
+ process_and_stopcm_grp(fplog, &vcm, *mdatoms, state->x.rvec_array(), state->v.rvec_array());
+ inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
+ }
wallcycle_start(wcycle, ewcUPDATE);
}
/* temperature scaling and pressure scaling to produce the extended variables at t+dt */
/* We need the kinetic energy at minus the half step for determining
* the full step kinetic energy and possibly for T-coupling.*/
/* This may not be quite working correctly yet . . . . */
- compute_globals(fplog, gstat, cr, ir, fr, ekind,
+ compute_globals(gstat, cr, ir, fr, ekind,
state->x.rvec_array(), state->v.rvec_array(), state->box, state->lambda[efptVDW],
mdatoms, nrnb, &vcm,
wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
{
/* erase F_EKIN and F_TEMP here? */
/* just compute the kinetic energy at the half step to perform a trotter step */
- compute_globals(fplog, gstat, cr, ir, fr, ekind,
+ compute_globals(gstat, cr, ir, fr, ekind,
state->x.rvec_array(), state->v.rvec_array(), state->box, state->lambda[efptVDW],
mdatoms, nrnb, &vcm,
wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
bool doIntraSimSignal = true;
SimulationSignaller signaller(&signals, cr, ms, doInterSimSignal, doIntraSimSignal);
- compute_globals(fplog, gstat, cr, ir, fr, ekind,
+ compute_globals(gstat, cr, ir, fr, ekind,
state->x.rvec_array(), state->v.rvec_array(), state->box, state->lambda[efptVDW],
mdatoms, nrnb, &vcm,
wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
top_global, &top, state->x.rvec_array(), state->box,
&shouldCheckNumberOfBondedInteractions);
+ if (!EI_VV(ir->eI) && bStopCM)
+ {
+ process_and_stopcm_grp(fplog, &vcm, *mdatoms, state->x.rvec_array(), state->v.rvec_array());
+ inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
+ }
}
}
}
{
- int cglo_flags = (CGLO_INITIALIZATION | CGLO_GSTAT |
+ int cglo_flags = (CGLO_GSTAT |
(shouldCheckNumberOfBondedInteractions ?
CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0));
bool bSumEkinhOld = false;
t_vcm *vcm = nullptr;
- compute_globals(fplog, gstat, cr, ir, fr, ekind,
+ compute_globals(gstat, cr, ir, fr, ekind,
state->x.rvec_array(), state->v.rvec_array(), state->box, state->lambda[efptVDW],
mdatoms, nrnb, vcm,
nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
t_vcm *vcm = nullptr;
SimulationSignaller signaller(&signals, cr, ms, doInterSimSignal, doIntraSimSignal);
- compute_globals(fplog, gstat, cr, ir, fr, ekind,
+ compute_globals(gstat, cr, ir, fr, ekind,
state->x.rvec_array(), state->v.rvec_array(), state->box, state->lambda[efptVDW],
mdatoms, nrnb, vcm,
wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
}
{
- int cglo_flags = (CGLO_INITIALIZATION | CGLO_GSTAT |
+ int cglo_flags = (CGLO_GSTAT |
(shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0));
bool bSumEkinhOld = false;
t_vcm *vcm = nullptr;
- compute_globals(fplog, gstat, cr, ir, fr, ekind,
+ compute_globals(gstat, cr, ir, fr, ekind,
state->x.rvec_array(), state->v.rvec_array(), state->box, state->lambda[efptVDW],
mdatoms, nrnb, vcm,
nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
t_vcm *vcm = nullptr;
SimulationSignaller signaller(&signals, cr, ms, doInterSimSignal, doIntraSimSignal);
- compute_globals(fplog, gstat, cr, ir, fr, ekind,
+ compute_globals(gstat, cr, ir, fr, ekind,
state->x.rvec_array(), state->v.rvec_array(), state->box, state->lambda[efptVDW],
mdatoms, nrnb, vcm,
wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff