Dear All,
Thanks for your reply and sorry for being late due to the time difference :)
Ross, in gromacs i perform a rerun on the gro file. So it is just energy calculation and it doesn't move the system. In amber, I perform one step and I look at the energy of step zero.
I computed the BOND and ANGLE energies. But it is only increasing my confusion ... Look: I have a simple ace-ala-nme system. This is the bond information both for amber and gromacs ... but I wrote the gromacs topology since it is easier to look:
[ bonds ]
; ai aj funct r k
2 3 1 1.0900e-01 2.8451e+05
2 4 1 1.0900e-01 2.8451e+05
1 2 1 1.0900e-01 2.8451e+05
11 12 1 1.0900e-01 2.8451e+05
11 13 1 1.0900e-01 2.8451e+05
11 14 1 1.0900e-01 2.8451e+05
9 10 1 1.0900e-01 2.8451e+05
7 8 1 1.0100e-01 3.6317e+05
19 20 1 1.0900e-01 2.8451e+05
19 21 1 1.0900e-01 2.8451e+05
19 22 1 1.0900e-01 2.8451e+05
17 18 1 1.0100e-01 3.6317e+05
5 6 1 1.2290e-01 4.7698e+05
5 7 1 1.3350e-01 4.1003e+05
2 5 1 1.5220e-01 2.6527e+05
15 16 1 1.2290e-01 4.7698e+05
15 17 1 1.3350e-01 4.1003e+05
9 11 1 1.5260e-01 2.5941e+05
9 15 1 1.5220e-01 2.6527e+05
7 9 1 1.4490e-01 2.8200e+05
17 19 1 1.4490e-01 2.8200e+05
[ angles ]
; ai aj ak funct theta cth
5 7 8 1 1.2000e+02 4.1840e+02
4 2 5 1 1.0950e+02 4.1840e+02
3 2 4 1 1.0950e+02 2.9288e+02
3 2 5 1 1.0950e+02 4.1840e+02
1 2 3 1 1.0950e+02 2.9288e+02
1 2 4 1 1.0950e+02 2.9288e+02
1 2 5 1 1.0950e+02 4.1840e+02
15 17 18 1 1.2000e+02 4.1840e+02
13 11 14 1 1.0950e+02 2.9288e+02
12 11 13 1 1.0950e+02 2.9288e+02
12 11 14 1 1.0950e+02 2.9288e+02
10 9 11 1 1.0950e+02 4.1840e+02
10 9 15 1 1.0950e+02 4.1840e+02
9 11 12 1 1.0950e+02 4.1840e+02
9 11 13 1 1.0950e+02 4.1840e+02
9 11 14 1 1.0950e+02 4.1840e+02
8 7 9 1 1.1804e+02 4.1840e+02
7 9 10 1 1.0950e+02 4.1840e+02
21 19 22 1 1.0950e+02 2.9288e+02
20 19 21 1 1.0950e+02 2.9288e+02
20 19 22 1 1.0950e+02 2.9288e+02
18 17 19 1 1.1804e+02 4.1840e+02
17 19 20 1 1.0950e+02 4.1840e+02
17 19 21 1 1.0950e+02 4.1840e+02
17 19 22 1 1.0950e+02 4.1840e+02
6 5 7 1 1.2290e+02 6.6944e+02
5 7 9 1 1.2190e+02 4.1840e+02
2 5 6 1 1.2040e+02 6.6944e+02
2 5 7 1 1.1660e+02 5.8576e+02
16 15 17 1 1.2290e+02 6.6944e+02
15 17 19 1 1.2190e+02 4.1840e+02
11 9 15 1 1.1110e+02 5.2718e+02
9 15 16 1 1.2040e+02 6.6944e+02
9 15 17 1 1.1660e+02 5.8576e+02
7 9 11 1 1.0970e+02 6.6944e+02
7 9 15 1 1.1010e+02 5.2718e+02
The units in gromacs are nm and kj/mol ... I computed the Bond energy which is:
BOND energy= Sum (ki*(ri*ri)) =10040.876 kjoul/mol = 2402.12 kcal/mol ...
ANGLE energy= Sum (ki(theta_i*theta_i)) = 204462689.8 kj/mol =48867755.6 kcal/mol
mmmmmmm! I don't know?
Looking forward to hear back from you and greeting from Boston :)
Golshan
________________________________
From: Scott Le Grand <varelse2005.gmail.com>
To: AMBER Mailing List <amber.ambermd.org>
Cc: Golshan Hejazi <golshan.hejazi.yahoo.com>
Sent: Thursday, August 8, 2013 7:36 PM
Subject: Re: [AMBER] BOND, ANGLE and TORSION energies discrepancy between amber and gromacs
It's very interesting to me that one never sees comparative energies in
various papers comparing themselves to pmemd.cuda, only comparative
performance, and even then, the benchmarks have frequently been tweaked to
make pmemd.cuda look bad. For the record, DPDP mode in pmemd.cuda matches
CPU PMEMD energies to many decimal places...
On Thu, Aug 8, 2013 at 12:48 PM, Ross Walker <ross.rosswalker.co.uk> wrote:
> Hi Golshan,
>
> You are correct the internal terms should be identical - there is no
> excuse for them not to be. I would caution that Gromacs uses KJ/mol and
> not Kcal but the difference here is not 4.184 so it isn't a unit
> conversion issue. It is also too big to be simple round off. Given the
> structure is small you should be able to calculate it yourself by hand.
> Make sure you are really doing just a single point energy calculation -
> I.e. step 0 of a minimization for example and that you are not using shake
> or anything else that might mess with the potential.
>
> Then on a piece of paper you can iterate through with the coordinates you
> have and calculate what all the bonds, angles and dihedrals should sum to.
> It will likely take a couple of hours but you'll be able to convince
> yourself what the correct value is.
>
> This will tell you which code is correct. (I suspect AMBER). Suffice to
> say the fact that Gromacs doesn't get this identical would cause me to be
> VERY VERY suspicious. If one cannot get bonds, angles and dihedrals to
> match then there is no hope for the code correctly implementing the force
> field.
>
> Note for the implementation of the Charmm force field in AMBER we get
> these terms matching charmm to 10^-14 kcal/mol so there is no excuse.
>
> http://www.rosswalker.co.uk/papers/Int_J_Quant_Chem_2009_109_15_3767-3772_C
> hamber.pdf
>
>
> All the best
> Ross
>
>
>
> On 8/8/13 12:26 PM, "Golshan Hejazi" <golshan.hejazi.yahoo.com> wrote:
>
> >Hello everyone.
> >
> >I am doing some simple tests as following: I have an amber crd and top
> >file. I perform an energy evaluation and I estimate the value of BOND,
> >ANGLE and TORSION energies for a simple system (ace-ala-nme).
> >
> >Then I use amb2gmx.pl to convert the crd and top file to gromacs format
> >and again I compute the BOND, ANGLE and TORSION energies of the same
> >structure ...
> >
> >I expect to have identical values between gromacs and amber. However, the
> >energis are not comparable.
> >
> > amber ... BOND= 0.0206 (kcal/mol) ANGLE=0.3620 (kcal/mol)
> >TORSION=8.1071 (kcal/mol) gromacs .... BOND=0.14044 (kcal/mol)
> >ANGLE=0.3780 (kcal/mol) TORSION=9.74190 (kcal/mol)
> >
> >I checked how amb2gmx converted the values of amber to gmx format. It
> >seems to me that everything is correct. There is a factor of two
> >difference between the force constants in amber and gromacs which I think
> >in the energy calculation term ... it is already taken care of!
> >
> >So I see nothing wrong but ... energies dont match!
> >Could you please help me to understan this
> >
> >Thanks
> >Golshan
> >_______________________________________________
> >AMBER mailing list
> >AMBER.ambermd.org
> >http://lists.ambermd.org/mailman/listinfo/amber
>
>
>
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Fri Aug 09 2013 - 07:00:03 PDT
