[gmx-users] OPLS-AA Coulomb-14 energies incorrect?
Malcolm Gillies
malcolm.b.gillies at anu.edu.au
Wed Aug 6 09:01:01 CEST 2003
Hi,
I appear to have found an error in the calculation of the Coulomb 1-4
energies using the GROMACS OPLS-AA implementation. Have potential
energies calculated with this force field been validated against those
from some other implementation (e.g. BOSS, MCPRO, Impact, MacroModel
etc)?
The Coulomb (SR) component of the GROMACS energy agrees with the
corresponding energy from TINKER (i.e. as calculated by setting
chg-14-scale to 0.0), but the Coulomb-14 energy is out by a sizeable
percentage. i.e for the attached file GROMACS reports 224.739 kJ/mol
while TINKER reports 258.3687 kJ/mol.
I've attached copies of my GROMACS and TINKER input files. I had to
modify the TINKER 4.0 oplsaa.prm file to include an atom type for the
GLY zwitterion CA with the correct charge; I've included this
modification as a patch file.
Comparing GROMACS (OPLS-AA parameter files checked out of CVS about
a month ago) and TINKER 4.0 (with the supplied oplsaa.prm parameters) on
a couple of potential energy calculations, I found that they agree on
the bond stretch, bend and LJ energies. The improper dihedral energies
differ, though this may be due to a problem with TINKER. GROMACS uses
some updated proper dihedral parameters, absent from the TINKER
oplsaa.prm file, so for most structures the proper dihedral energies are
not comparable (though they seem to match for the glycine zwitterion).
cheers,
Malcolm
--
Malcolm Gillies <Malcolm.B.Gillies at anu.edu.au>
Postdoctoral Fellow, Computational Proteomics and Therapy Design Group,
John Curtin School of Medical Research, Australian National University
-------------- next part --------------
;
; File 'GLY.top' was generated
; By user: mxg (252)
; On host: bakunin
; At date: Wed Aug 6 15:11:56 2003
;
; This is your topology file
; "That Was Really Cool" (Butthead)
;
; Include forcefield parameters
#include "ffoplsaa.itp"
[ moleculetype ]
; Name nrexcl
Protein_A 3
[ atoms ]
; nr type resnr residue atom cgnr charge mass typeB chargeB massB
1 opls_287 1 GLY N 1 -0.3 14.0067 ; qtot -0.3
2 opls_290 1 GLY H1 1 0.33 1.008 ; qtot 0.03
3 opls_290 1 GLY H2 1 0.33 1.008 ; qtot 0.36
4 opls_290 1 GLY H3 1 0.33 1.008 ; qtot 0.69
5 opls_298 1 GLY CA 1 0.09 12.011 ; qtot 0.78
6 opls_140 1 GLY HA1 1 0.06 1.008 ; qtot 0.84
7 opls_140 1 GLY HA2 1 0.06 1.008 ; qtot 0.9
8 opls_271 1 GLY C 2 0.7 12.011 ; qtot 1.6
9 opls_272 1 GLY O1 2 -0.8 15.9994 ; qtot 0.8
10 opls_272 1 GLY O2 2 -0.8 15.9994 ; qtot 0
[ bonds ]
; ai aj funct c0 c1 c2 c3
1 2 1
1 3 1
1 4 1
1 5 1
5 6 1
5 7 1
5 8 1
8 9 1
8 10 1
[ pairs ]
; ai aj funct c0 c1 c2 c3
1 9 1
1 10 1
2 8 1
3 8 1
4 8 1
6 9 1
6 10 1
7 9 1
7 10 1
[ angles ]
; ai aj ak funct c0 c1 c2 c3
2 1 3 1
2 1 4 1
2 1 5 1
3 1 4 1
3 1 5 1
4 1 5 1
1 5 6 1
1 5 7 1
1 5 8 1
6 5 7 1
6 5 8 1
7 5 8 1
5 8 9 1
5 8 10 1
9 8 10 1
[ dihedrals ]
; ai aj ak al funct c0 c1 c2 c3 c4 c5
2 1 5 6 3
2 1 5 7 3
2 1 5 8 3
3 1 5 6 3
3 1 5 7 3
3 1 5 8 3
4 1 5 6 3
4 1 5 7 3
4 1 5 8 3
1 5 8 9 3
1 5 8 10 3
6 5 8 9 3
6 5 8 10 3
7 5 8 9 3
7 5 8 10 3
[ dihedrals ]
; ai aj ak al funct c0 c1 c2 c3
5 10 9 8 1 improper_O_C_X_Y
; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif
; Include water topology
#ifdef FLEX_SPC
#include "flexspc.itp"
#else
#include "spc.itp"
#endif
#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
; i funct fcx fcy fcz
1 1 1000 1000 1000
#endif
[ system ]
; Name
Protein
[ molecules ]
; Compound #mols
Protein_A 1
-------------- next part --------------
S C A M O R G
10
1GLY N 1 0.442 0.702 -0.555
1GLY H1 2 0.470 0.621 -0.607
1GLY H2 3 0.456 0.685 -0.457
1GLY H3 4 0.345 0.721 -0.572
1GLY CA 5 0.522 0.816 -0.596
1GLY HA1 6 0.628 0.792 -0.578
1GLY HA2 7 0.509 0.832 -0.704
1GLY C 8 0.490 0.946 -0.524
1GLY O1 9 0.391 0.942 -0.430
1GLY O2 10 0.563 1.056 -0.559
0.28265 0.43423 0.27368
-------------- next part --------------
title = molecular dynamics
cpp = /lib/cpp
constraints = none
integrator = md
dt = 0.001 ; ps !
nsteps = 1
comm_mode = none
nstcomm = 0
nstxout = 1
nstvout = 10000
nstfout = 0
pbc = no
nstlist = 0
ns_type = simple
rlist = 0.0
rcoulomb = 0.0
rvdw = 0.0
;coulombtype = cutoff
;rcoulomb_switch = 0.0
dispcorr = no
;epsilon_surface = 78.0
;fourierspacing = 0.1
;pmeorder = 6
;loptimizefft = no
; Nose-Hoover coupling is on in two groups
Tcoupl = no
;tau_t = 0.1 0.1
;tc-grps = protein sol
;ref_t = 300 300
; Pressure coupling is on
Pcoupl = no
;Pcoupltype = isotropic
;tau_p = 0.5
;compressibility = 4.5e-5
;ref_p = 1.0
; Generate velocites is on at 300 K.
gen_vel = yes
gen_temp = 300.0
gen_seed = 173529
-------------- next part --------------
10
1 N3 4.420000 7.020000 -5.550000 146 2 5 6 7
2 CT 5.220000 8.160000 -5.960000 211 1 3 8 9
3 C_3 4.900000 9.460000 -5.240000 111 2 4 10
4 O2 3.91 9.42 -4.30 112 3
5 H3 4.70 6.21 -6.07 151 1
6 H3 4.56 6.85 -4.57 151 1
7 H3 3.45 7.21 -5.72 151 1
8 HC 6.280000 7.920000 -5.780000 6 2
9 HC 5.090000 8.320000 -7.040000 6 2
10 O2 5.63 10.56 -5.59 112 3
-------------- next part --------------
A non-text attachment was scrubbed...
Name: oplsaa.prm-patch
Type: text/x-patch
Size: 751 bytes
Desc: not available
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20030806/822da57b/attachment.bin>
More information about the gromacs.org_gmx-users
mailing list