[gmx-users] Update Amber force fields, ff12SB abd ff14SB gromacs implementation testing

[Man Hoang Viet]

Hi Justin,

> Date: Mon, 28 Dec 2015 08:18:45 -0500
> From: Justin Lemkul <jalemkul at vt.edu>
> To: gmx-users at gromacs.org
> Subject: Re: [gmx-users] Update Amber force fields, ff12SB abd ff14SB
> 	gromacs implementation testing
> Message-ID: <568136B5.5090109 at vt.edu>
> Content-Type: text/plain; charset=windows-1252; format=flowed
>
>
>
> On 12/28/15 1:40 AM, Man Hoang Viet wrote:
>> Hi,
>>
>>> Date: Sun, 27 Dec 2015 16:03:50 -0500
>>> From: Justin Lemkul <jalemkul at vt.edu>
>>
>>> The energy of the same configuration in both AMBER and GROMACS should
>>> be
>>> the
>>> same if the force field implementation is equivalent.  You don't need
>>> to
>>> run
>>> simulations, extract conformations, etc. to get this.  Take one
>>> coordinate
>>> file,
>>> calculate its energy in each program, and report the result.  The Sorin
>>> group
>>> did a slightly more complex version of this, by testing several
>>> different
>>> coordinate perturbations (the reason for which is explained on their
>>> website),
>>> but the for the purpose of troubleshooting, a single configuration
>>> should
>>> do
>>> initially.  It's really quite straightforward.  If the same coordinates
>>> produce
>>> different energies, there is a fundamental problem in the force field
>>> implementation.  If the problem is shared among Gly and other amino
>>> acids,
>>> this
>>> implies a backbone problem, which is very troubling, but the fact that
>>> you
>>> can't
>>> reproduce the ff99 result means you're likely doing something
>>> incorrectly
>>> in
>>> your assessment.  Be sure you're doing single-point energy calculations
>>> correctly using mdrun -rerun.
>>>
>>> -Justin
>>
>> My simulation for energy calculation is steep minimization with
>> nsteps=0.
>
> This is the wrong way to calculate a single point energy.
>
> http://www.gromacs.org/Documentation/How-tos/Single-Point_Energy
>
> Note the last paragraph.
>
> -Justin
>
>> This simulation gives energy information of the initial structures which
>> are the same in both AMBER and GROMACS simulation.
>> I can't not produce exact conformations as The Sorin done group for
>> testing without their code, however I do test for the initial structures
>> generated by Leap in AMBER by the same way I did for my previous tests.
>> This test result is shown in below table:
>> Error(%) & ff99   & ff99SB  & ff12SB  & ff14SB
>> BOND     & 0.001  &  0.001  &  0.001  &  0.001
>> ANGLE    & 0.000  &  0.000  &  0.000  &  0.000
>> DIHEDRAL & 0.000  &  0.000  &  0.000  &  0.000
>> 1-4vdW   & 0.000  &  0.000  &  0.000  &  0.000
>> 1-4QQ    & 0.000  &  0.000  &  0.000  &  0.000
>> vdW      & 0.000  &  0.000  &  0.000  &  0.000
>> QQ       & 0.002  &  0.002  &  0.002  &  0.002
>> Epot     & 0.005  &  0.004  &  0.004  &  0.004
>>
>> In my view, it is not clear what cause the dihedral energy difference.
>> My
>> test assessment, the force field implementations (ff99, ff99SB, ff12SB
>> and
>> ff14SB) and the different  dihedral energy calculations between GROMACS
>> and AMBER are all possible for making the problem. I will try to figure
>> out this issue.

Thank you very much for the information. I have corrected my test as below
commands:
grompp_s  -f eqtest.mdp -p  v.top -c  test.pdb -o  test.tpr
mdrun_s -nt 1 -s test.tpr -rerun test.pdb -g test.log
Then the energy information were taken from test.log

eqtest.mdp:
##############
integrator          =  md
emtol               =  100.0
emstep              =  0.005
dt                  =  0.000001    ; ps !
nsteps              =  1  ; total 100 ps
; Parameters describing how to find the neighbors of each atom and how to
calculate the interactions
nstlist         = 10
ns_type         = grid
rlist           = 1.0
coulombtype     = cut-off
rcoulomb        = 1.0
vdwtype         = cut-off
rvdw            = 1.0
ewald_rtol      = 1e-5
###############

However, the result is still the same as what I reported. Could you please
tell me other way to directly calculate diheral engery by GROMACS from a
give structure (or dihedral)? Till now, I do not think my test assessment
causes the difference of the dihedral energies.


Viet Man,
Postdoctoral Research Scholar,
Department of Physics, NCSU