chair flip

User 870ab5b546

06-08-2006 23:04:26

Here's part of a jsp page that takes a chair cyclohexane and flips its conformation.





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