I am trying to do some force field parameterization for a ligand similar to
how GAFF was developed. As a test, I decided to first attempt to recreate
Figure 3 in Junmei Wang's 2004 paper, "Development and Testing of a General
Amber Force Field." This figure shows how well the the proper torsional
parameters help the dihedral profile of 2,3-butanedione match that
determined at MP4/6-311G(d,p)//MP2/6-31G*. I first scanned the c-c-c-c
dihedral angle at 30 degree intervals from -180 to 180 with Gaussian03:
# MP2/6-31G* Opt=ModRedundant
dihedral scan of 2,3-butanedione
0 1
C 0.432 0.806 1.762
...
H -2.858 0.571 -1.527
1 5 3 2 -180.0 S 12 30.0
Then the energies of the optimized MP2 geometries were calculated at
MP4/6-311G(d,p):
# MP4/6-311g(d,p) SP
Also, the ESP derived charges for the molecule were determined at HF/6-31G*
as mentioned in Wang et al.:
# HF/6-31G*//HF/6-31G* pop=mk scf=tight iop(6/33=2) iop(6/42=6) test
and fit to the atomic centers with the two-stage fit automated by
'antechamber'. The resulting charges are found in the following link to a
prepi file created for the small molecule:
0 0 2
This is a remark line
molecule.res
BDE INT 0
CORRECT OMIT DU BEG
0.0000
1 DUMM DU M 0 -1 -2 0.000 .0 .0 .00000
2 DUMM DU M 1 0 -1 1.449 .0 .0 .00000
3 DUMM DU M 2 1 0 1.522 111.1 .0 .00000
4 C1 c3 M 3 2 1 1.540 111.208 180.000 -0.43474
5 H1 hc E 4 3 2 1.084 90.423 -53.479 0.12751
6 H2 hc E 4 3 2 1.080 34.705 -180.000 0.12751
7 H3 hc E 4 3 2 1.084 90.423 53.437 0.12751
8 C4 c M 4 3 2 1.506 143.914 -180.000 0.56915
9 O2 o E 8 4 3 1.190 124.063 0.000 -0.51693
10 C3 c M 8 4 3 1.538 116.849 180.000 0.56915
11 O1 o E 10 8 4 1.190 119.088 0.000 -0.51693
12 C2 c3 M 10 8 4 1.506 116.849 -180.000 -0.43474
13 H4 hc E 12 10 8 1.084 110.176 -58.835 0.12751
14 H5 hc E 12 10 8 1.084 110.154 58.874 0.12751
15 H6 hc E 12 10 8 1.080 109.209 180.000 0.12751
LOOP
IMPROPER
C1 C3 C4 O2
C2 C4 C3 O1
DONE
STOP
Finally the dihedral profile of the 2,3-butanedione was created in Amber9
using GAFF and the above prepi file with nmropt=1 to control the dihedral
angle with strong harmonic restraints:
Minimization of the entire molecular system
&cntrl
imin=1, ncyc=9, maxcyc=500000, drms=0.0001,
ntpr=1, ntb=0, igb=0, cut=12,
scnb = 2.0, scee = 1.2,
nmropt=1
/
&wt type='END'
/
LISTOUT=POUT
DISANG=RST
The following link contains my results compared to Figure 3 in Wang et al.
In general, my dihedral profile produced in Amber with GAFF is fairly close
to that used in Wang et al., however my ab initio dihedral profile is quite
different. I calculated a profile that is much more flat in the range from
-60 to 60 and much more steep in the other regions.
https://docs.google.com/open?id=0BwkQMO2EgiyoZWY0MWRjYTEtNTFjMC00NjNlLTg0MTUtM2MxZTExNjc5ZTVh
I scanned the PK values for the V2 potential and found that based on my
calculation methods, a PK value of 2.2 resulted in the lowest RMSD between
the Amber and Gaussian calculated dihedral profiles (13 points from -180 to
180). With PK=2.2, my calculated RMSD was 0.432. At PK=1.2, where Wang et
al. report their lowest RMSD of 0.2394, i calculate a RMSD of 1.285. As
you'll notice in the figure, by optimizing the dihedral profile based on
RMSD with only a V2 potential, i obtain a decent RMSD; however, the profile
is poor, with a pronounced local minimum at 0 degrees instead of a maximum!
Can anyone help point out where i may be going wrong in my attempt to
reproduce this figure/methodology? Thanks in advance for our help!
-ryan
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Feb 07 2012 - 08:00:03 PST
