AMBER Archive (2008)

Subject: RE: AMBER: from dihe parameters to torsional energy

From: Ilyas Yildirim (yildirim_at_pas.rochester.edu)
Date: Wed Dec 17 2008 - 11:10:13 CST

> 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.

Ross,

I do not think that the 1st equation is used to do the energy
calculations of dihedrals. I did not check the source code for
this, but people published new dihedral parameters, which have phase
angles other than 0 or 180 (amber barcelone force field).

Here is a link to the frcmod file for the barcelone force field:

http://mmb.pcb.ub.es/PARMBSC0/frcmod.parmbsc0

If the above equation is used to do the energy calculations, then
this means that amber force field does not allow to have phase angles
other than 0 or 180.

The first equation is used in nmode/ephi.f file. I will be very surprised
if the amber torsional energies are calculated this way in sander/pmemd.
nmode is used for normal-mode analysis (entropy calculations). It looks
like that if the phase is something other than 0 or 180, nmode will not
calculate the right torsional energy.

On Wed, 17 Dec 2008, Ross Walker wrote:

> 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_at_scripps.edu [mailto:owner-amber_at_scripps.edu] On Behalf Of
> Jeffrey
> Sent: Wednesday, December 17, 2008 1:15 AM
> To: amber_at_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
>
>
>
>
>
> 

-- 
  Ilyas Yildirim, Ph.D.
  ---------------------------------------------------------------
  = Hutchison Hall B#10          - Department of Chemistry      =
  =                              - University of Rochester      =
  = 585-275-6766 (office)        -                              =
  = http://www.pas.rochester.edu/~yildirim/                     =
  ---------------------------------------------------------------

-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber_at_scripps.edu
To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
      to majordomo_at_scripps.edu