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

From: Golshan Hejazi <golshan.hejazi.yahoo.com>
Date: Mon, 12 Aug 2013 06:32:20 -0700 (PDT)

Thanks guys for the helps and suggestions.

Ross, I am using Metadynamics and it is faster on gromacs. I don't know how the PLUMED and other stuff could be lunched on AMBER. That forces me to use gromacs. I will consider your suggestion on how to hack the script for printing with higher precision.

Oliver, it seems this script is for the topology file. Sure it is better to use the latest version but as you see the problem that I was dealing with is related to the coordinate file and I think we have to take into account these kind of simple checks about the bond, angle and dihedrals while converting to gromacs.

G.


________________________________
 From: Oliver Grant <olivercgrant.gmail.com>
To: AMBER Mailing List <amber.ambermd.org>
Cc: Golshan Hejazi <golshan.hejazi.yahoo.com>
Sent: Saturday, August 10, 2013 5:26 AM
Subject: Re: [AMBER] BOND, ANGLE and TORSION energies discrepancy between amber and gromacs
 


Hi all,
acpype is probably what you want to be using for converting amber to gmx:

https://code.google.com/p/acpype/wiki/TutorialAcpype4Gromacs 
It's been updated a lot more than the perl script. 
I don't use it though so I'm not aware of any quirks or if it's able to handle ACE and NME correctly.




On Sat, Aug 10, 2013 at 12:42 AM, Ross Walker <ross.rosswalker.co.uk> wrote:

Hi Golshan,
>
>I don't know of any other converters one could use - or even if Gromacs will
>accept more than 3 decimal places for the input coordinates. If it does then
>I'd suggest trying to hack the conversion script to print more decimal
>places. If it doesn't then the alternative (for doing sanity checks that the
>script is working) would be to quickly write something that zeros the amber
>inpcrd coordinates after the second decimal place.
>
>The real place I'd be wary that things are working or not is in the non
>bonds, so VDW and EEL which might be harder to compare since you have to do
>quite a bit of work to get the exact same calculation setup - but at least
>checking the dihedrals is a good initial sanity check.
>
>And who 'forces' you to use Gromacs? - You should switch to AMBER, get a 4 x
>GPU machine and you'll be able to run WAYYYY more sampling than you can with
>Gromacs. E.g. a 4-way GTX-Titan machine from Exxact can run 4 JAC-NVE runs
>at 111ns/day each so you get almost 450ns/day of aggregate throughput.
>Gromacs won't even get you close to that on such a machine.
>
>All the best
>Ross
>
>From:  Golshan Hejazi <golshan.hejazi.yahoo.com>
>Reply-To:  Golshan Hejazi <golshan.hejazi.yahoo.com>
>Date:  Friday, August 9, 2013 3:32 PM
>To:  Ross  Walker <ross.rosswalker.co.uk>, AMBER Mailing List
><amber.ambermd.org>
>
>Subject:  Re: [AMBER] BOND, ANGLE and TORSION energies discrepancy between
>amber and gromacs
>
>
>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
>


-- 
Oliver Grant
PhD Student,
Computational Glycoscience Laboratory,
School of Chemistry,
National University of Ireland, Galway
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Mon Aug 12 2013 - 07:00:02 PDT
Custom Search