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
 */
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
	 */
	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.
	 */
	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. 
	 */
	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. 
	 */
	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
	 */
	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.
	 */
    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
	 */
    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);