Project

General

Profile

Bug #160

Discrepancy in mean constraint force between free energy and pull code

Added by Chris Fennell over 12 years ago. Updated about 12 years ago.

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

Description

We ran into an issue when calculating PMF curves between charged objects, where the curves were distorted such that we didn't believe them. As a check, we decided to follow the procedure for NaCl in:

B. Hess, C. Holm; N. van der Vegt, J. Chem. Phys., 124, 164509 (

examples.tar.gz (120 KB) examples.tar.gz Example test cases for the 2 methods Chris Fennell, 08/16/2007 07:31 PM

History

#1 Updated by Erik Lindahl over 12 years ago

Hi,

This report seems to have been truncated, so we need more info and a test case if you have it.

#2 Updated by Chris Fennell over 12 years ago

Created an attachment (id=228)
Example test cases for the 2 methods

#3 Updated by Chris Fennell over 12 years ago

(In reply to comment #1)

Hi,

This report seems to have been truncated, so we need more info and a test case
if you have it.

That's bizarre. It was actually a rather long message when I submitted it. I'll try to reproduce it here:

...

We ran into an issue when calculating PMF curves between charged objects, where
the curves were distorted such that we didn't believe them. As a check, we
decided to follow the procedure for NaCl in:

Hess, C. Holm; N. van der Vegt, J. Chem. Phys., 124, 164509 (2006)

In the test, we made an NaCl molecule with the following constraint:

1 2 2 0.2 1.2

and performed several simulations with init-lambda values ranging from 0 to 1. The series of simulations were performed with all the same settings and specifics in the above paper; however, our test system contained 864 SPC/E molecules (rather than 1000) and we used the OPLS Na+ and Cl- parameters (rather than the KBFF parameters). The resulting mean forces for the sequence of simulations are as follows:

  1. r(nm) mf
    .20 -42191.92
    0.22 -11105.74
    0.24 -2843.42
    0.26 -526.12
    0.28 97.51
    0.30 220.77
    0.32 182.80
    0.34 92.70
    0.36 7.55
    0.38 -63.06
    0.40 -101.02
    0.42 -109.66
    0.44 -81.03
    0.46 -49.21
    0.48 -23.55
    0.50 -9.46
    0.54 29.80
    0.58 23.74
    0.62 -14.04
    0.66 -10.40
    0.70 0.47
    0.74 6.36
    0.78 13.38
    0.82 10.71
    0.86 3.90
    0.90 9.40
    0.94 11.47
    0.98 13.15
    1.02 12.14
    1.06 10.09
    1.10 15.35
    1.14 26.67
    1.18 35.29

If we do the exact same set of simulations, this time with the pull code (turn free-energy off and use an appropriately assembled pull.ppa and index file), we obtain the following sequence of mean force values:

  1. r(nm) mf
    0.20 -43375.30
    0.22 -11716.36
    0.24 -2913.94
    0.26 -537.87
    0.28 103.29
    0.30 214.88
    0.32 176.32
    0.34 93.92
    0.36 12.25
    0.38 -67.32
    0.40 -106.00
    0.42 -109.93
    0.44 -87.99
    0.46 -46.95
    0.48 -23.41
    0.50 -16.17
    0.54 27.15
    0.58 19.06
    0.62 -18.13
    0.66 -16.98
    0.70 -9.76
    0.74 0.06
    0.78 3.69
    0.82 -8.62
    0.86 -10.71
    0.90 -5.79
    0.94 -4.48
    0.98 -0.73
    1.02 -2.73
    1.06 -1.72
    1.10 -5.79
    1.14 -2.07
    1.18 -1.11

The PMF generated from this list of mean force values is a good match to the PMF obtained by inverting the NaCl g(r) from a simulation of a single ion pair in 864 SPC/E waters (run for 100 ns for decent bin statistics). The PMF generated from the previous free-energy code list continually increases with ion separation. Just looking at the mean force values, it appears as though there is an added attractive interaction in the free-energy code results.

We've been using the pull code as an immediate solution; however, we figured it would be a good idea to notify you guys of this issue. I'm including a set of files for doing a single simulation (0.9 nm separation) with the free energy code and with the pull code so that you can verify the mean force output we're seeing above. If you need anything else, just let me know.

-Chris

#4 Updated by Berk Hess about 12 years ago

Did you properly take care of the purely entropic volume
term that occur when constraining the distance between
two particles in 3D?

Berk.

#5 Updated by Chris Fennell about 12 years ago

(In reply to comment #4)

Did you properly take care of the purely entropic volume
term that occur when constraining the distance between
two particles in 3D?

Berk.

Hi Berk,

Yes, I believe so. The above list are the pure mean constraint force values that come from the simulation. Niether of them have the 4\pi r^2 term removed from the increase in available configuration space as the separation distance increases. We typically treat this term upon integration to obtain the potential of mean force. Doing so with the pull code mean force values generates the expected potential of mean force along a line, but this is not the case for the values from the free energy code simulations.

I was under the impression that the mechanics of both implementations of constraints were the same (i.e., only the separation distance is constrained and both are free to tumble in space). If this is not the case, I can certainly see this entropic surface area term playing a role, since the previously posted mean force values would not be equivalent.

It's possible we've been setting something else up incorrectly with the simulations using the free energy code (feel free to look at the previously posted example files). In the meantime, we're just using the pull code as a fix.

Thanks!

-Chris

#6 Updated by Berk Hess about 12 years ago

I think the issue is that you put both ions in the same charge group.

Berk.

#7 Updated by Chris Fennell about 12 years ago

(In reply to comment #6)

I think the issue is that you put both ions in the same charge group.

Berk.

Ugh. Yeah, that would probably do it. Thanks for pointing it out. Sorry for taking up your time with this.

Thanks again,

-Chris

#8 Updated by Berk Hess about 12 years ago

This problem was not a bug, but caused by an incorrect charge group setting
in the topology.

Also available in: Atom PDF