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

GeometricMinimizer.java

package org.openscience.cdk.modeling.forcefield;

import org.openscience.cdk.interfaces.IAtomContainer;
import org.openscience.cdk.interfaces.IMolecule;
import org.openscience.cdk.modeling.builder3d.ForceFieldConfigurator;

import javax.vecmath.GMatrix;
import javax.vecmath.GVector;
import java.util.Hashtable;
//import org.openscience.cdk.tools.LoggingTool;


/**
 * Call the minimization methods. Check the convergence.
 *
 * @author      vlabarta
 * @cdk.module  forcefield
 *
 * @cdk.keyword geometry
 * @cdk.keyword 3D-coordinates
 */
00022 public class GeometricMinimizer {

      Hashtable PotentialParameterSet = null;
      int SDMaximumIteration = 1000;
      int CGMaximumIteration = 500;
      int NRMaximumIteration = 100;

      double CGconvergenceCriterion = 0.001;
      double SDconvergenceCriterion = 0.001;
      double NRconvergenceCriterion = 0.001;

    GVector kCoordinates = null;
      GVector kplus1Coordinates = null;
      GVector g0 = null;
      int atomNumbers = 1;
      
      double fxk = 0;
      double fxkplus1 = 0;

      GVector gradientk = null;
      GVector gradientkplus1 = null;

      GVector steepestDescentsMinimum = null;
      GVector conjugateGradientMinimum = null;
      GVector newtonRaphsonMinimum = null;
      
      double minimumFunctionValueCGM;

      int iterationNumberkplus1 = 0;
      int iterationNumberk = -1;
      boolean convergence = false;
      double RMSD = 0;
      double RMS = 0;
      double d = 0;
      boolean gradientSmallEnoughFlag = false;
      double infiniteNorm;
      
      NewtonRaphsonMethod nrm = new NewtonRaphsonMethod();
      //private LoggingTool logger;

      IMolecule molecule;


      /**
       *  Constructor for the GeometricMinimizer object
       */
00068       public GeometricMinimizer() {
            //logger = new LoggingTool(this);
      }


      public void setMolecule(IMolecule mol, boolean clone) throws Exception {

            if (clone) {
                  this.molecule = (IMolecule) mol.clone();
            } else {
                  this.molecule = mol;
            }
      }

      public IMolecule getMolecule() {
            return this.molecule;
      }


      /**
       *  Assign MMFF94 atom types to the molecule.
       *  
       *@param  molecule  The molecule like an AtomContainer object.
       */
00092       public void setMMFF94Tables(IAtomContainer molecule) throws Exception {
            //logger.debug("Start setMMFF94Tables");
            ForceFieldConfigurator ffc = new ForceFieldConfigurator();
            //logger.debug("setForceFieldConfigurator");
            ffc.setForceFieldConfigurator("mmff94");
            //logger.debug("assignAtomTyps");
            ffc.assignAtomTyps((IMolecule) molecule); // returns non-used RingSet
            //logger.debug("PotentialParameterSet");
            PotentialParameterSet = ffc.getParameterSet();
            //logger.debug("PotentialParameterSet = " + PotentialParameterSet);
      }
      
      public Hashtable getPotentialParameterSet() {
            return PotentialParameterSet;
      }

      public double[] getConvergenceParametersForSDM(){
            double[] parameters={SDMaximumIteration,SDconvergenceCriterion};
            return parameters;
      }
      
      public double[] getConvergenceParametersForCGM(){
            double[] parameters={CGMaximumIteration,CGconvergenceCriterion};
            return parameters;
      }
            
      public double[] getConvergenceParametersForNRM(){
            double[] parameters={NRMaximumIteration,NRconvergenceCriterion};
            return parameters;
      }
      
      public void initializeMinimizationParameters(GVector initialCoord) {
            //kCoordinates.setSize(dimension);
            //kCoordinates.set(initialCoord);
            kCoordinates=initialCoord;
            //logger.debug("Coordinates at iteration 1: X1 = " + kCoordinates);
            gradientk = new GVector(kCoordinates.getSize());
            kplus1Coordinates = new GVector(kCoordinates.getSize());
            gradientkplus1 = new GVector(kCoordinates.getSize());
            
            convergence = false;
            iterationNumberkplus1 = 0;
            iterationNumberk = -1;
            
            atomNumbers = initialCoord.getSize()/3;
            
            g0 = new GVector(kCoordinates.getSize());
            g0.zero();
      }


      /**
       *  Calculate the Root Mean Square gradient (RMS)
       *
       * @param  gradient  Gradient of the energy function in the new calculated point xk+1. 
       */
00148       public double rootMeanSquareGradient(GVector gradient) {
      
            RMS = 0;
            RMS = gradient.dot(gradient);
            RMS = RMS / gradient.getSize();
            RMS = Math.sqrt(RMS);
            //logger.debug("RMS = " + RMS);
            
            return RMS;
      }

      
      /**
       *  Analyse if the gradient is small enough, using the criteria || gradient || < 10-5 (1 + |function|)
       *
       * @param  function  energy function value. 
       * @param  gradient  Gradient of the energy function. 
       */
00166       public void gradientSmallEnough(double function, GVector gradient) {
            
            //logger.debug("Analyse if the gradient is small enough");

            infiniteNorm = 0;
            
            for (int i=0; i < gradient.getSize(); i++) {
            infiniteNorm = Math.max(infiniteNorm, Math.abs(gradient.getElement(i)));
            }
            
            //logger.debug("infiniteNorm = " + infiniteNorm);
            //logger.debug("0.00005  * (1 + Math.abs(function)) = " + (0.00005  * (1 + Math.abs(function))));
            
            if (infiniteNorm < 0.00005  * (1 + Math.abs(function))) {gradientSmallEnoughFlag = true;}
            
      }

      
      /**
       *  To check convergence
       */
00187       public void checkConvergence(double convergenceCriterion) {
            
            //logger.debug("Checking convergence : ");
            
            RMS = rootMeanSquareGradient(gradientkplus1);
            //logger.debug("RMS = " + RMS);
            
            if (RMS < convergenceCriterion) {
                  convergence = true;
                  //logger.debug("RMS convergence");
                  //logger.debug("RMS = " + RMS);
            }
            
            gradientSmallEnough(fxkplus1, gradientkplus1);
            
            if (gradientSmallEnoughFlag == true) {
                  convergence = true;
                  //logger.debug("Gradient Small Enough");
            }
            

            if (Math.abs(this.fxk - this.fxkplus1) < 0.0001) {
                  RMSD = 0;
                  for (int i = 0; i < atomNumbers; i++) {
                        d = ForceFieldTools.distanceBetweenTwoAtomFromTwo3xNCoordinates(kplus1Coordinates, kCoordinates, i, i);
                        RMSD = RMSD + Math.pow(d, 2);
                  }
                  RMSD = RMSD / kCoordinates.getSize();
                  RMSD = Math.sqrt(RMSD);
                  //logger.debug("RMSD = " + RMSD);

                  if (RMSD < 0.000001) {
                        convergence = true;
                        //logger.debug("RMSD convergence");
                  }
            }
      }


      public void setkplus1Coordinates(GVector direction, double stepSize) {
            kplus1Coordinates.set(direction);
            kplus1Coordinates.scale(stepSize);
            kplus1Coordinates.add(kCoordinates);
            return;
      }


      /**
       *  Set convergence parameters for Steepest Decents Method
       *
       * @param  changeSDMaximumIteration  Maximum number of iteration for steepest descents method. 
       * @param  changeSDConvergenceCriterion  Convergence criterion for steepest descents method.
       */
00240     public void setConvergenceParametersForSDM(int changeSDMaximumIteration, double changeSDConvergenceCriterion){
            SDMaximumIteration = changeSDMaximumIteration;
            SDconvergenceCriterion = changeSDConvergenceCriterion;
    }


      /**
       *  Minimize the potential energy function using steepest descents method
       *
       * @param  forceField         The potential function to be used
       */
00251     public void steepestDescentsMinimization(GVector initialCoordinates, IPotentialFunction forceField) {
            
            initializeMinimizationParameters(initialCoordinates);
            fxk = forceField.energyFunction(initialCoordinates);
            //logger.debug("STEEPESTDM: initial coords:"+initialCoordinates);
            //logger.debug("STEEPESTDM: kcoordinates coords:"+kCoordinates);
            //logger.debug("f(x0) = " + fxk);
                        
            forceField.setEnergyGradient(kCoordinates);
            gradientk.set(forceField.getEnergyGradient());
            //logger.debug("Initial gradient : g0 = " + gradientk);

            GVector sk = new GVector(gradientk.getSize());
            GVector skplus1 = new GVector(gradientk.getSize());
            
            double linearFunctionDerivativek =1;
            
            double alphaInitialStep = 0.5;

            if (gradientk.equals(g0)) {
                  convergence = true;
                  kplus1Coordinates.set(kCoordinates);
            }           

            //logger.debug("");
            //logger.debug("FORCEFIELDTESTS steepestDescentTest");
            
            SteepestDescentsMethod sdm = new SteepestDescentsMethod(kCoordinates);
            LineSearchForTheWolfeConditions ls = new LineSearchForTheWolfeConditions(forceField, "sdm");
                        
            while ((iterationNumberkplus1 < SDMaximumIteration) & (convergence == false)) {
                  
                  iterationNumberkplus1 += 1;
                  iterationNumberk += 1;
                  

//                if (iterationNumberkplus1 % 50 == 0 | iterationNumberkplus1 == 2) {
                        //logger.debug("");
                        //logger.debug("SD Iteration number: " + iterationNumberkplus1);
//                }
                  
                  //logger.debug("gm.steepestDescentsMinimisation, Energy Gradient:"+forceField.getEnergyGradient());
                  
                  if (iterationNumberkplus1 != 1) {
                        alphaInitialStep = ls.alphaOptimum * linearFunctionDerivativek;                                        
                        kCoordinates.set(kplus1Coordinates);
                        fxk = fxkplus1;
                        gradientk.set(gradientkplus1);
                        sk.set(skplus1);
                  }                 
                  //logger.debug("Search direction: ");
                  sdm.setSk(gradientk);
                  skplus1.set(sdm.sk);
                  linearFunctionDerivativek = gradientk.dot(skplus1);
                  if (iterationNumberkplus1 != 1) {
                        alphaInitialStep = alphaInitialStep / linearFunctionDerivativek;
                  }
                  
                  ls.initialize(kCoordinates, fxk, gradientk, sdm.sk, linearFunctionDerivativek, alphaInitialStep);
                  ls.lineSearchAlgorithm (5);
                  setkplus1Coordinates(sdm.sk, ls.alphaOptimum);
                  fxkplus1 = ls.linearFunctionInAlphaOptimum;                 
                  //logger.debug("x" + iterationNumberkplus1 + " = " + kplus1Coordinates + "    ");
                  //logger.debug("f(x" + iterationNumberkplus1 + ") = " + fxkplus1);

                  gradientkplus1.set(ls.dfOptimum);
                  
                  checkConvergence(SDconvergenceCriterion);

                  /*if (fxkplus1 <= fxk + 0.0001 * ls.alphaOptimum * linearFunctionDerivativek) {}
                  else {
                        System.out.println("Sufficient Decrease Condition not satisfied");
                        break;
                  }
                  if ((Math.abs(gradientkplus1.dot(sdm.sk)) <= -0.1 * linearFunctionDerivativek) | (ls.alphaOptimum == 5)) {}
                  else {
                        System.out.println("Curvature Condition not satisfied");
                        //logger.debug("linearFunctionDerivativekplus1 = " + gradientkplus1.dot(sdm.sk));
                        //logger.debug("linearFunctionDerivativek = " + linearFunctionDerivativek);
                  }*/
                              
                  //if (iterationNumberkplus1 % 50 == 0 | iterationNumberkplus1 == 1) {
                        //logger.debug("");
                        //logger.debug("f(x" + iterationNumberk + ") = " + fxk);
                        //logger.debug("f(x" + iterationNumberkplus1 + ") = " + fxkplus1);
                        //logger.debug("fxkplus1 - fxk = " + (fxkplus1 - fxk));
                        //logger.debug("gradientkplus1, gradientk angle = " + gradientkplus1.angle(gradientk));
                  //}

                  /*if (iterationNumberkplus1 != 1) {
                        logger.debug("sk+1.sk = " + skplus1.dot(sk));
                        logger.debug("gk+1.gk = " + gradientkplus1.dot(gradientk));
                  }*/
                  
            }
            steepestDescentsMinimum = kplus1Coordinates;
            //logger.debug("The minimum energy is " + fxkplus1);
            //logger.debug("SD Iteration number: " + iterationNumberkplus1);
            
            //forceField.setEnergyHessian(steepestDescentsMinimum);
            //NewtonRaphsonMethod nrm = new NewtonRaphsonMethod();
            //kCoordinates.set(kplus1Coordinates);
            //nrm.gradientPerInverseHessian(forceField.getEnergyGradient(),forceField.getEnergyHessian());
            //setkplus1Coordinates(nrm.getGradientPerInverseHessian(), -1);
            //logger.debug("x" + iterationNumberkplus1 + " = " + kplus1Coordinates);
            //logger.debug("f(x" + iterationNumberkplus1 + ") = " + fxkplus1);

            
            //logger.debug("The SD minimum energy is at: " + steepestDescentsMinimum);
            
            if (molecule != null){
                //logger.debug("STEEPESTDM: kplus1Coordinates:"+kplus1Coordinates);
                ForceFieldTools.assignCoordinatesToMolecule(kplus1Coordinates, molecule);
            }

            return;
      }


      /**
       *  Gets the steepestDescentsMinimum attribute of the GeometricMinimizer object
       *
       *@return    The minimumCoordinates value
       */
00375       public GVector getSteepestDescentsMinimum() {
            return steepestDescentsMinimum;
      }
      
      
      /**
       *  Set convergence parameters for Conjugate Gradient Method
       *
       * @param  changeCGMaximumIteration  Maximum number of iteration for conjugated gradient method
       * @param  changeCGConvergenceCriterion  Convergence criterion for conjugated gradient method
       */
00386       public void setConvergenceParametersForCGM(int changeCGMaximumIteration, double changeCGConvergenceCriterion){
            CGMaximumIteration = changeCGMaximumIteration;
            CGconvergenceCriterion = changeCGConvergenceCriterion;
      }


      /**
       *  Minimize the potential energy function using conjugate gradient method
       *
       * @param  forceField         The potential function to be used
       */
00397       public void conjugateGradientMinimization(GVector initialCoordinates, IPotentialFunction forceField) {
            //logger.debug("");
            //logger.debug("FORCEFIELDTESTS ConjugatedGradientTest");
            
            initializeMinimizationParameters(initialCoordinates);
            fxk = forceField.energyFunction(initialCoordinates);
            //logger.debug("f(x0) = " + fxk);
            
            forceField.setEnergyGradient(kCoordinates);
            gradientk = forceField.getEnergyGradient();
            //logger.debug("gradient at iteration 1 : g1 = " + gradientk);

            double linearFunctionDerivativek =1;
            
            double alphaInitialStep = 0.01;

            if (gradientk.equals(g0)) {
                  convergence = true;
                  kplus1Coordinates.set(kCoordinates);
            }           
            
            ConjugateGradientMethod cgm = new ConjugateGradientMethod();
            cgm.initialize(gradientk);
            LineSearchForTheWolfeConditions ls = new LineSearchForTheWolfeConditions(forceField, "cgm");
            
            gradientSmallEnough(fxk, gradientk);
            if (gradientSmallEnoughFlag == true) {
                  convergence = true;
                  //logger.debug("Gradient Small Enough");
            }
            
                   
            while ((iterationNumberkplus1 < CGMaximumIteration) & (convergence == false)) {
                  
                  iterationNumberkplus1 += 1;
                  iterationNumberk += 1;
                  //logger.debug("");
                  //logger.debug("");
                  //if (iterationNumberkplus1 % 50 == 0 | iterationNumberkplus1 == 2) {
                        //logger.debug("");
                        //logger.debug("CG Iteration number: " + iterationNumberkplus1);
                  //}
                  
                  if (iterationNumberkplus1 != 1) {
                        //logger.debug("ls.alphaOptimum = " + ls.alphaOptimum);
                        //logger.debug("linearFunctionDerivativek = " + linearFunctionDerivativek);
                        alphaInitialStep = ls.alphaOptimum * linearFunctionDerivativek;                                        
                        kCoordinates.set(kplus1Coordinates);
                        fxk = fxkplus1;
                        cgm.setDirection(gradientkplus1, gradientk);
                        gradientk.set(gradientkplus1);
                  }

                  //logger.debug("Search direction: ");
                  linearFunctionDerivativek = gradientk.dot(cgm.conjugatedGradientDirection);
                  if (iterationNumberkplus1 != 1) {
                        alphaInitialStep = alphaInitialStep / linearFunctionDerivativek;
                        //alphaInitialStep = Math.min(1.01 * alphaInitialStep,1);
                  }

                  //logger.debug("linearFunctionDerivativek = " + linearFunctionDerivativek);
                  //logger.debug("alphaInitialStep = " + alphaInitialStep);
                  ls.initialize(kCoordinates, fxk, gradientk, cgm.conjugatedGradientDirection, linearFunctionDerivativek, alphaInitialStep);
                  ls.lineSearchAlgorithm(5);
                  //if (ls.alphaOptimum == 0) {convergence = true;}
                  //logger.debug("ls.alphaOptimum = " + ls.alphaOptimum);
                  setkplus1Coordinates(cgm.conjugatedGradientDirection, ls.alphaOptimum);             
                  fxkplus1 = ls.linearFunctionInAlphaOptimum;                 
                  //logger.debug("x" + iterationNumberkplus1 + " = " + kplus1Coordinates + "    ");
                  //logger.debug("f(x" + iterationNumberkplus1 + ") = " + fxkplus1);
                  
                  gradientkplus1.set(ls.dfOptimum);
                  //logger.debug("gradientkplus1 = " + gradientkplus1);
                  
                  
                  checkConvergence(CGconvergenceCriterion);
                  /*if (convergence == true) {
                        forceField.setEnergyHessian(kplus1Coordinates);
                        NewtonRaphsonMethod nrm = new NewtonRaphsonMethod();
                        //nrm.determinat(gradientkplus1,forceField.getEnergyHessian());
                        nrm.hessianEigenValues(forceField.getForEnergyHessian(),kplus1Coordinates.getSize());
                  }*/
                  
                  //logger.debug(" ");
                  //logger.debug("convergence = " + convergence);
                  
                  if (fxkplus1 <= fxk + 0.0001 * ls.alphaOptimum * linearFunctionDerivativek) {}
                  else {
                        //logger.debug("SUFFICIENT DECREASE CONDITION NOT SATISFIED");
                        break;
                  }
                  if (Math.abs(gradientkplus1.dot(cgm.conjugatedGradientDirection)) <= -0.09 * linearFunctionDerivativek) {}
                  else {
                        //logger.debug("CURVATURE CONDITION NOT SATISFIED");
                        //logger.debug("linearFunctionDerivativekplus1 = " + gradientkplus1.dot(cgm.conjugatedGradientDirection));
                        //logger.debug("linearFunctionDerivativek = " + linearFunctionDerivativek);
                  }
                              
                  //if (iterationNumberkplus1 % 50 == 0 | iterationNumberkplus1 == 1) {
                        //logger.debug("");
                        //logger.debug("f(x" + iterationNumberk + ") = " + fxk);
                        //logger.debug("f(x" + iterationNumberkplus1 + ") = " + fxkplus1);
                        //logger.debug("fxkplus1 - fxk = " + (fxkplus1 - fxk));
                  //}
                  
                  /*if (iterationNumberkplus1 != 1) {
                        logger.debug("gk+1.gk = " + gradientkplus1.dot(gradientk));
                  }*/
                  //forceField.setEnergyHessian(kplus1Coordinates);
                  //nrm.determinat(gradientkplus1, forceField.getEnergyHessian());
                  //nrm.hessianEigenValues(forceField.getForEnergyHessian(), kplus1Coordinates.getSize());

             }
             
             conjugateGradientMinimum = kplus1Coordinates;
             minimumFunctionValueCGM = fxkplus1;
             //logger.debug("conjugateGradientMinimum, forceField.getEnergyGradient().norm() = " + forceField.getEnergyGradient().norm());
             //logger.debug("The CG minimum energy is at: " + conjugateGradientMinimum);
             //logger.debug("f(x" + iterationNumberk + ") = " + fxk);
             //logger.debug("f(minimum) = " + fxkplus1);
             //logger.debug("CG converge at iteration " + iterationNumberkplus1);
             //logger.debug("Energy function evaluation number : " + forceField.functionEvaluationNumber);
             
            if (molecule !=null){
                //logger.debug("CGM: kplus1Coordinates:"+kplus1Coordinates);
                ForceFieldTools.assignCoordinatesToMolecule(kplus1Coordinates, molecule);
            }

            return;
      }


      /**
       *  Gets the conjugatedGradientMinimum attribute of the GeometricMinimizer object
       *
       *@return    The minimumCoordinates value
       */
00534       public GVector getConjugateGradientMinimum() {
            return conjugateGradientMinimum;
      }


      /**
       *  Gets the conjugatedGradientMinimum attribute of the GeometricMinimizer object
       *
       *@return    The minimumCoordinates value
       */
00544       public double getMinimumFunctionValueCGM() {
            return minimumFunctionValueCGM;
      }


      /**
       *  Set convergence parameters for Newton-Raphson Method
       *
       * @param  changeNRMaximumIteration  Maximum number of iteration for Newton-Raphson method
       * @param  changeNRConvergenceCriterion  Convergence criterion for Newton-Raphson method
       */
00555       public void setConvergenceParametersForNRM(int changeNRMaximumIteration, double changeNRConvergenceCriterion){
            NRMaximumIteration = changeNRMaximumIteration;
            NRconvergenceCriterion = changeNRConvergenceCriterion;
            return;
      }


      /**
       *  Minimize the potential energy function using the Newton-Raphson method
       *
       * @param  forceField         The potential function to be used
       */
00567       public void newtonRaphsonMinimization(GVector initialCoordinates, IPotentialFunction forceField) {
            
            initializeMinimizationParameters(initialCoordinates);

            forceField.setEnergyGradient(kCoordinates);
            gradientk = forceField.getEnergyGradient();
            //logger.debug("gradient at iteration 1 : g1 = " + gradientk);

            newtonRaphsonMinimum = new GVector(kCoordinates);

            //logger.debug("");
            //logger.debug("FORCEFIELDTESTS NewtonRaphsonTest");
            
            GMatrix hessian = new GMatrix(initialCoordinates.getSize(),initialCoordinates.getSize());
            
            while ((iterationNumberkplus1 < NRMaximumIteration) & (convergence == false)) {
                  
                  iterationNumberkplus1 += 1;
                  iterationNumberk += 1;
                  //logger.debug("");
                  //logger.debug("NR Iteration number: " + iterationNumberkplus1);

                  if (iterationNumberkplus1 != 1) {
                        kCoordinates.set(kplus1Coordinates);
                        gradientk.set(gradientkplus1);
                  }
                  forceField.setEnergyHessian(kCoordinates);
                  hessian.set(forceField.getEnergyHessian());
                  //logger.debug("hessian = " + hessian);
                  nrm.gradientPerInverseHessian(gradientk,hessian);
                  //logger.debug("GradientPerInverseHessian = " + nrm.getGradientPerInverseHessian());
                  
                  setkplus1Coordinates(nrm.getGradientPerInverseHessian(), -1);
                  //logger.debug("x" + iterationNumberkplus1 + " = " + kplus1Coordinates);
                  //logger.debug("f(x" + iterationNumberkplus1 + ") = " + fxkplus1);

                  forceField.setEnergyGradient(kplus1Coordinates);
                  gradientkplus1.set(forceField.getEnergyGradient());

                  checkConvergence(NRconvergenceCriterion);
                  //logger.debug("convergence: " + convergence);

                  //logger.debug("");
                  //logger.debug("f(x" + iterationNumberk + ") = " + fxk);
                  //logger.debug("f(x" + iterationNumberkplus1 + ") = " + fxkplus1);
                  //logger.debug("gradientkplus1 = " + gradientkplus1);
            }
               newtonRaphsonMinimum.set(kplus1Coordinates);
               //logger.debug("The NR minimum energy is at: " + newtonRaphsonMinimum);
             
            return;
      }


      /**
       *  Gets the newtonRaphsonMinimum attribute of the GeometricMinimizer object
       *
       *@return    The newtonRaphsonMinimum value
       */
00626       public GVector getNewtonRaphsonMinimum() {
            return newtonRaphsonMinimum;
      }

}


Generated by  Doxygen 1.6.0   Back to index