From 3fbf6749766ad7e9143da87a84383f88b7a8bf17 Mon Sep 17 00:00:00 2001 From: Pascal Merz Date: Sun, 22 Nov 2020 20:30:00 -0700 Subject: [PATCH] Refactor temperature annealing This refactors temperature annealing by dividing it into two subfunctions, one which returns a new temperature, keeping the inputrec constant and being independent of the Update object, and one which updates the inputrec and the Update object. This separation allows most of the code to be reused by modular simulator, and it is also a first step towards a design which does not need to modify the input record (#3854) Refs #3417, #3854 --- src/gromacs/mdlib/coupling.cpp | 122 +++++++++++++++++---------------- src/gromacs/mdlib/coupling.h | 10 +++ 2 files changed, 74 insertions(+), 58 deletions(-) diff --git a/src/gromacs/mdlib/coupling.cpp b/src/gromacs/mdlib/coupling.cpp index e8c2badf1f..4c24f6350c 100644 --- a/src/gromacs/mdlib/coupling.cpp +++ b/src/gromacs/mdlib/coupling.cpp @@ -2142,74 +2142,80 @@ bool initSimulatedAnnealing(t_inputrec* ir, gmx::Update* upd) return doSimAnnealing; } -/* set target temperatures if we are annealing */ -void update_annealing_target_temp(t_inputrec* ir, real t, gmx::Update* upd) +real computeAnnealingTargetTemperature(const t_inputrec& inputrec, int temperatureGroup, real time) { - for (int temperatureGroup = 0; temperatureGroup < ir->opts.ngtc; temperatureGroup++) + GMX_RELEASE_ASSERT(temperatureGroup >= 0 && temperatureGroup < inputrec.opts.ngtc, + "Invalid temperature group."); + if (inputrec.opts.annealing[temperatureGroup] == SimulatedAnnealing::No) { - const auto& inputrec = *ir; - if (inputrec.opts.annealing[temperatureGroup] == SimulatedAnnealing::No) - { - continue; - } - GMX_RELEASE_ASSERT( - inputrec.opts.annealing[temperatureGroup] == SimulatedAnnealing::Single - || inputrec.opts.annealing[temperatureGroup] == SimulatedAnnealing::Periodic, - gmx::formatString("Unknown simulated annealing algorithm for temperature group %d", temperatureGroup) - .c_str()); - real thist = 0; - const auto npoints = inputrec.opts.anneal_npoints[temperatureGroup]; - if (inputrec.opts.annealing[temperatureGroup] == SimulatedAnnealing::Periodic) - { - /* calculate time modulo the period */ - const auto pert = inputrec.opts.anneal_time[temperatureGroup][npoints - 1]; - const auto n = static_cast(t / pert); - thist = t - n * pert; /* modulo time */ - /* Make sure rounding didn't get us outside the interval */ - if (std::fabs(thist - pert) < GMX_REAL_EPS * 100) - { - thist = 0; - } - } - else if (inputrec.opts.annealing[temperatureGroup] == SimulatedAnnealing::Single) - { - thist = t; - } - /* We are doing annealing for this group if we got here, - * and we have the (relative) time as thist. - * calculate target temp */ - int j = 0; - while ((j < npoints - 1) && (thist > (inputrec.opts.anneal_time[temperatureGroup][j + 1]))) + // No change of temperature, return current reference temperature + return inputrec.opts.ref_t[temperatureGroup]; + } + GMX_RELEASE_ASSERT( + inputrec.opts.annealing[temperatureGroup] == SimulatedAnnealing::Single + || inputrec.opts.annealing[temperatureGroup] == SimulatedAnnealing::Periodic, + gmx::formatString("Unknown simulated annealing algorithm for temperature group %d", temperatureGroup) + .c_str()); + real thist = 0; + const auto npoints = inputrec.opts.anneal_npoints[temperatureGroup]; + if (inputrec.opts.annealing[temperatureGroup] == SimulatedAnnealing::Periodic) + { + /* calculate time modulo the period */ + const auto pert = inputrec.opts.anneal_time[temperatureGroup][npoints - 1]; + const auto n = static_cast(time / pert); + thist = time - n * pert; /* modulo time */ + /* Make sure rounding didn't get us outside the interval */ + if (std::fabs(thist - pert) < GMX_REAL_EPS * 100) { - j++; + thist = 0; } - if (j < npoints - 1) + } + else if (inputrec.opts.annealing[temperatureGroup] == SimulatedAnnealing::Single) + { + thist = time; + } + /* We are doing annealing for this group if we got here, + * and we have the (relative) time as thist. + * calculate target temp */ + int j = 0; + while ((j < npoints - 1) && (thist > (inputrec.opts.anneal_time[temperatureGroup][j + 1]))) + { + j++; + } + if (j < npoints - 1) + { + /* Found our position between points j and j+1. + * Interpolate: x is the amount from j+1, (1-x) from point j + * First treat possible jumps in temperature as a special case. + */ + if ((inputrec.opts.anneal_time[temperatureGroup][j + 1] + - inputrec.opts.anneal_time[temperatureGroup][j]) + < GMX_REAL_EPS * 100) { - /* Found our position between points j and j+1. - * Interpolate: x is the amount from j+1, (1-x) from point j - * First treat possible jumps in temperature as a special case. - */ - if ((inputrec.opts.anneal_time[temperatureGroup][j + 1] - - inputrec.opts.anneal_time[temperatureGroup][j]) - < GMX_REAL_EPS * 100) - { - ir->opts.ref_t[temperatureGroup] = inputrec.opts.anneal_temp[temperatureGroup][j + 1]; - } - else - { - const real x = ((thist - inputrec.opts.anneal_time[temperatureGroup][j]) - / (inputrec.opts.anneal_time[temperatureGroup][j + 1] - - inputrec.opts.anneal_time[temperatureGroup][j])); - ir->opts.ref_t[temperatureGroup] = - x * inputrec.opts.anneal_temp[temperatureGroup][j + 1] - + (1 - x) * inputrec.opts.anneal_temp[temperatureGroup][j]; - } + return inputrec.opts.anneal_temp[temperatureGroup][j + 1]; } else { - ir->opts.ref_t[temperatureGroup] = inputrec.opts.anneal_temp[temperatureGroup][npoints - 1]; + const real x = ((thist - inputrec.opts.anneal_time[temperatureGroup][j]) + / (inputrec.opts.anneal_time[temperatureGroup][j + 1] + - inputrec.opts.anneal_time[temperatureGroup][j])); + return x * inputrec.opts.anneal_temp[temperatureGroup][j + 1] + + (1 - x) * inputrec.opts.anneal_temp[temperatureGroup][j]; } } + else + { + return inputrec.opts.anneal_temp[temperatureGroup][npoints - 1]; + } +} + +/* set target temperatures if we are annealing */ +void update_annealing_target_temp(t_inputrec* ir, real t, gmx::Update* upd) +{ + for (int temperatureGroup = 0; temperatureGroup < ir->opts.ngtc; temperatureGroup++) + { + ir->opts.ref_t[temperatureGroup] = computeAnnealingTargetTemperature(*ir, temperatureGroup, t); + } upd->update_temperature_constants(*ir); } diff --git a/src/gromacs/mdlib/coupling.h b/src/gromacs/mdlib/coupling.h index 15211f5a83..4e8857c3b0 100644 --- a/src/gromacs/mdlib/coupling.h +++ b/src/gromacs/mdlib/coupling.h @@ -175,6 +175,16 @@ void rescale_velocities(const gmx_ekindata_t* ekind, gmx::ArrayRef v); /* Rescale the velocities with the scaling factor in ekind */ +/*! + * \brief Compute the new annealing temperature for a temperature group + * + * \param inputrec The input record + * \param temperatureGroup The temperature group + * \param time The current time + * \return The new reference temperature for the group + */ +real computeAnnealingTargetTemperature(const t_inputrec& inputrec, int temperatureGroup, real time); + //! Check whether we do simulated annealing. bool doSimulatedAnnealing(const t_inputrec* ir); -- 2.22.0