Discrepancy in mean constraint force between free energy and pull code
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 (
#3 Updated by Chris Fennell almost 13 years ago
(In reply to comment #1)
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:
- r(nm) mf
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:
- r(nm) mf
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.
#5 Updated by Chris Fennell almost 13 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?
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.