*/
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);