Hello,
On Mon, Oct 4, 2010 at 3:34 PM, Ilyas Yildirim
<i-yildirim.northwestern.edu>wrote:
> Dear All,
>
> It is not so clear in the manual on how the potential energy is calculated
> - explicitely - for some terms (esp. van der Waals). There are parameters
> derived for some molecules in other force fields, and transferring
> those parameters for AMBER simulations needs a good understanding of the
> AMBER potential energy function. Whether those parameters can be
> transferred or not is out of the scope of this email. Amber 7 manual has
> Appendix C written by Allison Howard and Bill Ross, which talks about
> AMBER parameter development. I would like to go step by step on each term
> in the potential energy function. If you see any mistakes, I will
> appreciate to hear your response on them.
>
> 1. Bond and Angle: The equation used to calculate the bond and angular
> energies in AMBER are the same: K*(x-x_eq)^2, where x is the bond/angle
> measured, x_eq is the equilibrium bond distance or angle, and K is the
> force constant. In some force fields, they use (0.5*K_r) as force constant
> to describe the bond and angular energies, but this is not true in AMBER
> potential energy. For instance, in parm99.dat the following entries define
> the CT-CT and CT-CT-CT parameters:
>
> CT-CT 310.0 1.526 JCC,7,(1986),230; AA, SUGARS
> CT-CT-CT 40.0 109.50
>
> What this says is that the bond and angular energies for CT-CT bond and
> CT-CT-CT angle are calculated as follows:
>
> E_bond = (310.0)*(r-1.526)^2
>
> E_angle = (40.0)*(x-109.50)^2
>
> The units of the parameters are such that at the end the energy is in
> kcal/mol.
>
Correct. You can verify this in the code as well
($AMBERHOME/src/sander/ene.f). This is done, I believe, for runtime
efficiency. It doesn't make sense to multiply by 0.5 for every bond for
every step when the factor of 1/2 can be taken out beforehand. This is true
for a lot of the parameters in Amber. In several instances, parameters are
transformed inside leap before being written to the prmtop for runtime
efficiency as well, though this is largely transparent to the user (until
they try to dissect the prmtop file).
> 2. Torsions: From Amber 7 Manual (pp 290), the equation used to calculate
> the torsional energy is as follows:
>
> T_tor = (PK/IDIVF)*(1+cos(PN*phi-PHASE))
>
> For instance in parm99, the parameters for CT-CT-CT-CT is
>
> CT-CT-CT-CT 1 0.18 0.0 -3.
> CT-CT-CT-CT 1 0.25 180.0 -2.
> CT-CT-CT-CT 1 0.20 180.0 1.
>
> What this says is, the energy for this torsion is explicitly calculated as
>
> E_tor = (0.20/1)*(1+cos(phi-180)) + (0.25/1)*(1+cos(2*phi-180)) +
> (0.18/1)*(1+cos(3*phi-0.0))
>
> In some force fields, the torsional part of the potential energy function
> is calculated as
>
> E_tor = Sum [ (Vn/2)*(1+cos(n*phi-phase))]
>
> In AMBER, the energy barrier Vn is not directly used. For instance, PK is
> given/written explicitly in the parm99.dat file, which is equal to Vn/2.
>
This is true, too. However, the prmtop file has only 3 "sections" for
dihedral stuff. They store only the height (PK), periodicity (PN), and
phase (PHASE). PK is divided by IDIVF inside LEaP before it's written to
the prmtop, so there's no way of recovering IDIVF from the prmtop (not that
it matters).
> There is also the 1-4 terms in the energy calculations, which I am not
> going to talk here. The main difference between different force fields is
> the scaling factor used in the 1-4 calculations. Again, the scope of this
> email message is not to discuss whether parameters can be transferred from
> one force field to another.
>
Only 1-4 non-bonded (EEL/VDW) terms are scaled, which only affect how the
dihedral terms are derived, not how they're calculated in sander/pmemd/nab.
> 3. van der Waals: According to AMBER 7 manual, the vdw parameters given in
> .dat files are R* and e. R* is half the interatomic distance at the
> minimum energy point, and e is the potential well depth of the minimum
> enery point between 2 atoms of the same type. Mixing different atom types
> are done by taking the root of the product.
>
> In some force fields, the vdw energy term explicitly uses the minimum
> interactomic distance R. So, the vdw pot. energy term in those
> force fields are
>
> E_vdw = E0 *[ (R0/R)^12 - 2(R0/R)^6].
>
> where E0 is the potential well depth of the minimum energy point, and R0
> is the interatomic distance at the minimum energy point.
>
> I am not so clear how AMBER explicitly calculates the vdw energy. For
> instance, in parm99.dat file, the vdw parameters for CT is
>
> CT 1.9080 0.1094 Spellmeyer
>
> where R/2=1.9080 and e=0.1094. If it is going to calculate the vdw energy
> of two CT type atom, how is it done using these parameter set? Is it
>
> E_vdw = 0.1994 * [ ((2*1.9080)/R )^12 - 2 * ((2*1.9080)/R)^6 ] ?
>
> PS: I am excluding the 1-4 term here.
>
The PDF linked by Bill is helpful here and explains exactly how they're
calculated (
http://ambermd.org/vdwequation.pdf). Again, the well depth and
radius are combined by leap into ACOEF and BCOEF, so the potential is
ACOEF/R^12 - BCOEF/R^6, which is cheaper than all of the other algebraic
manipulations. The end result is the same, though, as long as those are the
combining rules you're after.
> 4. electrostatic: This part uses the resp charges in AMBER, so a direct
> charge transfer from a force field to AMBER is not possible. It needs to
> be calculated.
>
This is unrelated to your statement (which is true), but the charges are
multiplied by 18.2223 before being printed to the prmtop, which is the
square root of "k" (~332) that gives electrostatic energies in kcal/mol.
>
> I am really hesitant to transfer parameters from one force field to
> another, but sometimes it needs to be done. I will appreciate any comments
> on this issue. Thanks in advance.
Agreed. Not trivial. What I would look for is how the different parameters
are derived (and, obviously whether or not the potential functions are the
same/similar) in addition to how they're calculated, since force field terms
are not orthogonal (1-4 scaling is proof of that).
Good luck!
Jason
> Best regards,
>
> Ilyas Yildirim, Ph.D.
> -----------------------------------------------------------
> = Department of Chemistry - 2145 Sheridan Road =
> = Northwestern University - Evanston, IL 60208 =
> = Ryan Hall #4035 (Nano Building) - Ph.: (847)467-4986 =
> = http://www.pas.rochester.edu/~yildirim/ =
> -----------------------------------------------------------
>
>
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>
--
Jason M. Swails
Quantum Theory Project,
University of Florida
Ph.D. Graduate Student
352-392-4032
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Mon Oct 04 2010 - 15:00:04 PDT
