From 3f970c00ffaf60b47c140003e72f941b429d275b Mon Sep 17 00:00:00 2001 From: Mark Abraham Date: Thu, 20 Jul 2017 09:16:37 +0200 Subject: [PATCH] Fix value of Ewald shift In all the short-ranged kernel flavours, sh_ewald is subtracted from rinv, which have inconsistent dimensions. Fortunately, rcutoff is often close to 1, and the inconsistent shifts often cancel in practice, and energy differences computed on neighbour lists of the same size will have the error cancel. The difference doesn't even show up in the regressiontests, but would if we had a unit test of a single interaction. Fixes #2215 Change-Id: Ia2ea57f3bd9d521879783b207353d9d6f4ccb4a8 --- src/gromacs/ewald/pme-load-balancing.cpp | 3 ++- src/gromacs/mdlib/forcerec.cpp | 8 +++++--- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/gromacs/ewald/pme-load-balancing.cpp b/src/gromacs/ewald/pme-load-balancing.cpp index c1200b2790..b7d4565788 100644 --- a/src/gromacs/ewald/pme-load-balancing.cpp +++ b/src/gromacs/ewald/pme-load-balancing.cpp @@ -820,7 +820,8 @@ pme_load_balance(pme_load_balancing_t *pme_lb, /* TODO: centralize the code that sets the potentials shifts */ if (ic->coulomb_modifier == eintmodPOTSHIFT) { - ic->sh_ewald = gmx_erfc(ic->ewaldcoeff_q*ic->rcoulomb); + GMX_RELEASE_ASSERT(ic->rcoulomb != 0, "Cutoff radius cannot be zero"); + ic->sh_ewald = gmx_erfc(ic->ewaldcoeff_q*ic->rcoulomb) / ic->rcoulomb; } if (EVDW_PME(ic->vdwtype)) { diff --git a/src/gromacs/mdlib/forcerec.cpp b/src/gromacs/mdlib/forcerec.cpp index e63381eaf1..dcc0c73cc3 100644 --- a/src/gromacs/mdlib/forcerec.cpp +++ b/src/gromacs/mdlib/forcerec.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014,2015,2016, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015,2016,2017, 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. @@ -83,6 +83,7 @@ #include "gromacs/topology/mtop_util.h" #include "gromacs/utility/exceptions.h" #include "gromacs/utility/fatalerror.h" +#include "gromacs/utility/gmxassert.h" #include "gromacs/utility/smalloc.h" #include "gromacs/utility/stringutil.h" @@ -2095,9 +2096,10 @@ init_interaction_const(FILE *fp, ic->epsfac = fr->epsfac; ic->ewaldcoeff_q = fr->ewaldcoeff_q; - if (fr->coulomb_modifier == eintmodPOTSHIFT) + if (EEL_PME_EWALD(ic->eeltype) && ic->coulomb_modifier == eintmodPOTSHIFT) { - ic->sh_ewald = gmx_erfc(ic->ewaldcoeff_q*ic->rcoulomb); + GMX_RELEASE_ASSERT(ic->rcoulomb != 0, "Cutoff radius cannot be zero"); + ic->sh_ewald = gmx_erfc(ic->ewaldcoeff_q*ic->rcoulomb) / ic->rcoulomb; } else { -- 2.11.4.GIT