Project

General

Profile

Bug #1463

Potential modifier potential-switch gives incorrect forces for perturbed pairs

Added by Berk Hess over 3 years ago. Updated about 3 years ago.

Status:
Closed
Priority:
High
Assignee:
Category:
mdrun
Target version:
Affected version - extra info:
Affected version:
Difficulty:
uncategorized
Close

Description

Potential modifier potential-switch gives for both Coulomb and VdW:
F*sw + V*sw*dsw/r
instead of:
F*sw - V*dsw
Since the second term is very small compared to the first for LJ, the relative error in the forces is only 10^-4, which made this hard to notice, if someone even ever used free-energy with this potential modifier.

tabulated.mdp - Old tabulated switched potentials (11 KB) Erik Lindahl, 04/18/2014 02:20 PM

modifiers.mdp - New analytical modifier switched potentials (11 KB) Erik Lindahl, 04/18/2014 02:20 PM

nonanol_group_original_rerun.0.mdp (5.55 KB) Michael Shirts, 04/18/2014 02:49 PM

nonanol.gro - Nonanol Gro (226 KB) Michael Shirts, 05/22/2014 07:41 PM

nonanol.top - Nonanol free energy topology (16 KB) Michael Shirts, 05/22/2014 07:41 PM

Associated revisions

Revision c97f2ec6 (diff)
Added by Roland Schulz about 3 years ago

Updated warnings for potential switch w/o defining switch regions

Related to #1463

Change-Id: I76da771b9960b770028bc75a260072588d0699af

Revision 53bfa7dd (diff)
Added by Michael Shirts about 3 years ago

Issue warnings for potential switch w/o defining switch regions.

If either older style switches or new potential-switch modifiers are
used, issue a warning if the PME regions is too long, resulting in
inaccurate energies. Also issue a warning if rvdw_switch is 0, which
occurs if no value is specified.

Discussion is in redmine #1463

Change-Id: I7a2c29f87ceb04712aab5076ac97f5f22f573671

Revision ba3dcd57 (diff)
Added by Erik Lindahl about 3 years ago

Updated reference with fixed potential-shift dispcorr

New reference energy files (single and double) that use the
proper potential-shift dispersion correction. It also fixes
a typo in the reference grompp output.

Refs. #1463.

Change-Id: If422a4326faede476ccb32be5a9da332e92f2b71

Revision dced970a (diff)
Added by Berk Hess about 3 years ago

Fixed shift and switch modifiers, particularly for free-energy

When using tabulated interactions (historically with PME-Switch), the
previous free-energy kernels used tabulated interactions which gave
correct results. However, as we have moved to using the new
interaction modifiers, Ewald short-ranged interactions are computed
analytically. To extend the range over which we apply the soft-core
interaction, the free-energy kernels evaluated interactions by
subtracting the reciprocal-space component, and then applying the
free-energy evaluation to the Coulomb (1/r) short-range
interaction. This works fine for vanilla PME, but led to problems when
combined with a switch modifier, since we are switching a different
function compared to the non-free-energy kernels. This could lead to
large artefacts where the free energy was 100x off if we were applying
the cutoff to r while the switch was applied to the scaled soft-core
radius.

This patch modifies the free-energy kernel so that the vanilla, shift,
and exact-cutoff versions still use the compensation trick, while the
switch modifier always operates on the traditional short-range Ewald
functional form.

The (very small) Ewald shift has also been added when computing free
energy in combination with Ewald summation and potential-shift
modifiers. As the perturbation goes to zero, the interaction will also
approach the non-free-energy interactions. Tested to match the
non-free-energy kernel to with 1e-8 in the fully coupled state, it
conserves energy, and produces reasonable free energies for ethanol in
water.

This also modifies table-generation, table-usage, and
dispersion-correction code to use shift/switch forms (and correctly),
when that has been selected in the interaction modifiers. This
provides much more accurate results for our new shifted interactions.

Correct (unmodified) tables are now generated for 1-4 interactions
in a few corner cases in the presence of modifiers for non-bonded
interactions.

Code paths for using exact cutoffs now work correctly when
rcoulomb-switch != rvdw-switch, or if only one kind of switch is
active.

Free-energy calculations using a plain Coulomb interaction now
incorporate a potential shift if one exists.

The GMX_NB_GENERIC environment variable can now be used to specify the
use of the generic kernel even with shifts or switches active.

Fixes #1463.

Change-Id: Ia63a1ed7d6c9cdf9cd9e6209b6326a49043060ec

History

#1 Updated by Gerrit Code Review Bot over 3 years ago

Gerrit received a related patchset '1' for Issue #1463.
Uploader: Berk Hess ()
Change-Id: Ia63a1ed7d6c9cdf9cd9e6209b6326a49043060ec
Gerrit URL: https://gerrit.gromacs.org/3275

#2 Updated by Gerrit Code Review Bot over 3 years ago

Gerrit received a related patchset '1' for Issue #1463.
Uploader: Erik Lindahl ()
Change-Id: Iced34f2ecfe5b34a0d1f5d51e572c9b3ba8b4d93
Gerrit URL: https://gerrit.gromacs.org/3290

#3 Updated by Gerrit Code Review Bot about 3 years ago

Gerrit received a related patchset '1' for Issue #1463.
Uploader: Mark Abraham ()
Change-Id: Ia63a1ed7d6c9cdf9cd9e6209b6326a49043060ec
Gerrit URL: https://gerrit.gromacs.org/3304

#4 Updated by Michael Shirts about 3 years ago

Ran some extensive calculations for: the old style group cutoffs + potential-shift + potential-switch with nonanol solvation. These were calculated by performing one set of free energy calculations with a reference potential, then rerunning with the new potential energy function. This allows differences in free energy differences to be calculated very accurately. I did need to tweak the dhdl output to use potential, not total energy.

Note: these calculations involve 20x1 ns direct calculation ONCE, and the rest are reruns. VERY fast. All this data generated on a laptop this evening.
To validate the uncertainty, I did one calculations with 3X the data, independent calculation (cut_11_long). Statistical Errors are consistent. VERY accurate method.

dG = free energy of solvation
ddG_Ref = difference of free energy change from our standard reference.
dG(lam=0) = difference of free energies at lambda=0 to the reference state at lambda=0 (the ENTIRE system, not just the molecule.
Start (lam=0)
EPot = potential energy - From the same single configuration at lam=0
Dhdlcoul = coulomb dH/dL - From the same single configuration at lam=0
dhdlVDW = vdw dH/dL - From the same single configuration at lam=0
Start (lam=1)
EPot = potential energy - From the same single configuration at lam=1
Dhdlcoul = coulomb dH/dL - From the same single configuration at lam=1
dhdlVDW = vdw dH/dL - From the same single configuration at lam=1

All numbers in kJ/mol

                                                                                                   Start (lam=0)                                                          End(lam=1)
                 dG                        ddG_Ref             dG(lam=0)             EPot   Dhdlcoul    DhdlVDW                                          EPot   Dhdlcoul    DhdlVDW

 cut_12          -4.514  +/-  0.587     -0.012  +/- 0.011    0.688  +/- 0.006    -44415.191 53.009216 -27.807858  -43861.848 362.63086 13.810789
 cut_11          -4.505  +/-  0.587     -0.021  +/- 0.012    0.502  +/- 0.009    -44415.461 53.010178 -27.808994  -43862.027 362.63229 13.820795 
 cut_10          -4.526  +/-  0.587     -0.000  +/- 0.000    0.000  +/- 0.000    -44415.766 52.974979 -27.839355  -43862.621 362.61575 13.811733
 cut_9           -4.478  +/-  0.588     -0.048  +/- 0.021   -1.699  +/- 0.017    -44417.836 52.969749 -27.807266  -43864.285 362.64990 13.779721 

 cut_11_long     -4.591  +/-  0.265     -0.010  +/- 0.005    0.513  +/- 0.005    -44415.461 53.010178 -27.808994  -43862.621 362.61575 13.811733
 (check -- should be statistically indistinguishable from cut_11)

 shift_12        -4.490  +/-  0.587     -0.036  +/- 0.010    0.188  +/- 0.008    -44415.625 53.011524 -27.813204  -43862.391 364.49475 13.784118
 shift_11        -4.465  +/-  0.587     -0.061  +/- 0.009   -0.464  +/- 0.008    -44416.363 52.998901 -27.812801  -43862.922 364.48834 13.766539 
 shift_10        -4.424  +/-  0.587     -0.102  +/- 0.008   -1.811  +/- 0.007    -44417.629 52.973351 -27.799637  -43864.363 364.48276 13.737906
 shift_9         -4.291  +/-  0.587     -0.235  +/- 0.015   -5.389  +/- 0.012    -44421.430 52.970177 -27.783443  -43867.871 364.51086 13.702047

 switch_12       -4.682  +/-  0.587      0.156  +/-  0.010   3.648  +/-  0.008   -44412.266 53.009220 -27.992271  -43859.059 364.49622 13.624015
 switch_11       -4.759  +/-  0.587      0.233  +/-  0.011   4.590  +/-  0.009   -44411.430 52.998150 -28.071299  -43858.121 364.48489 13.555407
 switch_10       -4.875  +/-  0.587      0.349  +/-  0.006   5.923  +/-  0.005   -44409.922 52.975014 -28.224947  -43856.988 364.48099 13.421916
 switch_9        -5.002  +/-  0.587      0.476  +/-  0.021   7.268  +/-  0.017   -44408.855 52.969612 -28.396847  -43855.719 364.51532 13.184030

Switch is much better than previously bad test. Still not nearly as good as the old switched cutoffs, though. Shift MAY be accurate enough to use for precision free energy calculations, but switch probably isn't. Error does appear to primarily be in the vdw part (these are done with linear charge, then vdw after, sequentially). dvhdl at lambda=0 seems to be particularly different. All methods have little radial dependence on Coulombic energy.

Parameters:
cut_10 (reference)
rlist = 1.25
cutoff-scheme = Group
coulombtype = PME-Switch
rcoulomb-switch = 0.98
rcoulomb = 1.0
vdw-type = Switch
rvdw-switch = 0.95
rvdw = 1.0
DispCorr = AllEnerPres

Cut_9
rlist = 1.25
cutoff-scheme = Group
coulombtype = PME-Switch
rcoulomb-switch = 0.88
rcoulomb = 0.9
vdw-type = Switch
rvdw-switch = 0.85
rvdw = 0.9
DispCorr = AllEnerPres

Cut_11
rlist = 1.35
cutoff-scheme = Group
coulombtype = PME-Switch
rcoulomb-switch = 1.08
rcoulomb = 1.1
vdw-type = Switch
rvdw-switch = 1.05
rvdw = 1.1
DispCorr = AllEnerPres

Cut_12
rlist = 1.45
cutoff-scheme = Group
coulombtype = PME-Switch
rcoulomb-switch = 1.18
rcoulomb = 1.2
vdw-type = Switch
rvdw-switch = 1.15
rvdw = 1.2
DispCorr = AllEnerPres

Shift_9
rlist = 1.25
cutoff-scheme = Group
coulomb_type = pme
rcoulomb = 0.9
vdw_type = cut-off
rvdw = 0.9
coulomb-modifier = potential-shift
vdw-modifier = potential-shift
dispcorr = enerpres

Shift_10
rlist = 1.25
cutoff-scheme = Group
coulomb_type = pme
rcoulomb = 1.0
vdw_type = cut-off
rvdw = 1.0
coulomb-modifier = potential-shift
vdw-modifier = potential-shift
dispcorr = enerpres

Shift_11
rlist = 1.35
cutoff-scheme = Group
coulomb_type = pme
rcoulomb = 1.1
vdw_type = cut-off
rvdw = 1.1
coulomb-modifier = potential-shift
vdw-modifier = potential-shift
dispcorr = enerpres

Shift_12
rlist = 1.45
cutoff-scheme = Group
coulomb_type = pme
rcoulomb = 1.2
vdw_type = cut-off
rvdw = 1.2
coulomb-modifier = potential-shift
vdw-modifier = potential-shift
dispcorr = enerpres

Switch_9
rlist = 1.25
cutoff-scheme = Group
coulomb_type = pme
rcoulomb = 0.9
coulomb-modifier = potential-switch
rcouloumb-switch = 0.88
vdw_type = cut-off
rvdw = 0.9
vdw-modifier = potential-switch
rvdw-switch = 0.85
dispcorr = enerpres

Switch_10
rlist = 1.25
cutoff-scheme = Group
coulomb_type = pme
rcoulomb = 1.0
coulomb-modifier = potential-switch
rcouloumb-switch = 0.98
vdw_type = cut-off
rvdw = 1.0
vdw-modifier = potential-switch
rvdw-switch = 0.95
dispcorr = enerpres

Switch_11
rlist = 1.35
cutoff-scheme = Group
coulomb_type = pme
rcoulomb = 1.1
coulomb-modifier = potential-switch
rcouloumb-switch = 1.08
vdw_type = cut-off
rvdw = 1.1
vdw-modifier = potential-switch
rvdw-switch = 1.05
dispcorr = enerpres

Switch_10
rlist = 1.45
cutoff-scheme = Group
coulomb_type = pme
rcoulomb = 1.2
coulomb-modifier = potential-switch
rcouloumb-switch = 1.18
vdw_type = cut-off
rvdw = 1.2
vdw-modifier = potential-switch
rvdw-switch = 1.15
dispcorr = enerpres

#5 Updated by Mark Abraham about 3 years ago

Great. Am assuming this was done with fdda24990 i.e. patch 9 of https://gerrit.gromacs.org/#/c/3275/.

Typo for rvdw=2.0 in cut_12?

rvdw-switch and rcoulomb-switch are parameters for the potential-switch modifiers, so we'd need to a) choose, and b) report that.

#6 Updated by Michael Shirts about 3 years ago

Mark Abraham wrote:

Great. Am assuming this was done with fdda24990 i.e. patch 9 of https://gerrit.gromacs.org/#/c/3275/.

yes:

Typo for rvdw=2.0 in cut_12?

yes:

rvdw-switch and rcoulomb-switch are parameters for the potential-switch modifiers, so we'd need to a) choose, and b) report that.

Aha, I wasn't keeping track of all the differences. Checking the mdout, they were both set to 0 as defaults. I thought we had fixed this behavior in 4.6.5 already?

Will rerun switch with proper modifications (should be done soon!).

#7 Updated by Mark Abraham about 3 years ago

Michael Shirts wrote:

Typo for rvdw=2.0 in cut_12?

yes:

The result for cut_12 looked perhaps like an outlier, so if the typo was in the original .mdp, a fix and rerun is in order?

rvdw-switch and rcoulomb-switch are parameters for the potential-switch modifiers, so we'd need to a) choose, and b) report that.

Aha, I wasn't keeping track of all the differences. Checking the mdout, they were both set to 0 as defaults. I thought we had fixed this behavior in 4.6.5 already?

f78f0dc83f55d588f0fbc049af667519d9cf868e made a warning for excessively large switch regions only for PME-Switch. Perhaps we should now extend that warning to also cover PME + modifier = potential-switch? (And later, in release-5-0, to LJPME? Probably immaterial, though.)

#8 Updated by Michael Shirts about 3 years ago

Switch information updated. Much better, though still noticebly worse that the older switch format and potential-shift.

The result for cut_12 looked perhaps like an outlier, so if the typo was in the original .mdp, a fix and rerun is in order?

No, it was a double typo in writing out the parameters in the issue tracker + in copying the data. It is consistent now below. I've also added info for energy and derivatives at lambda=1, where some trends are clearer.

f78f0dc83f55d588f0fbc049af667519d9cf868e made a warning for excessively large switch regions only for PME-Switch. Perhaps we should now extend that warning to also cover PME + modifier = potential-switch? (And later, in release-5-0, to LJPME? Probably immaterial, though.)

Anytime rX-switch and rX, where X is vdw or coul are different by too much should be flagged. rvdw-switch was getting set to 0. Someone was burned by this before -- I'll look to see where it was reported.

#9 Updated by Mark Abraham about 3 years ago

Michael Shirts wrote:

Anytime rX-switch and rX, where X is vdw or coul are different by too much should be flagged.

Right, but is it clear that warning should be given for non-PME?

#10 Updated by Michael Shirts about 3 years ago

Mark Abraham wrote:

Michael Shirts wrote:

Anytime rX-switch and rX, where X is vdw or coul are different by too much should be flagged.

Right, but is it clear that warning should be given for non-PME?

Yes, rvdw-switch is being set to zero when potential-switch is set. This is bad.

#11 Updated by Mark Abraham about 3 years ago

Michael Shirts wrote:

Mark Abraham wrote:

Michael Shirts wrote:

Anytime rX-switch and rX, where X is vdw or coul are different by too much should be flagged.

Right, but is it clear that warning should be given for non-PME?

Yes, rvdw-switch is being set to zero when potential-switch is set. This is bad.

It's not "being set to zero," it defaults to zero and is never assigned to again. Reasonable default behaviour would be to be somewhat smaller than rvdw, but that behaviour is not implemented that I can see. You are remembering that we added a warning in the commit that I linked. I do not think we did anything about it being set.

We won't change the default behaviour in release-4-6, but we could still make release-5-0 default to rX - 0.02, or something similar we think sound. That's why we're doing these kinds of testing.

Giving a warning when the switched range is too large is a different thing, and mitigates damage from the current default. The existing warning is sensible for PME, and PME plus switch modifier, but it is not clear to me that it is sensible for non-PME coulomb, which is what I was asking about. It does seem sensible to warn for rvdw=cutoff + rvdw-modifier=potential-switch (and to do so in release-4-6).

#12 Updated by Michael Shirts about 3 years ago

It's not "being set to zero," it defaults to zero and is never assigned to again. Reasonable default behaviour would be to be somewhat smaller than rvdw, but that behaviour is not implemented that I can see. You are remembering that we added a warning in the commit that I linked. I do not think we did anything about it being set.

Well, sure -- that's semantics. At creation, its value is set to zero by mdrun, and potential-switch uses that value, whether the user intended or not.

We won't change the default behaviour in release-4-6, but we could still make release-5-0 default to rX - 0.02, or something similar we think sound. That's why we're doing these kinds of testing.

Agreed. We have data to suggest defaults, though it might be good to issue a note if the defaults are used (again, that can be in 5.0).

Giving a warning when the switched range is too large is a different thing, and mitigates damage from the current default. The existing warning is sensible for PME, and PME plus switch modifier, but it is not clear to me that it is sensible for non-PME coulomb, which is what I was asking about. It does seem sensible to warn for rvdw=cutoff + rvdw-modifier=potential-switch (and to do so in release-4-6).

I'd say a switch distance any bigger than 0.1 nm should be a warning for vdw, and 0.05 for coul. Can cite http://pubs.acs.org/doi/pdf/10.1021/ct4005068 for people to investigate further. HOWEVER that's with the old switching framework, and might not be as applicable with the potential-switch framework.

I can add these warnings by the weekend.

#13 Updated by Mark Abraham about 3 years ago

Michael Shirts wrote:

It's not "being set to zero," it defaults to zero and is never assigned to again. Reasonable default behaviour would be to be somewhat smaller than rvdw, but that behaviour is not implemented that I can see. You are remembering that we added a warning in the commit that I linked. I do not think we did anything about it being set.

Well, sure -- that's semantics. At creation, its value is set to zero by mdrun, and potential-switch uses that value, whether the user intended or not.

Right - we would like to default to things that are sensible. There are other places where grompp changes things (the rlistlong historical #@&% comes to mind), and this causes other problems with things not always being what the user intends (e.g. #1413).

We won't change the default behaviour in release-4-6, but we could still make release-5-0 default to rX - 0.02, or something similar we think sound. That's why we're doing these kinds of testing.

Agreed. We have data to suggest defaults, though it might be good to issue a note if the defaults are used (again, that can be in 5.0).

OK, we can write a note in release-4-6 when rX-switch is zero, regardless of source (since someone using a default mdp will have it "set" to zero also) and regardless of other warnings.

Giving a warning when the switched range is too large is a different thing, and mitigates damage from the current default. The existing warning is sensible for PME, and PME plus switch modifier, but it is not clear to me that it is sensible for non-PME coulomb, which is what I was asking about. It does seem sensible to warn for rvdw=cutoff + rvdw-modifier=potential-switch (and to do so in release-4-6).

I'd say a switch distance any bigger than 0.1 nm should be a warning for vdw, and 0.05 for coul. Can cite http://pubs.acs.org/doi/pdf/10.1021/ct4005068 for people to investigate further. HOWEVER that's with the old switching framework, and might not be as applicable with the potential-switch framework.

I can add these warnings by the weekend.

OK great.

#14 Updated by Gerrit Code Review Bot about 3 years ago

Gerrit received a related patchset '1' for Issue #1463.
Uploader: Michael Shirts ()
Change-Id: I7a2c29f87ceb04712aab5076ac97f5f22f573671
Gerrit URL: https://gerrit.gromacs.org/3377

#15 Updated by Erik Lindahl about 3 years ago

#16 Updated by Erik Lindahl about 3 years ago

Attached examples files used to validate the new free energy kernels. Both these mdp files start the switch from 0.6 nm, and set the interaction to zero from 0.8 nm. By using rlistlong=1.0 we get another 0.2nm buffer zone. Using double precision and the latest switch-patch, I got forces and dhdl.xvg to match perfectly for old & new setups.

#17 Updated by Erik Lindahl about 3 years ago

Michael: To faciliate comparisons, feel free to put a copy of your files somewhere we can access them, or just upload them here.

#18 Updated by Michael Shirts about 3 years ago

I've attached a the file corresponding to switch_11, lambda=0; all the deltas from this file are outlined below.

Let me know if you identify any lines that would cause problems so we can write warnings for those settings combinations.

I'll try running with your .mdp's to test cutoff dependence. I'll try double precision to see if that has an effect as well. It might take until tomorrow to run because of meetings today.

#19 Updated by Michael Shirts about 3 years ago

I noticed that in these files that DispCorr = 0. This would suggest that the dispersion correction is not handling the potential modifiers correctly -- at least the first thing that comes to mind. I'll look into this when I run things within a day or so, but an idea if somebody has a chance to look before I do.

Just so I understand, how exactly is rlistlong used now? it seems like the new style is to set rlist = rcoulomb, and rlistlong = long enough to give a buffer. Sorry if this was explained somewhere else . . .

#20 Updated by Michael Shirts about 3 years ago

OK, some confirmation. From a single configuration check.

calculation                                             vdw_dhdl                       LJ  short       disp corr.           
lam=0 lam=1 (l=0)-(l=1) (l=0)-(l=1)
potential-switch LJ, with dispersion corr 13.421916 -28.224947 29.83 3.40
switched LJ, with dispersion corr 13.811733 -27.839355 29.83 3.15
potential-switch LJ, without dispersion corr 8.4929857 -33.108677
switched LJ, without dispersion corr 8.4930153 -33.109299

So the potential-switch LJ does not handle the dispersion LJ the same way as switched. I suspect the problem will be in the potential-switch modifier, since the switched potential has a (properly) lower dependence of the free energy on the cutoff. I'll try to dig into the code and figure out what is going on. LJ-PME will hopefully be the long term solution!

Now, one question is if it's a free energy + dispersion issue, or issue with dispersion independent of free energy calculations. If we look at the total LJ energy for lambda=0, we get:

Cutoff                              Potential-switch                                 Old vdw switch
LJ(short) disp Total LJ (short) disp Total
9 6702.367676 -263.353546 = 6439.01413 6716.306152 -286.568787 = 6429.737365
10 6627.729492 -191.984726 = 6435.744766 6636.926758 -207.172409 = 6429.754349
11 6577.887695 -144.240952 = 6433.646743 6584.094727 -154.579010 = 6429.515717
12 6543.599609 -111.102257 = 6432.497352 6547.950195 -118.377060 = 6429.573135

Clearly the old style switches have a smaller radial dependence on the total energy, which is what we would expect from a switch. So this is likely a general dispersion issue, and not a dispersion + free energy issue.

Looking at the potential shift modifier

Cutoff                              Potential-shift                                 Old vdw switch
LJ(short) disp Total LJ (short) disp Total
9 6953.062012 -526.707092 = 6426.35492 6716.306152 -286.568787 = 6429.737365
10 6811.983887 -383.969452 = 6428.014435 6636.926758 -207.172409 = 6429.754349
11 6717.224121 -288.481903 = 6428.742218 6584.094727 -154.579010 = 6429.515717
12 6651.304688 -222.204514 = 6429.100174 6547.950195 -118.377060 = 6429.573135

Better, but not nearly as good independence of the energy on the cutoff, as we would expect from the free energy calculations results

Ah, good times. My first contribution to gromacs was fixing the dispersion correction so that it would properly incorporate the switching function instead of assuming a strict cutoff.

I'll see what I can trace down in the code.

#21 Updated by Rossen Apostolov about 3 years ago

  • Status changed from New to Accepted

#22 Updated by Erik Lindahl about 3 years ago

OK, I think I found the issue, and if that's the case it's trivial. The old dispersion correction code (obviously) explicitly checks for the "shift" setting, not the potential modifiers. I'm boarding a US flight in 15 minutes, so I won't be able to fix it before then, but hopefully tonight.

#23 Updated by Michael Shirts about 3 years ago

Erik Lindahl wrote:

OK, I think I found the issue, and if that's the case it's trivial. The old dispersion correction code (obviously) explicitly checks for the "shift" setting, not the potential modifiers. I'm boarding a US flight in 15 minutes, so I won't be able to fix it before then, but hopefully tonight.

OK, I'll check over the solution when posted! I'll then move to validating the Verlet cutoffs in 5.0 using the same methods.

#24 Updated by Rossen Apostolov about 3 years ago

  • Assignee changed from Berk Hess to Erik Lindahl

#25 Updated by Erik Lindahl about 3 years ago

  • Status changed from Accepted to Fix uploaded

Switch/shift modifiers should now update both tables and the dispersion correction evaluation, and I'm getting identical dispersion correction results for the first step. It's pretty late here, so I haven't had time to check a longer runs yet.

#26 Updated by Michael Shirts about 3 years ago

Erik Lindahl wrote:

Switch/shift modifiers should now update both tables and the dispersion correction evaluation, and I'm getting identical dispersion correction results for the first step. It's pretty late here, so I haven't had time to check a longer runs yet.

Great! I'm in the middle of running a conference and have a talk tomorrow, so I will test this as soon as I can, but that might be Thursday.

I should figure out if there is someone there I can pass off these free energy sensitivity test to, as it could be useful for a number of things (for example, sensitivity analysis of LJ-PME). I'll run it for this test now, though!

#27 Updated by Erik Lindahl about 3 years ago

@Michael: As I wrote in response #17, it would be good if you could put a copy of your files in this issue so we are testing on the same thing. This is probably the only thing holding up the release though, so unless somebody else does testing right away, Mark & I will work on it until we see exact reproduction of results, and then be happy with it for 5.0 - it will take too long to wait a few days (which easily turns into weeks) for each step of patch/testing.

#28 Updated by Michael Shirts about 3 years ago

Erik Lindahl wrote:

@Michael: As I wrote in response #17, it would be good if you could put a copy of your files in this issue so we are testing on the same thing.

The main mdp I was using was posted in #18. There were about 10 different mdp files, and the deltas from the root file were all described. The analysis I was performing included running rerun on a large number of .tpr files, which I could not post here, and then running a reweighting calculation on top of that to generate free energies. As I said, it would be good to get it into shape that someone there could run the sorts of analysis I have been doing, but that will take a few hours as well. Since a bunch of files need to be pre-generated, it will probably take just as long for me to do it, which I can start late tonight after conference close + travel, finishing tomorrow.

This is probably the only thing holding up the release though, so unless somebody else does testing right away, Mark & I will work on it until we see > exact reproduction of results, and then be happy with it for 5.0 - it will take too long to wait a few days (which easily turns into weeks) for each step > of patch/testing.

I agree, though this does go both ways -- it took a few weeks for the last fix to be put in because you didn't have time when I was able to test, and I have time today. I would agree that over time I've held up things more because of staffing issues here, so I'm not throwing stones, just reporting what I can do when.

#29 Updated by Erik Lindahl about 3 years ago

Oh, I certainly meant that it's just as much my fault. I will have today+tomorrow, but then be busy, so that's why I figured I should try to push it in for Mark ;-)

Could you put the nonanol coordinates and topology somewhere? Then we're at least working on the very same system (and parameters), so we should be able to compare energies directly.

I think I've traced down the problems with shift too now. Apparently, the shift dispersion correction wasn't applied correctly with verlet kernels (since it was just specified by the potential modifiers). Now I've generated tables for this case too (rather than re-implementing an analytical switch/shift dispersion correction), and just need to track down one final NaN for the shifted verlet kernels - then I'll start comparing energies.

#30 Updated by Gerrit Code Review Bot about 3 years ago

Gerrit received a related patchset '1' for Issue #1463.
Uploader: Erik Lindahl ()
Change-Id: I052071328d35506d1d899d44687f712011d79e8c
Gerrit URL: https://gerrit.gromacs.org/3476

#31 Updated by Erik Lindahl about 3 years ago

Alright, I think we're getting close to mission accomplished. With the latest patch now get the following for an octanol system, similar to the single frame tests Michael did above:

SWITCH 
calculation w. dispcorr          dhdl                  LJ  short       disp corr.     LJ+dispcorr      
                          lam=0           lam=1       (l=0)-(l=1)      (l=0)-(l=1)
 switched LJ           -98.545613       27248.331     -666.421241      4.457706       -661.963535
 potential-switch LJ   -98.545613       27248.342     -666.421241      4.457706       -661.963535

#32 Updated by Erik Lindahl about 3 years ago

Comparison between old and new switch functions for different cutoffs:

Cutoff                              Potential-switch                         Old vdw switch
                    LJ(short)     disp          Total            LJ (short)      disp              Total                
       9           10023.169465  -401.063987  = 9622.105479     10023.169451  -401.063987    = 9622.105464   
      10            9911.060081  -289.945613  = 9621.114468      9911.060066  -289.945613    = 9621.114454   
      11            9837.648580  -216.339118  = 9621.309463      9837.648565  -216.339118    = 9621.309448   
      12            9787.013538  -165.673158  = 9621.340380      9787.013523  -165.673158    = 9621.340365   

#33 Updated by Erik Lindahl about 3 years ago

Using our new potential shift; note that shift matches switch (see above) pretty closely now:

SHIFT
calculation w. dispcorr          dhdl                  LJ  short       disp corr.    LJ+dispcorr         
                          lam=0           lam=1       (l=0)-(l=1)      (l=0)-(l=1)   (l=0)-(l=1)
 potential-shift LJ    -98.630685       27248.263     -670.127399      8.246379      -661.88102

... and finally the cutoff dependence (or rather, non-dependence) for potential shift - here too we now have quite a good match with switch:

Cutoff                              Potential-shift       
                    LJ(short)     disp          Total      
       9           10354.402077  -735.028229  = 9619.373848
      10           10156.750084  -536.383840  = 9620.366244
      11           10023.924171  -403.228251  = 9620.695921
      12            9931.662616  -310.699437  = 9620.963179

The current patch in gerrit will fail regression tests since complex/nbnxn_rf and complex/nbnxn_pme need to be updated with the new values for potential-shift dispersion correction.

#34 Updated by Michael Shirts about 3 years ago

Some updates:

cutoff type             dG           ddG_Ref       dG(lam=0)    EPot(l=0) dhdlC (l=0) dhdlVDW (l=0)  EPot (l=1) dhdlC (l=1) dhdlVDW (l=1)
switch_10_v2 -4.529+ 0.587 0.003+-0.005 -0.708+-0.004 -44416.379 52.975113 -27.836592 -43863.340 362.61572 13.813869
switch_9_v2 -4.470+-0.588 -0.056+-0.021 -2.459+-0.017 -44418.578 52.970840 -27.800959 -43865.031 362.64972 13.785430
shift_10_v2 -4.436+-0.587 -0.091+-0.008 -1.326+-0.006 -44417.156 52.973640 -27.811489 -43863.898 364.48270 13.726041
shift_9_v2 -4.312+-0.588 -0.215+-0.014 -4.499+-0.012 -44420.562 52.970573 -27.805754 -43866.977 364.51108 13.679724
cut_10 -4.526+-0.587 -0.000+-0.000 0.000+-0.000 -44415.766 52.974979 -27.839355 -43862.621 362.61575 13.811733

Reference is cutoff with 10. Switch now appears to be identical to the old switch (old switch at 9 was - 0.048 +/- 0.021 from ref).

shift is a little different. I'll keep examining this and look at the code. From what Erik has put, this looks like the new shift is indeed implemented preserving the functionality, it may just be that correcting the shifts as a function of lambda are problematic because of changing dispersion contributions. I.e. this is not a code error, this is a methodology issue.

There was a presentation at the Workshop on Free Energy Methods in Drug Design by a researcher at D.E. Shaw where they examined the errors that result from a shifted instead of a switched potential, because of the changing value of the cutoff with lambda. I'll see if I can get the notes to forward on.

I am OK with pushing the code fix through for 4.6.6, as new keywords now have no functionality loss from the old key words.

It appears that potential-switch is preferable for free energy calculations for now, as shift has a small-but-still-noticeable dependence on distance. I (happy for someone else to as well, of course) will investigate for 5.1 to see it it can be improved.

I will try to get the analysis method packaged up nicely in the next day or so other people can run them easily. After an initial investment of a few hours/trr for the trrs, the calculation runs in about 15-20 min per parameter set. A little rough for automated regression, but something done fairly easily in a separate run.

#35 Updated by Erik Lindahl about 3 years ago

Just formatting Michael's table for my reading:

cutoff type dG ddG_Ref dG(lam=0) EPot(l=0) Dhdlcoul(l=0) dhdlVDW(l=0) EPot(l=1) Dhdlcoul(l=1) dhdlVDW(l=1)
switch_10_v2 -4.529+ 0.587 0.003+-0.005 -0.708+-0.004 -44416.379 52.975113 -27.836592 -43863.340 362.61572 13.813869
switch_9_v2 -4.470+-0.588 -0.056+-0.021 -2.459+-0.017 -44418.578 52.970840 -27.800959 -43865.031 362.64972 13.785430
shift_10_v2 -4.436+-0.587 -0.091+-0.008 -1.326+-0.006 -44417.156 52.973640 -27.811489 -43863.898 364.48270 13.726041
shift_9_v2 -4.312+-0.588 -0.215+-0.014 -4.499+-0.012 -44420.562 52.970573 -27.805754 -43866.977 364.51108 13.679724
cut_10 -4.526+-0.587 -0.000+-0.000 0.000+-0.000 -44415.766 52.974979 -27.839355 -43862.621 362.61575 13.811733

#36 Updated by Erik Lindahl about 3 years ago

never mind, that' didn't help much - I'll simply trust Michael ;-)

#37 Updated by Michael Shirts about 3 years ago

Yeah, I tried for a while to get that to work. This seems to show up better.

#38 Updated by Mark Abraham about 3 years ago

There's a "pre" button for preformatted text or just

<pre>
Do it like this
</pre>

#39 Updated by Berk Hess about 3 years ago

  • Status changed from Fix uploaded to Resolved
  • % Done changed from 0 to 100

#40 Updated by Erik Lindahl about 3 years ago

  • Status changed from Resolved to Closed

Also available in: Atom PDF