Created
July 22, 2010 18:45
-
-
Save gilleain/486403 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
package org.openscience.cdk.tools; | |
import javax.vecmath.Point3d; | |
import javax.vecmath.Vector3d; | |
import org.openscience.cdk.interfaces.IAtom; | |
/** | |
* Methods 'borrowed' (stolen) from Jmol to determine or check the stereo | |
* class of a set of atoms. | |
* | |
* @author maclean | |
* | |
*/ | |
public class StereoTool { | |
/** | |
* Currently unused, but intended for the StereoTool to indicate what it | |
* 'means' by an assignment of some atoms to a class. | |
* | |
*/ | |
public enum StereoClass { TETRAHEDRAL, SQUARE_PLANAR, | |
TRIGONAL_BIPYRAMIDAL, OCTAHEDRAL } | |
/** | |
* The handedness of a tetrahedron, in terms of the point-plane distance | |
* of three of the corners, compared to the fourth. | |
* | |
* PLUS indices a positive point-plane distance, | |
* MINUS is a negative point-plane distance. | |
*/ | |
public enum TetrahedralSign { PLUS, MINUS } | |
/** | |
* The shape that four atoms take in a plane. | |
*/ | |
public enum SquarePlanarShape { U_SHAPE, FOUR_SHAPE, Z_SHAPE } | |
/** | |
* The maximum angle in radians for two lines to be 'diaxial'. | |
* Where 0.95 is about 172 degrees. | |
*/ | |
public static final double MAX_AXIS_ANGLE = 0.95; | |
/** | |
* Checks these four atoms for square planarity. | |
* | |
* @param atomA | |
* @param atomB | |
* @param atomC | |
* @param atomD | |
* @return | |
*/ | |
public static boolean isSquarePlanar( | |
IAtom atomA, IAtom atomB, IAtom atomC, IAtom atomD) { | |
Vector3d vNorm1 = new Vector3d(); | |
Vector3d vNorm2 = new Vector3d(); | |
Vector3d vNorm3 = new Vector3d(); | |
StereoTool.getPlaneNormals( | |
atomA, atomB, atomC, atomD, vNorm1, vNorm2, vNorm3); | |
// vNorm1 vNorm2 vNorm3 are right-hand normals for the given | |
// triangles | |
// 1-2-3, 2-3-4, 3-4-1 | |
// sp1 up up up U-shaped | |
// sp2 up up DOWN 4-shaped | |
// sp3 up DOWN DOWN Z-shaped | |
double norm12 = vNorm1.dot(vNorm2); | |
double norm23 = vNorm2.dot(vNorm3); | |
// XXX this is certainly wrong... (as in, a bad translation from Jmol) | |
return norm12 < 0 || norm23 < 0; | |
} | |
/** | |
* Checks these 7 atoms to see if they are at the points of an octahedron. | |
* | |
* @param atomA | |
* @param atomB | |
* @param atomC | |
* @param atomD | |
* @param atomE | |
* @param atomF | |
* @param atomG | |
* @return | |
*/ | |
public static boolean isOctahedral(IAtom atomA, IAtom atomB, IAtom atomC, | |
IAtom atomD, IAtom atomE, IAtom atomF, IAtom atomG) { | |
boolean isDiaxialAAGB = | |
isDiaxial(atomA, atomA, atomG, atomB, StereoTool.MAX_AXIS_ANGLE); | |
if (!isDiaxialAAGB) return false; // XXX ? | |
Vector3d vNorm1 = new Vector3d(); | |
Vector3d vNorm2 = new Vector3d(); | |
Vector3d vNorm3 = new Vector3d(); | |
StereoTool.getPlaneNormals( | |
atomC, atomD, atomE, atomF, vNorm1, vNorm2, vNorm3); | |
// XXX ? | |
if ((vNorm1.dot(vNorm2) < 0) || vNorm2.dot(vNorm3) < 0) return false; | |
// now check rotation in relation to the first atom | |
vNorm2.set(atomA.getPoint3d()); | |
vNorm2.sub(atomB.getPoint3d()); | |
return vNorm1.dot(vNorm2) < 0; // XXX ? | |
} | |
/** | |
* Checks these 6 atoms to see if they form a triangular-bipyramidal shape. | |
* | |
* @param atomA | |
* @param atomB | |
* @param atomC | |
* @param atomD | |
* @param atomE | |
* @param atomF | |
* @return | |
*/ | |
public static boolean isTrigonalBipyramidal(IAtom atomA, IAtom atomB, | |
IAtom atomC, IAtom atomD, IAtom atomE, IAtom atomF) { | |
boolean isDiaxialAAFB = StereoTool.isDiaxial( | |
atomA, atomA, atomF, atomB, StereoTool.MAX_AXIS_ANGLE); | |
if (isDiaxialAAFB) { | |
TetrahedralSign handednessCDEB = | |
StereoTool.getHandedness(atomC, atomD, atomE, atomB); | |
return handednessCDEB == TetrahedralSign.PLUS; // XXX ? not sure | |
} else { | |
return false; | |
} | |
} | |
/** | |
* Gets the tetrahedral handedness of four atoms. | |
* | |
* @param atomA | |
* @param atomB | |
* @param atomC | |
* @param atomD | |
* @return | |
*/ | |
public static TetrahedralSign getHandedness( | |
IAtom atomA, IAtom atomB, IAtom atomC, IAtom atomD) { | |
Point3d pointA = atomA.getPoint3d(); | |
Point3d pointB = atomB.getPoint3d(); | |
Point3d pointC = atomC.getPoint3d(); | |
Point3d pointD = atomD.getPoint3d(); | |
// these are re-used in the Jmol code by storing them | |
// in a struct-like inner class : 'SmilesSearch.VTemp' | |
Vector3d normal = new Vector3d(); | |
Vector3d vAB = new Vector3d(); | |
Vector3d vAC = new Vector3d(); | |
// XXX - is this name right? the dot product of two vectors | |
// is the angle between them, isn't it? | |
double normalAngle = StereoTool.getNormalThroughPoints( | |
pointA, pointB, pointC, normal, vAB, vAC); | |
double distance = distanceToPlane(normal, normalAngle, pointD); | |
// the point-plane distance is the absolute value, | |
// the sign of the distance gives the side of the plane the point is on | |
if (distance > 0) { | |
return TetrahedralSign.PLUS; | |
} else { | |
return TetrahedralSign.MINUS; | |
} | |
} | |
private static boolean isDiaxial( | |
IAtom atomA, IAtom atomB, IAtom atomC, IAtom atomD, double maxAngle) { | |
Point3d pointA = atomA.getPoint3d(); | |
Point3d pointB = atomB.getPoint3d(); | |
Point3d pointC = atomC.getPoint3d(); | |
Point3d pointD = atomD.getPoint3d(); | |
return isDiaxial(pointA, pointB, pointC, pointD, maxAngle); | |
} | |
// side-effect free version \o/ | |
private static boolean isDiaxial(Point3d pointA, Point3d pointB, | |
Point3d pointC, Point3d pointD, double maxAngle) { | |
// we don't care about these vectors... | |
Vector3d vAC = new Vector3d(); | |
Vector3d vBD = new Vector3d(); | |
return isDiaxial(pointA, pointB, pointC, pointD, vAC, vBD, maxAngle); | |
} | |
// NOTE : taken from Jmol class o.j.smiles.SmilesSearch | |
// I think that it is finding the angle between the lines A-C and B-D, and | |
// checking that angle is less than f... | |
// XXX the side-effect (boo!) of the method is to create | |
// the vectors corresponding to the lines | |
private static boolean isDiaxial(Point3d pointA, Point3d pointB, | |
Point3d pointC, Point3d pointD, Vector3d vAC, Vector3d vBD, | |
double maxAngle) { | |
vAC.set(pointA); | |
vBD.set(pointB); | |
vAC.sub(pointC); | |
vBD.sub(pointD); | |
// -0.95f about 172 degrees | |
return vAC.dot(vBD) < maxAngle; | |
} | |
// NOTE : taken from Jmol class o.j.smiles.SmilesAromatic | |
private static double distanceToPlane( | |
Vector3d normal, double w, Point3d point) { | |
if (normal == null) return Double.NaN; | |
double np = (normal.x * point.x) | |
+ (normal.y * point.y) | |
+ (normal.z * point.z) + w; | |
double nSQ = Math.sqrt((normal.x * normal.x) | |
+ (normal.y * normal.y) + (normal.z * normal.z)); | |
return np / nSQ; | |
} | |
// NOTE : taken from Jmol class o.j.smiles.SmilesAromatic | |
private static double getNormalThroughPoints( | |
Point3d pointA, Point3d pointB, Point3d pointC, | |
Vector3d vNorm, Vector3d vAB, Vector3d vAC) { | |
vAB.sub(pointB, pointA); | |
vAC.sub(pointC, pointA); | |
vNorm.cross(vAB, vAC); | |
vNorm.normalize(); | |
// COMMENT FROM JMOL | |
// ax + by + cz + d = 0 | |
// so if a point is in the plane, then N dot X = -d | |
vAB.set(pointA); | |
return -vAB.dot(vNorm); | |
} | |
private static void getPlaneNormals(IAtom atomA, IAtom atomB, IAtom atomC, | |
IAtom atomD, Vector3d vectorA, Vector3d vectorB, Vector3d vectorC) { | |
Point3d pointA = atomA.getPoint3d(); | |
Point3d pointB = atomB.getPoint3d(); | |
Point3d pointC = atomC.getPoint3d(); | |
Point3d pointD = atomD.getPoint3d(); | |
// these are temporary vectors | |
Vector3d vectorD = new Vector3d(); | |
Vector3d vectorE = new Vector3d(); | |
StereoTool.getNormalThroughPoints( | |
pointA, pointB, pointC, vectorA, vectorD, vectorE); | |
StereoTool.getNormalThroughPoints( | |
pointB, pointC, pointD, vectorB, vectorD, vectorE); | |
StereoTool.getNormalThroughPoints( | |
pointC, pointD, pointA, vectorC, vectorD, vectorE); | |
} | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Now in its own branch. Will add copyright header - but it is changing more and more from the original...