Hi Ryan,
I can basically reproduce your MM results and the Wang et al QM results. However, I took a slightly different approach for the quantum calculations and the partial atomic charges. Note that due to symmetry, one only needs to do a heavy atom rotation from 0-180 degrees (I also did this in 30 degree increments). I wasn't able to do the single-point calculation at the MP4 level, instead I opted for a larger basis set at the MP2 level.
For the QM constraint optimization about the C-C-C-C torsion angle I used MP2//6-31++G(d,p). I performed single-point calculations on these geometries at MP2/aug-cc-pVDZ (i.e. MP2/aug-cc-pVDZ//MP2//6-31++G(d,p)). For partial atomic charges I used R.E.D., with HF/6-31G(d) electrostatic potential using the Connolly surface algorithm, and a 2 stage RESP fit using weighting factors qwt=.0005/.001. (3 molecular orientations were used in the R.E.D. calculations.) All of this was done using GAMESS.
My partial atomic charges are slightly different than yours, but ultimately I don't think it matters too much for the internal heavy atom rotation of this molecule. Regardless, here are the charges:
1 C1 -0.7300 0.2420 0.0000 c 1 MOL 0.5239
2 C2 0.7300 -0.2420 0.0000 c 1 MOL 0.5239
3 C3 1.7980 0.8190 0.0000 c3 1 MOL -0.2739
4 H4 1.6930 1.4560 0.8710 hc 1 MOL 0.0861
5 H5 1.6930 1.4560 -0.8710 hc 1 MOL 0.0861
6 H6 2.7700 0.3460 0.0000 hc 1 MOL 0.0861
7 C7 -1.7980 -0.8190 0.0000 c3 1 MOL -0.2739
8 H8 -1.6930 -1.4560 0.8710 hc 1 MOL 0.0861
9 H9 -1.6930 -1.4560 -0.8710 hc 1 MOL 0.0861
10 H10 -2.7700 -0.3460 0.0000 hc 1 MOL 0.0861
11 O11 -0.9520 1.4110 -0.0000 o 1 MOL -0.5082
12 O12 0.9520 -1.4110 -0.0000 o 1 MOL -0.5082
Attached you can see the resulting potential energy curves. I can essentially reproduce the gaff curve and your curve. What I noticed is that the resulting molecular mechanics curve is extremely sensitive to the V2 parameter! By changing the parameter from 1.2 to 2.2 the resulting curve at 0 degrees goes from a transition state to a high energy minima.
Concerning the QM calculations, both MP2 calcuations have no local minima at 30 degrees, which align well with Wang et al MP4 curve. It is clear that the shape of this potential is dependent upon the basis set, and most likely the theory level too. I would double check your 30 degree conformation (perhaps there is an issue with the improper torsion geometry about the ketone functionality).
Concerning the RMSD. I obtain a lower RMSD (averaged over all 7 conformations, not just the minima) for gaff than for your parameter, which is opposite of what you find. For comparison to my MP2//6-31++G(d,p) geometries, I obtain an RMSD value of 0.056 for gaff, and 0.058 when using your V2.
I don't think this answers your question(s) directly, but I hope that helps some :) .
Cheers,
Karl
----- Original Message -----
From: "Ryan Pavlovicz" <pavlovicz.7.osu.edu>
To: "AMBER Mailing List" ambermd.org>
Sent: Tuesday, February 7, 2012 4:40:05 PM
Subject: [AMBER] Force Field Parameterization -- Torsion Potentials
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
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Wed Feb 08 2012 - 05:30:04 PST
