package name.mjw.jamber.util;
import org.apache.log4j.Logger;
import javax.vecmath.Point3d;
import javax.vecmath.Vector3d;
/**
* Given the Cartesian coordinates of four points i, j, k and l, this method
* calculates cosine and sine of the torsion angle, phi, as well as gradient of
* phi with respect to Cartesian coordinates of i,j,k,l.
* <p>
* The dihedral angle is the dot product of the normal of the plane formed by
* points i,j,k and the normal of the plane formed by the points j,k,l.
* <p>
* This is an implementation of the method presented by Hendrik Bekker in
* Chapter 5 of his thesis:
* {@link http://dissertations.ub.rug.nl/FILES/faculties/science/1996/h.bekker/thesis.pdf}
* <p>
* This is the dot product definition of theta
* <p>
* An alternative form (cross product form) can be found starting from equation
* 2.39 in
* {@link ftp://ftp.dl.ac.uk/ccp5/DL_POLY/DL_POLY_CLASSIC/DOCUMENTS/USRMAN2.17.pdf}
*/
public class Torsion {
final Logger LOG = Logger.getLogger(Torsion.class);
/**
* Returns the angle of the Torsion.
* @param i Position of atom i
* @param j Position of atom j
* @param k Position of atom k
* @param l Position of atom l
* @return angle Torsion angle is in degrees
*/
public double getAngle(Point3d i, Point3d j, Point3d k, Point3d l) {
double cosphi;
double sinphi;
double phi;
Vector3d rij = new Vector3d();
Vector3d rkj = new Vector3d();
Vector3d rjk = new Vector3d();
Vector3d rlk = new Vector3d();
Vector3d rkl = new Vector3d();
Vector3d R = new Vector3d();
Vector3d S = new Vector3d();
Vector3d n = new Vector3d();
rij.sub(i, j);
LOG.debug(" rij is " + rij);
LOG.debug(" rij length is " + rij.length());
rkj.sub(k, j);
LOG.debug(" rkj is " + rkj);
LOG.debug(" rkj length is " + rkj.length());
// equ. 5.3b and 5.3c, only this vector needs to be normalised
rkj.normalize();
rjk.sub(j, k);
LOG.debug(" rjk is " + rjk);
LOG.debug(" rjk length is " + rjk.length());
rlk.sub(l, k);
LOG.debug(" rlk is " + rlk);
LOG.debug(" rlk length is " + rlk.length());
rkl.sub(k, l);
LOG.debug(" rkl is " + rkl);
LOG.debug(" rkl length is " + rkl.length());
LOG.debug(" ");
// equ. 5.3b
rkj.scale(rij.dot(rkj));
R.sub(rij, rkj);
LOG.debug("R is " + R);
// Refresh vector
rkj.sub(k, j);
rkj.normalize();
// equ 5.3c
rkj.scale(rlk.dot(rkj));
S.sub(rlk, rkj);
LOG.debug("S is " + S);
// Normalise both these vectors
// equ. 5.3a
R.normalize();
S.normalize();
//TODO roll this into a private function?
cosphi = R.dot(S);
LOG.debug("cosphi is " + cosphi);
// When calculating phi, arccosine only returns values between 0 and 180
// degrees. If cosphi is 0.0 then phi could be either -90 degrees or +90
// degrees and we need to know the correct one when calculating the
// derivatives.
// This ambiguity is removed since we also know sinphi -
// hence if cosphi is 0.0 and sinphi is 1 then phi is +90 degrees.
// Conversely, if cosphi is 0.0 and sinphi is -1 then phi is -90
// degrees.
// equ. 5.2c ?
n.cross(R, S);
LOG.debug("n is " + n);
// equ. 5.28
sinphi = n.dot(rkj);
LOG.debug("sinphi is " + sinphi);
//TODO hack to prevent signum returning zero
phi = Math.acos(cosphi) * ((int) Math.signum(sinphi + 0.000000001) );
LOG.debug("phi is " + phi);
// Convert and return in degrees
LOG.debug("phi is " + Math.toDegrees(phi) + " degrees");
return Math.toDegrees(phi);
}
/**
* Returns a TorsionGradients
*
* @param i Position of atom i
* @param j Position of atom j
* @param k Position of atom k
* @param l Position of atom l
* @return TorsionGradients in units of kcal/(A mol)
*/
public TorsionGradients getGradients(Point3d i, Point3d j, Point3d k,
Point3d l) {
// Figure 5.1
Vector3d rij = new Vector3d();
Vector3d rkj = new Vector3d();
Vector3d rlk = new Vector3d();
Vector3d rkl = new Vector3d();
// Normal to plane ijk
Vector3d m = new Vector3d();
// Normal to plane jkl
Vector3d n = new Vector3d();
// Equ. 5.3a
Vector3d R = new Vector3d();
Vector3d S = new Vector3d();
rij.sub(i, j);
LOG.debug(" rij is " + rij);
LOG.debug(" rij length is " + rij.length());
rkj.sub(k, j);
LOG.debug(" rkj is " + rkj);
LOG.debug(" rkj length is " + rkj.length());
// Equ. 5.3b and 5.3c, only this vector needs to be normalised
rkj.normalize();
// Do we need this?
rlk.sub(l, k);
LOG.debug(" rlk is " + rlk);
LOG.debug(" rlk length is " + rlk.length());
// Equ. 5.3b
rkj.scale(rij.dot(rkj));
R.sub(rij, rkj);
LOG.debug("");
LOG.debug("R is " + R);
// Refresh vector since it has been rescaled in calculating R
rkj.sub(k, j);
rkj.normalize();
// Equ 5.3c
rkj.scale(rlk.dot(rkj));
S.sub(rlk, rkj);
LOG.debug("S is " + S);
// Normalise both these vectors
// Equ. 5.3a
R.normalize();
S.normalize();
// Now, the gradients
Vector3d gradPhiI = new Vector3d();
Vector3d gradPhiJ = new Vector3d();
Vector3d gradPhiK = new Vector3d();
Vector3d gradPhiL = new Vector3d();
Vector3d tempGradPhiL = new Vector3d();
TorsionGradients gradients = new TorsionGradients();
// Refresh vars
rij.sub(i, j);
rkj.sub(k, j);
rlk.sub(l, k);
rkl.sub(k, l);
// Equ. 5.2b; m is the normal to the plane ijk
m.cross(rij, rkj);
// m.normalize();
LOG.debug("");
LOG.debug("m (norm to ijk) is " + m);
LOG.debug("m.lengthSquared is " + m.lengthSquared());
// Equ. 5.2c; n is the normal to the plane jkl
n.cross(rlk, rkj);
// n.normalize();
LOG.debug("n (norm to jkl) is " + n);
LOG.debug("n.lengthSquared is " + n.lengthSquared());
// Equ. 5.11
gradPhiI.set(m);
LOG.debug("1.0 / m.length() is " + 1.0 / m.length());
gradPhiI.scale(1.0 / m.lengthSquared());
gradPhiI.scale(rkj.length());
gradients.setI(gradPhiI);
// Equ. 5.12
gradPhiL.set(n);
LOG.debug("1.0 / n.length() is " + 1.0 / n.length());
gradPhiL.scale(1.0 / n.lengthSquared());
gradPhiL.scale(rkj.length());
gradPhiL.negate();
gradients.setL(gradPhiL);
// Generate unknown vector, S
// Equ. 5.20
S.set(gradPhiI);
S.scale(rij.dot(rkj));
tempGradPhiL.set(gradPhiL);
tempGradPhiL.scale(rkl.dot(rkj));
S.sub(tempGradPhiL);
S.scale(1.0 / rkj.lengthSquared());
LOG.debug("S (unknown) is " + S);
// Equ. 5.13
gradPhiJ.set(gradients.getI());
gradPhiJ.negate();
gradPhiJ.add(S);
gradients.setJ(gradPhiJ);
// Equ. 5.14
gradPhiK.set(gradients.getL());
gradPhiK.negate();
gradPhiK.sub(S);
gradients.setK(gradPhiK);
return gradients;
}
}