stretchbendinteractions.java

化学图形处理软件

			}
			//else {System.out.println("deltarkj[" + i + "] was no recalculated");}
		}
		/*if 	(changedCoordinates == changeAtomCoordinates.length) {
			for (int m = 0; m < ab.angleNumber; m++) {
				System.out.println("phi[" + m + "] = " + Math.toDegrees(phi[m]));
			}
		}
		*/
		moleculeCurrentCoordinates.set(coords3d);
	}


	
	/**
	 *  Set the MMFF94 stretch-bend interaction term given the atoms cartesian
	 *  coordinates.
	 *
	 *@param  coords3d  Current molecule coordinates.
	 */
	public void setFunctionMMFF94SumEBA(GVector coords3d) {
		//ab.setAngleBendingFlag(false);
		if (currentCoordinates.equals(coords3d)) {
		} else {
			setDeltarijAndDeltarkj(coords3d);
			ab.setDeltav(coords3d);
			mmff94SumEBA = 0;
			for (int j = 0; j < ab.angleNumber; j++) {
				//logger.debug("kbaIJK[" + j + "] = " + kbaIJK[j]);
				//logger.debug("kbaKJI[" + j + "] = " + kbaKJI[j]);
				//logger.debug("deltarij[" + j + "] = " + deltarij[j]);
				//logger.debug("deltarkj[" + j + "] = " + deltarkj[j]);
				//logger.debug("ab.deltav[" + j + "] = " + ab.deltav[j]);
				mmff94SumEBA = mmff94SumEBA + 2.51210 * (kbaIJK[j] * deltarij[j] + kbaKJI[j] * deltarkj[j]) * ab.deltav[j];
				//logger.debug("mmff94SumEBA = " + mmff94SumEBA);
			}
			//mmff94SumEBA = Math.abs(mmff94SumEBA);
			//logger.debug("mmff94SumEBA = " + mmff94SumEBA);
			currentCoordinates.set(coords3d);
		}
	}


	/**
	 *  Get the MMFF94 stretch-bend interaction term.
	 *
	 *@return    MMFF94 stretch-bend interaction term value.
	 */
	public double getFunctionMMFF94SumEBA() {
		return mmff94SumEBA;
	}


	/**
	 *  Evaluate the gradient of the stretch-bend interaction term.
	 *
	 *@param  coords3d  Current molecule coordinates.
	 */
	public void setGradientMMFF94SumEBA(GVector coords3d) {

		if (currentCoordinates.equals(coords3d)) {
		} else {
			setFunctionMMFF94SumEBA(coords3d);
		}

		bs.setBondLengthsFirstDerivative(coords3d, deltarij, bondijAtomPosition);
		dDeltarij = bs.getBondLengthsFirstDerivative();
		bs.setBondLengthsFirstDerivative(coords3d, deltarkj, bondkjAtomPosition);
		dDeltarkj = bs.getBondLengthsFirstDerivative();
		ab.setAngleBending2ndOrderErrorApproximateGradient(coords3d);
		dDeltav = ab.getAngleBending2ndOrderErrorApproximateGradient();

		if (dDeltav == null) {logger.debug("setGradient: dDeltav null");} 
		double sumGradientEBA;
		for (int i = 0; i < gradientMMFF94SumEBA.getSize(); i++) {
			sumGradientEBA = 0;
			for (int j = 0; j < ab.angleNumber; j++) {
				sumGradientEBA = sumGradientEBA + (kbaIJK[j] * dDeltarij[i][j] + kbaKJI[j] * dDeltarkj[i][j]) * ab.deltav[j]
						 + (kbaIJK[j] * deltarij[j] + kbaKJI[j] * deltarkj[j]) * ab.angleBendingOrder2ndErrorApproximateGradient[i][j];
			}
			sumGradientEBA = sumGradientEBA * 2.51210;

			gradientMMFF94SumEBA.setElement(i, sumGradientEBA);
			gradientCurrentCoordinates.set(coords3d);
		}
		//logger.debug("gradientMMFF94SumEBA = " + gradientMMFF94SumEBA);
	}


	/**
	 *  Get the gradient of the stretch-bend interaction term.
	 *
	 *@return    stretch-bend interaction gradient value.
	 */
	public GVector getGradientMMFF94SumEBA() {
		return gradientMMFF94SumEBA;
	}


	/**
	 *  Evaluate a 2nd order approximation of the gradient for the stretch-bend interaction term, 
	 *  given the atoms coordinates.
	 *
	 *@param  coord3d  Current molecule coordinates.
	 */
	public void set2ndOrderErrorApproximateGradientMMFF94SumEBA(GVector coord3d) {
		order2ndErrorApproximateGradientMMFF94SumEBA = new GVector(coord3d.getSize());
		double sigma = Math.pow(0.000000000000001,0.33);
		GVector xplusSigma = new GVector(coord3d.getSize());
		GVector xminusSigma = new GVector(coord3d.getSize());
		double fInXplusSigma = 0;
		double fInXminusSigma = 0;

		for (int m = 0; m < order2ndErrorApproximateGradientMMFF94SumEBA.getSize(); m++) {
			xplusSigma.set(coord3d);
			xplusSigma.setElement(m,coord3d.getElement(m) + sigma);
			this.setFunctionMMFF94SumEBA(xplusSigma);
			fInXplusSigma = this.getFunctionMMFF94SumEBA();
			xminusSigma.set(coord3d);
			xminusSigma.setElement(m,coord3d.getElement(m) - sigma);
			this.setFunctionMMFF94SumEBA(xminusSigma);
			fInXminusSigma = this.getFunctionMMFF94SumEBA();
			order2ndErrorApproximateGradientMMFF94SumEBA.setElement(m,(fInXplusSigma - fInXminusSigma) / (2 * sigma));
		}
			
		//logger.debug("order2ndErrorApproximateGradientMMFF94SumEBA : " + order2ndErrorApproximateGradientMMFF94SumEBA);
	}


	/**
	 *  Get the 2nd order error approximate gradient for the stretch-bend term.
	 *
	 *
	 *@return           Stretch-bend interaction 2nd order error approximate gradient value.
	 */
	public GVector get2ndOrderErrorApproximateGradientMMFF94SumEBA() {
		return order2ndErrorApproximateGradientMMFF94SumEBA;
	}


	/**
	 *  Evaluate a 5th order error approximation of the gradient, of the stretch-bend interaction term, for a given atoms
	 *  coordinates
	 *
	 *@param  coord3d  Current molecule coordinates.
	 */
	public void set5thOrderErrorApproximateGradientMMFF94SumEBA(GVector coord3d) {
		order5thErrorApproximateGradientMMFF94SumEBA = new GVector(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());
		double fInXplusSigma = 0;
		double fInXminusSigma = 0;
		double fInXplus2Sigma = 0;
		double fInXminus2Sigma = 0;
		
		for (int m=0; m < order5thErrorApproximateGradientMMFF94SumEBA.getSize(); m++) {
			xplusSigma.set(coord3d);
			xplusSigma.setElement(m,coord3d.getElement(m) + sigma);
			this.setFunctionMMFF94SumEBA(xplusSigma);
			fInXplusSigma = this.getFunctionMMFF94SumEBA();
			xminusSigma.set(coord3d);
			xminusSigma.setElement(m,coord3d.getElement(m) - sigma);
			this.setFunctionMMFF94SumEBA(xminusSigma);
			fInXminusSigma = this.getFunctionMMFF94SumEBA();
			xplus2Sigma.set(coord3d);
			xplus2Sigma.setElement(m,coord3d.getElement(m) + 2 * sigma);
			this.setFunctionMMFF94SumEBA(xplus2Sigma);
			fInXplus2Sigma = this.getFunctionMMFF94SumEBA();
			xminus2Sigma.set(coord3d);
			xminus2Sigma.setElement(m,coord3d.getElement(m) - 2 * sigma);
			this.setFunctionMMFF94SumEBA(xminus2Sigma);
			fInXminus2Sigma = this.getFunctionMMFF94SumEBA();
			order5thErrorApproximateGradientMMFF94SumEBA.setElement(m, (8 * (fInXplusSigma - fInXminusSigma) - (fInXplus2Sigma - fInXminus2Sigma)) / (12 * sigma));
		}
			
		//logger.debug("order5thErrorApproximateGradientMMFF94SumEBA : " + order5thErrorApproximateGradientMMFF94SumEBA);
	}


	/**
	 *  Get the 5 order approximate gradient of the stretch-bend interaction term.
	 *
	 *@return        stretch-bend interaction 5 order approximate gradient value.
	 */
	public GVector get5thOrderErrorApproximateGradientMMFF94SumEBA() {
		return order5thErrorApproximateGradientMMFF94SumEBA;
	}


	/**
	 *  Evaluate the hessian of the stretch-bend interaction.
	 *
	 *@param  coords3d  Current molecule coordinates.
	 */
	public void setHessianMMFF94SumEBA(GVector coords3d) {

		double[] forHessian = new double[coords3d.getSize() * coords3d.getSize()];

		if (currentCoordinates.equals(coords3d)) {
		} else {
			setFunctionMMFF94SumEBA(coords3d);
		}

		if (dDeltarij == null) {logger.debug("dDeltarij null");} 
		if (gradientCurrentCoordinates.equals(coords3d) == false) {
			bs.setBondLengthsFirstDerivative(coords3d, deltarij, bondijAtomPosition);
			dDeltarij = bs.getBondLengthsFirstDerivative();
		}
		if (dDeltarkj == null) {logger.debug("dDeltarkj null");} 
		if (gradientCurrentCoordinates.equals(coords3d) == false) {
			bs.setBondLengthsFirstDerivative(coords3d, deltarkj, bondkjAtomPosition);
			dDeltarkj = bs.getBondLengthsFirstDerivative();
		}
		if (dDeltav == null) {logger.debug("setHessian: dDeltav null");} 
		if (gradientCurrentCoordinates.equals(coords3d) == false) {
			//logger.debug("ab.setAngleBending2ndOrderErrorApproximateGradient()");
			ab.setAngleBending2ndOrderErrorApproximateGradient(coords3d);
			dDeltav = ab.getAngleBending2ndOrderErrorApproximateGradient();
		}
		
		ab.setAngleBending2ndOrderErrorApproximateHessian(coords3d);
		double[][][] ddDeltav = ab.getAngleBending2ndOrderErrorApproximateHessian();
		bs.setBondLengthsSecondDerivative(coords3d, deltarij, bondijAtomPosition);
		double[][][] ddDeltarij = bs.getBondLengthsSecondDerivative();
		bs.setBondLengthsSecondDerivative(coords3d, deltarkj, bondkjAtomPosition);
		double[][][] ddDeltarkj = bs.getBondLengthsSecondDerivative();
		
		if (dDeltav == null) {logger.debug("setHessian: dDeltav null");} 
		//logger.debug("ab.angleNumber = " + ab.angleNumber);
		double sumHessianEBA;
		int forHessianIndex;
		for (int i = 0; i < coords3d.getSize(); i++) {
			for (int j = 0; j < coords3d.getSize(); j++) {
				forHessianIndex = i*coords3d.getSize()+j;
				sumHessianEBA = 0;
				for (int k = 0; k < ab.angleNumber; k++) {
					sumHessianEBA = sumHessianEBA + (kbaIJK[k] * ddDeltarij[i][j][k] + kbaKJI[k] * ddDeltarkj[i][j][k]) * ab.deltav[k]
						 + (kbaIJK[k] * dDeltarij[j][k] + kbaKJI[k] * dDeltarkj[j][k]) * dDeltav[j][k]
						 + (kbaIJK[k] * dDeltarij[j][k] + kbaKJI[k] * dDeltarkj[j][k]) * dDeltav[j][k]
						 + (kbaIJK[k] * deltarij[k] + kbaKJI[k] * deltarkj[k]) * ddDeltav[i][j][k];
				}
				forHessian[forHessianIndex] = sumHessianEBA;
			}
		}

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


	/**
	 *  Get the hessian of the stretch-bend interaction.
	 *
	 *@return    Hessian value of the stretch-bend interaction term.
	 */
	public GMatrix getHessianMMFF94SumEBA() {
		return hessianMMFF94SumEBA;
	}

}