Logo Search packages:      
Sourcecode: cdk version File versions  Download package

MMFF94EnergyFunction.java

package org.openscience.cdk.modeling.forcefield;

import java.util.Hashtable;

import javax.vecmath.GMatrix;
import javax.vecmath.GVector;

import org.openscience.cdk.interfaces.IAtomContainer;
//import org.openscience.cdk.tools.LoggingTool;


/**
 * MMFF94 energy function.
 *
 * @author     vlabarta
 * @cdk.created    March 14, 2005
 * @cdk.module     forcefield
 *
 * @cdk.keyword MMFF94
 * @cdk.keyword force field
 */
00022 public class MMFF94EnergyFunction implements IPotentialFunction {
      String energyFunctionShape = " MMFF94 energy ";
      double energy = 0;
      GVector energyGradient = null;
      GVector order2ndErrorApproximateEnergyGradient = new GVector(3);
      GVector order5thErrorApproximateEnergyGradient = new GVector(3);
      GMatrix energyHessian = null;
      double[] forHessian = null;
      
      //private LoggingTool logger;
      
      BondStretching bs = new BondStretching();
      AngleBending ab = new AngleBending();
      StretchBendInteractions sbi = new StretchBendInteractions();
      Torsions t =new Torsions();
      VanDerWaalsInteractions vdwi = new VanDerWaalsInteractions();
      ElectrostaticInteractions ei = new ElectrostaticInteractions();
      int functionEvaluationNumber = 0;

      /**
       *  Constructor for the MMFF94EnergyFunction object
       *
       */
00045       public MMFF94EnergyFunction(IAtomContainer molecule, Hashtable mmff94Tables) throws Exception {
            //logger.debug("MMFF94EnergyFunction Constructor");
            bs.setMMFF94BondStretchingParameters(molecule, mmff94Tables);
            ab.setMMFF94AngleBendingParameters(molecule, mmff94Tables,true);
            sbi.setMMFF94StretchBendParameters(molecule, mmff94Tables,false);
            t.setMMFF94TorsionsParameters(molecule, mmff94Tables);
            vdwi.setMMFF94VanDerWaalsParameters(molecule, mmff94Tables);        
            ei.setMMFF94ElectrostaticParameters(molecule, mmff94Tables);
            
            //logger = new LoggingTool(this);
      }


      /**
       *  Evaluate the MMFF94 energy function for a given molecule
       *
       *@param  molecule  Current molecule.
       *@return        MMFF94 energy function value.
       */
00064       public double energyFunctionOfAMolecule(IAtomContainer molecule) {
            
            GVector coords3d = ForceFieldTools.getCoordinates3xNVector(molecule);

            vdwi.setFunctionMMFF94SumEvdW(coords3d);
            sbi.setFunctionMMFF94SumEBA(coords3d);
            
            energy = bs.functionMMFF94SumEB(coords3d)
                  + ab.functionMMFF94SumEA(coords3d)
                  + sbi.getFunctionMMFF94SumEBA()
                  + t.functionMMFF94SumET(coords3d)
                  + vdwi.getFunctionMMFF94SumEvdW()
                  + ei.functionMMFF94SumEQ(coords3d);
            
            return energy;
      }


      /**
       *  Evaluate the MMFF94 energy function for a given 3xN point
       *
       *@param  coords3d  Current molecule 3xN coordinates.
       *@return        MMFF94 energy function value.
       */
00088       public double energyFunction(GVector coords3d) {

            vdwi.setFunctionMMFF94SumEvdW(coords3d);
            sbi.setFunctionMMFF94SumEBA(coords3d);
            
            energy = bs.functionMMFF94SumEB(coords3d)
                  + ab.functionMMFF94SumEA(coords3d)
                  + sbi.getFunctionMMFF94SumEBA()
                  + t.functionMMFF94SumET(coords3d)
                  + vdwi.getFunctionMMFF94SumEvdW()
                  + ei.functionMMFF94SumEQ(coords3d);
            
            //logger.debug("energy = " + bs.mmff94SumEB + "  " + ab.mmff94SumEA + "  " + sbi.getFunctionMMFF94SumEBA() + "  " + t.mmff94SumET + "  " + vdwi.getFunctionMMFF94SumEvdW() + "  " + ei.mmff94SumEQ);
            
            return energy;
      }


      /**
       *  Evaluate the gradient for the MMFF94 energy function in a given 3xN point
       *
       *@param  coords3d  Current molecule coordinates.
       */
00111       public void setEnergyGradient(GVector coords3d) {
            //setOrder2ndErrorApproximateEnergyGradient(coords3d);
            
            //logger.debug("coords3d : " + coords3d);
            energyGradient = new GVector(coords3d.getSize());
            
            bs.setGradientMMFF94SumEB(coords3d);
            ab.set2ndOrderErrorApproximateGradientMMFF94SumEA(coords3d);
            //ab.set5thOrderErrorApproximateGradientMMFF94SumEA(coords3d);
            //sbi.setGradientMMFF94SumEBA(coords3d);
            sbi.set2ndOrderErrorApproximateGradientMMFF94SumEBA(coords3d);
            t.set2ndOrderErrorApproximateGradientMMFF94SumET(coords3d);
            //t.set5thOrderErrorApproximateGradientMMFF94SumET(coords3d);
            vdwi.setGradientMMFF94SumEvdW(coords3d);
            ei.setGradientMMFF94SumEQ(coords3d);
            
            //logger.debug("bs.getGradientMMFF94SumEB() = " + bs.getGradientMMFF94SumEB());
            //logger.debug("ab.getGradientMMFF94SumEA() = " + ab.getGradientMMFF94SumEA());
            
            for (int i=0; i < energyGradient.getSize(); i++) {
                  energyGradient.setElement(i, 
                        bs.getGradientMMFF94SumEB().getElement(i) 
                        + ab.get2ndOrderErrorApproximateGradientMMFF94SumEA().getElement(i)
                        //+ sbi.getGradientMMFF94SumEBA().getElement(i)
                        + sbi.get2ndOrderErrorApproximateGradientMMFF94SumEBA().getElement(i)
                        + t.get2ndOrderErrorApproximateGradientMMFF94SumET().getElement(i)
                        + vdwi.getGradientMMFF94SumEvdW().getElement(i)
                        + ei.getGradientMMFF94SumEQ().getElement(i)
                        );
            }
      }


      /**
       *  Get the gradient for the MMFF94 energy function in a given 3xN point
       *
       *@return        MMFF94 energy gradient value
       */
00149       public GVector getEnergyGradient() {
            return energyGradient;
            //return order2ndErrorApproximateEnergyGradient;
      }


      /**
       *  Evaluate the 2nd order error approximate gradient for the MMFF94 energy function in a given 3xN point.
       *
       *@param  coords3d  Current molecule coordinates.
       */
00160       public void setOrder2ndErrorApproximateEnergyGradient(GVector coords3d) {
            //logger.debug("coords3d : " + coords3d);
            order2ndErrorApproximateEnergyGradient.setSize(coords3d.getSize());
            double sigma = Math.pow(0.0000000000000001,0.3333);
            GVector xplusSigma = new GVector(coords3d.getSize());
            GVector xminusSigma = new GVector(coords3d.getSize());
            
            for (int i=0; i < order2ndErrorApproximateEnergyGradient.getSize(); i++) {
                  xplusSigma.set(coords3d);
                  xplusSigma.setElement(i,coords3d.getElement(i) + sigma);
                  xminusSigma.set(coords3d);
                  xminusSigma.setElement(i,coords3d.getElement(i) - sigma);
                  order2ndErrorApproximateEnergyGradient.setElement(i, (energyFunction(xplusSigma) - energyFunction(xminusSigma)) / (2 * sigma));
            }
            //logger.debug("order2ndErrorApproximateEnergyGradient : " + order2ndErrorApproximateEnergyGradient);
      }


      /**
       *  Get the 2nd order error approximate gradient for the MMFF94 energy function.
       *
       *@return        2nd order error approximate MMFF94 energy gradient value
       */
00183       public GVector getOrder2ndErrorApproximateEnergyGradient() {
            return order2ndErrorApproximateEnergyGradient;
      }


      /**
       *  Evaluate the 5th order error approximate gradient for the MMFF94 energy function, given a 3xN point
       *
       *@param  coords3d  Current molecule coordinates.
       */
00193       public void setOrder5thErrorApproximateEnergyGradient(GVector coords3d) {
            //logger.debug("coords3d : " + coords3d);
            order5thErrorApproximateEnergyGradient.setSize(coords3d.getSize());
            double sigma = Math.pow(0.0000000000000001,0.2);
            GVector xplusSigma = new GVector(coords3d.getSize());
            GVector xminusSigma = new GVector(coords3d.getSize());
            GVector xplus2Sigma = new GVector(coords3d.getSize());
            GVector xminus2Sigma = new GVector(coords3d.getSize());
            
            for (int i=0; i < order5thErrorApproximateEnergyGradient.getSize(); i++) {
                  xplusSigma.set(coords3d);
                  xplusSigma.setElement(i,coords3d.getElement(i) + sigma);
                  xminusSigma.set(coords3d);
                  xminusSigma.setElement(i,coords3d.getElement(i) - sigma);
                  xplus2Sigma.set(coords3d);
                  xplus2Sigma.setElement(i,coords3d.getElement(i) + 2 * sigma);
                  xminus2Sigma.set(coords3d);
                  xminus2Sigma.setElement(i,coords3d.getElement(i) - 2 * sigma);
                  order5thErrorApproximateEnergyGradient.setElement(i, (8 * (energyFunction(xplusSigma) - energyFunction(xminusSigma)) - (energyFunction(xplus2Sigma) - energyFunction(xminus2Sigma))) / (12 * sigma));
            }
                  
            //logger.debug("order5thErrorApproximateEnergyGradient : " + order5thErrorApproximateEnergyGradient);
      }


      /**
       *  Get the 5th order error approximate gradient for the MMFF94 energy function.
       *
       *@return        5th order error approximate MMFF94 energy gradient value.
       */
00223       public GVector getOrder5thErrorApproximateEnergyGradient() {
            return order5thErrorApproximateEnergyGradient;
      }


      /**
       *  Evaluate the hessian for the MMFF94 energy function in a given 3xN point
       *
       *@param  coords3d  Current molecule coordinates.
       */
00233       public void setEnergyHessian(GVector coords3d) {
            
            forHessian = new double[coords3d.getSize() * coords3d.getSize()];
            
            bs.setHessianMMFF94SumEB(coords3d);
            //bs.set2ndOrderErrorApproximateHessianMMFF94SumEB(coords3d);
            //ab.setHessianMMFF94SumEA(coords3d);
            ab.set2ndOrderErrorApproximateHessianMMFF94SumEA(coords3d);
            sbi.setHessianMMFF94SumEBA(coords3d);
            //t.setHessianMMFF94SumET(coords3d);
            t.set2ndOrderErrorApproximateHessianMMFF94SumET(coords3d);
            vdwi.setHessianMMFF94SumEvdW(coords3d);
            ei.setHessianMMFF94SumEQ(coords3d);
            
            for (int i = 0; i < coords3d.getSize(); i++) {
                  for (int j = 0; j < coords3d.getSize(); j++) {
                        forHessian[i*coords3d.getSize()+j] = 
                              bs.getHessianMMFF94SumEB().getElement(i,j)
                              //bs.get2ndOrderErrorApproximateHessianMMFF94SumEB().getElement(i,j);
                              //+ ab.getHessianMMFF94SumEA().getElement(i,j)  //not working
                              + ab.get2ndOrderErrorApproximateHessianMMFF94SumEA().getElement(i,j)
                              + sbi.getHessianMMFF94SumEBA().getElement(i,j)
                              //+ t.getHessianMMFF94SumET().getElement(i,j)   //not working
                              + t.get2ndOrderErrorApproximateHessianMMFF94SumET().getElement(i,j)
                              + vdwi.getHessianMMFF94SumEvdW().getElement(i,j)
                              + ei.getHessianMMFF94SumEQ().getElement(i,j)
                              ;
                  }           
            }
            energyHessian = new GMatrix(coords3d.getSize(), coords3d.getSize(), forHessian);
      }


      /**
       *  Get the hessian for the MMFF94 energy function.
       *
       *@return        MMFF94 energy hessian value
       */
00271       public GMatrix getEnergyHessian() {
            return energyHessian;
      }


      /**
       *  Get the hessian for the MMFF94 energy function.
       *
       *@return        MMFF94 energy hessian value
       */
00281       public double[] getForEnergyHessian() {
            return forHessian;
      }

}

Generated by  Doxygen 1.6.0   Back to index