KUBERNETES_CPU_REQUEST: 1
KUBERNETES_MEMORY_REQUEST: 2Gi
script:
- # TODO (issue #3272) `master` is not appropriate for use on release-xxxx branches, how should we handle that?
- - REV=$(git fetch -q https://gitlab.com/gromacs/gromacs.git release-2021 && git show -s --pretty=format:"%h" `git merge-base FETCH_HEAD HEAD`)
+ - REV=$(git fetch -q https://gitlab.com/gromacs/gromacs.git master && git show -s --pretty=format:"%h" `git merge-base FETCH_HEAD HEAD`)
- HEAD_REV=$(git show -s --pretty=format:"%h" HEAD)
- if [[ "$REV" == "$HEAD_REV" ]] ; then
REV="HEAD~1" ;
#
# This file is part of the GROMACS molecular simulation package.
#
--# Copyright (c) 2020, by the GROMACS development team, led by
++# Copyright (c) 2020,2021, 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.
# Set up to use the libstdc++ from that g++. Note that we checked
# the existing contents of CMAKE_CXX_FLAGS* variables earlier, so
# we will not override any user settings here.
- if (CMAKE_CXX_COMPILER_ID MATCHES "Clang")
+ if (CMAKE_CXX_COMPILER_ID MATCHES "Clang" OR CMAKE_CXX_COMPILER_ID MATCHES "IntelLLVM")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} --gcc-toolchain=${GMX_GPLUSPLUS_PATH}")
- else() #Intel
- set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -gcc-name=${GMX_GPLUSPLUS_PATH}")
endif()
endif()
# all their stuff. It's not easy if you only want some of their
# stuff...
set(MKL_MANUALLY FALSE)
- if (GMX_FFT_LIBRARY STREQUAL "MKL" AND NOT GMX_ICC_NEXTGEN)
-if (GMX_FFT_LIBRARY STREQUAL "MKL" AND
- NOT ((CMAKE_C_COMPILER_ID MATCHES "Intel" AND CMAKE_C_COMPILER_VERSION VERSION_GREATER "11")
- OR GMX_INTEL_LLVM))
++if (GMX_FFT_LIBRARY STREQUAL "MKL" AND NOT GMX_INTEL_LLVM)
# The user will have to provide the set of magic libraries in
# MKL_LIBRARIES (see below), which we cache (non-advanced), so that they
# don't have to keep specifying it, and can easily see that
how-to/visualize.rst
install-guide/index.rst
release-notes/index.rst
+ release-notes/2022/major/highlights.rst
+ release-notes/2022/major/features.rst
+ release-notes/2022/major/performance.rst
+ release-notes/2022/major/tools.rst
+ release-notes/2022/major/bugs-fixed.rst
+ release-notes/2022/major/removed-functionality.rst
+ release-notes/2022/major/deprecated-functionality.rst
+ release-notes/2022/major/portability.rst
+ release-notes/2022/major/miscellaneous.rst
+ release-notes/2022/major/api.rst
+ release-notes/2021/2021.3.rst
+ release-notes/2021/2021.2.rst
release-notes/2021/2021.1.rst
release-notes/2021/major/highlights.rst
release-notes/2021/major/features.rst
* To help us fund GROMACS development, we humbly ask that you cite
* the research papers on the package. Check out http://www.gromacs.org.
*/
++
#include "gmxpre.h"
#include "read_params.h"
++#include <algorithm>
++
#include "gromacs/applied_forces/awh/awh.h"
#include "gromacs/fileio/readinp.h"
#include "gromacs/fileio/warninp.h"
Awh::registerAwhWithPull(*awhParams, pull_work);
}
- for (int k = 0; k < awhParams.numBias(); k++)
+void checkAwhParams(const AwhParams& awhParams, const t_inputrec& ir, warninp_t wi)
+{
+ std::string opt;
+ checkMtsConsistency(ir, wi);
+
+ opt = "awh-nstout";
+ if (awhParams.nstout() <= 0)
+ {
+ auto message = formatString("Not writing AWH output with AWH (%s = %d) does not make sense",
+ opt.c_str(),
+ awhParams.nstout());
+ warning_error(wi, message);
+ }
+ /* This restriction can be removed by changing a flag of print_ebin() */
+ if (ir.nstenergy == 0 || awhParams.nstout() % ir.nstenergy != 0)
+ {
+ auto message = formatString(
+ "%s (%d) should be a multiple of nstenergy (%d)", opt.c_str(), awhParams.nstout(), ir.nstenergy);
+ warning_error(wi, message);
+ }
+
+ opt = "awh-nsamples-update";
+ if (awhParams.numSamplesUpdateFreeEnergy() <= 0)
+ {
+ warning_error(wi, opt + " needs to be an integer > 0");
+ }
+
- checkBiasParams(awhParams.awhBiasParams()[k], prefixawh, ir, wi);
++ bool haveFepLambdaDim = false;
++ const auto awhBiasParams = awhParams.awhBiasParams();
++ for (int k = 0; k < awhParams.numBias() && !haveFepLambdaDim; k++)
+ {
+ std::string prefixawh = formatString("awh%d", k + 1);
++ checkBiasParams(awhBiasParams[k], prefixawh, ir, wi);
++ /* Check if there is a FEP lambda dimension. */
++ const auto dimParams = awhBiasParams[k].dimParams();
++ haveFepLambdaDim = std::any_of(dimParams.begin(), dimParams.end(), [](const auto& dimParam) {
++ return dimParam.coordinateProvider() == AwhCoordinateProviderType::FreeEnergyLambda;
++ });
++ }
++
++ if (haveFepLambdaDim)
++ {
++ if (awhParams.nstSampleCoord() % ir.nstcalcenergy != 0)
++ {
++ opt = "awh-nstsample";
++ auto message = formatString(
++ "%s (%d) should be a multiple of nstcalcenergy (%d) when using AWH for "
++ "sampling an FEP lambda dimension",
++ opt.c_str(),
++ awhParams.nstSampleCoord(),
++ ir.nstcalcenergy);
++ warning_error(wi, message);
++ }
++ if (awhParams.potential() != AwhPotentialType::Umbrella)
++ {
++ opt = "awh-potential";
++ auto message = formatString(
++ "%s (%s) must be set to %s when using AWH for sampling an FEP lambda dimension",
++ opt.c_str(),
++ enumValueToString(awhParams.potential()),
++ enumValueToString(AwhPotentialType::Umbrella));
++ warning_error(wi, message);
++ }
+ }
+
+ if (ir.init_step != 0)
+ {
+ warning_error(wi, "With AWH init-step should be 0");
+ }
+}
+
} // namespace gmx
*/
INSTANTIATE_TEST_CASE_P(WithParameters,
BiasFepLambdaStateTest,
- ::testing::Combine(::testing::Values(eawhgrowthLINEAR, eawhgrowthEXP_LINEAR),
- ::testing::Values(eawhpotentialUMBRELLA),
+ ::testing::Combine(::testing::Values(AwhHistogramGrowthType::Linear,
+ AwhHistogramGrowthType::ExponentialLinear),
- ::testing::Values(AwhPotentialType::Umbrella,
- AwhPotentialType::Convolved),
++ ::testing::Values(AwhPotentialType::Umbrella),
::testing::Values(BiasParams::DisableUpdateSkips::yes,
BiasParams::DisableUpdateSkips::no)));
// Test that we detect coverings and exit the initial stage at the correct step
TEST(BiasFepLambdaStateTest, DetectsCovering)
{
- const AwhFepLambdaStateTestParameters params =
- getAwhFepLambdaTestParameters(eawhgrowthEXP_LINEAR, eawhpotentialUMBRELLA);
+ constexpr AwhCoordinateProviderType coordinateProvider = AwhCoordinateProviderType::FreeEnergyLambda;
+ constexpr int coordIndex = 0;
+ constexpr double origin = 0;
+ constexpr double end = c_numLambdaStates - 1;
+ constexpr double period = 0;
+ constexpr double diffusion = 1e-4 / (0.12927243028700 * 2);
+ auto awhDimBuffer =
+ awhDimParamSerialized(coordinateProvider, coordIndex, origin, end, period, diffusion);
+ auto awhDimArrayRef = gmx::arrayRefFromArray(&awhDimBuffer, 1);
+ const AwhTestParameters params(getAwhTestParameters(AwhHistogramGrowthType::ExponentialLinear,
- AwhPotentialType::Convolved,
++ AwhPotentialType::Umbrella,
+ awhDimArrayRef,
+ false,
+ 0.4,
+ true,
+ 1.0,
+ c_numLambdaStates));
const double mdTimeStep = 0.1;
InteractionOfType::InteractionOfType(gmx::ArrayRef<const int> atoms,
gmx::ArrayRef<const real> params,
const std::string& name) :
-- atoms_(atoms.begin(), atoms.end()),
-- interactionTypeName_(name)
++ atoms_(atoms.begin(), atoms.end()), interactionTypeName_(name)
{
GMX_RELEASE_ASSERT(
params.size() <= forceParam_.size(),
/* check masses */
check_mol(&sys, wi);
+ if (haveFepPerturbedMassesInSettles(sys))
+ {
+ warning_error(wi,
+ "SETTLE is not implemented for atoms whose mass is perturbed. "
+ "You might instead use normal constraints.");
+ }
+
checkForUnboundAtoms(&sys, bVerbose, wi, logger);
- if (EI_DYNAMICS(ir->eI) && ir->eI != eiBD)
+ if (EI_DYNAMICS(ir->eI) && ir->eI != IntegrationAlgorithm::BD)
{
check_bonds_timestep(&sys, ir->delta_t, wi);
}
copy_mat(ir->compress, compressibility);
}
setStateDependentAwhParams(
- ir->awhParams, *ir->pull, pull, state.box, ir->pbcType, compressibility, &ir->opts,
- ir->efep != efepNO ? ir->fepvals->all_lambda[efptFEP][ir->fepvals->init_fep_state] : 0,
- sys, wi);
+ ir->awhParams.get(),
+ *ir->pull,
+ pull,
+ state.box,
+ ir->pbcType,
+ compressibility,
+ &ir->opts,
- ir->efep != FreeEnergyPerturbationType::No
- ? ir->fepvals->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Fep)]
- [ir->fepvals->init_fep_state]
- : 0,
++ ir->efep != FreeEnergyPerturbationType::No ? ir->fepvals->all_lambda[static_cast<int>(
++ FreeEnergyPerturbationCouplingType::Fep)][ir->fepvals->init_fep_state]
++ : 0,
+ sys,
+ wi);
}
if (ir->bPull)
struct RtpRename
{
RtpRename(const char* newGmx, const char* newMain, const char* newNter, const char* newCter, const char* newBter) :
-- gmx(newGmx),
-- main(newMain),
-- nter(newNter),
-- cter(newCter),
-- bter(newBter)
++ gmx(newGmx), main(newMain), nter(newNter), cter(newCter), bter(newBter)
{
}
std::string gmx;
j = 0;
while ((j < nr)
&& gmx_strcasecmp(
-- names[2 * i].c_str(),
-- *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][j]])))
++ names[2 * i].c_str(),
++ *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][j]])))
{
j++;
}
k = 0;
while ((k < nr)
&& gmx_strcasecmp(
-- names[2 * i + 1].c_str(),
-- *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][k]])))
++ names[2 * i + 1].c_str(),
++ *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][k]])))
{
k++;
}
}
}
- /* Check for position restraints */
- for (const auto ilist : IListRange(sys))
+ return absRef;
+ }
+
+ //! Returns whether position restraints are used for dimensions
+ static BasicVector<bool> havePositionRestraints(const gmx_mtop_t& sys)
+ {
+ BasicVector<bool> havePosres = { false, false, false };
+
- gmx_mtop_ilistloop_t iloop = gmx_mtop_ilistloop_init(&sys);
- int nmol;
- while (const InteractionLists* ilist = gmx_mtop_ilistloop_next(iloop, &nmol))
++ for (const auto ilists : IListRange(sys))
{
- if (ilist.nmol() > 0 && (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
- if (nmol > 0 && (!havePosres[XX] || !havePosres[YY] || !havePosres[ZZ]))
++ const auto& posResList = ilists.list()[F_POSRES];
++ const auto& fbPosResList = ilists.list()[F_FBPOSRES];
++ if (ilists.nmol() > 0 && (!havePosres[XX] || !havePosres[YY] || !havePosres[ZZ]))
{
- for (int i = 0; i < ilist.list()[F_POSRES].size(); i += 2)
- for (int i = 0; i < (*ilist)[F_POSRES].size(); i += 2)
++ for (int i = 0; i < posResList.size(); i += 2)
{
- pr = &sys.ffparams.iparams[ilist.list()[F_POSRES].iatoms[i]];
- for (d = 0; d < DIM; d++)
- const t_iparams& pr = sys.ffparams.iparams[(*ilist)[F_POSRES].iatoms[i]];
++ const t_iparams& pr = sys.ffparams.iparams[posResList.iatoms[i]];
+ for (int d = 0; d < DIM; d++)
{
- if (pr->posres.fcA[d] != 0)
+ if (pr.posres.fcA[d] != 0)
{
- AbsRef[d] = 1;
+ havePosres[d] = true;
}
}
}
- for (i = 0; i < ilist.list()[F_FBPOSRES].size(); i += 2)
- for (int i = 0; i < (*ilist)[F_FBPOSRES].size(); i += 2)
++ for (int i = 0; i < fbPosResList.size(); i += 2)
{
/* Check for flat-bottom posres */
- pr = &sys.ffparams.iparams[ilist.list()[F_FBPOSRES].iatoms[i]];
- if (pr->fbposres.k != 0)
- const t_iparams& pr = sys.ffparams.iparams[(*ilist)[F_FBPOSRES].iatoms[i]];
++ const t_iparams& pr = sys.ffparams.iparams[fbPosResList.iatoms[i]];
+ if (pr.fbposres.k != 0)
{
- switch (pr->fbposres.geom)
+ switch (pr.fbposres.geom)
{
- case efbposresSPHERE: AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1; break;
- case efbposresCYLINDERX: AbsRef[YY] = AbsRef[ZZ] = 1; break;
- case efbposresCYLINDERY: AbsRef[XX] = AbsRef[ZZ] = 1; break;
+ case efbposresSPHERE: havePosres = { true, true, true }; break;
+ case efbposresCYLINDERX: havePosres[YY] = havePosres[ZZ] = true; break;
+ case efbposresCYLINDERY: havePosres[XX] = havePosres[ZZ] = true; break;
case efbposresCYLINDER:
/* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
- case efbposresCYLINDERZ: AbsRef[XX] = AbsRef[YY] = 1; break;
+ case efbposresCYLINDERZ: havePosres[XX] = havePosres[YY] = true; break;
case efbposresX: /* d=XX */
case efbposresY: /* d=YY */
case efbposresZ: /* d=ZZ */
break;
default:
gmx_fatal(FARGS,
- " Invalid geometry for flat-bottom position restraint.\n"
+ "Invalid geometry for flat-bottom position restraint.\n"
"Expected nr between 1 and %d. Found %d\n",
- efbposresNR - 1, pr.fbposres.geom);
+ efbposresNR - 1,
- pr->fbposres.geom);
++ pr.fbposres.geom);
}
}
}
}
}
- return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
+ return havePosres;
}
-static void check_combination_rule_differences(const gmx_mtop_t* mtop,
+static void check_combination_rule_differences(const gmx_mtop_t& mtop,
int state,
bool* bC6ParametersWorkWithGeometricRules,
bool* bC6ParametersWorkWithLBRules,
ptr = getenv("GMX_LJCOMB_TOL");
if (ptr != nullptr)
{
-- double dbl;
++ double dbl;
double gmx_unused canary;
if (sscanf(ptr, "%lf%lf", &dbl, &canary) != 1)
char err_buf[STRLEN];
int i, m, c, nmol;
- bool bCharge, bAcc;
- real * mgrp, mt;
- rvec acc;
+ bool bCharge;
gmx_mtop_atomloop_block_t aloopb;
- ivec AbsRef;
char warn_buf[STRLEN];
set_warning_line(wi, mdparin, -1);
}
}
- if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD && ir->comm_mode == ecmNO
+ if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != IntegrationAlgorithm::BD
+ && ir->comm_mode == ComRemovalAlgorithm::No
- && !(absolute_reference(ir, *sys, FALSE, AbsRef) || ir->nsteps <= 10) && !ETC_ANDERSEN(ir->etc))
+ && !(allTrue(haveAbsoluteReference(*ir)) || allTrue(havePositionRestraints(*sys)) || ir->nsteps <= 10)
+ && !ETC_ANDERSEN(ir->etc))
{
warning(wi,
"You are not using center of mass motion removal (mdp option comm-mode), numerical "
}
/* Check for pressure coupling with absolute position restraints */
- if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
+ if (ir->epc != PressureCoupling::No && ir->refcoord_scaling == RefCoordScaling::No)
{
- absolute_reference(ir, *sys, TRUE, AbsRef);
+ const BasicVector<bool> havePosres = havePositionRestraints(*sys);
{
for (m = 0; m < DIM; m++)
{
#include "gromacs/mdtypes/commrec.h"
#include "gromacs/mdtypes/enerdata.h"
#include "gromacs/mdtypes/inputrec.h"
++#include "gromacs/utility/enumerationhelpers.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/smalloc.h"
ForeignLambdaTerms::ForeignLambdaTerms(int numLambdas) :
-- numLambdas_(numLambdas),
-- energies_(1 + numLambdas),
-- dhdl_(1 + numLambdas)
++ numLambdas_(numLambdas), energies_(1 + numLambdas), dhdl_(1 + numLambdas)
{
}
}
gmx_enerdata_t::gmx_enerdata_t(int numEnergyGroups, int numFepLambdas) :
-- grpp(numEnergyGroups),
-- foreignLambdaTerms(numFepLambdas),
-- foreign_grpp(numEnergyGroups)
++ grpp(numEnergyGroups), foreignLambdaTerms(numFepLambdas), foreign_grpp(numEnergyGroups)
{
}
}
}
- // Adds computed dV/dlambda contributions to the requested outputs
- static void set_dvdl_output(gmx_enerdata_t* enerd, const t_lambda& fepvals)
+ // Adds computed dH/dlambda contribution i to the requested output
-static void set_dhdl_output(gmx_enerdata_t* enerd, int i, const t_lambda& fepvals)
++static void set_dhdl_output(gmx_enerdata_t* enerd, FreeEnergyPerturbationCouplingType i, const t_lambda& fepvals)
{
- enerd->term[F_DVDL] = 0.0;
- for (auto i : keysOf(fepvals.separate_dvdl))
+ if (fepvals.separate_dvdl[i])
{
- if (fepvals.separate_dvdl[i])
+ /* Translate free-energy term indices to idef term indices.
+ * Could this be done more readably/compactly?
+ */
+ int index;
+ switch (i)
{
- /* Translate free-energy term indices to idef term indices.
- * Could this be done more readably/compactly?
- */
- int index;
- switch (i)
- {
- case (FreeEnergyPerturbationCouplingType::Mass): index = F_DKDL; break;
- case (FreeEnergyPerturbationCouplingType::Coul): index = F_DVDL_COUL; break;
- case (FreeEnergyPerturbationCouplingType::Vdw): index = F_DVDL_VDW; break;
- case (FreeEnergyPerturbationCouplingType::Bonded): index = F_DVDL_BONDED; break;
- case (FreeEnergyPerturbationCouplingType::Restraint):
- index = F_DVDL_RESTRAINT;
- break;
- default: index = F_DVDL; break;
- }
- enerd->term[index] = enerd->dvdl_lin[i] + enerd->dvdl_nonlin[i];
- if (debug)
- {
- fprintf(debug,
- "dvdl-%s[%2d]: %f: non-linear %f + linear %f\n",
- enumValueToString(i),
- static_cast<int>(i),
- enerd->term[index],
- enerd->dvdl_nonlin[i],
- enerd->dvdl_lin[i]);
- }
- case (efptMASS): index = F_DKDL; break;
- case (efptCOUL): index = F_DVDL_COUL; break;
- case (efptVDW): index = F_DVDL_VDW; break;
- case (efptBONDED): index = F_DVDL_BONDED; break;
- case (efptRESTRAINT): index = F_DVDL_RESTRAINT; break;
++ case (FreeEnergyPerturbationCouplingType::Mass): index = F_DKDL; break;
++ case (FreeEnergyPerturbationCouplingType::Coul): index = F_DVDL_COUL; break;
++ case (FreeEnergyPerturbationCouplingType::Vdw): index = F_DVDL_VDW; break;
++ case (FreeEnergyPerturbationCouplingType::Bonded): index = F_DVDL_BONDED; break;
++ case (FreeEnergyPerturbationCouplingType::Restraint): index = F_DVDL_RESTRAINT; break;
+ default: index = F_DVDL; break;
}
- else
+ enerd->term[index] = enerd->dvdl_lin[i] + enerd->dvdl_nonlin[i];
+ if (debug)
{
- enerd->term[F_DVDL] += enerd->dvdl_lin[i] + enerd->dvdl_nonlin[i];
- if (debug)
- {
- fprintf(debug,
- "dvd-%sl[%2d]: %f: non-linear %f + linear %f\n",
- enumValueToString(FreeEnergyPerturbationCouplingType::Fep),
- static_cast<int>(i),
- enerd->term[F_DVDL],
- enerd->dvdl_nonlin[i],
- enerd->dvdl_lin[i]);
- }
- fprintf(debug, "dvdl-%s[%2d]: %f: non-linear %f + linear %f\n", efpt_names[i], i,
- enerd->term[index], enerd->dvdl_nonlin[i], enerd->dvdl_lin[i]);
++ fprintf(debug,
++ "dvdl-%s[%2d]: %f: non-linear %f + linear %f\n",
++ enumValueToString(i),
++ static_cast<int>(i),
++ enerd->term[index],
++ enerd->dvdl_nonlin[i],
++ enerd->dvdl_lin[i]);
+ }
+ }
+ else
+ {
+ enerd->term[F_DVDL] += enerd->dvdl_lin[i] + enerd->dvdl_nonlin[i];
+ if (debug)
+ {
- fprintf(debug, "dvd-%sl[%2d]: %f: non-linear %f + linear %f\n", efpt_names[0], i,
- enerd->term[F_DVDL], enerd->dvdl_nonlin[i], enerd->dvdl_lin[i]);
++ fprintf(debug,
++ "dvd-%sl[%2d]: %f: non-linear %f + linear %f\n",
++ enumValueToString(FreeEnergyPerturbationCouplingType::Fep),
++ static_cast<int>(i),
++ enerd->term[F_DVDL],
++ enerd->dvdl_nonlin[i],
++ enerd->dvdl_lin[i]);
}
}
}
if (fepvals)
{
- set_dvdl_output(enerd, *fepvals);
+ enerd->term[F_DVDL] = 0.0;
- for (int i = 0; i < efptNR; i++)
++ for (auto i : gmx::EnumerationWrapper<FreeEnergyPerturbationCouplingType>{})
+ {
+ // Skip kinetic terms here, as those are not available here yet
- if (i != efptMASS)
++ if (i != FreeEnergyPerturbationCouplingType::Mass)
+ {
+ set_dhdl_output(enerd, i, *fepvals);
+ }
+ }
enerd->foreignLambdaTerms.finalizePotentialContributions(enerd->dvdl_lin, lambda, *fepvals);
}
enerd->term[F_DVDL] += enerd->term[F_DVDL_CONSTR];
}
- set_dhdl_output(enerd, efptMASS, fepvals);
+ // Add computed mass dH/dlambda contribution to the requested output
- enerd->foreignLambdaTerms.finalizeKineticContributions(enerd->term, enerd->dvdl_lin[efptMASS],
- lambda, fepvals);
++ set_dhdl_output(enerd, FreeEnergyPerturbationCouplingType::Mass, fepvals);
+
+ enerd->foreignLambdaTerms.finalizeKineticContributions(
+ enerd->term, enerd->dvdl_lin[FreeEnergyPerturbationCouplingType::Mass], lambda, fepvals);
/* The constrain contribution is now included in other terms, so clear it */
enerd->term[F_DVDL_CONSTR] = 0;
f_global = of->f_global;
if (mdof_flags & MDOF_F)
{
- auto globalFRef =
- MASTER(cr) ? gmx::arrayRefFromArray(reinterpret_cast<gmx::RVec*>(of->f_global),
- of->natoms_global)
- : gmx::ArrayRef<gmx::RVec>();
- dd_collect_vec(
- cr->dd, state_local->ddp_count, state_local->ddp_count_cg_gl, state_local->cg_gl, f_local,
- gmx::arrayRefFromArray(reinterpret_cast<gmx::RVec*>(of->f_global), f_local.size()));
++ auto globalFRef = MASTER(cr) ? gmx::arrayRefFromArray(
++ reinterpret_cast<gmx::RVec*>(of->f_global), of->natoms_global)
++ : gmx::ArrayRef<gmx::RVec>();
+ dd_collect_vec(cr->dd,
+ state_local->ddp_count,
+ state_local->ddp_count_cg_gl,
+ state_local->cg_gl,
+ f_local,
+ globalFRef);
}
}
else
int64_t count)
{
-- t_state *s1, *s2;
-- int start, end;
-- real dvdl_constr;
++ t_state * s1, *s2;
++ int start, end;
++ real dvdl_constr;
int nthreads gmx_unused;
bool validStep = true;
const bool do_or = ir->nstorireout != 0;
EnergyOutput::printAnnealingTemperatures(do_log ? fplog : nullptr, groups, &(ir->opts));
- energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or,
- do_log ? fplog : nullptr, step, t, fr->fcdata.get(), awh);
+ energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf),
+ do_ene,
+ do_dr,
+ do_or,
+ do_log ? fplog : nullptr,
+ step,
+ t,
+ fr->fcdata.get(),
+ awh);
+ if (ir->bPull)
+ {
+ pull_print_output(pull_work, step, t);
+ }
+
if (do_per_step(step, ir->nstlog))
{
if (fflush(fplog) != 0)
isInputCompatible =
isInputCompatible
&& conditionalAssert(
-- !inputrec->useMts,
-- "Multiple time stepping is not supported by the modular simulator.");
++ !inputrec->useMts,
++ "Multiple time stepping is not supported by the modular simulator.");
isInputCompatible =
isInputCompatible
&& conditionalAssert(!doRerun, "Rerun is not supported by the modular simulator.");
isInputCompatible =
isInputCompatible
&& conditionalAssert(
- inputrec->epc == PressureCoupling::No
- || inputrec->epc == PressureCoupling::ParrinelloRahman,
- inputrec->epc == epcNO || inputrec->epc == epcPARRINELLORAHMAN,
-- "Only Parrinello-Rahman barostat is supported by the modular simulator.");
++ inputrec->epc == PressureCoupling::No || inputrec->epc == PressureCoupling::ParrinelloRahman,
++ "Only Parrinello-Rahman barostat is supported by the modular simulator.");
isInputCompatible =
isInputCompatible
&& conditionalAssert(
-- !(inputrecNptTrotter(inputrec) || inputrecNphTrotter(inputrec)
-- || inputrecNvtTrotter(inputrec)),
-- "Legacy Trotter decomposition is not supported by the modular simulator.");
- isInputCompatible = isInputCompatible
- && conditionalAssert(inputrec->efep == efepNO || inputrec->efep == efepYES
- || inputrec->efep == efepSLOWGROWTH,
- "Expanded ensemble free energy calculation is not "
- "supported by the modular simulator.");
++ !(inputrecNptTrotter(inputrec) || inputrecNphTrotter(inputrec)
++ || inputrecNvtTrotter(inputrec)),
++ "Legacy Trotter decomposition is not supported by the modular simulator.");
+ isInputCompatible =
+ isInputCompatible
+ && conditionalAssert(inputrec->efep == FreeEnergyPerturbationType::No
+ || inputrec->efep == FreeEnergyPerturbationType::Yes
+ || inputrec->efep == FreeEnergyPerturbationType::SlowGrowth,
+ "Expanded ensemble free energy calculation is not "
+ "supported by the modular simulator.");
isInputCompatible = isInputCompatible
&& conditionalAssert(!inputrec->bPull,
"Pulling is not supported by the modular simulator.");
isInputCompatible =
isInputCompatible
&& conditionalAssert(
-- inputrec->deform[XX][XX] == 0.0 && inputrec->deform[XX][YY] == 0.0
-- && inputrec->deform[XX][ZZ] == 0.0 && inputrec->deform[YY][XX] == 0.0
-- && inputrec->deform[YY][YY] == 0.0 && inputrec->deform[YY][ZZ] == 0.0
-- && inputrec->deform[ZZ][XX] == 0.0 && inputrec->deform[ZZ][YY] == 0.0
-- && inputrec->deform[ZZ][ZZ] == 0.0,
-- "Deformation is not supported by the modular simulator.");
++ inputrec->deform[XX][XX] == 0.0 && inputrec->deform[XX][YY] == 0.0
++ && inputrec->deform[XX][ZZ] == 0.0 && inputrec->deform[YY][XX] == 0.0
++ && inputrec->deform[YY][YY] == 0.0 && inputrec->deform[YY][ZZ] == 0.0
++ && inputrec->deform[ZZ][XX] == 0.0 && inputrec->deform[ZZ][YY] == 0.0
++ && inputrec->deform[ZZ][ZZ] == 0.0,
++ "Deformation is not supported by the modular simulator.");
isInputCompatible =
isInputCompatible
&& conditionalAssert(gmx_mtop_interaction_count(globalTopology, IF_VSITE) == 0,
isInputCompatible =
isInputCompatible
&& conditionalAssert(
-- gmx_mtop_ftype_count(globalTopology, F_ORIRES) == 0,
-- "Orientation restraints are not supported by the modular simulator.");
++ gmx_mtop_ftype_count(globalTopology, F_ORIRES) == 0,
++ "Orientation restraints are not supported by the modular simulator.");
isInputCompatible =
isInputCompatible
&& conditionalAssert(ms == nullptr,
}
else
{
-- auto distantRestraintEnsembleEnvVar = getenv("GMX_DISRE_ENSEMBLE_SIZE");
++ auto* distantRestraintEnsembleEnvVar = getenv("GMX_DISRE_ENSEMBLE_SIZE");
numEnsembleRestraintSystems =
(ms != nullptr && distantRestraintEnsembleEnvVar != nullptr)
? static_cast<int>(strtol(distantRestraintEnsembleEnvVar, nullptr, 10))
isInputCompatible =
isInputCompatible
&& conditionalAssert(
-- getenv("GMX_FORCE_UPDATE_DEFAULT_GPU") == nullptr,
-- "Integration on the GPU is not supported by the modular simulator.");
++ getenv("GMX_FORCE_UPDATE_DEFAULT_GPU") == nullptr,
++ "Integration on the GPU is not supported by the modular simulator.");
// Modular simulator is centered around NS updates
// TODO: think how to handle nstlist == 0
isInputCompatible = isInputCompatible
isInputCompatible = isInputCompatible
&& conditionalAssert(!GMX_FAHCORE,
"GMX_FAHCORE not supported by the modular simulator.");
- GMX_RELEASE_ASSERT(
- isInputCompatible
- || !(inputrec->eI == IntegrationAlgorithm::VV
- && inputrec->epc == PressureCoupling::ParrinelloRahman),
- "Requested Parrinello-Rahman barostat with md-vv, but other options are not compatible "
- "with the modular simulator. The Parrinello-Rahman barostat is not implemented for "
- "md-vv in the legacy simulator. Use a different pressure control algorithm.");
-
- if (!isInputCompatible && (inputrec->eI == eiVV && inputrec->epc == epcPARRINELLORAHMAN))
++ if (!isInputCompatible
++ && (inputrec->eI == IntegrationAlgorithm::VV && inputrec->epc == PressureCoupling::ParrinelloRahman))
+ {
+ gmx_fatal(FARGS,
+ "Requested Parrinello-Rahman barostat with md-vv. This combination is only "
+ "available in the modular simulator. Some other selected options are, however, "
+ "only available in the legacy simulator. Use a different pressure control "
+ "algorithm.");
+ }
-
return isInputCompatible;
}
return TRUE;
}
-typedef struct gmx_mtop_ilistloop
+IListIterator::IListIterator(const gmx_mtop_t& mtop, size_t molblock) :
- mtop_(&mtop),
- mblock_(molblock)
++ mtop_(&mtop), mblock_(molblock)
{
- const gmx_mtop_t* mtop;
- int mblock;
-} t_gmx_mtop_ilist;
+}
-gmx_mtop_ilistloop_t gmx_mtop_ilistloop_init(const gmx_mtop_t* mtop)
+IListIterator& IListIterator::operator++()
{
- struct gmx_mtop_ilistloop* iloop;
-
- snew(iloop, 1);
-
- iloop->mtop = mtop;
- iloop->mblock = -1;
-
- return iloop;
+ mblock_++;
+ return *this;
}
-gmx_mtop_ilistloop_t gmx_mtop_ilistloop_init(const gmx_mtop_t& mtop)
+bool IListIterator::operator==(const IListIterator& o) const
{
- return gmx_mtop_ilistloop_init(&mtop);
+ return mtop_ == o.mtop_ && mblock_ == o.mblock_;
}
-static void gmx_mtop_ilistloop_destroy(gmx_mtop_ilistloop_t iloop)
+const InteractionLists& IListProxy::list() const
{
- sfree(iloop);
+ // one past the end means we want to take the
+ // intermolecular list instead.
+ if (it_->mblock_ == it_->mtop_->molblock.size())
+ {
+ return *it_->mtop_->intermolecular_ilist;
+ }
+ else
+ {
+ return it_->mtop_->moltype[it_->mtop_->molblock[it_->mblock_].type].ilist;
+ }
}
-const InteractionLists* gmx_mtop_ilistloop_next(gmx_mtop_ilistloop_t iloop, int* nmol)
+int IListProxy::nmol() const
{
- if (iloop == nullptr)
+ // one past the end means we want to take the
+ // intermolecular list instead.
+ if (it_->mblock_ == it_->mtop_->molblock.size())
{
- gmx_incons("gmx_mtop_ilistloop_next called without calling gmx_mtop_ilistloop_init");
+ return 1;
}
-
- iloop->mblock++;
- if (iloop->mblock >= gmx::ssize(iloop->mtop->molblock))
+ else
{
- if (iloop->mblock == gmx::ssize(iloop->mtop->molblock) && iloop->mtop->bIntermolecularInteractions)
- {
- *nmol = 1;
- return iloop->mtop->intermolecular_ilist.get();
- }
-
- gmx_mtop_ilistloop_destroy(iloop);
- return nullptr;
+ return it_->mtop_->molblock[it_->mblock_].nmol;
}
-
- *nmol = iloop->mtop->molblock[iloop->mblock].nmol;
-
- return &iloop->mtop->moltype[iloop->mtop->molblock[iloop->mblock].type].ilist;
}
-typedef struct gmx_mtop_ilistloop_all
-{
- const gmx_mtop_t* mtop;
- size_t mblock;
- int mol;
- int a_offset;
-} t_gmx_mtop_ilist_all;
-int gmx_mtop_ftype_count(const gmx_mtop_t* mtop, int ftype)
+IListRange::IListRange(const gmx_mtop_t& mtop) : begin_(mtop), end_(mtop, mtop.molblock.size())
{
- gmx_mtop_ilistloop_t iloop;
- int n, nmol;
+ if (mtop.bIntermolecularInteractions)
+ {
+ end_ = IListIterator(mtop, mtop.molblock.size() + 1);
+ }
+}
- n = 0;
+int gmx_mtop_ftype_count(const gmx_mtop_t& mtop, int ftype)
+{
+ int n = 0;
- iloop = gmx_mtop_ilistloop_init(mtop);
- while (const InteractionLists* il = gmx_mtop_ilistloop_next(iloop, &nmol))
+ for (const IListProxy il : IListRange(mtop))
{
- n += nmol * (*il)[ftype].size() / (1 + NRAL(ftype));
+ n += il.nmol() * il.list()[ftype].size() / (1 + NRAL(ftype));
}
- if (mtop->bIntermolecularInteractions)
+ if (mtop.bIntermolecularInteractions)
{
- n += (*mtop->intermolecular_ilist)[ftype].size() / (1 + NRAL(ftype));
+ n += (*mtop.intermolecular_ilist)[ftype].size() / (1 + NRAL(ftype));
}
return n;
* KeyValueTreeObjectBuilder enables writing of module internal data to
* .tpr files.
*/
- BuildMDModulesNotifier<EnergyCalculationFrequencyErrors*, IndexGroupsAndNames, KeyValueTreeObjectBuilder>::type preProcessingNotifier_;
- registerMdModuleNotification<EnergyCalculationFrequencyErrors*, const IndexGroupsAndNames&, KeyValueTreeObjectBuilder>::type preProcessingNotifications_;
++ BuildMDModulesNotifier<EnergyCalculationFrequencyErrors*, const IndexGroupsAndNames&, KeyValueTreeObjectBuilder>::type preProcessingNotifier_;
- /*! \brief Checkpointing callback functions.
+ /*! \brief Handles subscribing and calling checkpointing callback functions.
*
- * MdModulesCheckpointReadingDataOnMaster provides modules with their
+ * MDModulesCheckpointReadingDataOnMaster provides modules with their
* checkpointed data on the master
* node and checkpoint file version
- * MdModulesCheckpointReadingBroadcast provides modules with a communicator
+ * MDModulesCheckpointReadingBroadcast provides modules with a communicator
* and the checkpoint file version to
* distribute their data
- * MdModulesWriteCheckpointData provides the modules with a key-value-tree
+ * MDModulesWriteCheckpointData provides the modules with a key-value-tree
* builder to store their checkpoint data and
* the checkpoint file version
*/
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff