New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Reproducibility of Charmm results in original paper #4
Comments
Hi, I have launched several NAMD simulations of DPPC Here are some results on pure DPPC, in a simulation similar to the one published in I'll give some more results about other simulations settings soon. Pure DPPC (128 lipids and 8000 water) at 323K , 1 atm isotropic NPT , fixed Lx/Ly, beta1 -3 0.06259640550975085 0.1777529934747562 sn2sn1 15 0.09253427959845012 0.06864963229660581 0.11180686965525233 0.06694917560353464 |
Thanks! I did calculate the acyl chain order parameters from this data: http://dx.doi.org/10.5281/zenodo.15549 For area per lipid I get 60.2 Å^2 One detail attracted my attention; you have used 10-12 cutoff and vdwForceSwitching ON. In CHARMM36 paper (Klauda et al. 2010) it is written that: |
It had struck me that this effect can be easily produced by a barostat if misconfigured or if it has some glitches or implementation inaccuracies... The Areas per lipid are different and they might be the source of the difference in the OPs. Right now, I'm running pure POPC (Charmm36) with openMM7.1 and Gromacs5.1.2 to see, whether there are differences -- and I'll probably try exchanging the barostat in Gromacs to see whether there's an effect... |
Indeed, the barostat may play a role. Also, sometimes in gromacs I think you can activate the calculation of the correction to pressure due to the LJ cutoffs. This has no sense in membrane since the assumption of a homogeneous medium does not hold, but maybe there is a spurious default somewhere ? If such correction would be inserted in the barostat, this may change the area/lipid, maybe. Claire Claire LOISON Theoretical Physical Chemistry Group
De : Josef Melcr notifications@github.com It had struck me that this effect can be easily produced by a barostat if misconfigured or if it has some glitches or implementation inaccuracies... The Areas per lipid are different and they might be the source of the difference in the OPs. Right now, I'm running pure POPC (Charmm36) with openMM7.1 and Gromacs5.1.2 to see, whether there are differences -- and I'll probably try exchanging the barostat in Gromacs to see whether there's an effect... You are receiving this because you commented. |
I think that the barostat issue is unlikely. There is one thing to note about previous discussion: |
Exactly, pure POPC running now with Charmm36. the mentioned Gromacs5 setting for C36 is used exactly (btw. standard output of Charmm-gui agrees with it). |
Matti Javanainen also uploaded yesterday some CHARMM36 data of POPC/cholesterol mixtures in Zenodo: http://dx.doi.org/10.5281/zenodo.61649 I ran the analysis for pure POPC data (analysis for systems with cholesterol yet to be done) and added the results in figures: |
Concerning the reproducibility of the original C36 results (Klauda 2010) Despite small differences, I consider that the NPAT SCD results at 64AA /lipid in Figure 3 (bottom) of Klauda's article are very close to ours. Note that our settings are slightly On the opposite, forgetting the vdwForceSwitching and simulating in NPT ps: The impact of LJ cutoffs and timestep is also discussed in section 3.1.4 of this article : |
Nice plot, thank you! |
Where is the experimental data by Douliez actually taken from? By looking the table 1 in http://dx.doi.org/10.1016%2FS0006-3495(95)80350-4, it seems to me that the values at T=323K would be smaller than in the current plot. |
For the (no vdwFS ) simulation, the NAMD parameter vdwForceSwitching is turned off. This parameters changes the way the switch of LJ interactions between the two cutoffs is performed. As far as I understood, the default NAMD is a LJ energy switch between the two cutoffs, but the force is not continuous. In CHARMM, both LJ energy and forces are continuous at both cutoffs. Here is the explanation form NAMD UG : "If both switching and vdwForceSwitching are set to on, then CHARMM force see also http://www.ks.uiuc.edu/Research/namd/mailing_list/namd-l.2011-2012/2609.html The reason for studying this is twice for us :
Since the area per lipid is different with different switching functions, it is interesting so investigate the impact of different area per lipid on results.
they also had vdwForceSwitching? turned off for their optimisation of the parameters. De : Josef Melcr notifications@github.com Nice plot, thank you! You are receiving this because you commented. |
Hi, thanks for your attentive remark. For Douliez et al. I have taken the data directly from the article. Douliez JP, Léonard A, Dufourc EJ. Restatement of order parameters in biomembranes: calculation of C-C bond order parameters from C-D quadrupolar splittings. Biophysical journal. 1995;68(5):1727-1739. doi:10.1016/S0006-3495(95)80350-4. Except for mistakes, I took the value of Table 3 at 50 degres. Another remark considering this article : maybe one could also try to see whether the even-odd effect is obtained in the simulations at various temperatures ? SCC is as easy as the SCD to calculate from simulations. Since the even-odd effect correspond to the derivative of the SCC It seems even more difficult to reproduce than the general trend of SCD. Even more at various temperatures... ? De : ohsOllila notifications@github.com Where is the experimental data by Douliez actually taken from? By looking the table 1 in http://dx.doi.org/10.1016%2FS0006-3495(95)80350-4, it seems to me that the values at T=323K would be smaller than in the current plot. You are receiving this because you commented. |
Ok. My mistake. I was looking the table from wrong lipid. |
Thanks! Could you share also the numerical values so that we can put the Gromacs results in the same plot? The numerical values from Gromacs are here: |
Here is the data for our NAMD-C36 results : the values are followed by their variance, the error estimate is about sqrt(variance)/26. 72 POPC + 2200 watersNPT, 303K, 1atm, Lx/Ly fixedNAMDbeta1 -3 0.07640677098708852 0.1747629068218557 beta2 -3 0.07508321598444838 0.17484878758022854 alpha1 -2 -0.03807238141291847 0.20956378114476149 alpha2 -2 -0.0257067386401897 0.20767275880799413 g31 -1 0.2388885768775786 0.08607526591923258 g32 -1 0.2530631511040515 0.08004482593269276 g2 0 0.20660125026805878 0.09925535984770309 g11 1 0.18019627359407536 0.12318823224085698 g12 1 0.04018602628337821 0.14488592088378308 sn2sn1 2 0.22525701866525205 0.044788047891244534 0.0913459913206161 0.06298796645580687 sn2sn1 3 0.17072383587417905 0.05983885256685818 0.20537798031009685 0.054061884452644586 sn2sn1 4 0.20233165522851435 0.05475314491851802 0.20034825282091387 0.05440962814752503 sn2sn1 5 0.20231784487335827 0.054147959989756224 0.20976385371144332 0.052973475529369055 sn2sn1 6 0.21541941789042587 0.05201180330379692 0.18644196996350526 0.05766554736771268 sn2sn1 7 0.20515775030881273 0.05377655463402382 0.17649410132895962 0.05985269125024771 sn2sn1 8 0.20555136158534454 0.05382470579698103 0.09802525883090739 0.07196403989847149 sn2sn1 9 0.18762370378771778 0.05709155452568148 0.05536004669025257 0.18209310216161567 sn2sn1 10 0.1799601505917554 0.0581497470636555 0.044096302426066006 0.1851044843406602 sn2sn1 11 0.16026646313264598 0.06125112428528533 0.07007914618895074 0.0762407815353155 sn2sn1 12 0.14813536524386606 0.0629708392759908 0.1142848323421472 0.06684090802165985 sn2sn1 13 0.12713875204681396 0.06532279602467393 0.12123498946244926 0.0656729125376915 sn2sn1 14 0.11131470465324475 0.06677228205642875 0.1263016863989463 0.06497669721425287 sn2sn1 15 0.08248560786365253 0.06862334500419658 0.11467703440384294 0.06636734644952041 sn2sn1 16 0.0 0.0 0.10380040692872766 0.06699027824107469 sn2sn1 17 0.0 0.0 0.08121245520460374 0.06855480483810632 Here is the CHARMM simulation data I extracted from the paper by Klauda et al. using the figure 9 and http://arohatgi.info/WebPlotDigitizer/app/ POPC from the C36 paper,DOI: 10.1021/jp101759qfigure 9 , upper picturesgreen circles (sn-2)
11.0 0.07207 12.0 0.1145 13.0 0.1200 14.0 0.1214 15.0 0.1102 16.0 0.1008 17.0 0.08027? POPC from the C36 paper,DOI: 10.1021/jp101759qfigure 9 , upper picturesblue diamonds (sn-1)
Claire LOISON Theoretical Physical Chemistry Group
De : ohsOllila notifications@github.com Thanks! Could you share also the numerical values so that we can put the Gromacs results in the same plot? The numerical values from Gromacs are here: You are receiving this because you commented. |
Great, I'll bring up sim-data from openMM and Gmx512 later today. |
Updated the figure - added data from Claire's post (NAMD simulation and Charmm36 paper data), only inverted the values.
According to the obtained values, there is close to no difference in between current simulation softwares -- no secret hidden magic in obscure simulation settings and VdW cutoffs. Once the simulations use same setting, the results begin to match within some (in)accuracy -- as they should. Howgh. I'll upload the obtained trajectories to Zenodo for all of you to check and use. In addition, we will have a set of trajectories of the same system, which were generated by 3 different simulation engines and results that match well. |
The curves look more similar than the ones in the Charmm-GUI paper (DOI: 10.1021/acs.jctc.5b00935), where the problem with pre-5.0 era Gromacs simulations with Charmm36 is also documented. I wonder what's the reason behind the better agreement observed here. Any ideas? I provided links to some Charmm36 POPC trajectories with and without cholesterol in the blog (the ones Samuli refers to above); perhaps, using this data, we should check whether Gromacs 5.0 and Gromacs 5.1 give identical results. Anyway, we should clearly stop using any Charmm36 data produced with pre-5.0 era Gromacs. |
They claim they used Gromacs 5.0 in the Charmm-GUI paper. In release notes to version 5.0.2 there is:
So it also depends what settings was used in the paper (non-GPU/GPU). I didn't find it in the Charmm-GUI paper.
The differnce between force/potential-based switching can also be expected, since they change the potential form of the classical model -> interactions are different -> ensembles are different. And this difference is pronounced in aggregates, in which VdW dispersion plays a large role -- lipid bilayers.
Yes, use the best available software for the task =]
My guess is that it would be the easiest for you to compare -- significant differences should be visible after short runs. Also depends what 5.0.x version you used and whether it included this tunepme - or some other bug =] |
Well, now we have both Gromacs 5.1.2 and Gromacs 5.0.7 trajectories available for comparison, Samuli could you take a look at this? My arguments for staying in Gromacs (>5.0) are:
|
I will take a closer look at this, but I am currently pretty busy with applications. Based on this discussion it seems to me that with Gromacs 5 with correct cut-off treatment we get results in reasonable agreement with NAMD and literature results. For NMRlipids III this is enough and we already have the simulations of POPC cholesterol mixtures with Gromacs 5 waiting in Zenodo for the analysis. About 5.1.2 vs. 5.0.7 comparison, we indeed have both trajectories but I think that they are with little bit different temperatures. Anyway, based on general picture it seems that the results are consistent enough for NMRlipids III. Just for curiosity, how large is the performace difference between NAMD, openMM and Gromacs? |
On our desktops Gromacs is ~3–5(!) times faster than NAMD depending on whether GPU is used. For this test we used the files as provided by Charmm-GUI. This drastic difference in performance is also evident here: https://research.csc.fi/sisu-scalability-tests I have no experience with openMM. However in my opinion, due to the lack of the correct LJ modifiers (this is stated in the Charmm-GUI paper), openMM should not really be considered for the Charmm36 simulations, just like Gromacs versions <5.0. |
Hi, The order parameters are too high for the GROMACS 4.5 simulations due to the use of potential rather than force switching. This can be seen in both the above NAMD simulations where the no vdwFS simulation uses potential switching (it is probably also worth noting that this simulation is in agreement with a previous NAMD simulation I did with this setting; see table S3 of http://pubs.acs.org/doi/suppl/10.1021/ct3003157) and also, as already mentioned, seen in the simulations reported by Lee et al. using both force and potential switching in GROMACS 5. If you actually wanted to get force switching in GROMACS 4.5 you would have needed to have used a shift cut-off. I've had a look at the mdp's for DPPC/POPC and they don't use this. I didn't realise when we published our work that a shift cut-off could achieve force switching in GROMACS. As for the order parameters in GROMACS 5, it is very interesting that Joseph's simulations show close agreement between GROMACS 5 and other softwares in contrast to Lee et al. who claimed that the membranes were consistently slightly too ordered. I guess we could also run a simulation using GROMACS 5.0 as per Lee et al. (with/without GPU) to test Joseph's hypothesis (and seeing if the cut-off is changed). That said, I imagine it would probably be easier to contact the corresponding author of the paper and ask if they can look at the log files to see if this was the case (and to also confirm the details, such as did they use 5.0 or another of the 5.0.x series). In any case, I've just set running some DPPC/POPC membranes in 5.0.6 to confirm what Joseph's simulations show, so we can at least be sure that the problem is something with the Lee et al. simulations. Cheers Tom |
Tom, if you did not notice, I already posted some POPC simulations run with Gromacs 5.0.7 to the blog some time ago, see: https://doi.org/10.5281/zenodo.61649 |
Great, thanks. I hadn't noticed. I've actually found a GROMACS 5.1.2 simulation I had for POPC too (500 ns, 0.8/1.2 nm cut-off). I'm just having a look at the order parameters, etc. now. |
It is also good to note that the numerical results from Matti's simulation are here: |
So despite the simulations of both Matti and Joseph having slightly more disordered membranes (making GROMACS in closer agreement with other softwares), the simulation I mentioned above performed using GROMACS 5.1.2 is in better agreement with the results of Lee et al. The differences aren't huge but then again the differences reported by Lee et al (using the force switching) weren't huge either. I should also mention again that this simulation is 500 ns in length and is using a 0.8/1.2 force switch cut-off (which should likely make the membrane slightly more disordered than using 1.0/1.2 as was used in the simulations of Joseph/Matti). The area per lipid from the simulation is shown in the plot below (the average is 0.640 nm^2 using the whole trajectory), matching closely to the Lee et al. GROMACS/POPC results of 0.641 (albeit with a 1.0/1.2 nm cut-off). Below are the order parameters for the sn-1 chain (as calculated using the MATCH scripts and I also checked this with a VMD tcl script (calc_op.tcl) that averages over the hydrogens; Claire I guess you probably used this for your NAMD sims?): label Order_Parameter_1 Order_Parameter_2 And plotted against Joseph's values reported from his GROMACS simulation: And also plotted against Claire's NAMD simulations (for both sn-1 and sn-2 chains here): While the differences aren't huge, there are consistent differences in the order parameters. The differences between GROMACS and NAMD are also of a similar order and nature to those reported by Lee et al. (the POPC order parameter plots are in the SI of the paper, although it is hard to work out the exact numbers given the large scales and symbols used in these plots), who also said the increases were consistent but not significant. As to why I am seeing differences to Joseph and Matti, I am not sure. I am running a series of simulations/tests to examine starting structure, cut-off, where the force field/lipid parameters are from and simulation length so as to hopefully figure this out. I should also say there are repeats of everything running too, so as we can hopefully be confident with the results. I will report back and can also share these simulations once completed. The differences really aren't huge anyway but it still would be nice to work out why this consistent small increase in order of the membranes is (well, if it is!) happening. That said, I really don't think it is anything to get overly concerned about. |
I'm afraid this is not generally valid. In addition, one should use the cutoff scheme belonging to the Force field, otherwise the potential form is being changed (and hence the results).
Bearing in mind that you have a modified potential form for LJ interaction (through different switching), the agreement of the reported curves is OK. Notice that in Fig 2 in charmm-gui paper the difference is larger than the one shown here. Could you do the simulation with the same settings as we did? The results must match within convergence-precision. @mattijavanainen (choice of software)
|
Hi Josef (PS, sorry for the incorrect spelling of your name before!). Yes, good to chat again, although my memories of the conference are somewhat hazy given the free booze on offer! As for the cut-off settings with CHARMM, I probably should have been a bit more explicit in my previous message. I completely agree with you that random changes to the cut-off for a force field should not be done. Indeed some of my previous work was the first (I think) to highlight how both the change from force to potential switch and also which variant of TIP3P water model you use, can both massively impact upon the membrane properties with CHARMM36. Indeed, both of these (and other) settings should be classified as much as part of the force field as any of the atomtypes, etc. However, the one thing with the CHARMM force fields that has been a little more variable over time is the point at which the van der Waals interactions are switched off. In the older CHARMM force fields this value was 0.8 nm and indeed this was the value primarily used in the original CHARMM36 lipid force field parameterisation of Klauda et al. (http://pubs.acs.org/doi/abs/10.1021/jp101759q). However, Klauda also tested using a switching point of 1.1 nm due to the fact that this had been frequently used in the NAMD software. This resulted in a slightly more ordered membrane (e.g. lower APL, higher order parameters, thicker membrane, etc.) as can be understood through inclusion of some more of the longer range, attractive, van der Waals interactions. Therefore, it could be argued that either one of these values could be valid for CHARMM36 membrane simulations but based upon the membrane properties the most reasonable value to use would be 0.8 nm. Moving forwards from the original parameterisation of the CHARMM36 lipids, I believe the newer CHARMM36 protein force field was subsequently parameterised with a point to begin switching off the van der Waals interactions of 1.0 nm (IIRC this is true for other more recent parts of the CHARMM force field too but please don't take my word as gospel on this without looking). This is why in the paper of Lee et al., they tested both the 0.8 and 1.0 nm switching points with a complete cut-off at 1.2 nm for the van der Waals. While the former actually in general provides better membrane properties (and, given the fact that it was used in the vast majority of the original CHARMM36 lipid parameterisation, should probably be used for pure CHARMM36 membrane simulations IMHO) they recommended the latter (i.e. 1.0/1.2) to be consistent with other parts of the CHARMM36 force field parameters. As to why they also tested 0.8/1.0 in GROMACS, I am not completely sure. Perhaps to try and optimise membrane properties for pure membrane simulations I guess. Anyway, so when I previously talked about simulations using different cut-off's I was simply referring to simulations with 0.8/1.2 and 1.0/1.2 nm for the van der Waals interactions. Both of these are valid settings and, as I mentioned above, the former (IMHO) is more appropriate for pure membrane systems. All that said, yes I also agree that simulations reproducing what has already been done by yourself and Matti is the easiest way to work out what is going on. This is why the vast majority of the simulations I have set going are using the 1.0/1.2 nm cut-off and indeed I have 4 simulations running using the 200 ns structure you provided as input (my runs are cpu only). However, given what is known from the different settings (0.8/1.2 versus 1.0/1.2) and the fact the GROMACS 5.1.2 simulation I had to hand used the 0.8/1.2 settings that should have resulted in lower order parameters than your simulations, not higher, I am not convinced that there is agreement between our simulations. The simulations that are now running will show whether this is the case or not! Cheers Tom |
We should perhaps try to organize the results with different software, versions, cut offs, temperatures etc. in a table format as I find it quite hard to follow the results both here and in the blog. @jmelcr: Good to hear that the proper LJ options are available, and that openMM performs so well on a GPU. Have you tried it on a CPU – or more importantly – on a supercomputer/cluster? How is the performance there? Actually I am not sure whether all MD codes are expected to provide the exact same result as these depend heavily on implementation details, which are not bugs. Therefore the obvious question is whether a force field should be used on multiple MD codes at all, or whether one should stick to the one used for the parameterization. Or in other words, should the implementation of algorithms be considered to be a fixed part of the force field? Unfortunately our exercise with Charmm36 seems to support the former view. |
FYI, I'm running now POPC+30%Chol (using files from Matti) to see, whether we get same/different results. #HashTagDoubleCheck =] edit: results regarding DOPs match well. |
I plotted the results from this discussion into a figure below. It seems to me now that with CHARMM, NAMD, Gromacs 5.x and openMM we are able to get results for sn-1 chain which are not quite the same but still pretty close to each others. It seems that this accuracy is currently enough for NMRlipids III. Thus understanding the fine details and origins of these differences is not currently necessary for this project, however, it may relevant for some other reasons. The results from all program packages give slightly larger order parameters from CHARMM36 than the experimental values from Ferreira et al. The results from Gromacs 4.5 are larger than from other programs. Based on this discussion and the literature origin of the reason for this seems quite clear. In conclusion, we should not use simulations from Gromacs 4.5 for CHARMM36 anymore. There is, however, a new problem: The results for sn-2 chain from Gromacs 5.0. simulation (http://doi.org/10.5281/zenodo.61649) are significantly different than the reported results from other simulations (see figure below). Now we need to find out what is the reason for this and if we can use already available simulations from Gromacs 5.0 with cholesterol in NMRlipids III. I will double check the order parameter calculation. The simulation running by Melcr is useful for this as well. @TomPiggot |
The temperature was 300 K in this one simulation. The sn-2 order parameters are: label Order_Parameter_1 Order_Parameter_2 The series of simulations I mentioned above have just recently finished so I will report back with the analysis soon on these. |
…te the results in correct temperatures
I added also the results by Tom Piggot in the figure: |
* Add figure based on results at NMRLipids#4 * Add results by Tom Piggot from NMRLipids#4
I opened a new issue about the CHARMM36 results from Gromacs 5.0 (http://dx.doi.org/10.5281/zenodo.61649, see comment #4 (comment) above), because this does not seem to be a version issue: #9 |
Results for POC+30%Chol presented in issue #9. |
Firstly apologies for my slowness on this, I've been rushed off my feet with work. The series of simulations I mentioned above (GROMACS 5.0.6, 4 x 500 ns POPC started from the 200 ns structure of Josef and also 2 x 500 ns POPC from a different starting structure and different CHARMM/GROMACS parameter conversion I already had) are in pretty close agreement with the reported results of Lee et al. in the CHARMM-GUI paper, i.e. that GROMACS does have a slightly more ordered membrane (lower APL, higher order parameters) than those reported for other softwares. Below are the APL stats (first 100 ns not included): From Josef's structure (simulation version, mean and standard deviation as outputted by g_analyze) From my structure And all the results from the order parameters (again for the final 400 ns of the simulations), firstly the simulations started from Josef's 200 ns structure (sn-1 followed by sn-2) v1 2 -0.069326 -0.10744 v2 2 -0.0702182 -0.108638 v3 2 -0.070868 -0.108084 v4 2 -0.0742355 -0.107933 And also from the simulations of the structure and force field parameters I already had: v1 2 -0.0743307 -0.1071 v2 2 -0.0727035 -0.107453 These results reported above (sorry for all the numbers) are from simulations performed using the 1.0/1.2 nm vdW cut-off's. The same set of simulations performed with 0.8/1.2 nm vdW cut-off's demonstrate how changing from a 1.0 nm to 0.8 nm switching point (to match the original values of Klauda et al.) can slightly decrease the order: From Josef's structure From my structure Hopefully this puts some of the doubts around this issue to bed. Yes, there is a slight differences in the order of POPC membranes when using GROMACS 5 compared to reported values when simulating in other MD packages (e.g. NAMD) but they are fairly minor and also they can be reproduced with different membrane structures and different origins of the GROMACS force field parameters. As to why they occur, I am still unsure. I also think this demonstrates my point that, if you are simulating a membrane system only, you should use the 0.8/1.2 nm cut-off scheme for the van der Waals interactions. I'll upload these simulations and analysis to Zenodo when I get the chance. Finally it's worth noting that I ran some DPPC simulations too that also agree with these points above on matching with Lee et al. and that 0.8/1.2 nm switching is better in GROMACS (I didn't include these results in this post as there are quite enough numbers already!). |
The POPC simulation files are on Zenodo now: https://doi.org/10.5281/zenodo.164206 The trajectories have been processed with trjconv -skip 10 to make the uploads a reasonable size but the analysis reported above was all on the original trajectories. |
I think that these issues are nowadays known and discussed in several publications after this discussion. However, I am not really sure if origins of all differences are understood. Anyway, I will close this issue now. |
The order parameters for POPC acyl chain with CHARMM36 model are larger than in experiments in figure https://github.com/NMRLipids/NmrLipidsCholXray/blob/master/FIGS/OrderParametersCHOL.jpg.
The results look similar as reported by Piggot et al. [Fig. 5 in dx.doi.org/10.1021/ct3003157 | J. Chem. Theory Comput. 2012, 8, 4593−4609].
However, in the original CHARMM36 paper a better agreement with experiments is reported [Fig. 9 inhttp://dx.doi.org/10.1021/jp101759q, J. Phys. Chem. B, Vol. 114, No. 23, 2010]. This is similar difference as in area per molecule.
I think that it is necessary to find out what is the reason for this. Does CHARMM36 really have too condensed bilayer, or if there is something wrong in the Gromacs conversion?
The text was updated successfully, but these errors were encountered: