# This file is part of the GROMACS molecular simulation package.
#
# Copyright (c) 2014,2015,2016,2017,2018 by the GROMACS development team.
-# Copyright (c) 2019,2020, by the GROMACS development team, led by
+# Copyright (c) 2019,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.
release-notes/2020/2020.3.rst
release-notes/2020/2020.4.rst
release-notes/2020/2020.5.rst
+ release-notes/2020/2020.6.rst
release-notes/2020/major/highlights.rst
release-notes/2020/major/features.rst
release-notes/2020/major/performance.rst
# This file is part of the GROMACS molecular simulation package.
#
# Copyright (c) 2015,2016,2017,2018,2019 by the GROMACS development team.
-# 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.
.. |Gromacs| replace:: GROMACS
.. _gmx-manual: manual-{gmx_version_string}.pdf
.. _gmx-manual-parent-dir: ../manual-{gmx_version_string}.pdf
-.. |gmx-source-package-ftp| replace:: As ftp ftp://ftp.gromacs.org/pub/gromacs/gromacs-{gmx_version_string}.tar.gz
-.. |gmx-source-package-http| replace:: As http http://ftp.gromacs.org/pub/gromacs/gromacs-{gmx_version_string}.tar.gz
-.. |gmx-regressiontests-package| replace:: http://ftp.gromacs.org/pub/regressiontests/regressiontests-{regressiontest_version}.tar.gz
-.. _up-to-date installation instructions: http://manual.gromacs.org/documentation/current/install-guide/index.html
+.. |gmx-source-package-ftp| replace:: As ftp ftp://ftp.gromacs.org/gromacs/gromacs-{gmx_version_string}.tar.gz
+.. |gmx-source-package-http| replace:: As https https://ftp.gromacs.org/gromacs/gromacs-{gmx_version_string}.tar.gz
+.. |gmx-regressiontests-package| replace:: https://ftp.gromacs.org/regressiontests/regressiontests-{regressiontest_version}.tar.gz
+.. _up-to-date installation instructions: https://manual.gromacs.org/documentation/current/install-guide/index.html
.. _CUDA: http://www.nvidia.com/object/cuda_home_new.html
.. _OpenCL: https://www.khronos.org/opencl/
.. _OpenMPI: http://www.open-mpi.org
.. _VMD: http://www.ks.uiuc.edu/Research/vmd/
.. _PyMOL: http://www.pymol.org
.. _webpage: http://www.gromacs.org
-.. _ftp site: ftp://ftp.gromacs.org/pub/gromacs/
-.. _tutorials: http://www.gromacs.org/Documentation/Tutorials
+.. _ftp site: ftp://ftp.gromacs.org/gromacs/
+.. _tutorials: http://www.mdtutorials.com/gmx/
.. _issue tracker: https://gitlab.com/gromacs/gromacs/-/issues/
.. _gitlab: https://gitlab.com/gromacs/gromacs/
.. _download: ../download.html
GROMACS 2020.5 release notes
----------------------------
-This version was released on TODO, 2020. These release notes
+This version was released on January 6th, 2021. These release notes
document the changes that have taken place in GROMACS since the
previous 2020.4 version, to fix known issues. It also incorporates all
fixes made in version 2019.6 and earlier, which you can find described
:issue:`3750`
+Fix incorrect AWH free-energies when multiple walkers share a bias
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+The AWH free-energy output was incorrect when multiple walkers shared
+an AWH bias. The error went up quadratically with the free-energy update
+interval, as well as with the number of walkers. The error decreases as
+update size decreases with time. This meant that with default AWH settings
+the error was negligible. With a free-energy update interval of 2 ps,
+we observed an error about equal to the statistical error with 32 walkers
+for a rather fast reaction coordinate. For slower coordinates the error
+will be smaller than the statistical error.
+
+:issue:`3828`
+
Fixed conserved energy for MTTK
"""""""""""""""""""""""""""""""
:issue:`3796`
+Fixed conserved energy for Nose-Hoover
+""""""""""""""""""""""""""""""""""""""
+
+When using `tcoupl=nose-hoover` and one or more temperature groups with
+non-integer number of degrees of freedom, the calculated conserved
+energy was incorrect due to an error dating back to GROMACS 2018.
+Reported conserved energies using Nose-Hoover temperature coupling and
+non-integer number of degrees of freedom since GROMACS 2018 are likely to
+be slightly off. Note that this error does not impact the dynamics, as the
+conserved energy is only reported, but never used in calculations. Also note
+that this will only be noticeable when using small temperature groups or
+small systems.
+
+:issue:`3831`
+
+Fixed kinetic energy and temperature reporting for MTTK
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+When using `pcoupl=MTTK` and `tcoupl=nose-hoover`, the reported kinetic
+energy and temperature were very slightly off. The integration of the
+temperature coupling trailed the reporting by half a time step. Note that
+these errors did not impact the dynamics, as the quantities were correctly
+integrated and only wrongly reported. Also note that the difference is so
+small that it is unlikely to have been significant for any application
+except for rigorous algorithm validation. Finally, note that this bug
+only affects this exact combination of temperature / pressure coupling
+algorithms.
+
+:issue:`3832`
+
+Fix pull error message with angles and dihedrals
+""""""""""""""""""""""""""""""""""""""""""""""""
+
+The COM pull code could print incorrect pull group indices when mdrun exited
+with an error about a too long pull distance in angle and dihedral geometries.
+
+:issue:`3613`
+
+Fix numerical issues in expanded ensemble
+"""""""""""""""""""""""""""""""""""""""""
+
+When performing simulated tempering or expanded ensemble simulations
+with changes in the Hamiltonian that were too large, then Monte Carlo
+proposals to states that were sufficiently unlikely would underflow,
+causing division by zero errors. This was fixed by numerically
+hardening the logical flow so that such proposals would be rejected
+instead.
+
+:issue:`3304`
+
+Fix incorrect electric field strength with applied electric field
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+The electric field generated by the electric field module would be incorrect when
+used together with domain decomposition due to an error with indexing the field
+to all atoms instead of just those on the current domain.
+
+In overlap regions between domains, which have the thickness of the pairlist
+cut-off distance, the electric field would be doubled (or more with 2D or
+3D domain decomposition).
+
+To validate if a simulation has been affected by the issue, users should calculate
+the actual potential across the simulation box using the Poisson equation.
+If this potential agrees with the one provided as the input, a simulation was not affected.
+
+:issue:`3800`
Fixes for ``gmx`` tools
^^^^^^^^^^^^^^^^^^^^^^^
:issue:`3820`
+Fix pull group index handling
+"""""""""""""""""""""""""""""
+
+The pull code would not validate its index groups correctly, leading
+to infinite loops or assertions being triggered at grompp time.
+
+:issue:`3810`
+
Fixes that affect portability
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
--- /dev/null
+GROMACS 2020.6 release notes
+----------------------------
+
+This version was released on TODO, 2021. These release notes
+document the changes that have taken place in GROMACS since the
+previous 2020.5 version, to fix known issues.
+
+.. Note to developers!
+ Please use """"""" to underline the individual entries for fixed issues in the subfolders,
+ otherwise the formatting on the webpage is messed up.
+ Also, please use the syntax :issue:`number` to reference issues on redmine, without the
+ a space between the colon and number!
+
+Fixes where mdrun could behave incorrectly
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Fixes for ``gmx`` tools
+^^^^^^^^^^^^^^^^^^^^^^^
+
+Fixes that affect portability
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Miscellaneous
+^^^^^^^^^^^^^
.. toctree::
:maxdepth: 1
+ 2020/2020.6
2020/2020.5
2020/2020.4
2020/2020.3
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2015,2016,2017,2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2015,2016,2017,2018,2019, The GROMACS development team.
+ * 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.
/* Need to temporarily exponentiate the log weights to sum over simulations */
for (size_t i = 0; i < buffer.size(); i++)
{
- buffer[i] = pointState[i].inTargetRegion() ? std::exp(-pointState[i].logPmfSum()) : 0;
+ buffer[i] = pointState[i].inTargetRegion() ? std::exp(pointState[i].logPmfSum()) : 0;
}
sumOverSimulations(gmx::ArrayRef<double>(buffer), commRecord, multiSimComm);
{
if (pointState[i].inTargetRegion())
{
- pointState[i].setLogPmfSum(-std::log(buffer[i] * normFac));
+ pointState[i].setLogPmfSum(std::log(buffer[i] * normFac));
}
}
}
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2015,2016,2017,2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2015,2016,2017,2018,2019, The GROMACS development team.
+ * 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.
if (fieldStrength != 0)
{
// TODO: Check parallellism
- for (index i = 0; i != ssize(f); ++i)
+ for (int i = 0; i < mdatoms.homenr; ++i)
{
// NOTE: Not correct with perturbed charges
f[i][m] += mdatoms.chargeA[i] * fieldStrength;
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
* Copyright (c) 2013,2014,2015,2016,2017, The GROMACS development team.
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,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.
#include "gromacs/topology/mtop_util.h"
#include "gromacs/topology/symtab.h"
#include "gromacs/topology/topology.h"
+#include "gromacs/utility/arrayref.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/exceptions.h"
#include "gromacs/utility/fatalerror.h"
s);
}
+static void atomGroupRangeValidation(int natoms, int groupIndex, const t_blocka& block)
+{
+ /* Now go over the atoms in the group */
+ for (int j = block.index[groupIndex]; (j < block.index[groupIndex + 1]); j++)
+ {
+ int aj = block.a[j];
+
+ /* Range checking */
+ if ((aj < 0) || (aj >= natoms))
+ {
+ gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj + 1);
+ }
+ }
+}
+
static void do_numbering(int natoms,
SimulationGroups* groups,
gmx::ArrayRef<std::string> groupsFromMdpFile,
{
unsigned short* cbuf;
AtomGroupIndices* grps = &(groups->groups[gtype]);
- int j, gid, aj, ognr, ntot = 0;
+ int ntot = 0;
const char* title;
char warn_buf[STRLEN];
for (int i = 0; i != groupsFromMdpFile.ssize(); ++i)
{
/* Lookup the group name in the block structure */
- gid = search_string(groupsFromMdpFile[i].c_str(), block->nr, gnames);
+ const int gid = search_string(groupsFromMdpFile[i].c_str(), block->nr, gnames);
if ((grptp != egrptpONE) || (i == 0))
{
grps->emplace_back(gid);
}
-
+ GMX_ASSERT(block, "Can't have a nullptr block");
+ atomGroupRangeValidation(natoms, gid, *block);
/* Now go over the atoms in the group */
- for (j = block->index[gid]; (j < block->index[gid + 1]); j++)
+ for (int j = block->index[gid]; (j < block->index[gid + 1]); j++)
{
-
- aj = block->a[j];
-
- /* Range checking */
- if ((aj < 0) || (aj >= natoms))
- {
- gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj + 1);
- }
+ const int aj = block->a[j];
/* Lookup up the old group number */
- ognr = cbuf[aj];
+ const int ognr = cbuf[aj];
if (ognr != NOGID)
{
gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)", aj + 1, title,
warning_note(wi, warn_buf);
}
/* Assign all atoms currently unassigned to a rest group */
- for (j = 0; (j < natoms); j++)
+ for (int j = 0; (j < natoms); j++)
{
if (cbuf[j] == NOGID)
{
grps->emplace_back(restnm);
/* Assign the rest name to all atoms not currently assigned to a group */
- for (j = 0; (j < natoms); j++)
+ for (int j = 0; (j < natoms); j++)
{
if (cbuf[j] == NOGID)
{
if (ir->bPull)
{
+ for (int i = 1; i < ir->pull->ngroup; i++)
+ {
+ const int gid = search_string(inputrecStrings->pullGroupNames[i].c_str(),
+ defaultIndexGroups->nr, gnames);
+ GMX_ASSERT(defaultIndexGroups, "Must have initialized default index groups");
+ atomGroupRangeValidation(natoms, gid, *defaultIndexGroups);
+ }
+
process_pull_groups(ir->pull->group, inputrecStrings->pullGroupNames, defaultIndexGroups, gnames);
checkPullCoords(ir->pull->group, ir->pull->coord);
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
* Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,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.
const double* ivxi = &state->nosehoover_vxi[i * nh];
const double* iQinv = &(MassQ->Qinv[i * nh]);
- int nd = static_cast<int>(ir->opts.nrdf[i]);
+ real nd = ir->opts.nrdf[i];
real reft = std::max<real>(ir->opts.ref_t[i], 0);
real kT = BOLTZ * reft;
{
energy += 0.5 * gmx::square(ivxi[j]) / iQinv[j];
/* contribution from the thermal variable of the NH chain */
- int ndj;
+ real ndj = 0;
if (j == 0)
{
ndj = nd;
* This file is part of the GROMACS molecular simulation package.
*
* Copyright (c) 2012-2018, The GROMACS development team.
- * Copyright (c) 2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2019,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.
#include "gromacs/utility/gmxmpi.h"
#include "gromacs/utility/smalloc.h"
+#include "expanded_internal.h"
+
static void init_df_history_weights(df_history_t* dfhist, const t_expanded* expand, int nlim)
{
int i;
int64_t step)
{
gmx_bool bSufficientSamples;
+ real acceptanceWeight;
int i;
- int n0, np1, nm1, nval, min_nvalm, min_nvalp, maxc;
- real omega_m1_0, omega_p1_0, clam_osum;
- real de, de_function;
- real cnval, zero_sum_weights;
+ int min_nvalm, min_nvalp, maxc;
+ real omega_m1_0, omega_p1_0;
+ real zero_sum_weights;
real *omegam_array, *weightsm_array, *omegap_array, *weightsp_array, *varm_array, *varp_array,
*dwp_array, *dwm_array;
- real clam_varm, clam_varp, clam_weightsm, clam_weightsp, clam_minvar;
+ real clam_varm, clam_varp, clam_osum, clam_weightsm, clam_weightsp, clam_minvar;
real * lam_variance, *lam_dg;
double* p_k;
double pks = 0;
- real chi_m1_0, chi_p1_0, chi_m2_0, chi_p2_0, chi_p1_m1, chi_p2_m1, chi_m1_p1, chi_m2_p1;
- /* if we have equilibrated the weights, exit now */
+ /* Future potential todos for this function (see #3848):
+ * - Update the names in the dhist structure to be clearer. Not done for now since this
+ * a bugfix update and we are mininizing other code changes.
+ * - Modularize the code some more.
+ * - potentially merge with accelerated weight histogram functionality, since it's very similar.
+ */
+ /* if we have equilibrated the expanded ensemble weights, we are not updating them, so exit now */
if (dfhist->bEquil)
{
return FALSE;
if (EWL(expand->elamstats))
{
- if (expand->elamstats == elamstatsWL) /* Standard Wang-Landau */
+ if (expand->elamstats == elamstatsWL) /* Using standard Wang-Landau for weight updates */
{
dfhist->sum_weights[fep_state] -= dfhist->wl_delta;
dfhist->wl_histo[fep_state] += 1.0;
}
- else if (expand->elamstats == elamstatsWWL) /* Weighted Wang-Landau */
+ else if (expand->elamstats == elamstatsWWL)
+ /* Using weighted Wang-Landau for weight updates.
+ * Very closly equivalent to accelerated weight histogram approach
+ * applied to expanded ensemble. */
{
snew(p_k, nlim);
if (expand->elamstats == elamstatsBARKER || expand->elamstats == elamstatsMETROPOLIS
|| expand->elamstats == elamstatsMINVAR)
{
-
- de_function = 0; /* to get rid of warnings, but this value will not be used because of the logic */
maxc = 2 * expand->c_range + 1;
snew(lam_dg, nlim);
snew(varm_array, maxc);
snew(dwm_array, maxc);
- /* unpack the current lambdas -- we will only update 2 of these */
+ /* unpack the values of the free energy differences and the
+ * variance in their estimates between nearby lambdas. We will
+ * only actually update 2 of these, the state we are currently
+ * at and the one we end up moving to
+ */
for (i = 0; i < nlim - 1; i++)
{ /* only through the second to last */
gmx::square(dfhist->sum_variance[i + 1]) - gmx::square(dfhist->sum_variance[i]);
}
- /* accumulate running averages */
- for (nval = 0; nval < maxc; nval++)
+ /* accumulate running averages of thermodynamic averages for Bennett Acceptance Ratio-based
+ * estimates of the free energy .
+ * Rather than peforming self-consistent estimation of the free energies at each step,
+ * we keep track of an array of possible different free energies (cnvals),
+ * and we self-consistently choose the best one. The one that leads to a free energy estimate
+ * that is closest to itself is the best estimate of the free energy. It is essentially a
+ * parallellized version of self-consistent iteration. maxc is the number of these constants. */
+
+ for (int nval = 0; nval < maxc; nval++)
{
- /* constants for later use */
- cnval = static_cast<real>(nval - expand->c_range);
- /* actually, should be able to rewrite it w/o exponential, for better numerical stability */
+ const real cnval = static_cast<real>(nval - expand->c_range);
+
+ /* Compute acceptance criterion weight to the state below this one for use in averages.
+ * Note we do not have to have just moved from that state to use this free energy
+ * estimate; these are essentially "virtual" moves. */
+
if (fep_state > 0)
{
- de = std::exp(cnval - (scaled_lamee[fep_state] - scaled_lamee[fep_state - 1]));
- if (expand->elamstats == elamstatsBARKER || expand->elamstats == elamstatsMINVAR)
- {
- de_function = 1.0 / (1.0 + de);
- }
- else if (expand->elamstats == elamstatsMETROPOLIS)
- {
- if (de < 1.0)
- {
- de_function = 1.0;
- }
- else
- {
- de_function = 1.0 / de;
- }
- }
- dfhist->accum_m[fep_state][nval] += de_function;
- dfhist->accum_m2[fep_state][nval] += de_function * de_function;
+ const auto lambdaEnergyDifference =
+ cnval - (scaled_lamee[fep_state] - scaled_lamee[fep_state - 1]);
+ acceptanceWeight =
+ gmx::calculateAcceptanceWeight(expand->elamstats, lambdaEnergyDifference);
+ dfhist->accum_m[fep_state][nval] += acceptanceWeight;
+ dfhist->accum_m2[fep_state][nval] += acceptanceWeight * acceptanceWeight;
}
+ // Compute acceptance criterion weight to transition to the next state
if (fep_state < nlim - 1)
{
- de = std::exp(-cnval + (scaled_lamee[fep_state + 1] - scaled_lamee[fep_state]));
- if (expand->elamstats == elamstatsBARKER || expand->elamstats == elamstatsMINVAR)
- {
- de_function = 1.0 / (1.0 + de);
- }
- else if (expand->elamstats == elamstatsMETROPOLIS)
- {
- if (de < 1.0)
- {
- de_function = 1.0;
- }
- else
- {
- de_function = 1.0 / de;
- }
- }
- dfhist->accum_p[fep_state][nval] += de_function;
- dfhist->accum_p2[fep_state][nval] += de_function * de_function;
+ const auto lambdaEnergyDifference =
+ -cnval + (scaled_lamee[fep_state + 1] - scaled_lamee[fep_state]);
+ acceptanceWeight =
+ gmx::calculateAcceptanceWeight(expand->elamstats, lambdaEnergyDifference);
+ dfhist->accum_p[fep_state][nval] += acceptanceWeight;
+ dfhist->accum_p2[fep_state][nval] += acceptanceWeight * acceptanceWeight;
}
- /* Metropolis transition and Barker transition (unoptimized Bennett) acceptance weight determination */
+ /* Determination of Metropolis transition and Barker transition weights */
- n0 = dfhist->n_at_lam[fep_state];
+ int numObservationsCurrentState = dfhist->n_at_lam[fep_state];
+ /* determine the number of observations above and below the current state */
+ int numObservationsLowerState = 0;
if (fep_state > 0)
{
- nm1 = dfhist->n_at_lam[fep_state - 1];
- }
- else
- {
- nm1 = 0;
+ numObservationsLowerState = dfhist->n_at_lam[fep_state - 1];
}
+ int numObservationsHigherState = 0;
if (fep_state < nlim - 1)
{
- np1 = dfhist->n_at_lam[fep_state + 1];
- }
- else
- {
- np1 = 0;
- }
-
- /* logic SHOULD keep these all set correctly whatever the logic, but apparently it can't figure it out. */
- chi_m1_0 = chi_p1_0 = chi_m2_0 = chi_p2_0 = chi_p1_m1 = chi_p2_m1 = chi_m1_p1 = chi_m2_p1 = 0;
-
- if (n0 > 0)
- {
- chi_m1_0 = dfhist->accum_m[fep_state][nval] / n0;
- chi_p1_0 = dfhist->accum_p[fep_state][nval] / n0;
- chi_m2_0 = dfhist->accum_m2[fep_state][nval] / n0;
- chi_p2_0 = dfhist->accum_p2[fep_state][nval] / n0;
+ numObservationsHigherState = dfhist->n_at_lam[fep_state + 1];
}
- if ((fep_state > 0) && (nm1 > 0))
- {
- chi_p1_m1 = dfhist->accum_p[fep_state - 1][nval] / nm1;
- chi_p2_m1 = dfhist->accum_p2[fep_state - 1][nval] / nm1;
- }
-
- if ((fep_state < nlim - 1) && (np1 > 0))
- {
- chi_m1_p1 = dfhist->accum_m[fep_state + 1][nval] / np1;
- chi_m2_p1 = dfhist->accum_m2[fep_state + 1][nval] / np1;
- }
+ /* Calculate the biases for each expanded ensemble state that minimize the total
+ * variance, as implemented in Martinez-Veracoechea and Escobedo,
+ * J. Phys. Chem. B 2008, 112, 8120-8128
+ *
+ * The variance associated with the free energy estimate between two states i and j
+ * is calculated as
+ * Var(i,j) = {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1} / numObservations(i->j)
+ * + {avg[xi(j->i)^2] / avg[xi(j->i)]^2 - 1} / numObservations(j->i)
+ * where xi(i->j) is the acceptance factor / weight associated with moving from state i to j
+ * As we are calculating the acceptance factor to the neighbors every time we're visiting
+ * a state, numObservations(i->j) == numObservations(i) and numObservations(j->i) == numObservations(j)
+ */
- omega_m1_0 = 0;
- omega_p1_0 = 0;
- clam_weightsm = 0;
- clam_weightsp = 0;
- clam_varm = 0;
- clam_varp = 0;
+ /* Accumulation of acceptance weight averages between the current state and the
+ * states +1 (p1) and -1 (m1), averaged at current state (0)
+ */
+ real avgAcceptanceCurrentToLower = 0;
+ real avgAcceptanceCurrentToHigher = 0;
+ /* Accumulation of acceptance weight averages quantities between states 0
+ * and states +1 and -1, squared
+ */
+ real avgAcceptanceCurrentToLowerSquared = 0;
+ real avgAcceptanceCurrentToHigherSquared = 0;
+ /* Accumulation of free energy quantities from lower state (m1) to current state (0) and squared */
+ real avgAcceptanceLowerToCurrent = 0;
+ real avgAcceptanceLowerToCurrentSquared = 0;
+ /* Accumulation of free energy quantities from upper state (p1) to current state (0) and squared */
+ real avgAcceptanceHigherToCurrent = 0;
+ real avgAcceptanceHigherToCurrentSquared = 0;
+
+ if (numObservationsCurrentState > 0)
+ {
+ avgAcceptanceCurrentToLower = dfhist->accum_m[fep_state][nval] / numObservationsCurrentState;
+ avgAcceptanceCurrentToHigher =
+ dfhist->accum_p[fep_state][nval] / numObservationsCurrentState;
+ avgAcceptanceCurrentToLowerSquared =
+ dfhist->accum_m2[fep_state][nval] / numObservationsCurrentState;
+ avgAcceptanceCurrentToHigherSquared =
+ dfhist->accum_p2[fep_state][nval] / numObservationsCurrentState;
+ }
+
+ if ((fep_state > 0) && (numObservationsLowerState > 0))
+ {
+ avgAcceptanceLowerToCurrent =
+ dfhist->accum_p[fep_state - 1][nval] / numObservationsLowerState;
+ avgAcceptanceLowerToCurrentSquared =
+ dfhist->accum_p2[fep_state - 1][nval] / numObservationsLowerState;
+ }
+
+ if ((fep_state < nlim - 1) && (numObservationsHigherState > 0))
+ {
+ avgAcceptanceHigherToCurrent =
+ dfhist->accum_m[fep_state + 1][nval] / numObservationsHigherState;
+ avgAcceptanceHigherToCurrentSquared =
+ dfhist->accum_m2[fep_state + 1][nval] / numObservationsHigherState;
+ }
+ /* These are accumulation of positive values (see definition of acceptance functions
+ * above), or of squares of positive values.
+ * We're taking this for granted in the following calculation, so make sure
+ * here that nothing weird happened. Although technically all values should be positive,
+ * because of floating point precisions, they might be numerically zero. */
+ GMX_RELEASE_ASSERT(
+ avgAcceptanceCurrentToLower >= 0 && avgAcceptanceCurrentToLowerSquared >= 0
+ && avgAcceptanceCurrentToHigher >= 0
+ && avgAcceptanceCurrentToHigherSquared >= 0 && avgAcceptanceLowerToCurrent >= 0
+ && avgAcceptanceLowerToCurrentSquared >= 0 && avgAcceptanceHigherToCurrent >= 0
+ && avgAcceptanceHigherToCurrentSquared >= 0,
+ "By definition, the acceptance factors should all be nonnegative.");
+
+ real varianceCurrentToLower = 0;
+ real varianceCurrentToHigher = 0;
+ real weightDifferenceToLower = 0;
+ real weightDifferenceToHigher = 0;
+ real varianceToLower = 0;
+ real varianceToHigher = 0;
if (fep_state > 0)
{
- if (n0 > 0)
+ if (numObservationsCurrentState > 0)
{
- omega_m1_0 = chi_m2_0 / (chi_m1_0 * chi_m1_0) - 1.0;
- if (nm1 > 0)
+ /* Calculate {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1}
+ *
+ * Note that if avg[xi(i->j)] == 0, also avg[xi(i->j)^2] == 0 (since the
+ * acceptances are all positive!), and hence
+ * {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1} -> 0 for avg[xi(i->j)] -> 0
+ * We're catching that case explicitly to avoid numerical
+ * problems dividing by zero when the overlap between states is small (#3304)
+ */
+ if (avgAcceptanceCurrentToLower > 0)
{
- real omega_p1_m1 = chi_p2_m1 / (chi_p1_m1 * chi_p1_m1) - 1.0;
- clam_weightsm = (std::log(chi_m1_0) - std::log(chi_p1_m1)) + cnval;
- clam_varm = (1.0 / n0) * (omega_m1_0) + (1.0 / nm1) * (omega_p1_m1);
+ varianceCurrentToLower =
+ avgAcceptanceCurrentToLowerSquared
+ / (avgAcceptanceCurrentToLower * avgAcceptanceCurrentToLower)
+ - 1.0;
+ }
+ if (numObservationsLowerState > 0)
+ {
+ /* Calculate {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1}
+ *
+ * Note that if avg[xi(i->j)] == 0, also avg[xi(i->j)^2] == 0 (since the
+ * acceptances are all positive!), and hence
+ * {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1} -> 0 for avg[xi(i->j)] -> 0
+ * We're catching that case explicitly to avoid numerical
+ * problems dividing by zero when the overlap between states is small (#3304)
+ */
+ real varianceLowerToCurrent = 0;
+ if (avgAcceptanceLowerToCurrent > 0)
+ {
+ varianceLowerToCurrent =
+ avgAcceptanceLowerToCurrentSquared
+ / (avgAcceptanceLowerToCurrent * avgAcceptanceLowerToCurrent)
+ - 1.0;
+ }
+ /* Free energy difference to the state one state lower */
+ /* if these either of these quantities are zero, the energies are */
+ /* way too large for the dynamic range. We need an alternate guesstimate */
+ if ((avgAcceptanceCurrentToLower == 0) || (avgAcceptanceLowerToCurrent == 0))
+ {
+ weightDifferenceToLower =
+ (scaled_lamee[fep_state] - scaled_lamee[fep_state - 1]);
+ }
+ else
+ {
+ weightDifferenceToLower = (std::log(avgAcceptanceCurrentToLower)
+ - std::log(avgAcceptanceLowerToCurrent))
+ + cnval;
+ }
+ /* Variance of the free energy difference to the one state lower */
+ varianceToLower =
+ (1.0 / numObservationsCurrentState) * (varianceCurrentToLower)
+ + (1.0 / numObservationsLowerState) * (varianceLowerToCurrent);
}
}
}
if (fep_state < nlim - 1)
{
- if (n0 > 0)
+ if (numObservationsCurrentState > 0)
{
- omega_p1_0 = chi_p2_0 / (chi_p1_0 * chi_p1_0) - 1.0;
- if (np1 > 0)
+ /* Calculate {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1}
+ *
+ * Note that if avg[xi(i->j)] == 0, also avg[xi(i->j)^2] == 0 (since the
+ * acceptances are all positive!), and hence
+ * {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1} -> 0 for avg[xi(i->j)] -> 0
+ * We're catching that case explicitly to avoid numerical
+ * problems dividing by zero when the overlap between states is small (#3304)
+ */
+
+ if (avgAcceptanceCurrentToHigher < 0)
+ {
+ varianceCurrentToHigher =
+ avgAcceptanceCurrentToHigherSquared
+ / (avgAcceptanceCurrentToHigher * avgAcceptanceCurrentToHigher)
+ - 1.0;
+ }
+ if (numObservationsHigherState > 0)
{
- real omega_m1_p1 = chi_m2_p1 / (chi_m1_p1 * chi_m1_p1) - 1.0;
- clam_weightsp = (std::log(chi_m1_p1) - std::log(chi_p1_0)) + cnval;
- clam_varp = (1.0 / np1) * (omega_m1_p1) + (1.0 / n0) * (omega_p1_0);
+ /* Calculate {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1}
+ *
+ * Note that if avg[xi(i->j)] == 0, also avg[xi(i->j)^2] == 0 (since the
+ * acceptances are all positive!), and hence
+ * {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1} -> 0 for avg[xi(i->j)] -> 0
+ * We're catching that case explicitly to avoid numerical
+ * problems dividing by zero when the overlap between states is small (#3304)
+ */
+ real varianceHigherToCurrent = 0;
+ if (avgAcceptanceHigherToCurrent > 0)
+ {
+ varianceHigherToCurrent =
+ avgAcceptanceHigherToCurrentSquared
+ / (avgAcceptanceHigherToCurrent * avgAcceptanceHigherToCurrent)
+ - 1.0;
+ }
+ /* Free energy difference to the state one state higher */
+ /* if these either of these quantities are zero, the energies are */
+ /* way too large for the dynamic range. We need an alternate guesstimate */
+ if ((avgAcceptanceHigherToCurrent == 0) || (avgAcceptanceCurrentToHigher == 0))
+ {
+ weightDifferenceToHigher =
+ (scaled_lamee[fep_state + 1] - scaled_lamee[fep_state]);
+ }
+ else
+ {
+ weightDifferenceToHigher = (std::log(avgAcceptanceHigherToCurrent)
+ - std::log(avgAcceptanceCurrentToHigher))
+ + cnval;
+ }
+ /* Variance of the free energy difference to the one state higher */
+ varianceToHigher =
+ (1.0 / numObservationsHigherState) * (varianceHigherToCurrent)
+ + (1.0 / numObservationsCurrentState) * (varianceCurrentToHigher);
}
}
}
- if (n0 > 0)
+ if (numObservationsCurrentState > 0)
{
- omegam_array[nval] = omega_m1_0;
+ omegam_array[nval] = varianceCurrentToLower;
}
else
{
omegam_array[nval] = 0;
}
- weightsm_array[nval] = clam_weightsm;
- varm_array[nval] = clam_varm;
- if (nm1 > 0)
+ weightsm_array[nval] = weightDifferenceToLower;
+ varm_array[nval] = varianceToLower;
+ if (numObservationsLowerState > 0)
{
- dwm_array[nval] = fabs((cnval + std::log((1.0 * n0) / nm1)) - lam_dg[fep_state - 1]);
+ dwm_array[nval] =
+ fabs((cnval + std::log((1.0 * numObservationsCurrentState) / numObservationsLowerState))
+ - lam_dg[fep_state - 1]);
}
else
{
dwm_array[nval] = std::fabs(cnval - lam_dg[fep_state - 1]);
}
- if (n0 > 0)
+ if (numObservationsCurrentState > 0)
{
- omegap_array[nval] = omega_p1_0;
+ omegap_array[nval] = varianceCurrentToHigher;
}
else
{
omegap_array[nval] = 0;
}
- weightsp_array[nval] = clam_weightsp;
- varp_array[nval] = clam_varp;
- if ((np1 > 0) && (n0 > 0))
+ weightsp_array[nval] = weightDifferenceToHigher;
+ varp_array[nval] = varianceToHigher;
+ if ((numObservationsHigherState > 0) && (numObservationsCurrentState > 0))
{
- dwp_array[nval] = fabs((cnval + std::log((1.0 * np1) / n0)) - lam_dg[fep_state]);
+ dwp_array[nval] =
+ fabs((cnval + std::log((1.0 * numObservationsHigherState) / numObservationsCurrentState))
+ - lam_dg[fep_state]);
}
else
{
}
}
- /* find the C's closest to the old weights value */
+ /* find the free energy estimate closest to the guessed weight's value */
min_nvalm = FindMinimum(dwm_array, maxc);
omega_m1_0 = omegam_array[min_nvalm];
if (expand->elamstats == elamstatsMINVAR)
{
bSufficientSamples = TRUE;
- /* make sure they are all past a threshold */
+ /* make sure the number of samples in each state are all
+ * past a user-specified threshold
+ */
for (i = 0; i < nlim; i++)
{
if (dfhist->n_at_lam[i] < expand->minvarmin)
de = weighted_lamee[lamtrial] - weighted_lamee[fep_state];
if (expand->elmcmove == elmcmoveMETROPOLIS)
{
- tprob = 1.0;
- trialprob = std::exp(de);
- if (trialprob < tprob)
+ tprob = 1.0;
+ if (de < 0)
{
- tprob = trialprob;
+ tprob = std::exp(de);
}
propose[fep_state] = 0;
propose[lamtrial] = 1.0; /* note that this overwrites the above line if fep_state = ntrial, which only occurs at the ends */
}
else if (expand->elmcmove == elmcmoveBARKER)
{
- tprob = 1.0 / (1.0 + std::exp(-de));
-
+ if (de > 0) /* Numerically stable version */
+ {
+ tprob = 1.0 / (1.0 + std::exp(-de));
+ }
+ else if (de < 0)
+ {
+ tprob = std::exp(de) / (std::exp(de) + 1.0);
+ }
propose[fep_state] = (1 - tprob);
propose[lamtrial] +=
tprob; /* we add, to account for the fact that at the end, they might be the same point */
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * 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.
+ *
+ * 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 internal functionality for expanded ensemble
+ *
+ * \author Pascal Merz <pascal.merz@me.com>
+ * \author Michael Shirts <michael.shirts@colorado.edu>
+ * \ingroup module_mdlib
+ */
+#include "gmxpre.h"
+
+#include "expanded_internal.h"
+
+#include <cmath>
+
+#include "gromacs/mdtypes/md_enums.h"
+#include "gromacs/utility/exceptions.h"
+
+namespace gmx
+{
+real calculateAcceptanceWeight(int calculationMode, real lambdaEnergyDifference)
+{
+ if (calculationMode == elamstatsBARKER || calculationMode == elamstatsMINVAR)
+ {
+ /* Barker acceptance rule forumula is used for accumulation of probability for
+ * both the Barker variant of the weight accumulation algorithm and the
+ * minimum variance variant of the weight accumulation algorithm.
+ *
+ * Barker acceptance rule for a jump from state i -> j is defined as
+ * exp(-E_i)/exp(-Ei)+exp(-Ej) = 1 / (1 + exp(dE_ij))
+ * where dE_ij is the potential energy difference between the two states
+ * plus a constant offset that can be removed at the end for numerical stability.
+ * dE_ij = FE_j - FE_i + offset
+ * Numerically, this computation can be unstable if dE gets large. (See #3304)
+ * To avoid numerical instability, we're calculating it as
+ * 1 / (1 + exp(dE_ij)) (if dE < 0)
+ * exp(-dE_ij) / (exp(-dE_ij) + 1) (if dE > 0)
+ */
+ if (lambdaEnergyDifference < 0)
+ {
+ return 1.0 / (1.0 + std::exp(lambdaEnergyDifference));
+ }
+ else
+ {
+ return std::exp(-lambdaEnergyDifference) / (1.0 + std::exp(-lambdaEnergyDifference));
+ }
+ }
+ else if (calculationMode == elamstatsMETROPOLIS)
+ {
+ /* Metropolis acceptance rule for a jump from state i -> j is defined as
+ * 1 (if dE_ij < 0)
+ * exp(-dE_ij) (if dE_ij >= 0)
+ * where dE_ij is the potential energy difference between the two states
+ * plus a free energy offset that can be subtracted off later:
+ * dE_ij = FE_j - FE_i + offset
+ */
+ if (lambdaEnergyDifference < 0)
+ {
+ return 1.0;
+ }
+ else
+ {
+ return std::exp(-lambdaEnergyDifference);
+ }
+ }
+
+ GMX_THROW(NotImplementedError("Unknown acceptance calculation mode"));
+}
+
+} // namespace gmx
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * 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.
+ *
+ * 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 internal functionality for expanded ensemble
+ *
+ * This file is only used by expanded.cpp and tests/expanded.cpp.
+ *
+ * \author Pascal Merz <pascal.merz@me.com>
+ * \author Michael Shirts <michael.shirts@colorado.edu>
+ * \ingroup module_mdlib
+ */
+#ifndef GMX_MDLIB_EXPANDEDINTERNAL_H
+#define GMX_MDLIB_EXPANDEDINTERNAL_H
+
+#include "gromacs/utility/real.h"
+
+namespace gmx
+{
+/*! \brief Calculates the acceptance weight for a lambda state transition
+ *
+ * \param[in] calculationMode How the lambda weights are calculated
+ * \param[in] lambdaEnergyDifference The difference in energy between the two states
+ * \return The acceptance weight
+ */
+real calculateAcceptanceWeight(int calculationMode, real lambdaEnergyDifference);
+} // namespace gmx
+
+#endif // GMX_MDLIB_EXPANDEDINTERNAL_H
#
# This file is part of the GROMACS molecular simulation package.
#
-# Copyright (c) 2014,2016,2017,2018,2019,2020, by the GROMACS development team, led by
+# Copyright (c) 2014,2016,2017,2018,2019,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.
ebin.cpp
energydrifttracker.cpp
energyoutput.cpp
+ expanded.cpp
freeenergyparameters.cpp
leapfrog.cpp
leapfrogtestdata.cpp
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * 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.
+ *
+ * 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 expanded ensemble
+ *
+ * This file contains unit tests for functions used by the expanded
+ * ensemble.
+ *
+ * \todo Add more tests as the expanded ensemble implementation
+ * gets more modular (#3848).
+ *
+ * \author Pascal Merz <pascal.merz@me.com>
+ * \author Michael Shirts <michael.shirts@colorado.edu>
+ * \ingroup module_mdlib
+ */
+#include "gmxpre.h"
+
+#include <cmath>
+
+#include <gtest/gtest.h>
+
+#include "gromacs/mdlib/expanded_internal.h"
+#include "gromacs/mdtypes/md_enums.h"
+
+#include "testutils/testasserts.h"
+
+namespace gmx
+{
+namespace test
+{
+namespace
+{
+
+//! Test fixture accepting a value to pass into calculateAcceptanceWeight
+class CalculateAcceptanceWeightSimple : public ::testing::Test, public ::testing::WithParamInterface<real>
+{
+};
+// Check that unimplemented calculation modes throw
+TEST_P(CalculateAcceptanceWeightSimple, UnknownCalculationModeThrows)
+{
+ for (auto calculationMode = 0; calculationMode < elamstatsNR; ++calculationMode)
+ {
+ if (calculationMode != elamstatsBARKER && calculationMode != elamstatsMINVAR
+ && calculationMode != elamstatsMETROPOLIS)
+ {
+ EXPECT_THROW_GMX(calculateAcceptanceWeight(calculationMode, GetParam()), NotImplementedError);
+ }
+ }
+}
+// Check that implemented calculation modes don't throw
+TEST_P(CalculateAcceptanceWeightSimple, KnownCalculationModeDoesNotThrow)
+{
+ EXPECT_NO_THROW(calculateAcceptanceWeight(elamstatsMETROPOLIS, GetParam()));
+ EXPECT_NO_THROW(calculateAcceptanceWeight(elamstatsBARKER, GetParam()));
+ EXPECT_NO_THROW(calculateAcceptanceWeight(elamstatsMINVAR, GetParam()));
+}
+// Barker and MinVar are expected to be equal
+TEST_P(CalculateAcceptanceWeightSimple, BarkerAndMinVarAreIdentical)
+{
+ EXPECT_EQ(calculateAcceptanceWeight(elamstatsBARKER, GetParam()),
+ calculateAcceptanceWeight(elamstatsMINVAR, GetParam()));
+}
+
+/*! \brief Test fixture accepting a calculation mode and an input value for
+ * calculateAcceptanceWeight as well as the expected output value
+ */
+using RegressionTuple = std::tuple<int, real, real>;
+class CalculateAcceptanceWeightRangeRegression :
+ public ::testing::Test,
+ public ::testing::WithParamInterface<RegressionTuple>
+{
+};
+// Check that output is as expected
+TEST_P(CalculateAcceptanceWeightRangeRegression, ValuesMatch)
+{
+ const auto calculationMode = std::get<0>(GetParam());
+ const auto inputValue = std::get<1>(GetParam());
+ const auto expectedOutput = std::get<2>(GetParam());
+
+ EXPECT_REAL_EQ(expectedOutput, calculateAcceptanceWeight(calculationMode, inputValue));
+}
+
+INSTANTIATE_TEST_CASE_P(
+ SimpleTests,
+ CalculateAcceptanceWeightSimple,
+ ::testing::Values(1., -1., 0., GMX_REAL_NEGZERO, GMX_REAL_EPS, -GMX_REAL_EPS, GMX_REAL_MAX, -GMX_REAL_MAX));
+INSTANTIATE_TEST_CASE_P(
+ RegressionTests,
+ CalculateAcceptanceWeightRangeRegression,
+ ::testing::Values(RegressionTuple{ elamstatsMETROPOLIS, 0.0, 1.0 },
+ RegressionTuple{ elamstatsMETROPOLIS, GMX_REAL_NEGZERO, 1.0 },
+ RegressionTuple{ elamstatsMETROPOLIS, GMX_REAL_EPS, 1.0 },
+ RegressionTuple{ elamstatsMETROPOLIS, -1.0, 1.0 },
+ RegressionTuple{ elamstatsMETROPOLIS, -GMX_REAL_MAX, 1.0 },
+ RegressionTuple{ elamstatsMETROPOLIS, 1.0, std::exp(-1.0) },
+ RegressionTuple{ elamstatsMETROPOLIS, GMX_REAL_MAX, 0.0 },
+ RegressionTuple{ elamstatsBARKER, 0.0, 0.5 },
+ RegressionTuple{ elamstatsBARKER, GMX_REAL_NEGZERO, 0.5 },
+ RegressionTuple{ elamstatsBARKER, GMX_REAL_EPS, 0.5 },
+ RegressionTuple{ elamstatsBARKER, -1.0, 1.0 / (1.0 + std::exp(-1.0)) },
+ RegressionTuple{ elamstatsBARKER, -GMX_REAL_MAX, 1.0 },
+ RegressionTuple{ elamstatsBARKER, 1.0, 1.0 / (1.0 + std::exp(1.0)) },
+ RegressionTuple{ elamstatsBARKER, GMX_REAL_MAX, 0.0 }));
+
+} // namespace
+} // namespace test
+} // namespace gmx
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2011-2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2011-2019,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.
copy_mat(shake_vir, state->svir_prev);
copy_mat(force_vir, state->fvir_prev);
}
- if (inputrecNvtTrotter(ir) && ir->eI == eiVV)
+ if ((inputrecNptTrotter(ir) || inputrecNvtTrotter(ir)) && ir->eI == eiVV)
{
/* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
enerd->term[F_TEMP] =
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
* Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,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.
/* This function returns the distance based on coordinates xg and xref.
* Note that the pull coordinate struct pcrd is not modified.
+ *
+ * \param[in] pull The pull struct
+ * \param[in] pcrd The pull coordinate to compute a distance for
+ * \param[in] pbc The periodic boundary conditions
+ * \param[in] xg The coordinate of group 1
+ * \param[in] xref The coordinate of group 0
+ * \param[in] groupIndex0 The index of group 0 in the pcrd->params.group
+ * \param[in] groupIndex1 The index of group 1 in the pcrd->params.group
+ * \param[in] max_dist2 The maximum distance squared
+ * \param[out] dr The distance vector
*/
static void low_get_pull_coord_dr(const struct pull_t* pull,
const pull_coord_work_t* pcrd,
const t_pbc* pbc,
- dvec xg,
+ const dvec xg,
dvec xref,
- double max_dist2,
+ const int groupIndex0,
+ const int groupIndex1,
+ const double max_dist2,
dvec dr)
{
const pull_group_work_t* pgrp0 = &pull->group[pcrd->params.group[0]];
gmx_fatal(FARGS,
"Distance between pull groups %d and %d (%f nm) is larger than 0.49 times the "
"box size (%f).\n%s",
- pcrd->params.group[0], pcrd->params.group[1], sqrt(dr2),
+ pcrd->params.group[groupIndex0], pcrd->params.group[groupIndex1], sqrt(dr2),
sqrt(0.98 * 0.98 * max_dist2), pcrd->params.eGeom == epullgDIR ? "You might want to consider using \"pull-geometry = direction-periodic\" instead.\n" : "");
}
pull_group_work_t* pgrp1 = &pull->group[pcrd->params.group[1]];
low_get_pull_coord_dr(pull, pcrd, pbc, pgrp1->x,
- pcrd->params.eGeom == epullgCYL ? pull->dyna[coord_ind].x : pgrp0->x, md2,
- spatialData.dr01);
+ pcrd->params.eGeom == epullgCYL ? pull->dyna[coord_ind].x : pgrp0->x, 0,
+ 1, md2, spatialData.dr01);
if (pcrd->params.ngroup >= 4)
{
pgrp2 = &pull->group[pcrd->params.group[2]];
pgrp3 = &pull->group[pcrd->params.group[3]];
- low_get_pull_coord_dr(pull, pcrd, pbc, pgrp3->x, pgrp2->x, md2, spatialData.dr23);
+ low_get_pull_coord_dr(pull, pcrd, pbc, pgrp3->x, pgrp2->x, 2, 3, md2, spatialData.dr23);
}
if (pcrd->params.ngroup >= 6)
{
pgrp4 = &pull->group[pcrd->params.group[4]];
pgrp5 = &pull->group[pcrd->params.group[5]];
- low_get_pull_coord_dr(pull, pcrd, pbc, pgrp5->x, pgrp4->x, md2, spatialData.dr45);
+ low_get_pull_coord_dr(pull, pcrd, pbc, pgrp5->x, pgrp4->x, 4, 5, md2, spatialData.dr45);
}
}
/* Get the current difference vector */
low_get_pull_coord_dr(pull, pcrd, pbc, rnew[pcrd->params.group[1]],
- rnew[pcrd->params.group[0]], -1, unc_ij);
+ rnew[pcrd->params.group[0]], 0, 1, -1, unc_ij);
if (debug)
{
g0 = pcrd->params.group[0];
g1 = pcrd->params.group[1];
- low_get_pull_coord_dr(pull, pcrd, pbc, rnew[g1], rnew[g0], -1, tmp);
- low_get_pull_coord_dr(pull, pcrd, pbc, dr1, dr0, -1, tmp3);
+ low_get_pull_coord_dr(pull, pcrd, pbc, rnew[g1], rnew[g0], 0, 1, -1, tmp);
+ low_get_pull_coord_dr(pull, pcrd, pbc, dr1, dr0, 0, 1, -1, tmp3);
fprintf(debug, "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n", rnew[g0][0],
rnew[g0][1], rnew[g0][2], rnew[g1][0], rnew[g1][1], rnew[g1][2], dnorm(tmp));
fprintf(debug, "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f\n", "", "", "", "", "",
}
low_get_pull_coord_dr(pull, &coord, pbc, rnew[coord.params.group[1]],
- rnew[coord.params.group[0]], -1, unc_ij);
+ rnew[coord.params.group[0]], 0, 1, -1, unc_ij);
switch (coord.params.eGeom)
{
# This file is part of the GROMACS molecular simulation package.
#
# Copyright (c) 2012,2013,2014,2015,2016, The GROMACS development team.
-# Copyright (c) 2017,2018,2019,2020, by the GROMACS development team, led by
+# Copyright (c) 2017,2018,2019,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(REGRESSIONTEST_URL https://gitlab.com/gromacs/gromacs-regressiontests/-/archive/${REGRESSIONTEST_BRANCH}/gromacs-regressiontests-${REGRESSIONTEST_BRANCH}.tar.gz)
# REGRESSIONTEST_PATH for dev trees is set later based on the dirname found in the tar
else()
- set(REGRESSIONTEST_URL http://ftp.gromacs.org/regressiontests/regressiontests-${REGRESSIONTEST_VERSION}.tar.gz)
+ set(REGRESSIONTEST_URL https://ftp.gromacs.org/regressiontests/regressiontests-${REGRESSIONTEST_VERSION}.tar.gz)
set(REGRESSIONTEST_PATH
"${CMAKE_CURRENT_BINARY_DIR}/regressiontests-${REGRESSIONTEST_VERSION}"
CACHE PATH "Path to auto-downloaded regressiontests" FORCE)
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff