Hi,
As Carlos mentioned, due to the way coordinates are printed in Amber,
the energies that are printed in the MDOUT file correspond to the
coordinates of the step _just prior_ to the coordinates that are
written to MDCRD (i.e. you get energies of i and coordinates of i+1).
The reasons are historical and have to do with getting kinetic energy
from the Leapfrog scheme.
Second, unless you have a version of cpptraj where the 'energy' action
has the 'nobondstoh' keyword (>= 6.12.0, available on GitHub) the
'energy' bond calculation is including all bonds with hydrogen, which
sander will ignore when ntf=2 (typical when SHAKE is on for H, ntc=2).
When I run SANDER on your GAFF Ala2 system with ntc=ntf=1 (one step
minimization) I get perfect agreement between cpptraj and sander
(Amber 22):
My mdin:
minimization
&cntrl
imin = 1, ntmin = 2, maxcyc = 1,
ntx=1, irest = 0, ig = -1,
ntwx = 1, ioutfm = 1, ntpr = 1, ntwr = 10000,
ntc = 1, ntf = 1, ntb = 0, cut = 9999.0,
igb = 6,
&end
sander:
NSTEP ENERGY RMS GMAX NAME NUMBER
1 5.0924E+01 1.4008E+01 5.6101E+01 C 21
BOND = 2.8260 ANGLE = 7.0458 DIHED = 9.4700
VDWAALS = 1.9340 EEL = -274.5862 EGB = 0.0000
1-4 VDW = 5.5213 1-4 EEL = 298.7130 RESTRAINT = 0.0000
cpptraj input:
parm ala2_gaff.prmtop
trajin ala2_gaff.rst7
energy out cpptraj.energy.dat
cpptraj:
#Frame ENE_00001[bond] ENE_00001[angle] ENE_00001[dih]
ENE_00001[vdw14] ENE_00001[elec14] ENE_00001[vdw] ENE_00001[elec]
ENE_00001[total]
1 2.8260 7.0458 9.4700
5.5213 298.7130 1.9340 -274.5862
50.9239
Third, I don't know how you're doing your MD calculation (i.e. what
your MDIN is etc), but I can't reproduce the numbers you have in
ala2_gaff.info file even with ntc=ntf=2. Your numbers:
NSTEP = 0 TIME(PS) = 0.000 TEMP(K) = 0.00 PRESS = 0.0
Etot = 50.3558 EKtot = 0.0000 EPtot = 50.3558
BOND = 2.2579 ANGLE = 7.0458 DIHED = 9.4700
1-4 NB = 5.5213 1-4 EEL = 298.7130 VDWAALS = 1.9340
EELEC = -274.5862 EHBOND = 0.0000 RESTRAINT = 0.0000
My sander numbers (same mdin as above but ntc=ntf=2):
NSTEP ENERGY RMS GMAX NAME NUMBER
1 5.0161E+01 1.4929E+01 5.6167E+01 C 21
BOND = 2.3174 ANGLE = 7.0012 DIHED = 9.4710
VDWAALS = 1.9844 EEL = -274.5609 EGB = 0.0000
1-4 VDW = 5.5033 1-4 EEL = 298.4448 RESTRAINT = 0.0000
Are you sure you were using the right coordinates in your calculation?
-Dan
On Mon, Aug 29, 2022 at 8:32 AM Ruihan Zhou via AMBER <amber.ambermd.org> wrote:
>
> Dear Amber Users,
>
> I'm working on converting amber force fields to another MD engine LAMMPS.
> During the conversion I've noticed a bond energy discrepancy between the
> results from amber simulations and cpptraj results when using GAFF and I
> want to know why this happens.
>
> Take dialanine as an example, we've run amber simulations without periodic
> boundary conditions under NVE with a cutoff effectively covering the whole
> molecule using FF14SB, FF19SB and GAFF. With cpptraj on the same molecule,
> we take in the initial configuration and output the energy.
>
> With FF14SB, FF19SB, all molecular energies have matched, but for GAFF the
> bond energy was different. This observation holds true for some other
> molecules we've tested like n-methylacetamide and trialanine. (As a side
> note, the bond energies from cpptraj matched all the LAMMPS test we've ran.)
>
> The bond energy comparisons for ala2 are listed here:
> ff14sb 0.7569 cpptraj 0.7569
> ff19sb 0.7846 cpptraj 0.7846
> gaff 2.2579 cpptraj 2.8260
>
> For ala3:
> ff14sb 0.7572 cpptraj 0.7572
> ff19sb 0.7723 cpptraj 0.7723
> gaff 3.2580 cpptraj 3.8895
>
> The NVE.in script I've used is like this:
> Production NVE
> &cntrl
> imin=0,
> ntx=1,
> irest=0,
> nstlim=1,
> dt=0.001,
> ntf=2,
> ntc=2,
> ntpr=10,
> ntwx=10,
> cut=999.0,
> ntb=0,
> ntp=0,
> ntt=0,
> ig=-1,
> tol = 0.0000001,
> nscm=0,
> /
>
> Here's the cpptraj log:
> # Args: -p ala2_gaff.prmtop -c ala2_gaff.rst7
> # Loaded topologies:
> # ala2_gaff.prmtop
> energy system out sys.agr
> run
> quit
>
> In the attachment I have the prmtop and rst7 files for ala2 of all the
> force fields tested. All info files are from amber simulation and all agr
> files are energy outputs from cpptraj.
>
> I would appreciate any help on this issue! Thanks in advance!
>
>
> Best,
> Ruihan
> _______________________________________________
> 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 Tue Aug 30 2022 - 07:30:02 PDT
