set dihedral angle about a bond

This code allows you to specify the dihedral angle about any bond in a molecule.





The code comes from this jsp page, hence the getParameter() commands at the top, but you can easily adapt it to a Java page.





Requires a Molecule object or a molecule in a standard string format, four atom indices separated by commas or colons, and an angle in degrees.

Code:

String moleculeStr = request.getParameter("molecule");    String atomString = request.getParameter("atomString");    String dihedralSetStr = request.getParameter("dihedralSet");    boolean doCalculation = true;    Molecule molecule = null;    if (moleculeStr != null)       molecule = MolImporter.importMol(moleculeStr);    else doCalculation = false;    int[] atomIndex = new int[4];    if (atomString != null && atomString != "") {       String[] atoms = atomString.split(":");       if (atoms.length == 1) atoms = atomString.split(",");       if ((atoms.length == 4) && doCalculation) {          int maxIndex = molecule.getAtomCount() - 1;          for (int j=0; j<4; j++) {             atomIndex[j] = Integer.parseInt(atoms[j]) - 1;             if (atomIndex[j] > maxIndex) doCalculation = false;          }       } else doCalculation = false;    } else doCalculation = false; // no atoms entered           int dihedralSetDegrees = (dihedralSetStr == null || dihedralSetStr == "")       ? 0 : Integer.parseInt(dihedralSetStr);    double dihedralSetRadians = dihedralSetDegrees*Math.PI/180;    if (doCalculation) {       MolAtom originAtom = molecule.getAtom(atomIndex[1]);       MolAtom pivotAtom = molecule.getAtom(atomIndex[2]);       DPoint3 originLoc = originAtom.getLocation();       DPoint3 pivotLoc = pivotAtom.getLocation();       CTransform3D translate = new CTransform3D();       translate.setTranslation(          -originLoc.x,          -originLoc.y,          -originLoc.z          );       System.out.println("Before translation to put origin at bond midpoint: ");       System.out.println("   " +          originAtom.getSymbol()          + (atomIndex[1]+1) + " has coords " + originLoc.x + ","          + originLoc.y + "," + originLoc.z + ";"          );       System.out.println("   " +          pivotAtom.getSymbol()          + (atomIndex[2]+1) + " has coords " + pivotLoc.x + ","          + pivotLoc.y + "," + pivotLoc.z + "."          );       molecule.transform(translate);       originLoc = originAtom.getLocation();       pivotLoc = pivotAtom.getLocation();       System.out.println("After translation: ");       System.out.println("   " +          originAtom.getSymbol()          + (atomIndex[1]+1) + " has coords " + originLoc.x + ","          + originLoc.y + "," + originLoc.z + ";"          );       System.out.println("   " +          pivotAtom.getSymbol()          + (atomIndex[2]+1) + " has coords " + pivotLoc.x + ","          + pivotLoc.y + "," + pivotLoc.z + "."          );       // now rotate entire molecule       double lengthPivotOrigin = pivotLoc.distance(originLoc);       // make rotation vector average of pivot atom and vector along       // z axis with length same as pivot atom; then rotation angle of 180 degrees       // puts pivot atom on Z axis       DPoint3 vectorRotate = new DPoint3(          (0 + pivotLoc.x)/2,          (0 + pivotLoc.y)/2,          (lengthPivotOrigin + pivotLoc.z)/2          );       CTransform3D rotate = new CTransform3D();       rotate.setRotationCenter(originLoc);       rotate.setRotation(vectorRotate.x, vectorRotate.y, vectorRotate.z, Math.PI);       molecule.transform(rotate);       System.out.println("After rotation of whole molecule: ");       originLoc = originAtom.getLocation();       pivotLoc = pivotAtom.getLocation();       System.out.println("   " +          originAtom.getSymbol()          + (atomIndex[1]+1) + " has coords " + originLoc.x + ","          + originLoc.y + "," + originLoc.z + ";"          );       System.out.println("   " +          pivotAtom.getSymbol()          + (atomIndex[2]+1) + " has coords " + pivotLoc.x + ","          + pivotLoc.y + "," + pivotLoc.z + "."          );       DPoint3 atom4Loc = molecule.getAtom(atomIndex[3]).getLocation();       DPoint3 atom1Loc = molecule.getAtom(atomIndex[0]).getLocation();       // get current dihedral angle       double atom4angleOrig = Math.atan(atom4Loc.y/atom4Loc.x);       if (atom4Loc.x < 0) atom4angleOrig += Math.PI;       double atom1angleOrig = Math.atan(atom1Loc.y/atom1Loc.x);       if (atom1Loc.x < 0) atom1angleOrig += Math.PI;       double angleDiffOrig = atom4angleOrig - atom1angleOrig;       double moveAngle = dihedralSetRadians - angleDiffOrig;       for (int i=0; i<pivotAtom.getBondCount(); i++) {          MolAtom ligand = (MolAtom) pivotAtom.getLigand(i);          if (ligand != originAtom) {             DPoint3 ligLoc = ligand.getLocation();             System.out.println("   Before bond rotation, " +                ligand.getSymbol()                + (molecule.indexOf(ligand)+1) + " has coords " + ligLoc.x + ","                + ligLoc.y + "," + ligLoc.z + "."                );             double ligAngle = Math.atan(ligLoc.y/ligLoc.x);             double ligLength = Math.sqrt( Math.pow(ligLoc.x, 2) + Math.pow(ligLoc.y, 2) );             if (ligLoc.x < 0) ligAngle += Math.PI;             ligAngle += moveAngle;             DPoint3 newLoc = new DPoint3(                Math.cos(ligAngle)*ligLength,                Math.sin(ligAngle)*ligLength,                ligLoc.z                );             ligand.setLocation(newLoc);             ligLoc = ligand.getLocation();             System.out.println("   After bond rotation, " +                ligand.getSymbol()                + (molecule.indexOf(ligand)+1) + " has coords " + ligLoc.x + ","                + ligLoc.y + "," + ligLoc.z + "."                );          } // if not the origin       } // for each ligand of pivot atom       // now rotate entire molecule back       molecule.transform(rotate);       System.out.println("After rotation of whole molecule back: ");       originLoc = originAtom.getLocation();       pivotLoc = pivotAtom.getLocation();       System.out.println("   " +          originAtom.getSymbol()          + (atomIndex[1]+1) + " has coords " + originLoc.x + ","          + originLoc.y + "," + originLoc.z + ";"          );       System.out.println("   " +          pivotAtom.getSymbol()          + (atomIndex[2]+1) + " has coords " + pivotLoc.x + ","          + pivotLoc.y + "," + pivotLoc.z + "."          );       moleculeStr = molecule.toFormat("mrv");    }

The following code is better. The previous code worked well only when the fourth atom was not attached to anything other than the third atom. The following code works well for any structure.

Code:

String moleculeStr = request.getParameter("molecule");    String atomString = request.getParameter("atomString");    String dihedralSetStr = request.getParameter("dihedralSet");    boolean doCalculation = true;    Molecule molecule = null;    if (moleculeStr != null)       molecule = MolImporter.importMol(moleculeStr);    else doCalculation = false;    int[] atomIndex = new int[4];    if (atomString != null && atomString != "") {       String[] atoms = atomString.split(":");       if (atoms.length == 1) atoms = atomString.split(",");       if ((atoms.length == 4) && doCalculation) {          int maxIndex = molecule.getAtomCount() - 1;          for (int j=0; j<4; j++) {             atomIndex[j] = Integer.parseInt(atoms[j]) - 1;             if (atomIndex[j] > maxIndex) doCalculation = false;          }       } else doCalculation = false;    } else doCalculation = false; // no atoms entered           int dihedralSetDegrees = (dihedralSetStr == null || dihedralSetStr == "")       ? 0 : Integer.parseInt(dihedralSetStr);    double dihedralSetRadians = dihedralSetDegrees*Math.PI/180;    if (doCalculation) {       MolAtom originAtom = molecule.getAtom(atomIndex[1]);       MolAtom pivotAtom = molecule.getAtom(atomIndex[2]);       DPoint3 originLoc = originAtom.getLocation();       DPoint3 pivotLoc = pivotAtom.getLocation();       // to make the calculations easier, do translation       // and rotation to place 2-3 bond along the z axis       // set atom2 to origin       CTransform3D translate = new CTransform3D();       translate.setTranslation(          -originLoc.x,          -originLoc.y,          -originLoc.z          );       molecule.transform(translate);       // now rotate entire molecule;       // first three values of CTransform3D correspond to a vector around       // which the entire molecule rotates, like the axis of a top;       // fourth value is the angle by which the top is turned.       originLoc = originAtom.getLocation();       pivotLoc = pivotAtom.getLocation();       double lengthPivotOrigin = pivotLoc.distance(originLoc);       // make rotation vector average of pivot atom and vector along       // z axis with length same as pivot atom; then rotation angle of 180 degrees       // puts pivot atom on Z axis       DPoint3 vectorMolRotate = new DPoint3(          (0 + pivotLoc.x)/2,          (0 + pivotLoc.y)/2,          (lengthPivotOrigin + pivotLoc.z)/2          );       CTransform3D molRotate = new CTransform3D();       molRotate.setRotationCenter(originLoc);       molRotate.setRotation(          vectorMolRotate.x,          vectorMolRotate.y,          vectorMolRotate.z,          Math.PI);       molecule.transform(molRotate);       // calculate angle to rotate around bond       originLoc = originAtom.getLocation();       pivotLoc = pivotAtom.getLocation();       DPoint3 atom4Loc = molecule.getAtom(atomIndex[3]).getLocation();       DPoint3 atom1Loc = molecule.getAtom(atomIndex[0]).getLocation();       // get current dihedral angle       double atom4angleOrig = Math.atan(atom4Loc.y/atom4Loc.x);       if (atom4Loc.x < 0) atom4angleOrig += Math.PI;       double atom1angleOrig = Math.atan(atom1Loc.y/atom1Loc.x);       if (atom1Loc.x < 0) atom1angleOrig += Math.PI;       double angleDiffOrig = atom4angleOrig - atom1angleOrig;       double moveAngle = dihedralSetRadians - angleDiffOrig;       // set bond rotation values       originAtom = molecule.getAtom(atomIndex[1]);       pivotAtom = molecule.getAtom(atomIndex[2]);       originLoc = originAtom.getLocation();       pivotLoc = pivotAtom.getLocation();       // bond rotation vector is the 2-3 bond itself       DPoint3 vectorFragRotate = new DPoint3(          pivotLoc.x - originLoc.x,          pivotLoc.y - originLoc.y,          pivotLoc.z - originLoc.z          );       CTransform3D fragRotate = new CTransform3D();       fragRotate.setRotationCenter(originLoc);       fragRotate.setRotation(          vectorFragRotate.x,          vectorFragRotate.y,          vectorFragRotate.z,          moveAngle);       // remove rotatable bond to allow selection of one fragment;       // will restore soon       MolBond bond3To2 = (MolBond) pivotAtom.getEdgeTo(originAtom);       int bond3To2Type = bond3To2.getType();       molecule.removeEdge(bond3To2);       // find (newly) disconnected fragments       SelectionMolecule[] fragments = molecule.findFrags();       // restore the deleted bond       MolBond newBond3To2 = new MolBond(pivotAtom, originAtom, bond3To2Type);       molecule.add(newBond3To2);       // find the fragment containing atom3       CGraph pivotGraph = new CGraph();       pivotGraph.add(pivotAtom);       int i=0;       for (i=0; i<fragments.length; i++)          if (fragments[i].contains(pivotGraph)) break;       // now rotate the selected fragment by the calculated angle       fragments[i].transform(fragRotate);       // now rotate the whole molecule back to its original orientation;       // because original rotation was 180 degrees, a second application       // brings it right back where it started       molecule.transform(molRotate);       moleculeStr = molecule.toFormat("mrv");    }

This code may do more transformations than necessary, but the transformations make the trigonometry easier.

I didn't realize that Miklos had posted a fix to his dihedral angle calculator that accounts for the direction of the rotation. With this fix, we can write a much simpler algorithm to rotate about a bond. Here the direction of rotation is obtained by putting the fixed bond in the back and the rotating bond in the front. I also added a parameter to allow the user to set the dihedral angle or simply to rotate by a specified amount.

Code:

String moleculeStr = request.getParameter("molecule");    String atomString = request.getParameter("atomString");    String userAngleStr = request.getParameter("userAngle");    String setNotRotateStr = request.getParameter("setNotRotate");    boolean doCalculation = true;    boolean setNotRotate = true;    if (setNotRotateStr != null)       setNotRotate = setNotRotateStr.equals("0");    else doCalculation = false;    Molecule molecule = null;    if (moleculeStr != null)       molecule = MolImporter.importMol(moleculeStr);    else doCalculation = false;    int[] atomIndex = new int[4];    if (atomString != null && atomString != "") {       String[] atoms = atomString.split("-");       if (atoms.length == 1) atoms = atomString.split(":");       if (atoms.length == 1) atoms = atomString.split(",");       if ((atoms.length == 4) && doCalculation) {          int maxIndex = molecule.getAtomCount() - 1;          for (int j=0; j<4; j++) {             atomIndex[j] = Integer.parseInt(atoms[j]) - 1;             if (atomIndex[j] > maxIndex)  // bad atom index                doCalculation = false;          }  // for j       } else doCalculation = false;  // four atoms not entered    } else doCalculation = false; // no atoms entered           int userAngleDegrees = (userAngleStr == null || userAngleStr == "")       ? 0 : Integer.parseInt(userAngleStr);    double userAngleRadians = userAngleDegrees*Math.PI/180;    if (doCalculation) {       MolAtom atom2 = molecule.getAtom(atomIndex[1]);       MolAtom atom3 = molecule.getAtom(atomIndex[2]);       double moveAngle = 0;       if (setNotRotate) {          Dihedral calc = new Dihedral();          calc.setMolecule(molecule);          MolAtom atom1 = molecule.getAtom(atomIndex[0]);          MolAtom atom4 = molecule.getAtom(atomIndex[3]);          double initDihedralDeg =             calc.calcDihedral(atom1, atom2, atom3, atom4);          System.out.println("   current dihedral = "             + initDihedralDeg + " degrees.");          moveAngle = userAngleRadians - initDihedralDeg*Math.PI/180;          System.out.println("   moveAngle = "             + moveAngle*180/Math.PI + " degrees.");       } else moveAngle = userAngleRadians;       // set bond rotation values       // rotation vector is 3->2 bond       DPoint3 atom2Loc = atom2.getLocation();       DPoint3 atom3Loc = atom3.getLocation();       CTransform3D fragRotate = new CTransform3D();       fragRotate.setRotationCenter(atom2Loc);       fragRotate.setRotation(          atom2Loc.x - atom3Loc.x,          atom2Loc.y - atom3Loc.y,          atom2Loc.z - atom3Loc.z,          moveAngle          );       // remove rotatable bond to allow selection of one fragment; will restore soon       MolBond bond3To2 = (MolBond) atom3.getEdgeTo(atom2);       int bond3To2Type = bond3To2.getType();       molecule.removeEdge(bond3To2);       // find (newly) disconnected fragments       SelectionMolecule[] fragments = molecule.findFrags();       // restore the deleted bond       MolBond newBond3To2 = new MolBond(atom3, atom2, bond3To2Type);       molecule.add(newBond3To2);       // find the fragment containing atom 3 and rotate it       CGraph atom3Graph = new CGraph();       atom3Graph.add(atom3);       for (int i=0; i<fragments.length; i++)          if (fragments[i].contains(atom3Graph)) {             fragments[i].transform(fragRotate);             break;          } // if fragment contains atom3       moleculeStr = molecule.toFormat("mrv");    }  // if doCalculation

The lines,

Code:

fragRotate.setRotationCenter(atom2Loc);       fragRotate.setRotation(          atom2Loc.x - atom3Loc.x,          atom2Loc.y - atom3Loc.y,          atom2Loc.z - atom3Loc.z,          moveAngle          );

should be switched, or the rotation center may not be set properly, and bad results will be obtained.