Re: [AMBER] Force Field Parameterization -- Torsion Potentials

From: Karl N. Kirschner <kkirsch.scai.fraunhofer.de>
Date: Mon, 19 Mar 2012 13:50:46 +0100 (CET)

Hi Mohd,

  I understand your question better now. Do a minimization of the single molecule with Amber. You want to compare "apples to apples" so to speak.

  Here is a brief explanation why. Using quantum mechanics you performed a gas-phase (I'm assuming you didn't use an implicit solvent model) energetic minimization of a molecule until it reached the minimum of the local potential energy surface, and based on a certain convergence criteria. Unless you did something very special in the QM calculation to account for the effects of temperature and vibrations, your minimization represents the molecule's internal coordinates at 0 Kelvin. Thus, you want to perform a very similar minimization using molecular mechanics. (If you included an geometric constraint in the QM calculation, then you must include it in the MM calculation.) Molecular dynamics is a different theory, one that uses molecular mechanics, but it does not perform minimizations and it includes temperature effects. Therefore, directly comparing a QM minimized geometry to geometries generated during an MD simulation is comparing "apples to oranges."

Cheers,
Karl

----- Original Message -----
From: "Mohd F. Ismail" <farid.ou.edu>
To: "AMBER Mailing List" ambermd.org>
Sent: Friday, March 16, 2012 6:31:02 AM
Subject: Re: [AMBER] Force Field Parameterization -- Torsion Potentials

Dr. Kirschner,

I want to compare the structure i.e. the bonds and the angles first. And later the torsion energy.

So I computed the MP2 geometry optimization with GAMESS. Now I am not sure how to get the GAFF version of the bonds and angles. Do I do minimization of a single molecule with Amber? Do I run a single MD?

Sorry for using this thread, but I thought it is somewhat related.

*******************************
Mohd Farid Ismail
Graduate Student
Dept. of Chemistry/Biochemistry
University of Oklahoma
Norman 73019

________________________________________
From: Karl N. Kirschner [kkirsch.scai.fraunhofer.de]
Sent: Wednesday, March 14, 2012 6:53 AM
To: AMBER Mailing List
Subject: Re: [AMBER] Force Field Parameterization -- Torsion Potentials

Hello Mohd,

  Your question is still a bit unclear. Is your question regarding reproducing torsion potentials (as indicated in the subject line), or is it regarding molecular structure comparison of minima, or is it in regards to nonbonded parameters? What are the general goals of your study?

Cheers,
Karl

----- Original Message -----
From: "Mohd F. Ismail" <farid.ou.edu>
To: "AMBER Mailing List" ambermd.org>
Sent: Tuesday, March 13, 2012 11:14:08 PM
Subject: Re: [AMBER] Force Field Parameterization -- Torsion Potentials

Ahh, I apologize, that was too short a message. What I want to do is perform regular organic solvent simulation. Since some of the densities do not match, I want to compare the result of MP2 calculation to the GAFF-based calculation. I have perform the calculation for the geometry with GAMESS. However, how do I get the GAFF-based structure using Amber? Do I just perform the minimization of 1 molecule?

*******************************
Mohd Farid Ismail
Graduate Student
Dept. of Chemistry/Biochemistry
University of Oklahoma
Norman 73019
________________________________________
From: Ismail, Mohd F. [farid.ou.edu]
Sent: Tuesday, March 13, 2012 5:08 PM
To: AMBER Mailing List
Subject: Re: [AMBER] Force Field Parameterization -- Torsion Potentials

Dr. Kirschner,

Can I ask how did you calculate the MM part?

*******************************
Mohd Farid Ismail
Graduate Student
Dept. of Chemistry/Biochemistry
University of Oklahoma
Norman 73019
________________________________________
From: Karl N. Kirschner [kkirsch.scai.fraunhofer.de]
Sent: Thursday, February 09, 2012 3:25 AM
To: AMBER Mailing List
Subject: Re: [AMBER] Force Field Parameterization -- Torsion Potentials

Hi Ryan,

  In regards to the partial atomic charges. I would recommend reading the following paper and visiting the R.E.D. group's website (http://q4md-forcefieldtools.org/RED):

F.-Y. Dupradeau, A. Pigache, T. Zaffran, C. Savineau, R. Lelong, N. Grivel, D. Lelong, W. Rosanski & P. Cieplak, The R.E.D. tools: Advances in RESP and ESP charge derivation and force field library building, Phys. Chem. Chem. Phys. 2010, 12, 7821-7839

  The conformation I used was the minima at 180 degrees for the ESP QM calculations. Three ESP QM calculations were done, each using a different molecular orientation in x,y,z coordinate space of the 180 degree conformation. This is need due to slight differences that arise when the ESP places a grid around the molecule. The use of three orientation help make the partial atomic charges reproducible.

Cheers,
Karl

----- Original Message -----
From: "Ryan Pavlovicz" <pavlovicz.7.osu.edu>
To: "AMBER Mailing List" ambermd.org>
Sent: Wednesday, February 8, 2012 6:14:45 PM
Subject: Re: [AMBER] Force Field Parameterization -- Torsion Potentials

Hi Karl,

Thanks a lot for your help!

Although the MM profile seems to be the smaller of my two problems, I have
played around with the partial atomic charges that you obtained from your
three-conformation fit, and have found that using the charges you derived
increases my relative MM maximum at 0 degrees by 0.5 kcal/mol. This
ultimately helps lower that RMSD from the QM profile as show figure I link
to below:

https://docs.google.com/open?id=0BwkQMO2EgiyoZmU2NzYyYjctMmIxNS00MzM3LWI2MDEtMGU2N2RiYjJiNGM1

I am wondering which three 2,3-butanedione conformations you used for your
charge calculations.

Also, I am now playing around with regenerating my QM profile, which seems
to be the most important in recreating the results shown in Figure 3 of
Wang et al. I'll get back to you with the results when the computations
have completed. Thanks again, Karl!

-ryan

On Wed, Feb 8, 2012 at 8:11 AM, Karl N. Kirschner <
kkirsch.scai.fraunhofer.de> wrote:

> 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

_______________________________________________
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 Mon Mar 19 2012 - 06:00:04 PDT
Custom Search