Ross, thanks a lot. I really appreciate your effort to clarify the source of problem. This was so Nasty!!!!
I did exactly the same calculation, but I used the coordinate of gromacs ...!!!
In fact, I checked and luckily, the different energy of two structure are comparable in gromacs and amber ... regardless that the absolute value of energies are not matching !!!
So, at this point, what is your suggestion! I am forced to use gromacs ... but I will use amber force fields. Do you know of any better converter ...!
Any suggestion will help :)
Best
G.
________________________________
From: Ross Walker <ross.rosswalker.co.uk>
To: Golshan Hejazi <golshan.hejazi.yahoo.com>; AMBER Mailing List <amber.ambermd.org>
Sent: Friday, August 9, 2013 5:54 PM
Subject: Re: [AMBER] BOND, ANGLE and TORSION energies discrepancy between amber and gromacs
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 - 16:00:03 PDT
