Hi all,
I’m a little stuck on deriving an angular force constant for the amber forcefield.
Consider in amber, that the potential energy from the angle ABC is: U_{ABC} = K(theta - theta_eq)^2 The K is the force constant with the units kcal/(mol rad^2) which reduced to kcal/mol.
Using Gaussian, I am able to get out the Force Constant Matrix (Hessian) from the check point file which is an nxn matrix (n = number of atoms), second derivative of the potential energy relative to the cartesian displacement. The Hessian matrix is in using the Hartree/bohr-radius^2. This can be converted to kcal/(mol ang^2).
Now using the attached maths (see the file) I can pull out the angular force constant simply by taking two 3x3 subsets from the hessian of the displaced atom and the atom experiencing the force. But the problem I face is that I do not know how I can get the obvious displacement units, into radians. If you look at equation (9) there is even a distance_AB and distance_CB term, which is in Angstrom.
How can I get these calculations to reflect the same units at those in the amber force field for angular potential energy?
Thanks
Anthony
Anyone got any thoughts?
Can anyone explain how the distance term (see the image) in the final equations loses it's Angstrom^2 unit? Given that it must, if it were to fit into the Amber force field for a bond angle i.e., K (kcal/(mol rad^2)) * (Angle(rad) - Angle_eq(rad))^2
?
I am really looking for the same answer and I really couldn't understand why all "forcefield parametrize/development articles" don't explain step by step how to develop a forcefield from quantum mechanic calculations. In all articles, some points are missing and you can not even follow half of the instructions.
Hi Anthony Nash and Barış Kurt
I am going through the same process of deriving force fields from frequency analysis in Gaussian.
I am a bit confused about the angular force constant.
K_thetha in eq9 is in Hartree/radians^2 or Hartree/degree^2?
One way the distance term (Angstrom^2 or borh^2) in the denominator of eq9 is lost is that the sum term in the denominator has a Angstrom^-2 or borh^-2 component.
kcal/(mol ang^2).
Eq9 is then similar to:
1/hartree.rad^-2 = 1/borh^2*(hartree.rad^-2.borh^-2) + 1/borh^2*(hartree.rad^-2.borh^-2).
I don't know if that makes sense.
Vuyani Moses any idea?
@Bakary Ntji Diallo you should apply seminaro approach for bonds and angles by using eigen valuse in the gaussian (frequency noraman calculation) output file. ( J. M. Seminario (1996) Calculation of intramolecular force fields from second-derivative tensors. Internat. J. Quant. Chem. 60:1271-1277 )
For dihedrals you should do scan with gaussian
and for ESP you should do POP=MK calculation.
if you want to convert eigen values from mDyne/Å to kcalmol−1Å−2 use this equation: 1mDyne/Å =143,842 kcalmol−1Å−2 but you can use this for only bond force constans (http://archive.ambermd.org/200910/0146.html or http://perso.ibcp.fr/gilbert.deleage/Cours/modelisation.pdf in some articles some authors use 144 directly). After that you have to multiple the result by 0,5 for AMBER force field.
ALso you can use this information for angles: https://www.researchgate.net/post/How_to_Convert_Gaussian_Angular_Force_Constants_units_to_kcal_mol_radian2_units
..And again you should mutiple result by 0,5.
- Badges
- Science topic