package name.mjw.jamber;
import javax.vecmath.Vector3d;
import name.mjw.jamber.util.Torsion;
import name.mjw.jamber.util.TorsionGradients;
/**
* Dihedral interaction between four atoms. It is a function of the plane
* generated by the first three atoms and the plane generated by the last three
* atoms
*
* @author mjw
*
*/
public class Dihedral extends ForceFieldTerm {
/**
* i--j
* /
* k--l
*
* Atom i in dihedral term.
*/
private Atom i;
/**
* Atom j in dihedral term.
*/
private Atom j;
/**
* Atom k in dihedral term.
*/
private Atom k;
/**
* Atom l in dihedral term.
*/
private Atom l;
/**
* Barrier height, also known as PK
*/
private double barrierHeight = 0.0;
/**
* Phase, the phase shift angle in the torsional function, in degrees.
*/
private double phase = 0.0;
/**
* Number of repeats of the cosine function in the range -180 to 180 degrees
* Also known as PN
*/
private double periodicity = 0.0;
/**
* The factor by which the torsional barrier is divided. Also know as IDIVF
*/
private final double multiplicity = 1.0;
/**
* Is the dihedral made up of multiple terms?
*/
private boolean multiTermed = false;
/**
* Should the the 1-4 atom in this dihedral contribute to the Non Bonded Vdw
* and Electrostatic terms?
*
* This is to discriminate against some ring dihedrals that may be counted
* twice by the 14NB terms.
*/
private boolean contributeTo14NonBondedTerms = true;
/**
* Is this an improper dihedral?
*
*/
private boolean improper = false;
/**
* Dihedral (PK/IDIVF) * (1 + cos(PN*phi - PHASE))
*
* @param i
* Atom i in dihedral term.
* @param j
* Atom j in dihedral term.
* @param k
* Atom k in dihedral term.
* @param l
* Atom l in dihedral term.
* @param barrierHeight
* Barrier height, also known as PK
* @param phase
* Phase in degrees
* @param periodicity
* Number of repeats of the cosine function in the range -180 to
* 180 degrees. Also known as PN
*/
public Dihedral(Atom i, Atom j, Atom k, Atom l, Double barrierHeight,
Double phase, Double periodicity) {
setAtomI(i);
setAtomJ(j);
setAtomK(k);
setAtomL(l);
setBarrierHeight(barrierHeight);
setPhase(phase);
setPeriodicity(periodicity);
}
/**
* Sets atom i in the dihedral.
*/
private void setAtomI(Atom i) {
this.i = i;
}
/**
* Gets atom i in the dihedral.
*/
public Atom getAtomI() {
return i;
}
/**
* Sets atom j in the dihedral.
*/
private void setAtomJ(Atom j) {
this.j = j;
}
/**
* Gets atom j in the dihedral.
*/
public Atom getAtomJ() {
return j;
}
/**
* Sets atom k in the dihedral.
*/
private void setAtomK(Atom k) {
this.k = k;
}
/**
* Gets atom k in the dihedral.
*/
public Atom getAtomK() {
return k;
}
/**
* Sets atom l in the dihedral.
*/
private void setAtomL(Atom l) {
this.l = l;
}
/**
* Gets atom l in the dihedral.
*/
public Atom getAtomL() {
return l;
}
/**
* Sets the dihedral parameter barrier height (PK) in kcal/mol
*/
private void setBarrierHeight(double barrierHeight) {
this.barrierHeight = barrierHeight;
}
/**
* Returns the dihedral parameter barrier height (PK) in kcal/mol
*/
public double getBarrierHeight() {
return barrierHeight;
}
/**
* Sets the dihedral parameter phase (PHASE) in degrees
*/
private void setPhase(Double phase) {
this.phase = phase;
}
/**
* Returns the dihedral parameter phase (PHASE) in degrees
*
* @return dihedral parameter phase (PHASE) in degrees
*/
public double getPhase() {
return phase;
}
/**
* Sets the dihedral parameter periodicity (PN)
*/
private void setPeriodicity(Double periodicity) {
this.periodicity = periodicity;
}
/**
* Gets the dihedral parameter periodicity (PN)
*/
public double getPeriodicicity() {
return periodicity;
}
/**
* Gets the dihedral parameter multiplicity (IDIVF)
*/
public double getMultiplicity() {
return multiplicity;
}
/**
* Gets the dihedral angle (phi) in degrees
*/
public double getCurrentDihedralAngle() {
Torsion torsion = new Torsion();
return torsion.getAngle(i.getPosition(), j.getPosition(),
k.getPosition(), l.getPosition());
}
@Override
/**
* Returns the potential energy of the dihedral in kcal/mol
*/
public double getPotentialEnergy() {
// http://ambermd.org/doc8/amber8.pdf page 257
// Etors = ( PK / IDIVF ) * ( 1 + cos( PN * phi - PHASE) )
double phi = Math.toRadians(getCurrentDihedralAngle());
return (barrierHeight / multiplicity)
* (1 + Math.cos(periodicity * phi - phase));
}
/**
* Returns the current (scalar) gradient of the potential energy term. This
* is the analytical derivative of the potential energy at this point wrt to
* phi
*/
@Override
public double getAnalyticalGradient() {
double phi = Math.toRadians(getCurrentDihedralAngle());
return (barrierHeight / multiplicity)
* (periodicity * Math.sin(phase - (periodicity * phi)));
}
@Override
public void evaluateForce() {
Vector3d temp3d = new Vector3d();
temp3d.sub(j.getPosition(), k.getPosition());
double unitForce = getAnalyticalGradient() / temp3d.length();
Torsion gt = new Torsion();
TorsionGradients grad = gt.getGradients(i.getPosition(), j.getPosition(),
k.getPosition(), l.getPosition());
Vector3d tempI = grad.getI();
tempI.scale(unitForce);
Vector3d tempJ = grad.getJ();
tempJ.scale(unitForce);
Vector3d tempK = grad.getK();
tempK.scale(unitForce);
Vector3d tempL = grad.getL();
tempL.scale(unitForce);
i.subForce(tempI);
j.subForce(tempJ);
k.subForce(tempK);
l.subForce(tempL);
}
public String toString() {
StringBuilder result = new StringBuilder();
String NEW_LINE = System.getProperty("line.separator");
result.append(this.getClass().getName()).append(" Object {")
.append(NEW_LINE);
result.append(
String.format(" Atoms:\t\t\t\t\t\t%-2s-%-2s-%-2s-%-2s",
i.getName(), j.getName(), k.getName(), l.getName()))
.append(NEW_LINE);
result.append(
String.format(" AtomTypes:\t\t\t\t\t%-2s-%-2s-%-2s-%-2s",
i.getAMBERAtomType(), j.getAMBERAtomType(),
k.getAMBERAtomType(), l.getAMBERAtomType())).append(
NEW_LINE);
result.append(" barrierHeight: " + "\t\t\t\t")
.append(barrierHeight * 2).append(NEW_LINE);
result.append(" phase (radians): " + "\t\t\t\t").append(phase)
.append(NEW_LINE);
result.append(" periodicity: " + "\t\t\t\t\t").append(periodicity)
.append(NEW_LINE);
result.append(" current dihedral angle (degrees) " + "\t\t")
.append(this.getCurrentDihedralAngle()).append(NEW_LINE);
result.append(" current potential contribution: (kcal/mol) " + "\t")
.append(this.getPotentialEnergy()).append(NEW_LINE);
result.append(" contributeTo14NonBondedTerms:\t\t\t")
.append(this.isContributeTo14NonBondedTerms()).append(NEW_LINE);
result.append(" improper:\t\t\t\t\t").append(this.isImproper())
.append(NEW_LINE);
result.append("}");
return result.toString();
}
public void setMultiTermed() {
this.multiTermed = true;
}
public boolean isMultiTermed() {
return multiTermed;
}
public void setContributeTo14NonBondedTerms(
boolean contributeTo14NonBondedTerms) {
this.contributeTo14NonBondedTerms = contributeTo14NonBondedTerms;
}
public boolean isContributeTo14NonBondedTerms() {
return contributeTo14NonBondedTerms;
}
public void setImproper() {
this.improper = true;
}
public boolean isImproper() {
return improper;
}
}