# Modify class2 dihedral style of Eaat coupling term

Dear All,

Hi. I was trying to modify class2 dihedral energy. (http://lammps.sandia.gov/doc/dihedral_class2.html) The original format of Eaat is written as M * (theta_ijk - theta_1) * (theta_jkl - theta_2) * cos(phi). I would like to rewrite it as M * (cos(theta_ijk) - cos(theta_1)) * (cos(theta_jkl) - cos(theta_2)) * cos(phi).

The original source code for Eaat of dihedral_class2.cpp is expressed as
da1 = acos(costh12) - aat_theta0_1[type];
da2 = acos(costh23) - aat_theta0_2[type];
if (eflag) edihedral += aat_k[type]da1da2*cosphi:
for energy term and

// for (i = 0; i < 4; i++)
// for (j = 0; j < 3; j++)
// fabcd[i][j] -= aat_k[type] * (cosphi * (da2dthetadr[i][j] - da1dthetadr[i][j]) + sinphi * da1da2dphidr[i][j]);

for force term.

Therefore I rewrite this term as
double cos_da1, cos_da2;
cos_da1 = costh12 - cos(aat_theta0_1[type]);
cos_da2 = costh23 - cos(aat_theta0_2[type]);
if (eflag) edihedral += aat_k[type]cos_da1cos_da2*cosphi;
and

// for (i = 0; i < 4; i++)
// for (j = 0; j < 3; j++)
// fabcd[i][j] -= aat_k[type] *
// sinphi * cos_da1
cos_da2
dphidr[i][j]);

I am not sure (How to validate?) whether the force term is correct because the sign(direction of the force) for different atom(angle) is defined differently according to the source code. For instance, the angle/torsion coupling term, the energy is the same on angle1 and angle2 but the force is in opposite direction. Thanks in advance.

Cheers,
FC

I am not sure (How to validate?) whether the force term is correct because the sign(direction of the force) for different atom(angle) is defined differently >according to the source code. For instance, the angle/torsion coupling term, the energy is the same on angle1 and angle2 but the force is in opposite >direction.

Not clear what sign issue you are referring to. However, to
validate a modified potential, you should make
sure the energy and force are consistent. E.g. if you
take a single small molecule with one of these angles
and do dynamics, is energy conserved? Or if you
create 2 snapshots of the molecule with one atom
displaced by a small dx (or dy or dz), and do a “run 0”
on both snapshots to see the energy and force, do you get
Fx = -dE/dx, etc.

Steve