Re: [gmx-users] Gromos DPPC bilayer: different results for 4.0.7 and 4.6.5
Hi, I answer a long time after the last discussion of this thread, but I did some tests in between and discussed with some people. I would like to follow-up on this issue and sum-up what I found. After digging into the mailing list, I finally found that the same type of issue was discussed in 2011 between Antonio Baptista/Miguel Machuqueiro and gromacs devs. There a typical GROMOS setup (twin-range cutoff 0.8/1.4 + RF) was giving LINCS crashes (protein simulations in water) in 4.5 whereas it was stable in 4.0. It looks like that what we see with the Poger lipids has the same origin (the new reversible algorithm introduced in 4.5). We don't have crashes because there are no explicit H in the lipids. So to sum up the situation, if we want a stable simulation with 4.5 and the GROMOS typical simulation parameters, that give similar results to 4.0 there seem to be 2 possibilies (maybe there are others, but at least these 2 have been tested): 1) reduce nstlist to 1 (2 seems to work according to Miguel) 2) increase rlist up to rlist=rcoulomb (I tried Samuli's suggestion rlist=rcoulomb and recover a correct area on the DPPC system, see http://redmine.gromacs.org/issues/1400 for the plot). Both solutions increase a lot the computational cost. For example 2) gives 1.5 ns/h on my eight-core computer, whereas in 4.0 I had ~ 1.0 ns/h. In the discussion of 2011 Berk proposed some other possible solutions to speed up with 4.6 (https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-developers/2011-December/005506.html), but these were not tested to my knowledge. My concern is that a twin-range/RF setup in version = 4.0 gives a given result, but does'nt give the same result in = 4.5 and possibly crashes. So I shall say that it would be nice to warn the users that they will encounter problems when using twin-range/RF with versions = 4.5. Maybe with a grompp warning? Section 3.4.7 in the manual may not be the first place a user will look at. Anyhow, if for some reason someone wants to use the twin-range cut-off as it is done in the GROMOS software (with all versions of GROMOS force fields), one has to stick to 4.0.* or lower. Best, Patrick Le 18/12/2013 21:26, Ollila Samuli a écrit : Hi, On Dec 18, 2013, at 8:50 AM, Sabine Reisser sabine.reis...@kit.edu wrote: Up to Gromacs 4.0.7, the long rance forces, which are calculated at every nstlist'th step, were added constantly to the short range force in the subsequent steps. Starting from Gromacs 4.5., the long range force is added nstlist times to the short range force on the nstlist'th step, while there is nothing added to the short range force on the subsequent steps. This creates some kind of impulse force which we believe is responsible for the different behaviour of the membrane, although we haven't figured out the exact mechanism yet. That's interesting. Now I see that this change is documented in the manual (Section 3.4.7). Would that imply that this effect does not depend on the model for long-range electrostatics (Reaction Field or PME), but it comes from using twin range cutoffs? I also noticed the same section in the manual today. I think that what possibly happens is that with the reaction field the forces between rlist and rcoulumb are not used in the steps between neighbour list update. However, in these steps the reaction field force calculation assumes the dielectric envinronment beoynd rcoulomb, not beoynd rlist even though the forces only up to rlist are used. I am not sure if this would have any practical relevance, though? I am also not sure how this is handled in the pressure calculation? I think that the reason why Patrick did not get the difference in the PME results between versions might be that he probably used rcoulomb=rlist in those simulations. If I try to set rcoulombrlist with PME, I get error: With coulombtype = PME, rcoulomb must be equal to rlist If you want optimal energy conservation or exact integration use PME-Switch If it is true what I am saying, then the area per molecules should be the same between versions also using RF and setting rlist=rcoulomb. BR, Samuli Ollila I expect you to get the 'right' area per lipid if you use nstlist 1 as already suggested by Mirco, because in that case the multistep algorithms give the same answer. This is also our experience with G54A7. I started a run yesterday with nstlist=1, and so far (about 12 ns) the area per lipid indeed stays at the correct value of around 0.63 nm^2. Thanks, Lutz -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org. -- ___ Patrick
Re: [gmx-users] Gromos DPPC bilayer: different results for 4.0.7 and 4.6.5
On Sat, Dec 28, 2013 at 1:43 AM, Lutz Maibaum lutz.maib...@gmail.com wrote: On Dec 24, 2013, at 2:29 AM, Szilárd Páll pall.szil...@gmail.com wrote: Just wondering, has anybody done a comparison with the Verlet scheme? It could be useful to know whether it produces results consistent with the 4.6 group scheme implementation or exhibits different behavior. I updated the comparison chart with another simulation that uses the Verlet cutoff scheme: http://faculty.washington.edu/maibaum/dppc_comparison/comparison.png It might be too early to tell, but I’d say that there is no significant difference between using the Verlet and the group scheme in Gromacs 4.6.5 — they both give an area per lipid that too low (compared to results obtained with Gromacs 4.0.7). Using the Verlet scheme is about two times slower in performance than using the group cutoff scheme. Note that the Verlet scheme calculated the buffer based on the given drift, so when using it nstlist can be chosen freely. As it has far better energy conservation than the group scheme, it may not make much difference, but you can make the drift (originating from the non-bondeds) even lower by setting the verlet-buffer-drift mdp option. I do not expect this to change much though, as what you are observing is probably caused by something else. Cheers, -- Szilárd Best, Lutz -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Gromos DPPC bilayer: different results for 4.0.7 and 4.6.5
On Dec 24, 2013, at 2:29 AM, Szilárd Páll pall.szil...@gmail.com wrote: Just wondering, has anybody done a comparison with the Verlet scheme? It could be useful to know whether it produces results consistent with the 4.6 group scheme implementation or exhibits different behavior. I updated the comparison chart with another simulation that uses the Verlet cutoff scheme: http://faculty.washington.edu/maibaum/dppc_comparison/comparison.png It might be too early to tell, but I’d say that there is no significant difference between using the Verlet and the group scheme in Gromacs 4.6.5 — they both give an area per lipid that too low (compared to results obtained with Gromacs 4.0.7). Using the Verlet scheme is about two times slower in performance than using the group cutoff scheme. Best, Lutz -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Gromos DPPC bilayer: different results for 4.0.7 and 4.6.5
Just wondering, has anybody done a comparison with the Verlet scheme? It could be useful to know whether it produces results consistent with the 4.6 group scheme implementation or exhibits different behavior. Cheers, -- Szilárd On Sun, Dec 22, 2013 at 1:08 AM, Lutz Maibaum lutz.maib...@gmail.com wrote: On Dec 20, 2013, at 3:17 PM, Patrick Fuchs patrick.fu...@univ-paris-diderot.fr wrote: to follow up on this, the simulation with Reaction-Field-nec under 4.5.3 has completed. The final area is 0.60 nm^2 (see http://redmine.gromacs.org/issues/1400 for the plot). Lutz, you mentionned you tried a simulation with nstlist=1, could you give feedback on red mine? I uploaded another comparison to redmine: http://redmine.gromacs.org/attachments/download/1114/comparison.png It indeed looks like setting nstlist=1 allows one to recover the results form Gromacs 4.0.7 with nstlist=5 (but it is computationally not really feasible). Best, Lutz -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Gromos DPPC bilayer: different results for 4.0.7 and 4.6.5
On 21.12.2013 00:17, Patrick Fuchs wrote: Hi, to follow up on this, the simulation with Reaction-Field-nec under 4.5.3 has completed. The final area is 0.60 nm^2 (see http://redmine.gromacs.org/issues/1400 for the plot). Lutz, you mentionned you tried a simulation with nstlist=1, could you give feedback on redmine? nstlist=5 vs nstlist=1 (preliminary results): http://imgur.com/eQlC0oy Regards M. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Gromos DPPC bilayer: different results for 4.0.7 and 4.6.5
On Dec 20, 2013, at 3:17 PM, Patrick Fuchs patrick.fu...@univ-paris-diderot.fr wrote: to follow up on this, the simulation with Reaction-Field-nec under 4.5.3 has completed. The final area is 0.60 nm^2 (see http://redmine.gromacs.org/issues/1400 for the plot). Lutz, you mentionned you tried a simulation with nstlist=1, could you give feedback on red mine? I uploaded another comparison to redmine: http://redmine.gromacs.org/attachments/download/1114/comparison.png It indeed looks like setting nstlist=1 allows one to recover the results form Gromacs 4.0.7 with nstlist=5 (but it is computationally not really feasible). Best, Lutz -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Gromos DPPC bilayer: different results for 4.0.7 and 4.6.5
Hi, I ran into the same problem comparing the subsequent G54A7 lipid force field for DPPC. We could trace the relevant (trajectory-changing) differences in the code down to the transition Gromacs 4.0.7 - 4.5., and here to the change in the multistep update algorithm. Up to Gromacs 4.0.7, the long rance forces, which are calculated at every nstlist'th step, were added constantly to the short range force in the subsequent steps. Starting from Gromacs 4.5., the long range force is added nstlist times to the short range force on the nstlist'th step, while there is nothing added to the short range force on the subsequent steps. This creates some kind of impulse force which we believe is responsible for the different behaviour of the membrane, although we haven't figured out the exact mechanism yet. I expect you to get the 'right' area per lipid if you use nstlist 1 as already suggested by Mirco, because in that case the multistep algorithms give the same answer. This is also our experience with G54A7. Cheers Sabine On 12/17/2013 11:45 PM, Mirco Wahab wrote: On 17.12.2013 20:58, Lutz Maibaum wrote: What would the next step be trying to debug this issue? You could, just for debugging, give nstlist=1 a try (instead of nstlist=5)? I tried this today on a quadrupled version of your system (duplicated in x, and y direction, 512 lipids), which resulted at an area of around 0.605 nm² (still growing, incredible slow computation). The very same system at nstlist=5 ends up at around 0.575 nm². Maybe this is a hint for the developers? Regards M. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Gromos DPPC bilayer: different results for 4.0.7 and 4.6.5
On Dec 18, 2013, at 8:50 AM, Sabine Reisser sabine.reis...@kit.edu wrote: Up to Gromacs 4.0.7, the long rance forces, which are calculated at every nstlist'th step, were added constantly to the short range force in the subsequent steps. Starting from Gromacs 4.5., the long range force is added nstlist times to the short range force on the nstlist'th step, while there is nothing added to the short range force on the subsequent steps. This creates some kind of impulse force which we believe is responsible for the different behaviour of the membrane, although we haven't figured out the exact mechanism yet. That's interesting. Now I see that this change is documented in the manual (Section 3.4.7). Would that imply that this effect does not depend on the model for long-range electrostatics (Reaction Field or PME), but it comes from using twin range cutoffs? I expect you to get the 'right' area per lipid if you use nstlist 1 as already suggested by Mirco, because in that case the multistep algorithms give the same answer. This is also our experience with G54A7. I started a run yesterday with nstlist=1, and so far (about 12 ns) the area per lipid indeed stays at the correct value of around 0.63 nm^2. Thanks, Lutz -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Gromos DPPC bilayer: different results for 4.0.7 and 4.6.5
Thank you all for your comments! It's interesting that Patrick ran into the same issues just a few days earlier. Following your suggestions, I started some simulations with the Reaction-Field-nec coulombtype. So far those have run for only about 30 nanoseconds, but based on the attached graph I would assume that they will still result in an area per lipid that is lower than that obtained with Gromacs 4.0.7. I will update the bug report on redmine once the simulations have progressed a little more. What would the next step be trying to debug this issue? Is it possible to output only the various components of the forces (in particular, those induced by the reaction field)? Best, Lutz -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Gromos DPPC bilayer: different results for 4.0.7 and 4.6.5
On Wed, Dec 18, 2013 at 6:58 AM, Lutz Maibaum lutz.maib...@gmail.comwrote: Thank you all for your comments! It's interesting that Patrick ran into the same issues just a few days earlier. Following your suggestions, I started some simulations with the Reaction-Field-nec coulombtype. So far those have run for only about 30 nanoseconds, but based on the attached graph I would assume that they will still result in an area per lipid that is lower than that obtained with Gromacs 4.0.7. I will update the bug report on redmine once the simulations have progressed a little more. What would the next step be trying to debug this issue? I would suggest using mdrun -rerun magic-configuration.gro -s some.tpr on whatever combinations of .mdp settings and GROMACS versions make sense, to observe that there is a difference that can be reproduced, and which .edr fields show it. Then I'd suggest browsing the release notes for any obvious clues, once there's a more clear idea what to look for (e.g. 4.0.7 was fortuitously correct despite having a bug that was fixed). In extremis, one could use a git bisect to walk back through history and find the commit in which a difference arose, but doing that is reasonable only when one can write a script that executes reasonably fast that can serve to identify whether the behaviour is old or new (hence the mdrun -rerun suggestion; equilibrating a membrane is not on!). Once the code change is identified, the correctness of one version the code can be assessed. Is it possible to output only the various components of the forces (in particular, those induced by the reaction field)? Not directly. The non-bonded kernels for RF compute both 1/r and the correction terms before recording any data. Note that the RF terms need not be the primary problem - using the wrong form of long-range correction would be a possible candidate, and there's doubtless others. Mark Best, Lutz -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Gromos DPPC bilayer: different results for 4.0.7 and 4.6.5
On Dec 17, 2013, at 1:55 PM, Mark Abraham mark.j.abra...@gmail.com wrote: I would suggest using mdrun -rerun magic-configuration.gro -s some.tpr on whatever combinations of .mdp settings and GROMACS versions make sense, to observe that there is a difference that can be reproduced, and which .edr fields show it. I did that for the 53A6-L forcefield and the configuration posted at http://faculty.washington.edu/maibaum/dppc_comparison/sample.gro The various energies reported by Gromacs 4.0.7 and 4.6.5 are those: http://faculty.washington.edu/maibaum/dppc_comparison/g_energy_comparison.pdf It does look like the main differences come from the virial (and therefore pressure). I was surprised that even the temperatures are slightly different. These calculations are done in single precision. correction terms before recording any data. Note that the RF terms need not be the primary problem - using the wrong form of long-range correction would be a possible candidate, and there's doubtless others. Dispersion correction has been turned off, so hopefully that's not the culprit (but of course you are right that there are plenty others). Best, Lutz -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Gromos DPPC bilayer: different results for 4.0.7 and 4.6.5
Hi, Essentially similar results have been reported recently in redmine: http://redmine.gromacs.org/issues/1400 The issue seems to be related to the reaction field since with PME the results were the same between the versions. It might be useful to check if you get the similar pressure differences between versions also with PME? If not, then the pressure calculation with RF in 4.5. and 4.6 might be the issue. In this case also the different components of pressure might be interesting. Maybe it is better to report this and further findings by updating redmine issue 1400? BR, Samuli Ollila From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se [gromacs.org_gmx-users-boun...@maillist.sys.kth.se] on behalf of Lutz Maibaum [lutz.maib...@gmail.com] Sent: Monday, December 16, 2013 10:10 PM To: gromacs.org_gmx-users@maillist.sys.kth.se Subject: [gmx-users] Gromos DPPC bilayer: different results for 4.0.7 and 4.6.5 I am running MD simulations of DPPC bilayers with the Gromos force field, and I am seeing some differences between using Gromacs 4.0.7 and 4.5.5/4.6.5 that I do not understand. It would be great if someone had any insight into what's going on here. When I use the Gromos 53A6-L force field, which is the default 53A6 force field combined with improved lipid parameters that can be downloaded from http://compbio.chemistry.uq.edu.au/~david/research/lipids.htm, I obtain an average area per lipid of 0.627 nm^2, in good agreement with both the paper that describes these new parameters (0.629 nm^2, obtained with Gromacs 3.2.1, Ref. [1]) and another follow-up study (0.631 nm^2 and 0.623 nm^2, Gromacs 4.0.7, Ref. [2]). Now, if I use the same force field and mdp file, and the same initial configuration (which is a pre-equilibrated DPPC bilayer from http://compbio.biosci.uq.edu.au/atb/system_download.py?boxid=32 and randomly generated velocities), but use Gromacs 4.6.5 instread, I get a lower value of about 0.59 nm^2. I also tried the Gromos 54A7 force field, which is included with Gromacs 4.6.5 and that should be identical to 53A6-L (plus it has some other improvements over 53A6 that shouldn't be relevant here), I also get the lower area per lipid of ~0.59 nm^2. If have attached a plot of the area per lipid for these three simulations, each more than 100ns long. This looks to me like 4.6.5 give systematically lower area per lipid than 4.0.7. Running additional simulations with Gromacs 4.5.5 suggest that that also results in the lower area per lipid. Does anyone know why this might be? I have uploaded the relevant files in case that is helpful: http://faculty.washington.edu/maibaum/dppc_comparison/ To see if there are any differences between the energies that 4.0.7 and 4.6.5 compute, I picked a configuration, and used mdrun -rerun with the three different gromacs / force field combinations. Here is what I get for the sample.gro configuration (included in the link above): Gromacs 4.0.7 + Gromos53A6-L: Energies (kJ/mol) G96Bond G96AngleProper Dih. Improper Dih. LJ-14 1.76906e+021.29521e+049.57922e+034.97638e+02 -1.38126e+03 Coulomb-14LJ (SR)LJ (LR) Coulomb (SR) Coulomb (LR) 1.50791e+041.23533e+04 -7.49664e+03 -3.78489e+05 -8.97363e+02 RF excl. PotentialKinetic En. Total EnergyTemperature -5.25414e+04 -3.90168e+056.47185e+04 -3.25449e+052.87013e+02 Pressure (bar) 2.00740e+02 Gromacs 4.6.5 + Gromos53A6-L: Energies (kJ/mol) G96Bond G96AngleProper Dih. Improper Dih. LJ-14 1.76905e+021.29522e+049.57922e+034.97637e+02 -1.38126e+03 Coulomb-14LJ (SR)LJ (LR) Coulomb (SR) Coulomb (LR) 1.50791e+041.23533e+04 -7.49659e+03 -3.78489e+05 -8.97344e+02 RF excl. PotentialKinetic En. Total EnergyTemperature -5.25414e+04 -3.90168e+056.47458e+04 -3.25422e+052.87134e+02 Pressure (bar) 2.15012e+02 Gromacs 4.6.5 + Gromos54A7: Energies (kJ/mol) G96Bond G96AngleProper Dih. Improper Dih. LJ-14 1.76905e+021.30332e+049.57922e+034.97637e+02 -1.38126e+03 Coulomb-14LJ (SR)LJ (LR) Coulomb (SR) Coulomb (LR) 1.50791e+041.25510e+04 -7.32610e+03 -3.78489e+05 -8.97344e+02 RF excl. PotentialKinetic En. Total EnergyTemperature -5.25414e+04 -3.89718e+056.47631e+04 -3.24955e+052.87211e+02 Pressure (bar) 2.89622e+02 I don't see any significant difference between what 4.0.7 and what 4.6.5 compute. I don't know if the somewhat higher pressure with the 4.6.5/54A7 combination is meaningful. If anyone can shed any light on this, or has ideas for how to debug this further, I'd be most grateful. Best, Lutz [1] A new force field for simulating phosphatidylcholine bilayers.
Re: [gmx-users] Gromos DPPC bilayer: different results for 4.0.7 and 4.6.5
Hi, following Justin's suggestion on redmine (issue #1400), a simulation with Reaction-Field-nec using 4.5.3 is running. It hasn't completed yet but the mean area per lipid around 50 ns is 0.59 mean area 0.60 nm^2. I will update on redmine once the simulation is completed (probably in a couple of days). Best, Patrick Le 16/12/2013 22:30, Justin Lemkul a écrit : On Mon, Dec 16, 2013 at 3:10 PM, Lutz Maibaum lutz.maib...@gmail.comwrote: I am running MD simulations of DPPC bilayers with the Gromos force field, and I am seeing some differences between using Gromacs 4.0.7 and 4.5.5/4.6.5 that I do not understand. It would be great if someone had any insight into what's going on here. When I use the Gromos 53A6-L force field, which is the default 53A6 force field combined with improved lipid parameters that can be downloaded from http://compbio.chemistry.uq.edu.au/~david/research/lipids.htm, I obtain an average area per lipid of 0.627 nm^2, in good agreement with both the paper that describes these new parameters (0.629 nm^2, obtained with Gromacs 3.2.1, Ref. [1]) and another follow-up study (0.631 nm^2 and 0.623 nm^2, Gromacs 4.0.7, Ref. [2]). Now, if I use the same force field and mdp file, and the same initial configuration (which is a pre-equilibrated DPPC bilayer from http://compbio.biosci.uq.edu.au/atb/system_download.py?boxid=32 and randomly generated velocities), but use Gromacs 4.6.5 instread, I get a lower value of about 0.59 nm^2. I also tried the Gromos 54A7 force field, which is included with Gromacs 4.6.5 and that should be identical to 53A6-L (plus it has some other improvements over 53A6 that shouldn't be relevant here), I also get the lower area per lipid of ~0.59 nm^2. If have attached a plot of the area per lipid for these three simulations, each more than 100ns long. This looks to me like 4.6.5 give systematically lower area per lipid than 4.0.7. Running additional simulations with Gromacs 4.5.5 suggest that that also results in the lower area per lipid. Does anyone know why this might be? I have uploaded the relevant files in case that is helpful: http://faculty.washington.edu/maibaum/dppc_comparison/ To see if there are any differences between the energies that 4.0.7 and 4.6.5 compute, I picked a configuration, and used mdrun -rerun with the three different gromacs / force field combinations. Here is what I get for the sample.gro configuration (included in the link above): Gromacs 4.0.7 + Gromos53A6-L: Energies (kJ/mol) G96Bond G96AngleProper Dih. Improper Dih. LJ-14 1.76906e+021.29521e+049.57922e+034.97638e+02 -1.38126e+03 Coulomb-14LJ (SR)LJ (LR) Coulomb (SR) Coulomb (LR) 1.50791e+041.23533e+04 -7.49664e+03 -3.78489e+05 -8.97363e+02 RF excl. PotentialKinetic En. Total EnergyTemperature -5.25414e+04 -3.90168e+056.47185e+04 -3.25449e+052.87013e+02 Pressure (bar) 2.00740e+02 Gromacs 4.6.5 + Gromos53A6-L: Energies (kJ/mol) G96Bond G96AngleProper Dih. Improper Dih. LJ-14 1.76905e+021.29522e+049.57922e+034.97637e+02 -1.38126e+03 Coulomb-14LJ (SR)LJ (LR) Coulomb (SR) Coulomb (LR) 1.50791e+041.23533e+04 -7.49659e+03 -3.78489e+05 -8.97344e+02 RF excl. PotentialKinetic En. Total EnergyTemperature -5.25414e+04 -3.90168e+056.47458e+04 -3.25422e+052.87134e+02 Pressure (bar) 2.15012e+02 Gromacs 4.6.5 + Gromos54A7: Energies (kJ/mol) G96Bond G96AngleProper Dih. Improper Dih. LJ-14 1.76905e+021.30332e+049.57922e+034.97637e+02 -1.38126e+03 Coulomb-14LJ (SR)LJ (LR) Coulomb (SR) Coulomb (LR) 1.50791e+041.25510e+04 -7.32610e+03 -3.78489e+05 -8.97344e+02 RF excl. PotentialKinetic En. Total EnergyTemperature -5.25414e+04 -3.89718e+056.47631e+04 -3.24955e+052.87211e+02 Pressure (bar) 2.89622e+02 I don't see any significant difference between what 4.0.7 and what 4.6.5 compute. I don't know if the somewhat higher pressure with the 4.6.5/54A7 combination is meaningful. If anyone can shed any light on this, or has ideas for how to debug this further, I'd be most grateful. I suspect you're running into the same issue that has been reported here: http://redmine.gromacs.org/issues/1400. Can you check to see if Reaction-Field-nec fixes the issue? We're still waiting on feedback from the Redmine issue. The RF methods changed somewhere along the way, and we need to make sure that the appropriate algorithms are being tested for an apples-to-apples comparison. -Justin -- ___ Patrick FUCHS Dynamique des membranes et trafic intracellulaire Institut Jacques Monod, CNRS UMR 7592, Université Paris Diderot Bâtiment
Re: [gmx-users] Gromos DPPC bilayer: different results for 4.0.7 and 4.6.5
On 16.12.2013 21:10, Lutz Maibaum wrote: I am running MD simulations of DPPC bilayers with the Gromos force field, and I am seeing some differences between using Gromacs 4.0.7 and 4.5.5/4.6.5 that I do not understand. It would be great if someone had any insight into what's going on here. Lutz, from your provided files I can read you use twin-range cutoff RF-electrostatics for 54A7: rcoulomb = rvdw = 1.4 at rlist = 0.8 w/nstlist = 5 are these values intended? I remember a posting of Berk Hess some time ago where he mentioned that the newer Gromacs versions would get a better intergation scheme for this case (I could not find the source). This might be one of the reasons for such discrepancies? BTW - a good check would be a simple box of say 15,625 SPC water - at which water density would the combination of your parameters and the actual algorithm knock in. Regards M. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
