[gmx-users] Comparing Charmm 27 energies from GROMACS and NAMD

Reza resal81 at gmail.com
Mon May 13 21:57:24 CEST 2013


So I think I figured out what was causing the discrepancy of Charmm27 energies between gromacs and NAMD. It appears that it's related to the gromacs version: energies from 4.5.7 and NAMD match very well while 4.6.1 gives energies that are different from both 4.5.7 and NAMD.

Here are the results form my tests on POPC membrane:

===============================
  BOND        ANGLE        DIH         IMP       COUL        LJ          POT
gromacs-4.5.7/gromacs.log
  1539.120    3111.902    1250.540      16.284   -1705.306   -1219.345    2993.188

gromacs-4.6.1/gromacs.log
  1539.120    3111.902    1250.540      16.284   -5997.992     350.084     269.945

namd-2.9/namd.out
  1539.116    3111.920    1250.542      16.284   -1705.307   -1219.343    2993.211
===============================


I think this could be due to either 1) the simulation setup should be different between 4.5.7 and 4.6.1 when running single point calculations, or 2) my 4.6.1 compilation is not correct (although all regression tests passed), or 3) there is a bug somewhere.

I'll attach a script that automates the comparison between GROMACS and NAMD in a separate post (it will make this post > 50kB). If you run it and got similar/different results as mine, I'll be happy to know. Also any comment on the simulations setup is highly appreciated.

Reza



On May 3, 2013, at 11:57 AM, Reza <resal81 at gmail.com> wrote:

> 
> Actually pdb which is similarly doesn't have velocity information.  In this case however I'm mainly interested to see how the potential terms compare between the two packages with Charmm27 ff.
> 
> Reza
> 
> 
>> 
>> 
>> On 5/3/13 11:33 AM, Reza wrote:
>>> Thanks Mark!
>>> 
>>>> No, they weren't. See
>>>> http://www.gromacs.org/Documentation/How-tos/Single-Point_Energy
>>>> You cannot hope to reproduce accurately the energy of a configuration if
>>>> you let the coordinates be manipulated.
>>>> 
>>> 
>>> 
>>> The "rerun" option is interesting - in my case however the potential energy terms stayed identical but the kinetic term became zero.
>>> 
>> 
>> If you are using an .xtc file for the rerun, this makes sense since the .xtc does not store velocities and hence kinetic energy cannot be calculated.
>> 
>> -Justin
>> 
>>>> 
>>>> Something you think is equivalent is not :-) Move to testing a system with
>>>> two lipids. Inspect all the logfile outputs very carefully for clues.
>>> 
>>> 
>>> I totally agree :) So far I found out that for "no cut-off" simulation in Gromacs, rather that specifying a large cut-off, it needs rlist=rvdw=rcoulomb=0 and pbc=no along with ns_type=simple and nstlist=0 (according to the manual). I am running various tests and will update if I find out what is causing the discrepancy.
>>> 
>>> Reza
>>> 
>>> On May 1, 2013, at 5:46 PM, Mark Abraham <mark.j.abraham at gmail.com> wrote:
>>> 
>>>> On Wed, May 1, 2013 at 2:32 PM, Reza Salari <resal81 at gmail.com> wrote:
>>>> 
>>>>> Hi Justin,
>>>>> 
>>>>> I actually did :) but it ended up being bigger than 50 kb so it needed
>>>>> moderator approval to show up. I was hoping it would've been released by
>>>>> now. I'll attach a the details below.
>>>>> 
>>>>> Any help/hint is highly appreciated.
>>>>> 
>>>>> Reza
>>>>> 
>>>>> Details:
>>>>> 
>>>>> *1)* Both systems were prepared using VMD "membrane" package and then
>>>>> waters were removed. I used Gromacs 4.6.1 and NAMD 2.9.
>>>>> 
>>>>> *2)* Simulations were run in vacuum as a single-point energy calculations
>>>>> (0 step).
>>>> 
>>>> 
>>>> No, they weren't. See
>>>> http://www.gromacs.org/Documentation/How-tos/Single-Point_Energy
>>>> You cannot hope to reproduce accurately the energy of a configuration if
>>>> you let the coordinates be manipulated.
>>>> 
>>>> PME was not used.
>>>>> 
>>>>> *3) *For Gromacs runs, pdb2gmx was used to prepare the system and the
>>>>> output was saved as the pdb format. The mdp file:
>>>>> 
>>>>> integrator    = md
>>>>> nsteps        = 0
>>>>> nstlog      = 1
>>>>> nstlist        = 1
>>>>> ns_type        = grid
>>>>> rlist        = 100.0
>>>>> coulombtype    = cut-off
>>>>> rcoulomb    = 100.0
>>>>> rvdw        = 100.0
>>>>> pbc         = no
>>>>> 
>>>>> *4) *NAMD input file:
>>>>> 
>>>>> structure          ../0_prep/memb_nowat.psf
>>>>> paratypecharmm  on
>>>>> parameters        par_all27_prot_lipid.prm
>>>>> exclude         scaled1-4
>>>>> 1-4scaling         1.0
>>>>> switching     off
>>>>> switchdist     8
>>>>> cutoff         1000
>>>>> pairlistdist     1000
>>>>> margin             0
>>>>> timestep         1.0
>>>>> outputenergies     1
>>>>> outputtiming     1
>>>>> binaryoutput     no
>>>>> coordinates     ../0_prep/memb_nowat.pdb
>>>>> outputname         out
>>>>> dcdfreq         10
>>>>> temperature     300
>>>>> run 0
>>>>> 
>>>>> *5)* Energies:
>>>>> 
>>>>> For Single POPC  (kcal/mol)
>>>>> 
>>>>> 
>>>>> Gromacs NAMD Diff
>>>>> 
>>>>> Bond 43.0022 43.0015 -0.0007
>>>>> Angle 80.6568 80.6571 0.0003
>>>>> Dih 29.8083 29.8083 0.0000
>>>>> Imp 0.8452 0.8452 0.0000
>>>>> 
>>>>> Coul -17.2983 -17.2983 0.0000
>>>>> LJ -7.0798 -7.0798 0.0000
>>>>> 
>>>>> Pot 129.9343 129.9340 -0.0003
>>>>> 
>>>>> 
>>>> The intra-molecule terms look fine. Since this is a lipid, there are
>>>> non-bonded interactions that are intra-molecule, so the non-bondeds also
>>>> seem fine.
>>>> 
>>>> For POPC Memb (kcal/mol)
>>>>> 
>>>>> Gromacs NAMD Diff
>>>>> 
>>>>> Bond 1539.1181 1539.1162 -0.0019
>>>>> Angle 3111.9264 3111.9197 -0.0067
>>>>> Dih 1250.5425 1250.5421 -0.0004
>>>>> Imp 16.2837 16.2837 0.0000
>>>>> 
>>>>> Coul -1837.8585 -1705.3075 132.5510
>>>>> LJ -995.0311   -1219.3432 -224.3121
>>>>> 
>>>>> Pot 3084.9904 2993.2110 -91.7794
>>>>> 
>>>> 
>>>> Something you think is equivalent is not :-) Move to testing a system with
>>>> two lipids. Inspect all the logfile outputs very carefully for clues.
>>>> 
>>>> Mark
>>>> 
>>>> 
>>>>> 
>>>>> 
>>>>> 
>>>>> 
>>>>> 
>>>>> On Tue, Apr 30, 2013 at 8:33 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>>>>> 
>>>>>> 
>>>>>> 
>>>>>> On 4/30/13 4:19 PM, Reza Salari wrote:
>>>>>> 
>>>>>>> Hi
>>>>>>> 
>>>>>>> I have set up two small systems, one with a single POPC lipid, and
>>>>> another
>>>>>>> system with 23 POPC's arranged as a lipid bilayer. I am hoping to do a
>>>>>>> similar comparison as in Table 1 in the Par Bjelkmar et al paper
>>>>> (porting
>>>>>>> charmm ff to gromacs) for my systems. My main question is that for the
>>>>>>> single POPC, all the potential energy terms match very well, but for the
>>>>>>> membrane system the non-bonding terms differ significantly.
>>>>>>> 
>>>>>>> I am providing the full details below and greatly appreciate any hint
>>>>> for
>>>>>>> better comparison of the energies.
>>>>>>> 
>>>>>>> 
>>>>>>> Thanks,
>>>>>>> Reza Salari
>>>>>>> 
>>>>>>> Details:
>>>>>>> 
>>>>>>> 1) Both systems were prepared using VMD "membrane" package:
>>>>>>> 
>>>>>>> 2)
>>>>>>> 
>>>>>>> 
>>>>>> It appears you hit "send" too early.  Please provide the entirety of the
>>>>>> details you intended.  Complete .mdp files and actual quantification of
>>>>> the
>>>>>> differences you are observing are also very important.
>>>>>> 
>>>>>> -Justin
>>>>>> 
>>>>>> --
>>>>>> ==============================**==========
>>>>>> 
>>>>>> Justin A. Lemkul, Ph.D.
>>>>>> Research Scientist
>>>>>> Department of Biochemistry
>>>>>> Virginia Tech
>>>>>> Blacksburg, VA
>>>>>> jalemkul[at]vt.edu | (540) 231-9080
>>>>>> http://www.bevanlab.biochem.**vt.edu/Pages/Personal/justin<
>>>>> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin>
>>>>>> 
>>>>>> ==============================**==========
>>>>>> --
>>>>>> gmx-users mailing list    gmx-users at gromacs.org
>>>>>> http://lists.gromacs.org/**mailman/listinfo/gmx-users<
>>>>> http://lists.gromacs.org/mailman/listinfo/gmx-users>
>>>>>> * Please search the archive at http://www.gromacs.org/**
>>>>>> Support/Mailing_Lists/Search<
>>>>> http://www.gromacs.org/Support/Mailing_Lists/Search>before posting!
>>>>>> * Please don't post (un)subscribe requests to the list. Use the www
>>>>>> interface or send it to gmx-users-request at gromacs.org.
>>>>>> * Can't post? Read http://www.gromacs.org/**Support/Mailing_Lists<
>>>>> http://www.gromacs.org/Support/Mailing_Lists>
>>>>>> 
>>>>> --
>>>>> gmx-users mailing list    gmx-users at gromacs.org
>>>>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>>>>> * Please search the archive at
>>>>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>>>>> * Please don't post (un)subscribe requests to the list. Use the
>>>>> www interface or send it to gmx-users-request at gromacs.org.
>>>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>>> 
>>>> --
>>>> gmx-users mailing list    gmx-users at gromacs.org
>>>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>>>> * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>>>> * Please don't post (un)subscribe requests to the list. Use the
>>>> www interface or send it to gmx-users-request at gromacs.org.
>>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>> 
>> 
>> -- 
>> ========================================
>> 
>> Justin A. Lemkul, Ph.D.
>> Research Scientist
>> Department of Biochemistry
>> Virginia Tech
>> Blacksburg, VA
>> jalemkul[at]vt.edu | (540) 231-9080
>> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
>> 
>> ========================================
>> -- 
>> gmx-users mailing list    gmx-users at gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>> * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>> * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-request at gromacs.org.
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> 




More information about the gromacs.org_gmx-users mailing list