torsions.java

来自「化学图形处理软件」· Java 代码 · 共 489 行 · 第 1/2 页

		}
		*/
		moleculeCurrentCoordinates.set(coords3d);
		
	}


	/**
	 *  Evaluate the MMFF94 torsions term.
	 *
	 *@param  coords3d  Current molecule coordinates.
	 *@return        MMFF94 torsions term value.
	 */
	public double functionMMFF94SumET(GVector coords3d) {
		//logger.debug("SetPhi for torsion energy evaluation");
		setPhi(coords3d);
		mmff94SumET = 0;
		double torsionEnergy=0;
		for (int m = 0; m < torsionNumber; m++) {
			//logger.debug("phi[" + m + "] = " + Math.round(Math.toDegrees(phi[m])) + ",	cos(phi[" + m + "]) = " + Math.round(Math.cos(phi[m])) + ",	cos(2 * phi[" + m + "]) = " + Math.round(Math.cos(2 * phi[m])) + ",	cos(3 * phi[" + m + "]) = " + Math.round(Math.cos(3 * phi[m]))); 
			torsionEnergy = v1[m] * (1 + Math.cos(phi[m])) + v2[m] * (1 - Math.cos(2 * phi[m])) + v3[m] * (1 + Math.cos(3 * phi[m]));
			//logger.debug("phi[" + m + "] = " + Math.toDegrees(phi[m]) + ", cph" + Math.cos(phi[m]) + ", c2ph" + Math.cos(2 * phi[m]) + ", c3ph" + Math.cos(3 * phi[m]) + ", te=" + torsionEnergy);
			//if (torsionEnergy < 0) {
			//	torsionEnergy= (-1) * torsionEnergy;
			//}
			mmff94SumET = mmff94SumET + torsionEnergy;
			
			//mmff94SumET = mmff94SumET + v1[m] * (1 + phi[m]) + v2[m] * (1 - 2 * phi[m]) + v3[m] * (1 + 3 * phi[m]);
		}
		//logger.debug("mmff94SumET = " + mmff94SumET);
		return mmff94SumET;
	}


	/**
	 *  Evaluate the gradient of the torsions term.
	 *  
	 *
	 *@param  coords3d  Current molecule coordinates.
	 */
	public void setGradientMMFF94SumET(GVector coords3d) {

		gradientMMFF94SumET.setSize(coords3d.getSize());
		//logger.debug("Set phi for torsion energy gradient calculation");
		setPhi(coords3d);
		dPhi.setSize(coords3d.getSize());

		double sumGradientET;
		for (int i = 0; i < gradientMMFF94SumET.getSize(); i++) {

			sumGradientET = 0;
			dPhi.setElement(i,1);                 // dPhi : partial derivative of phi. To change in the future

			for (int m = 0; m < torsionNumber; m++) {

				sumGradientET = sumGradientET - v1[m] * Math.sin(phi[m]) * dPhi.getElement(i) + 
					v2[m] * Math.sin(2 * phi[m]) * 2 * dPhi.getElement(i) - 
					v3[m] * Math.sin(3 * phi[m]) * 3 * dPhi.getElement(i);
			}
			gradientMMFF94SumET.setElement(i, sumGradientET);
		}
		//logger.debug("gradientMMFF94SumET = " + gradientMMFF94SumET);
	}


	/**
	 *  Get the gradient of the torsions term. 
	 *  
	 *
	 *@return           torsions gradient value.
	 */
	public GVector getGradientMMFF94SumET() {
		return gradientMMFF94SumET;
	}


	/**
	 *  Evaluate a 2nd order error approximation of the gradient, for the torsion term, 
	 *  given the atoms coordinates.
	 *
	 *@param  coord3d  Current molecule coordinates.
	 */
	public void set2ndOrderErrorApproximateGradientMMFF94SumET(GVector coord3d) {
		//logger.debug("Set the approximative gradient of the torsion energy");
		order2ndErrorApproximateGradientMMFF94SumET.setSize(coord3d.getSize());
		xplusSigma = new GVector(coord3d.getSize());
		xminusSigma = new GVector(coord3d.getSize());
		
		for (int m = 0; m < order2ndErrorApproximateGradientMMFF94SumET.getSize(); m++) { 
			//logger.debug("m = " + m);
			xplusSigma.set(coord3d);
			xplusSigma.setElement(m,coord3d.getElement(m) + sigma);
			
			xminusSigma.set(coord3d);
			xminusSigma.setElement(m,coord3d.getElement(m) - sigma);
			order2ndErrorApproximateGradientMMFF94SumET.setElement(m,(functionMMFF94SumET(xplusSigma) - functionMMFF94SumET(xminusSigma)) / (2 * sigma));
		}
			
		//logger.debug("order2ndErrorApproximateGradientMMFF94SumET : " + order2ndErrorApproximateGradientMMFF94SumET);
	}


	/**
	 *  Get the 2nd order error approximate gradient of the torsion term.
	 *
	 *
	 *@return           torsion approximate gradient value
	 */
	public GVector get2ndOrderErrorApproximateGradientMMFF94SumET() {
		return order2ndErrorApproximateGradientMMFF94SumET;
	}


	/**
	 *  Evaluate an 5 order approximation of the gradient, of the torsion term, 
	 *  given the atoms coordinates
	 *
	 *@param  coord3d  Current molecule coordinates.
	 */
	public void set5thOrderApproximateGradientMMFF94SumET(GVector coord3d) {
		order5thErrorApproximateGradientMMFF94SumET.setSize(coord3d.getSize());
		double sigma = Math.pow(0.000000000000001,0.2);
		GVector xplusSigma = new GVector(coord3d.getSize());
		GVector xminusSigma = new GVector(coord3d.getSize());
		GVector xplus2Sigma = new GVector(coord3d.getSize());
		GVector xminus2Sigma = new GVector(coord3d.getSize());
		
		for (int m=0; m < order5thErrorApproximateGradientMMFF94SumET.getSize(); m++) {
			xplusSigma.set(coord3d);
			xplusSigma.setElement(m,coord3d.getElement(m) + sigma);
			xminusSigma.set(coord3d);
			xminusSigma.setElement(m,coord3d.getElement(m) - sigma);
			xplus2Sigma.set(coord3d);
			xplus2Sigma.setElement(m,coord3d.getElement(m) + 2 * sigma);
			xminus2Sigma.set(coord3d);
			xminus2Sigma.setElement(m,coord3d.getElement(m) - 2 * sigma);
			order5thErrorApproximateGradientMMFF94SumET.setElement(m, (8 * (functionMMFF94SumET(xplusSigma) - functionMMFF94SumET(xminusSigma)) - (functionMMFF94SumET(xplus2Sigma) - functionMMFF94SumET(xminus2Sigma))) / (12 * sigma));
		}
			
		//logger.debug("order5thErrorApproximateGradientMMFF94SumET : " + order5thErrorApproximateGradientMMFF94SumET);
	}


	/**
	 *  Get the 5th order error approximate gradient of the torsion term.
	 *
	 *@return        Torsion 5th order error approximate gradient value.
	 */
	public GVector get5thOrderErrorApproximateGradientMMFF94SumET() {
		return order5thErrorApproximateGradientMMFF94SumET;
	}


	/**
	 *  Evaluate the hessian of the torsions.
	 *
	 *@param  coords3d  Current molecule coordinates.
	 */
	public void setHessianMMFF94SumET(GVector coords3d) {

		double[] forHessian = new double[coords3d.getSize() * coords3d.getSize()];
		setPhi(coords3d);
		double[] ddPhi = new double[coords3d.getSize() * coords3d.getSize()];
		
		double sumHessianET = 0;
		for (int i = 0; i < coords3d.getSize(); i++) {
			for (int j = 0; j < dPhi.getSize(); j++) {
				ddPhi[i*j] = 0;
				for (int m = 0; m < torsionNumber; m++) {
					sumHessianET = sumHessianET - v1[m] * (Math.cos(phi[m]) * dPhi.getElement(i) * dPhi.getElement(j) + Math.sin(phi[m]) * ddPhi[i*j]) +
					2 * v2[m] * (Math.cos(2 * phi[m]) * 2 * dPhi.getElement(i) * dPhi.getElement(j) + Math.sin(2 * phi[m]) * ddPhi[i*j]) -
					3 * v3[m] * (Math.cos(3 * phi[m]) * 3 * dPhi.getElement(i) * dPhi.getElement(j) + Math.sin(3 * phi[m]) * ddPhi[i*j]);
				}
			}
			forHessian[i] = 0.5 * sumHessianET;
		}

		hessianMMFF94SumET.setSize(coords3d.getSize(), coords3d.getSize());
		hessianMMFF94SumET.set(forHessian); 
		//logger.debug("hessianMMFF94SumET : " + hessianMMFF94SumET);
	}


	/**
	 *  Get the hessian of the torsions.
	 *
	 *@return        Hessian value of the torsions term.
	 */
	public GMatrix getHessianMMFF94SumET() {
		return hessianMMFF94SumET;
	}


	/**
	 *  Evaluate a 2nd order approximation of the Hessian, for the torsion energy term,
	 *  given the atoms coordinates.
	 *
	 *@param  coord3d  Current molecule coordinates.
	 */
	public void set2ndOrderErrorApproximateHessianMMFF94SumET(GVector coord3d) {
		forOrder2ndErrorApproximateHessian = new double[coord3d.getSize() * coord3d.getSize()];
		
		double sigma = Math.pow(0.000000000000001,0.33);
		GVector xminusSigma = new GVector(coord3d.getSize());
		GVector xplusSigma = new GVector(coord3d.getSize());
		GVector gradientAtXminusSigma = new GVector(coord3d.getSize());
		GVector gradientAtXplusSigma = new GVector(coord3d.getSize());
		
		int forHessianIndex;
		for (int i = 0; i < coord3d.getSize(); i++) {
			xminusSigma.set(coord3d);
			xminusSigma.setElement(i,coord3d.getElement(i) - sigma);
			setGradientMMFF94SumET(xminusSigma);
			gradientAtXminusSigma.set(gradientMMFF94SumET);
			xplusSigma.set(coord3d);
			xplusSigma.setElement(i,coord3d.getElement(i) + sigma);
			setGradientMMFF94SumET(xplusSigma);
			gradientAtXplusSigma.set(gradientMMFF94SumET);
			for (int j = 0; j < coord3d.getSize(); j++) {
				forHessianIndex = i*coord3d.getSize()+j;
				forOrder2ndErrorApproximateHessian[forHessianIndex] = (gradientAtXplusSigma.getElement(j) - gradientAtXminusSigma.getElement(j)) / (2 * sigma);
				//(functionMMFF94SumET(xplusSigma) - 2 * fx + functionMMFF94SumET(xminusSigma)) / Math.pow(sigma,2);
			}
		}
		
		order2ndErrorApproximateHessianMMFF94SumET = new GMatrix(coord3d.getSize(), coord3d.getSize());
		order2ndErrorApproximateHessianMMFF94SumET.set(forOrder2ndErrorApproximateHessian);
		//logger.debug("order2ndErrorApproximateHessianMMFF94SumET : " + order2ndErrorApproximateHessianMMFF94SumET);
	}


	/**
	 *  Get the 2nd order error approximate Hessian for the torsion term.
	 *
	 *
	 *@return           Torsion 2nd order error approximate Hessian value.
	 */
	public GMatrix get2ndOrderErrorApproximateHessianMMFF94SumET() {
		return order2ndErrorApproximateHessianMMFF94SumET;
	}

}