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[0][i][j] - da1dthetadr[1][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] *
// (cosphi * (-cos_da2sqrt(1.0-costh12costh12)dthetadr[0][i][j] + cos_da1sqrt(1.0-costh23costh23)dthetadr[1][i][j]) +
// sinphi * cos_da1cos_da2dphidr[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