User 870ab5b546
06-08-2006 23:04:26
Here's part of a jsp page that takes a chair cyclohexane and flips its conformation.
Vector math procedures used by the code above are below. I stored them in the public class Dihedral along with slightly enhanced versions of Miklos' procedure for calculating a dihedral angle. (The enhancements allow the user to specify a 0 to 180 or -180 to 180 range and whether to return degrees or radians.)
| Code: |
| <%@ page language="java" %> <%@ page import="com.prenhall.epoch.chem.ChemUtils, com.prenhall.epoch.chem.Dihedral, chemaxon.struc.Molecule, chemaxon.struc.MolAtom, chemaxon.struc.MolBond, chemaxon.struc.SelectionMolecule, chemaxon.struc.DPoint3, chemaxon.struc.CTransform3D, chemaxon.struc.CGraph, chemaxon.sss.search.Search, chemaxon.sss.search.MolSearch, chemaxon.formats.MolImporter, java.math.*" %> <% String moleculeStr = request.getParameter("molecule"); if (moleculeStr != null && !moleculeStr.equals("")) { Molecule molecule = MolImporter.importMol(moleculeStr); MolSearch ourSearch = new MolSearch(); Molecule C6H12 = MolImporter.importMol("C1CCCCC1"); ourSearch.setTarget(molecule); ourSearch.setQuery(C6H12); ourSearch.setExactBondMatching(true); int[] ringIndices = ourSearch.findFirst(); if (ringIndices != null) { // a cyclohexane was found System.out.println("Found a cyclohexane."); Dihedral vector = new Dihedral(); vector.setMolecule(molecule); MolAtom[] ringAtoms = new MolAtom[6]; DPoint3[] ringAtomLocs = new DPoint3[6]; for (int i=0; i<6; i++) { ringAtoms[i] = molecule.getAtom(ringIndices[i]); ringAtomLocs[i] = ringAtoms[i].getLocation(); } // will rotate each ring atom and its ligands // store this info now before we start moving atoms DPoint3[] vecPrev2ToNext2 = new DPoint3[6]; double[] rotateAngle = new double[6]; for (int i=0; i<6; i++) { // rotation axis is along vector between 3rd and 5th atoms vecPrev2ToNext2[i] = vector.diff( ringAtomLocs[( (i+2>5) ? i-4 : i+2 )], ringAtomLocs[( (i-2<0) ? i+4 : i-2 )] ); // rotation angle is close to +/-60, but not exactly; // the amount that dihedral of ring bonds is short of 60 // is the amount that rotation should be more than 60; // sign of rotation is opposite to sign of dihedral rotateAngle[i] = vector.calcDihedral( ringAtoms[( (i-1<0) ? i+5 : i-1 )], ringAtoms[i], ringAtoms[( (i+1>5) ? i-5 : i+1 )], ringAtoms[( (i+2>5) ? i-4 : i+2 )], "radians", 360 ); rotateAngle[i] = ( (rotateAngle[i] > 0) ? (rotateAngle[i] - Math.PI)/2 : (rotateAngle[i] + Math.PI)/2 ); } // find vector connecting planes of alternating ring atoms // will translate alternating atoms up and down along this vector DPoint3 middleOdd = new DPoint3( (ringAtomLocs[1].x + ringAtomLocs[3].x + ringAtomLocs[5].x)/3, (ringAtomLocs[1].y + ringAtomLocs[3].y + ringAtomLocs[5].y)/3, (ringAtomLocs[1].z + ringAtomLocs[3].z + ringAtomLocs[5].z)/3 ); DPoint3 middleEven = new DPoint3( (ringAtomLocs[0].x + ringAtomLocs[2].x + ringAtomLocs[4].x)/3, (ringAtomLocs[0].y + ringAtomLocs[2].y + ringAtomLocs[4].y)/3, (ringAtomLocs[0].z + ringAtomLocs[2].z + ringAtomLocs[4].z)/3 ); DPoint3 vec1to3 = vector.diff(ringAtomLocs[3], ringAtomLocs[1]); DPoint3 vec1to5 = vector.diff(ringAtomLocs[5], ringAtomLocs[1]); DPoint3 normalOdd = vector.crossProd(vec1to3, vec1to5); DPoint3 vecMidEvenToMidOdd = vector.diff(middleEven, middleOdd); DPoint3 projMidEvenOnPlaneOdd = vector.proj1onPlaneNormalTo2(vecMidEvenToMidOdd, normalOdd); DPoint3 vecPlaneEvenToPlaneOdd = vector.diff(vecMidEvenToMidOdd, projMidEvenOnPlaneOdd); System.out.println("translate by " + vecPlaneEvenToPlaneOdd.x + ", " + vecPlaneEvenToPlaneOdd.y + ", " + vecPlaneEvenToPlaneOdd.z); CTransform3D moveRingAtom = new CTransform3D(); CTransform3D turnRingAtom = new CTransform3D(); for (int i=0; i<6; i++) { int prev = ( (i-1<0) ? i+5 : i-1); int next = ( (i+1>5) ? i-5 : i+1); System.out.println("processing atom " + (ringIndices[i]+1)); System.out.println("this atom location = " + ringAtomLocs[i].x + ", " + ringAtomLocs[i].y + ", " + ringAtomLocs[i].z); // remove ring bonds to allow selection of one fragment; // will restore soon MolBond prevBond = (MolBond) ringAtoms[prev].getEdgeTo(ringAtoms[i]); MolBond nextBond = (MolBond) ringAtoms[i].getEdgeTo(ringAtoms[next]); molecule.removeEdge(nextBond); molecule.removeEdge(prevBond); // find (newly) disconnected fragments SelectionMolecule[] fragments = molecule.findFrags(); // find the fragment j containing atom i CGraph atomGraph = new CGraph(); atomGraph.add(ringAtoms[i]); int j; for (j=0; j<fragments.length; j++) if (fragments[j].contains(atomGraph)) break; // rotate fragment j System.out.println("Turning by " + rotateAngle[i]*180/Math.PI + " degrees"); turnRingAtom.setRotation( vecPrev2ToNext2[i].x, vecPrev2ToNext2[i].y, vecPrev2ToNext2[i].z, rotateAngle[i]); turnRingAtom.setRotationCenter(ringAtomLocs[i]); fragments[j].transform(turnRingAtom); DPoint3 loc = ringAtoms[i].getLocation(); System.out.println(" location after rotation is " + loc.x + ", " + loc.y + ", " + loc.z); // translate fragment j DPoint3 vecTranslate = ( (i % 2 == 0) ? vector.scalarProd(vecPlaneEvenToPlaneOdd, -1) : vecPlaneEvenToPlaneOdd); moveRingAtom.setTranslation(vecTranslate); fragments[j].transform(moveRingAtom); loc = ringAtoms[i].getLocation(); System.out.println(" location after translation is " + loc.x + ", " + loc.y + ", " + loc.z); // restore the deleted bonds molecule.add(nextBond); molecule.add(prevBond); } // repeat for each ring atom moleculeStr = molecule.toFormat("mrv"); System.out.println(" flipped chair is:"); System.out.println(moleculeStr); } // a cyclohexane ring was found else moleculeStr = null; } // molecule is not null %> |
Vector math procedures used by the code above are below. I stored them in the public class Dihedral along with slightly enhanced versions of Miklos' procedure for calculating a dihedral angle. (The enhancements allow the user to specify a 0 to 180 or -180 to 180 range and whether to return degrees or radians.)
| Code: |
| public void setMolecule( Molecule m ) { mol = m; } // setMolecule /** * Calculates and returns the dihedral angle specified by the given atoms. * Atom indices refer to the <code>Molecule</code> object * set by the last call of <code>setMolecule()</code>. The angle is measured * around the bond <code>atom2</code>-<code>atom3</code> with respect to * <code>atom1</code> and <code>atom4</code>. Typically, <code>atom1</code> * is a neighbour of <code>atom2</code>, while <code>atom4</code> is a * neighbour of <code>atom3</code>. <br> * The angle calculated is the angle between planes <i>atom1 atom2 atom3</i> * and <i>atom2 atom3 atom4</i>. * @param atom1 neighbour of <code>atom2</code> * @param atom2 first atom of the bond for which dihedral is calculated * @param atom3 second atom of the bond for which dihedral is calculated * @param atom4 neighbour of <code>atom3</code> */ public double calcDihedral( MolAtom atom1, MolAtom atom2, MolAtom atom3, MolAtom atom4, String degOrRadian, int range ) { Simple3DVector p1 = new Simple3DVector( atom1.getX(), atom1.getY(), atom1.getZ() ); Simple3DVector p2 = new Simple3DVector( atom2.getX(), atom2.getY(), atom2.getZ() ); Simple3DVector p3 = new Simple3DVector( atom3.getX(), atom3.getY(), atom3.getZ() ); Simple3DVector p4 = new Simple3DVector( atom4.getX(), atom4.getY(), atom4.getZ() ); // p1-p2-p3-p4 and we need dihedral around the p2-p3 bond // this is the angle between two planes: (p1 p2 p3) and (p2 p3 p4) Simple3DVector v1 = p1.diff( p2 ); // points from p2 to p1 Simple3DVector v2 = p3.diff( p2 ); // points from p2 to p3 Simple3DVector n1 = v1.crossProd( v2 ); // perpendicular to (p1 p2 p3) plane Simple3DVector v3 = p2.diff( p3 ); // points from p3 to p2 Simple3DVector v4 = p4.diff( p3 ); // points from p3 to p4 Simple3DVector n2 = v3.crossProd( v4 ); // perpendicular to (p2 p3 p4) plane double da = n1.angle( n2 ); // angle between two vectors perpendicular to (p1 p2 p3) and (p2 p3 p4) // planes, respectively. angle between the two vectors is equal to the // angle between the two planes if (((Double) da).isNaN()) da = 180; // correction of NaN bug for coplanar planes da = ( ( (range == 360) && (n1.dotProd(v4) < 0) ) ? -da : da); // direction of angle return (degOrRadian.equals("radians") ? da*Math.PI/180 : da); // radians or degrees } // calcDihedral(4 atom indices, degrees/radians, 180/360) public double calcDihedral( int i1, int i2, int i3, int i4, String degOrRadian, int range ) { MolAtom atom1 = mol.getAtom( i1 ); MolAtom atom2 = mol.getAtom( i2 ); MolAtom atom3 = mol.getAtom( i3 ); MolAtom atom4 = mol.getAtom( i4 ); return calcDihedral( atom1, atom2, atom3, atom4, degOrRadian, range ); } // calcDihedral(4 atom indices, degrees/radians, 180/360) // if there's no fifth integer parameter, assume range of 0 to 180 is desired public double calcDihedral( MolAtom atom1, MolAtom atom2, MolAtom atom3, MolAtom atom4, String degOrRadian ) { return calcDihedral( atom1, atom2, atom3, atom4, degOrRadian, 180 ); } public double calcDihedral( int i1, int i2, int i3, int i4, String degOrRadian ) { return calcDihedral( i1, i2, i3, i4, degOrRadian, 180 ); } // if there's no string parameter, assume degrees is desired public double calcDihedral( MolAtom atom1, MolAtom atom2, MolAtom atom3, MolAtom atom4, int range ) { return calcDihedral( atom1, atom2, atom3, atom4, "degrees", range ); } public double calcDihedral( int i1, int i2, int i3, int i4, int range ) { return calcDihedral( i1, i2, i3, i4, "degrees", range ); } // if there's no string or fifth integer parameter, // assume degrees in range of 0 to 180 is desired public double calcDihedral( MolAtom atom1, MolAtom atom2, MolAtom atom3, MolAtom atom4 ) { return calcDihedral( atom1, atom2, atom3, atom4, "degrees", 180 ); } public double calcDihedral( int i1, int i2, int i3, int i4 ) { return calcDihedral( i1, i2, i3, i4, "degrees", 180 ); } /*** vector functions available for public use ***/ public DPoint3 crossProd ( DPoint3 vector1, DPoint3 vector2 ) { DPoint3 crossProduct = new DPoint3( vector1.y * vector2.z - vector1.z * vector2.y, -vector1.x * vector2.z + vector1.z * vector2.x, vector1.x * vector2.y - vector1.y * vector2.x ); return crossProduct; } // crossProd public double dotProd ( DPoint3 vector1, DPoint3 vector2 ) { return ( vector1.x * vector2.x + vector1.y * vector2.y + vector1.z * vector2.z ); } // dotProd public DPoint3 scalarProd ( DPoint3 vector, double scalar ) { DPoint3 scalarProduct = new DPoint3( vector.x * scalar, vector.y * scalar, vector.z * scalar ); return scalarProduct; } // scalarProd public DPoint3 scalarQuot ( DPoint3 vector, double scalar ) { DPoint3 scalarQuotient = new DPoint3( vector.x / scalar, vector.y / scalar, vector.z / scalar ); return scalarQuotient; } // scalarQuot public DPoint3 diff ( DPoint3 vector1, DPoint3 vector2 ) { DPoint3 vectorDifference = new DPoint3( vector1.x - vector2.x, vector1.y - vector2.y, vector1.z - vector2.z ); return vectorDifference; } // diff public double length ( DPoint3 vector ) { return Math.sqrt( Math.pow(vector.x, 2) + Math.pow(vector.y, 2) + Math.pow(vector.z, 2) ); } // length public double angle ( DPoint3 vector1, DPoint3 vector2 ) { return Math.acos( dotProd(vector1, vector2) / (length(vector1) * length(vector2)) ); } // angle // Generates a projection of vector1 on a plane perpendicular to // vector2. Formula from // http://www.euclideanspace.com/maths/geometry/elements/plane/lineOnPlane/index.htm public DPoint3 proj1onPlaneNormalTo2 ( DPoint3 vector1, DPoint3 vector2) { double length2 = length(vector2); DPoint3 unit2 = scalarQuot(vector2, length2); DPoint3 proj1 = scalarQuot(crossProd(vector2, crossProd(vector1, unit2)), length2); return proj1; } // proj1onPlaneNormalTo2 |