normal have problem to repeat PMF with martini polarizable water model

  • neofloat
  • neofloat's Avatar Topic Author
  • Offline
  • Fresh Boarder
More
11 years 1 month ago - 11 years 1 month ago #1480 by neofloat
Hi All,

This is the first time we use Martini force filed and CG polarizable water model. So we decide to repeat the result in the recent published paper.
Title: Polarizable Water Model for the Coarse-Grained MARTINI Force Field
link:
http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1000810
We use Gromacs version 4.6-beta3, single precision.
we want to repeat figure 9 in this paper, the polarizable martini sodium result.

We build a system containing 128 DPPC, 1 NA+, 1CL-, 2000 CG polarizable water, and we equilibrate the system for about 500ns. The area per lipid we got is around 0.635 nm2. it is the same value reported in the paper. The number density of different cg groups of DPPC and water also match the density reported in the paper.
Then we set up 40 windows to compute the pmf of translocation of NA+ through a DPPC bilayer, by using umbrella sampling method.
To set up each window, we place sodium ion at the ideal place and turn off the interaction between sodium and other beads in the system
, then we slowly grow in the vdw interactions between sodium and the the rest of the system in 10 ns, after that we grow in the vdwq interaction in the same way.



The results is shown in this figure.

This image is hidden for guests.
Please log in or register to see it.


http://udel.edu/~yuanhu/pmf-plos-repeat.png

see the mdp file, and log file of the center of bilayer window:

http://udel.edu/~yuanhu/tau-us.mdp
http://udel.edu/~yuanhu/us.mdp (grompp output mdp file)
http://udel.edu/~yuanhu/us.log


We add the umbrella sampling code to the end the example mdp file from this website ( md.chem.rug.nl/cgmartini/images/paramete...ini_v2.P_example.mdp ), it's also the same mdp file shown in the supporting information of this paper.
===============================
; Pull code
pull = umbrella
pull_geometry = distance
pull_dim = N N Y
pull_start = no
pull_ngroups = 1
pull_group0 = DPPC
pull_group1 = NA+
pull_init1 = 0.9
pull_rate1 = 0.0
pull_k1 = 400

ld_seed=-1
===============================

In order to match the PMF, we tried several test, but all of the tests give the similar results.



The systems we have tested are as follow:

a. T=320k, tau_p=1.0, tau-t=1.0, nstlist=10, rlist=1.2,compressibility=3e-4, tc-grps=DPPC PW_ION
h. T=323k, tau_p=5.0, tau-t=1.0, nstlist=10, rlist=1.2,compressibility=3e-4, tc-grps=DPPC PW_ION
i. T=323k, tau_p=5.0, tau-t=1.0, nstlist=5, rlist=2,compressibility=3e-4, tc-grps=DPPC PW_ION
j. T=323k, tau_p=3.0, tau-t=0.3, nstlist=5, rlist=2,compressibility=3e-4, tc-grps=DPPC PW_ION
k. T=323k, tau_p=3.0, tau-t=0.3, nstlist=5, rlist=2,compressibility=4.5e-5, tc-grps=DPPC PW_ION
l. T=323k, tau_p=3.0, tau-t=0.3, nstlist=5, rlist=2,compressibility=4.5e-5, pull_vec1= 0 0 1
o. T=323k, tau_p=3.0, tau-t=0.3, nstlist=5, rlist=2,compressibility=4.5e-5,comm-grps= DPPC PW_ION
p. T=323k, tau_p=3.0, tau-t=0.3, nstlist=5, rlist=2,compressibility=4.5e-5,tc-grps=DPPC_ION PW
q. T=323k, tau_p=3.0, tau-t=0.3, nstlist=5, rlist=2,compressibility=4.5e-5,comm-grps= DPPC PW_ION, tc-grps=DPPC PW ION
s. T=323k, tau_p=3.0, tau-t=0.3, compressibility=4.5e-5, comm-grps= DPPC PW ION, tc-grps=DPPC PW ION
r. T=325k, tau_p=3.0, tau-t=0.3, nstlist=5, rlist=2,compressibility=4.5e-5,comm-grps= DPPC PW_ION


we also tried to use gromacs 3.3.1 version (single precision), but the pmf didn't match either.
We don't know where is the problem. We run out of ideas now. Could someone help us? Thanks.

Please Log in or Create an account to join the conversation.

More
11 years 1 month ago #1481 by djurre
I saw there is a double post at the Gromacs mailing list, but I'll still try to answer it here.
First of all I don't think the difference of +/-3kJ/mol in the deepest well of you PMF is so much. It is most likely less the accuracy of the martini model. However, it seems to be significant and therefor should not be the case.

A few things that might be reason:
-The PMF simulations in the paper might also have been done with 4.0.2 or 4.0.5 (since those are also listed in the paper). There might be differences in the mdrun, but also in g_wham?
-There might be differences in counter ions. You write you have one Na+ and one Cl- ion in each simulation. In the paper the authors are not very clear, but it might be they have no, or more ions.

Before you spend more time on this, I send an e-mail to one of the authors (neither is in our lab any more). Let's wait for that.

Please Log in or Create an account to join the conversation.

  • neofloat
  • neofloat's Avatar Topic Author
  • Offline
  • Fresh Boarder
More
11 years 1 month ago #1482 by neofloat
Thanks for your reply. I really appreciate. We spend about a month to solve the problem. We don't know whether different version of Gromacs(or g_wham) will give different result. We tried 0.4M NaCl and no counterion systems. They didn't change the shape of the PMF.
Will the precision of the gromacs give different results?
We use single precision. We don't know which type of precision that the author used.

Please Log in or Create an account to join the conversation.

More
11 years 4 weeks ago #1487 by xavier
dear neofloat,

I agree with Djurre the difference you see might not be important and t judge this it would be good to see a plot in which you place all the PMFs you obtained with the different setups you tried.

neofloat wrote: Thanks for your reply. I really appreciate. We spend about a month to solve the problem. We don't know whether different version of Gromacs(or g_wham) will give different result. We tried 0.4M NaCl and no counterion systems. They didn't change the shape of the PMF.
Will the precision of the gromacs give different results?
We use single precision. We don't know which type of precision that the author used.

Please Log in or Create an account to join the conversation.

More
11 years 4 weeks ago #1491 by kmcallenberg
Hi everyone,

I see similar results to neofloat when I perform these calculations. In other words, the ion is slightly less stable in the headgroup region than in the PLoS paper. However until neofloat's post, I was happy to find them so close to the paper's results (haha). I am also using Gromacs 4.5.3. I will post my data soon.

Keith

Please Log in or Create an account to join the conversation.

More
11 years 4 weeks ago #1492 by xavier
Thanks for reporting.

I now realize you both use more recent gromacs versions than the originally used (3.3.1/4.0.2/4.0.5). If any of you could repeat the PMF using the 4.0.5 version that would be useful but first let's look at your data first if you can share it.

kmcallenberg wrote: Hi everyone,

I see similar results to neofloat when I perform these calculations. In other words, the ion is slightly less stable in the headgroup region than in the PLoS paper. However until neofloat's post, I was happy to find them so close to the paper's results (haha). I am also using Gromacs 4.5.3. I will post my data soon.

Keith

Please Log in or Create an account to join the conversation.

  • neofloat
  • neofloat's Avatar Topic Author
  • Offline
  • Fresh Boarder
More
11 years 4 weeks ago #1493 by neofloat
Hi xavier and kmcallenberg,

Thanks for your reply.
You can find the test result from the following link.
udel.edu/~yuanhu/pmf.pdf

We tested the version 3.3.1/4.0.2/4.0.5/4.6.beta3 (single precision), and 4.6.1 (double precision). The PMFs almost overlap with the result from 4.6.beta3 version.



This is the first time we use gromacs. We hope we didn't miss some important things.
Thanks for your help.

Please Log in or Create an account to join the conversation.

More
11 years 2 days ago #1548 by djurre
Dear Neofloat, Keith,
I've been asking around in the lab about this problem, but we can't really come up with a completely satisfactory explanation (that is why it took so long, sorry). I can't reach the person who did the original experiments. Possible explanations that came up are the way the systems have been prepared that might have led to not well equilibrated simulations or the way the data has been treated afterwards (WHAM, graph smoothing?).

However, I would like to stress, as I did in my original post, that the difference between the graphs is small, 3kJ/mol. 3kJ/mol does correspond to a relative concentration difference of about 3-fold. This seems a lot, but it is a local effect. Over the whole membrane the effect will be negligible.

For simulations with the model, one would normally not expect large differences. For systems where this very small difference does play a role, the Martini model is probably not your model of choice, it is in the end a "semi-quantitative" model.

Please Log in or Create an account to join the conversation.

Time to create page: 0.112 seconds