geometricminimizer.java

化学图形处理软件

			
			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
	 */
	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
	 */
	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
	 */
	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
	 */
	public GVector getConjugateGradientMinimum() {
		return conjugateGradientMinimum;
	}


	/**
	 *  Gets the conjugatedGradientMinimum attribute of the GeometricMinimizer object
	 *
	 *@return    The minimumCoordinates value
	 */
	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
	 */
	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
	 */
	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
	 */
	public GVector getNewtonRaphsonMinimum() {
		return newtonRaphsonMinimum;
	}

}