Re: [AMBER] BOND, ANGLE and TORSION energies discrepancy between amber and gromacs

From: Ross Walker <ross.rosswalker.co.uk>
Date: Fri, 09 Aug 2013 14:54:37 -0700

Ok, here's what I get for AMBER.

1) Using RDPARM to get the bond info ($AMBERHOME/bin/rdparm
ALA-amber.top bonds)

    Bond Kb Req atom names (numbers)
       1: 434.00 1.010 :1.N :1.H (1,2)
       2: 340.00 1.090 :1.CA :1.HA (3,4)
       3: 340.00 1.090 :1.CB :1.HB1 (5,6)
       4: 340.00 1.090 :1.CB :1.HB2 (5,7)
       5: 340.00 1.090 :1.CB :1.HB3 (5,8)

       6: 337.00 1.449 :1.N :1.CA (1,3)
       7: 310.00 1.526 :1.CA :1.CB (3,5)
       8: 317.00 1.522 :1.CA :1.C (3,9)
       9: 570.00 1.229 :1.C :1.O (9,10)




Then for the structure you provided going through each bond in turn:

sqrt([Xi-Xj]^2+[Yi-Yj]^2+[Zi-Zj]^2)

1 - 3.3257700 1.5479090 -0.0000016 3.9094070 0.7236110 -0.0000027
= 1.010
2 - 3.9700480 2.8457950 -0.0000001 3.6716630 3.4001290 -0.8898200
= 1.090
3 - 3.5769650 3.6538380 1.2321430 3.8774840 3.1157950 2.1311970 = 1.090
4 - 3.5769650 3.6538380 1.2321430 4.0750590 4.6230170 1.2057860 = 1.090
5 - 3.5769650 3.6538380 1.2321430 2.4969950 3.8010750 1.2413790 = 1.090
6 - 3.3257700 1.5479090 -0.0000016 3.9700480 2.8457950 -0.0000001 = 1.449
7 - 3.9700480 2.8457950 -0.0000001 3.5769650 3.6538380 1.2321430 = 1.525
8 - 3.9700480 2.8457950 -0.0000001 5.4855410 2.7052070 -0.0000044 = 1.522
9 - 5.4855410 2.7052070 -0.0000044 6.0088240 1.5931750 -0.0000084 = 1.229

So then we have

Ebond = kr(r - req)^2

the ONLY bond here that is not at it's equilibrium bond length is 7 which
is 1.525 vs 1.526

Hence the total bond energy is 310.0 * (1.525 - 1.526)^2 = 0.00031 Kcal/mol

Which is exactly what AMBER is reporting. So I'm not sure how you got the
values you did for the AMBER calculation.

Now, looking at your gromacs coordinates:

1ALA N 1 0.333 0.155 -0.000
    1ALA H 2 0.391 0.072 -0.000
    1ALA CA 3 0.397 0.285 -0.000
    1ALA HA 4 0.367 0.340 -0.089
    1ALA CB 5 0.358 0.365 0.123
    1ALA HB1 6 0.388 0.312 0.213
    1ALA HB2 7 0.408 0.462 0.121
    1ALA HB3 8 0.250 0.380 0.124
    1ALA C 9 0.549 0.271 -0.000
    1ALA O 10 0.601 0.159 -0.000

We first see that the coordinates are horribly truncated at 3 decimal
places. So I think what you are seeing is that the coordinate file being
spat out by the conversion script is extremely low precision and this is
what is causing your problems. Nasty.


All the best
Ross


On 8/9/13 10:02 AM, "Golshan Hejazi" <golshan.hejazi.yahoo.com> wrote:

>Dear Adrian,
>
>Thanks for the correction. I computed by hand!!!! for amber and gromacs
>the BOND energy of ALA according to amber topology.
>Here is the results:
>
>GROMACS:
>The number that I obtained: .440148 (kj/mol)
>The number in the gromacs log file: .220074 (kj.mol)
>
>There is a factor of 2 difference ... which I think it is correct ... and
>they take care of it in the energy file.
>
>AMBER:
>The number that I obtained: .052598574280 (kcal/mol)
>The number in the AMBER mdout file: 0.0003 (kcal/mol)
>
>
>What is really the correct value? I think .220074 (kj.mol) is the correct
>value ... because when I convert it to kcal/mol ... I arrive to the same
>number .052598574280 (kcal/mol) ... that I obtained by manually
>computing the BOND energy for amber.
>
>I attached some files if you wish to control it. There is crd and top
>file of amber, gro and top file of gromacs and another file(called test)
>which contains to make it easier for you to check:
>
>r_eq r(distance of the atom according to amber.coor) k(gromacs)
>k(amber)
>
>What do you think?
>
>G.
>
>
>
>
>________________________________
> From: Adrian Roitberg <roitberg.ufl.edu>
>To: Golshan Hejazi <golshan.hejazi.yahoo.com>; AMBER Mailing List
><amber.ambermd.org>
>Sent: Friday, August 9, 2013 11:05 AM
>Subject: Re: [AMBER] BOND, ANGLE and TORSION energies discrepancy between
>amber and gromacs
>
>
>Why would you do
>
>sum (ki*(r_i*r_i)) ??
>
>
>That is not the correct equation. What you need is k*(r_i - r_re_i)**2,
>which is the difference between the actual bond lengths minus the
>equilibrium value.
>
>
>On 8/9/13 10:53 AM, Golshan Hejazi wrote:
>> This is for ALA to avoid any complication!
>>
>>
>> [ bonds ]
>> ; ai aj funct r k
>> 1 2 1 1.0100e-01 3.6317e+05
>> 3 4 1 1.0900e-01 2.8451e+05
>> 5 6 1 1.0900e-01 2.8451e+05
>> 5 7 1 1.0900e-01 2.8451e+05
>> 5 8 1 1.0900e-01 2.8451e+05
>> 1 3 1 1.4490e-01 2.8200e+05
>> 3 5 1 1.5260e-01 2.5941e+05
>> 3 9 1 1.5220e-01 2.6527e+05
>> 9 10 1 1.2290e-01 4.7698e+05
>>
>> >From here, I computed manually the bond energy which is : sum
>>(ki*(r_i*r_i)) = 42536.87 (kj/mol)
>> however, when I perform a rerun on the gro file ... this is the bond
>>energy: 2.20074e-01 (kj/mol) in gromacs
>> and 0.0003 (kcal/mol) from amber
>>
>>
>>
>> ________________________________
>> From: David A Case <case.biomaps.rutgers.edu>
>> To: Golshan Hejazi <golshan.hejazi.yahoo.com>; AMBER Mailing List
>><amber.ambermd.org>
>> Sent: Friday, August 9, 2013 9:46 AM
>> Subject: Re: [AMBER] BOND, ANGLE and TORSION energies discrepancy
>>between amber and gromacs
>>
>>
>> On Fri, Aug 09, 2013, Golshan Hejazi wrote:
>>> I computed the BOND and ANGLE energies. But it is only increasing my
>>> confusion ... Look: I have a simple ace-ala-nme system.
>> How about going to just ace-nme? It would much easier to compare each
>> individual term, and you would still have bonds, angles and dihedrals.
>>
>> While I am thinking about it, be sure that the conversion of the
>>blocking
>> groups (ace and nme) between Amber and Gromacs is correct (i.e. not just
>> the amino acids themselves.)
>>
>> ...dac
>>
>>
>> _______________________________________________
>> 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
>
>--
> Dr. Adrian E. Roitberg
>
>Colonel Allan R. and Margaret G. Crow Term Professor.
>Quantum Theory Project, Department of Chemistry
>University of Florida
>roitberg.ufl.edu
>352-392-6972
>
>
>_______________________________________________
>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 - 15:00:02 PDT
Custom Search