Fixes issue with missing GMX_ASSERT as well.
constraining with LINCS and SHAKE still introduces significant drift,
which limits the system size to 100 to 200 nm).
+mdrun now reports energy drift
+""""""""""""""""""""""""""""""
+
+With conservative integrators, mdrun now reports the drift of the conserved
+energy quantity in the log file.
+
FEP using AWH
"""""""""""""
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, 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 Implements functions from the EnergyDriftTracker class.
+ *
+ * \author Berk Hess <hess@kth.se>
+ *
+ * \ingroup module_mdlib
+ */
+#include "gmxpre.h"
+
+#include "energydrifttracker.h"
+
+#include <cmath>
+
+#include <string>
+
+#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/real.h"
+#include "gromacs/utility/stringutil.h"
+
+namespace gmx
+{
+
+void EnergyDriftTracker::addPoint(double time, double energy)
+{
+ GMX_ASSERT(std::isfinite(energy), "Non-finite energy encountered!");
+
+ if (!storedFirst_)
+ {
+ firstTime_ = time;
+ firstEnergy_ = energy;
+ storedFirst_ = true;
+ }
+ lastTime_ = time;
+ lastEnergy_ = energy;
+}
+
+double EnergyDriftTracker::energyDrift() const
+{
+ if (timeInterval() > 0)
+ {
+ return (lastEnergy_ - firstEnergy_) / (timeInterval() * numAtoms_);
+ }
+ else
+ {
+ return 0;
+ }
+}
+
+std::string EnergyDriftTracker::energyDriftString(const std::string& partName) const
+{
+ std::string mesg;
+
+ if (timeInterval() > 0)
+ {
+ mesg = formatString("Energy conservation over %s of length %g ns, time %g to %g ns\n",
+ partName.c_str(), timeInterval(), firstTime_, lastTime_);
+ mesg += formatString(" Conserved energy drift: %.2e kJ/mol/ps per atom\n", energyDrift());
+ }
+ else
+ {
+ mesg = formatString(
+ "Time interval for measuring conserved energy has length 0, time %g to %g\n",
+ firstTime_, lastTime_);
+ }
+
+ return mesg;
+}
+
+} // namespace gmx
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, 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 Declares and defines the EnergyDriftTracker class.
+ *
+ * \author Berk Hess <hess@kth.se>
+ *
+ * \ingroup module_mdlib
+ */
+#ifndef GMX_MDLIB_ENERGYDRIFTTRACKER_H
+#define GMX_MDLIB_ENERGYDRIFTTRACKER_H
+
+#include <string>
+
+#include "gromacs/utility/real.h"
+
+namespace gmx
+{
+
+/*! \internal
+ * \brief Class for tracking and printing the drift in the conserved energy quantity
+ */
+class EnergyDriftTracker
+{
+public:
+ /*! \brief Constructor
+ *
+ * \param[in] numAtoms The total number of atoms in the system
+ */
+ EnergyDriftTracker(int numAtoms) : numAtoms_(numAtoms) {}
+
+ //! Add a point to the conserved energy tracking
+ void addPoint(double time, double energy);
+
+ //! Returns the time of the last point minus the time of the first point
+ double timeInterval() const { return lastTime_ - firstTime_; }
+
+ //! Returns the energy drift over the measured interval
+ double energyDrift() const;
+
+ /*! \brief Returns two-line string with the time interval and drift over the interval
+ *
+ * \param[in] partName A descriptive name for the period over which the tracking occured
+ */
+ std::string energyDriftString(const std::string& partName) const;
+
+private:
+ //! Whether we stored the first point
+ bool storedFirst_ = false;
+ //! The first time stored
+ double firstTime_ = 0;
+ //! The energy for the first time point
+ double firstEnergy_ = 0;
+ //! The last time stored
+ double lastTime_ = 0;
+ //! The energy for the last time point
+ double lastEnergy_ = 0;
+ //! The number of atoms in the system
+ int numAtoms_;
+};
+
+} // namespace gmx
+
+#endif
#include "gromacs/utility/smalloc.h"
#include "gromacs/utility/stringutil.h"
+#include "energydrifttracker.h"
+
//! Labels for energy file quantities
//! \{
static const char* conrmsd_nm[] = { "Constr. rmsd", "Constr.2 rmsd" };
FILE* fp_dhdl,
bool isRerun,
const StartingBehavior startingBehavior,
+ const bool simulationsShareState,
const MdModulesNotifier& mdModulesNotifier)
{
const char* ener_nm[F_NRE];
{
numTemperatures_ = 0;
}
+
+ if (EI_MD(ir->eI) && !simulationsShareState)
+ {
+ conservedEnergyTracker_ = std::make_unique<EnergyDriftTracker>(mtop->natoms);
+ }
}
EnergyOutput::~EnergyOutput()
store_dhdl, dE_ + fep->lambda_start_n, time);
}
}
+
+ if (conservedEnergyTracker_)
+ {
+ conservedEnergyTracker_->addPoint(
+ time, bEner_[F_ECONSERVED] ? enerd->term[F_ECONSERVED] : enerd->term[F_ETOT]);
+ }
}
void EnergyOutput::recordNonEnergyStep()
return ebin_->nener;
}
+void EnergyOutput::printEnergyConservation(FILE* fplog, int simulationPart, bool usingMdIntegrator) const
+{
+ if (fplog == nullptr)
+ {
+ return;
+ }
+
+ if (conservedEnergyTracker_)
+ {
+ std::string partName = formatString("simulation part #%d", simulationPart);
+ fprintf(fplog, "\n%s\n", conservedEnergyTracker_->energyDriftString(partName).c_str());
+ }
+ else if (usingMdIntegrator)
+ {
+ fprintf(fplog,
+ "\nCannot report drift of the conserved energy quantity because simulations share "
+ "state\n\n");
+ }
+}
+
} // namespace gmx
#include <cstdio>
+#include <memory>
+
#include "gromacs/mdtypes/enerdata.h"
class energyhistory_t;
namespace gmx
{
+class EnergyDriftTracker;
/*! \internal
* \brief Arrays connected to Pressure and Temperature coupling
* \param[in] fp_dhdl FEP file.
* \param[in] isRerun Is this is a rerun instead of the simulations.
* \param[in] startingBehavior Run starting behavior.
+ * \param[in] simulationsShareState Tells whether the physical state is shared over simulations
* \param[in] mdModulesNotifier Notifications to MD modules.
*/
EnergyOutput(ener_file* fp_ene,
FILE* fp_dhdl,
bool isRerun,
StartingBehavior startingBehavior,
+ bool simulationsShareState,
const MdModulesNotifier& mdModulesNotifier);
~EnergyOutput();
//! Print an output header to the log file.
static void printHeader(FILE* log, int64_t steps, double time);
+ /*! \brief Print conserved energy drift message to \p fplog
+ *
+ * Note that this is only over the current run (times are printed),
+ * this is not from the original start time for runs with continuation.
+ * This has the advantage that one can find if conservation issues are
+ * from the current run with the current settings on the current hardware.
+ */
+ void printEnergyConservation(FILE* fplog, int simulationPart, bool usingMdIntegrator) const;
+
private:
//! Timestep
double delta_t_ = 0;
real* temperatures_ = nullptr;
//! Number of temperatures actually saved
int numTemperatures_ = 0;
+
+ //! For tracking the conserved or total energy
+ std::unique_ptr<EnergyDriftTracker> conservedEnergyTracker_;
};
} // namespace gmx
constrtestdata.cpp
constrtestrunners.cpp
ebin.cpp
+ energydrifttracker.cpp
energyoutput.cpp
freeenergyparameters.cpp
leapfrog.cpp
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, 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 Tests for the EnergyDriftTacker class
+ *
+ * \author berk Hess <hess@kth.se>
+ * \ingroup module_mdlib
+ */
+#include "gmxpre.h"
+
+#include "gromacs/mdlib/energydrifttracker.h"
+
+#include <gtest/gtest.h>
+
+#include "testutils/testasserts.h"
+
+namespace gmx
+{
+
+namespace
+{
+
+TEST(EnergyDriftTracker, emptyWorks)
+{
+ EnergyDriftTracker tracker(1);
+
+ EXPECT_EQ(tracker.timeInterval(), 0);
+ EXPECT_EQ(tracker.energyDrift(), 0);
+}
+
+TEST(EnergyDriftTracker, onePointWorks)
+{
+ EnergyDriftTracker tracker(1);
+
+ tracker.addPoint(1.5, -3.5_real);
+ EXPECT_EQ(tracker.timeInterval(), 0);
+ EXPECT_EQ(tracker.energyDrift(), 0);
+}
+
+TEST(EnergyDriftTracker, manyPointsWorks)
+{
+ EnergyDriftTracker tracker(10);
+
+ tracker.addPoint(1.5, 2.5_real);
+ tracker.addPoint(3.5, 4.0_real);
+ tracker.addPoint(5.5, -5.5_real);
+ EXPECT_FLOAT_EQ(tracker.timeInterval(), 4.0_real);
+ EXPECT_FLOAT_EQ(tracker.energyDrift(), -0.2_real);
+}
+
+} // namespace
+
+} // namespace gmx
MdModulesNotifier mdModulesNotifier;
std::unique_ptr<EnergyOutput> energyOutput = std::make_unique<EnergyOutput>(
energyFile_, &mtop_, &inputrec_, nullptr, nullptr, parameters.isRerun,
- StartingBehavior::NewSimulation, mdModulesNotifier);
+ StartingBehavior::NewSimulation, false, mdModulesNotifier);
// Add synthetic data for a single step
double testValue = 10.0;
init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider, mdModulesNotifier, ir,
top_global, oenv, wcycle, startingBehavior, simulationsShareState, ms);
gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), top_global, ir, pull_work,
- mdoutf_get_fp_dhdl(outf), false, startingBehavior, mdModulesNotifier);
+ mdoutf_get_fp_dhdl(outf), false, startingBehavior,
+ simulationsShareState, mdModulesNotifier);
gstat = global_stat_init(ir);
{
if (ir->nstcalcenergy > 0)
{
+ energyOutput.printEnergyConservation(fplog, ir->simulation_part, EI_MD(ir->eI));
+
gmx::EnergyOutput::printAnnealingTemperatures(fplog, groups, &(ir->opts));
energyOutput.printAverages(fplog, groups);
}
StartingBehavior::NewSimulation, simulationsShareState, ms);
gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), top_global, ir, pull_work,
mdoutf_get_fp_dhdl(outf), true, StartingBehavior::NewSimulation,
- mdModulesNotifier);
+ simulationsShareState, mdModulesNotifier);
gstat = global_stat_init(ir);
gmx_mdoutf* outf = init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider,
mdModulesNotifier, inputrec, top_global, nullptr, wcycle,
StartingBehavior::NewSimulation, simulationsShareState, ms);
- gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), top_global, inputrec, pull_work, nullptr,
- false, StartingBehavior::NewSimulation, mdModulesNotifier);
+ gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), top_global, inputrec, pull_work,
+ nullptr, false, StartingBehavior::NewSimulation,
+ simulationsShareState, mdModulesNotifier);
/* Print to log file */
print_em_start(fplog, cr, walltime_accounting, wcycle, CG);
gmx_mdoutf* outf = init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider,
mdModulesNotifier, inputrec, top_global, nullptr, wcycle,
StartingBehavior::NewSimulation, simulationsShareState, ms);
- gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), top_global, inputrec, pull_work, nullptr,
- false, StartingBehavior::NewSimulation, mdModulesNotifier);
+ gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), top_global, inputrec, pull_work,
+ nullptr, false, StartingBehavior::NewSimulation,
+ simulationsShareState, mdModulesNotifier);
start = 0;
end = mdatoms->homenr;
gmx_mdoutf* outf = init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider,
mdModulesNotifier, inputrec, top_global, nullptr, wcycle,
StartingBehavior::NewSimulation, simulationsShareState, ms);
- gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), top_global, inputrec, pull_work, nullptr,
- false, StartingBehavior::NewSimulation, mdModulesNotifier);
+ gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), top_global, inputrec, pull_work,
+ nullptr, false, StartingBehavior::NewSimulation,
+ simulationsShareState, mdModulesNotifier);
/* Print to log file */
print_em_start(fplog, cr, walltime_accounting, wcycle, SD);
StartingBehavior::NewSimulation, simulationsShareState, ms);
gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), top_global, ir, pull_work,
mdoutf_get_fp_dhdl(outf), true, StartingBehavior::NewSimulation,
- mdModulesNotifier);
+ simulationsShareState, mdModulesNotifier);
gstat = global_stat_init(ir);
const MdModulesNotifier& mdModulesNotifier,
bool isMasterRank,
ObservablesHistory* observablesHistory,
- StartingBehavior startingBehavior) :
+ StartingBehavior startingBehavior,
+ bool simulationsShareState) :
element_(std::make_unique<Element>(this, isMasterRank)),
isMasterRank_(isMasterRank),
forceVirialStep_(-1),
fcd_(fcd),
mdModulesNotifier_(mdModulesNotifier),
groups_(&globalTopology->groups),
- observablesHistory_(observablesHistory)
+ observablesHistory_(observablesHistory),
+ simulationsShareState_(simulationsShareState)
{
clear_mat(forceVirial_);
clear_mat(shakeVirial_);
void EnergyData::setup(gmx_mdoutf* outf)
{
pull_t* pull_work = nullptr;
- energyOutput_ = std::make_unique<EnergyOutput>(mdoutf_get_fp_ene(outf), top_global_, inputrec_,
- pull_work, mdoutf_get_fp_dhdl(outf), false,
- startingBehavior_, mdModulesNotifier_);
+ energyOutput_ = std::make_unique<EnergyOutput>(
+ mdoutf_get_fp_ene(outf), top_global_, inputrec_, pull_work, mdoutf_get_fp_dhdl(outf),
+ false, startingBehavior_, simulationsShareState_, mdModulesNotifier_);
if (!isMasterRank_)
{
const MdModulesNotifier& mdModulesNotifier,
bool isMasterRank,
ObservablesHistory* observablesHistory,
- StartingBehavior startingBehavior);
+ StartingBehavior startingBehavior,
+ bool simulationsShareState);
/*! \brief Final output
*
const SimulationGroups* groups_;
//! History of simulation observables.
ObservablesHistory* observablesHistory_;
+ //! Whether simulations share the state
+ bool simulationsShareState_;
};
/*! \internal
opt2fn("-c", legacySimulatorData->nfile, legacySimulatorData->fnm), legacySimulatorData->inputrec,
legacySimulatorData->mdAtoms->mdatoms(), legacySimulatorData->top_global);
+ // Multi sim is turned off
+ const bool simulationsShareState = false;
+
energyData_ = std::make_unique<EnergyData>(
- statePropagatorData_.get(), freeEnergyPerturbationData_.get(),
- legacySimulatorData->top_global, legacySimulatorData->inputrec, legacySimulatorData->mdAtoms,
- legacySimulatorData->enerd, legacySimulatorData->ekind, legacySimulatorData->constr,
- legacySimulatorData->fplog, legacySimulatorData->fr->fcdata.get(),
- legacySimulatorData->mdModulesNotifier, MASTER(legacySimulatorData->cr),
- legacySimulatorData->observablesHistory, legacySimulatorData->startingBehavior);
+ statePropagatorData_.get(), freeEnergyPerturbationData_.get(), legacySimulatorData->top_global,
+ legacySimulatorData->inputrec, legacySimulatorData->mdAtoms, legacySimulatorData->enerd,
+ legacySimulatorData->ekind, legacySimulatorData->constr, legacySimulatorData->fplog,
+ legacySimulatorData->fr->fcdata.get(), legacySimulatorData->mdModulesNotifier,
+ MASTER(legacySimulatorData->cr), legacySimulatorData->observablesHistory,
+ legacySimulatorData->startingBehavior, simulationsShareState);
}
ModularSimulatorAlgorithm ModularSimulatorAlgorithmBuilder::build()
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff