Dear Jeffrey,
My understanding of this is that the two equations are the same thing given
certain conditions with respect to phase.
The 'real' equation as written in terms of the AMBER force field is:
e = pk(1.0+cos(pn*phi-phase))
If my trigonometry is correct (it is early in the morning so buyer beware),
The cos term can be expanded here as:
e = pk(1.0+cos(pn*phi)*cos(-phase))
You would not want to do this computationally since cos functions are
expensive, however, in the case of phase = 0 then cos(-0) = 1 and if phase =
180 then cos(-180) = -1
Hence the above equation becomes:
e = pk(1.0(+/-)cos(pn*phi))
The net result though is that it is the same function as the original
definition above. Hence it should yield the same results. I don't know for
sure but I suspect that it is written this way either to make the
derivatives easier to express or compute or that the standard way of writing
it has some kind of instability / accuracy issue in certain circumstances,
i.e. when the dihedral gets close to either 0, 90 or 180.
Hence the short answer is that it shouldn't matter which route through the
code gets taken here - in both cases you should end up with the same answer.
If you want a list of dihedrals from the prmtop then you can run:
$AMBERHOME/exe/rdparm prmtop
and type 'dihedrals'
This will list everything from which given a set of coordinates and the
original equation at the top above you should end up with the same energy as
is calculated using the expanded terms above.
I hope that helps,
All the best
Ross
From: owner-amber.scripps.edu [mailto:owner-amber.scripps.edu] On Behalf Of
Jeffrey
Sent: Wednesday, December 17, 2008 1:15 AM
To: amber.scripps.edu
Subject: AMBER: from dihe parameters to torsional energy
Dear all,
I would like to figure out how dihedral parameters are used to
calculate the torsional energy. After checking the subroutine subroutine
ephi (src/nmode/ephi.f), I found two expressions were adopted for the
dihedral energy:
1>. e = pk(ic) * (1.0+phase*cos(pn(ic)*phi)
where phase = 1.0 or -1.0, and pn = 1,2,3,4, or 6
2>. e = pk*( 1.0+cos(pn*phi-phase) ----- the old energy form,
if phase angle is other than 0 or pi then assume this angle.
Taking CT-CT-CT-CT in parm99.dat as an example,
CT-CT-CT-CT 1 0.18 0.0 -3. Junmei et al,
1999
CT-CT-CT-CT 1 0.25 180.0 -2. Junmei et al,
1999
CT-CT-CT-CT 1 0.20 180.0 1. Junmei et al,
1999
should the total torsional energy for this dihedral be one of the two forms:
a. (pk1*(1 + phase*Cos[pn1*phi ]) + pk2*(1 + phase*Cos[pn2*phi ]) + pk3*(1
+ phase*Cos[pn3*phi]))
b. (pk1*(1 + Cos[pn1*phi - phase1]) + pk2*(1 + Cos[pn2*phi - phase2]) +
pk3*(1 + Cos[pn3*phi - phase3]))
where pk1 = 0.18; pk2 = 0.25; pk3 = 0.20; phase1 = 0; phase2 = Pi; phase3 =
Pi; pn1 = 1; pn2 = 2; pn3 = 3;
phase=1 or -1
Another question is how to output each energy term in AMBER to check if I
have assigned a correct form of the expression?
Any suggestion is greatly appreciated.
Thanks very much for your time.
-------
Jeffrey
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
to majordomo.scripps.edu
Received on Fri Dec 19 2008 - 01:09:32 PST
