anglebending.java

化学图形处理软件

	}


	/**
	 *  Set a 2nd order approximation of the gradient, for the angle bending, given the atoms
	 *  coordinates
	 *
	 *@param  coord3d  Current molecule coordinates.
	 */
	public void setAngleBending2ndOrderErrorApproximateGradient(GVector coord3d) {
		angleBendingOrder2ndErrorApproximateGradient = new double[coord3d.getSize()][];
		double sigma = Math.pow(0.000000000000001,0.33);
		GVector xplusSigma = new GVector(coord3d.getSize());
		GVector xminusSigma = new GVector(coord3d.getSize());
		/*boolean[] sigmaChange = new boolean[coord3d.getSize()];
		for (int i=1; i < sigmaChange.length; i++) {sigmaChange[i] = false;}
		*/
		double[] deltavInXplusSigma = null;
		double[] deltavInXminusSigma = null;
		for (int m = 0; m < angleBendingOrder2ndErrorApproximateGradient.length; m++) {
			angleBendingOrder2ndErrorApproximateGradient[m] = new double[angleNumber];
			xplusSigma.set(coord3d);
			xplusSigma.setElement(m,coord3d.getElement(m) + sigma);
			setDeltav(xplusSigma);
			deltavInXplusSigma = getDeltav();
			xminusSigma.set(coord3d);
			xminusSigma.setElement(m,coord3d.getElement(m) - sigma);
			setDeltav(xminusSigma);
			deltavInXminusSigma = getDeltav();
			for (int l=0; l < angleNumber; l++) {
				angleBendingOrder2ndErrorApproximateGradient[m][l] = (deltavInXplusSigma[l] - deltavInXminusSigma[l]) / (2 * sigma);
			}
		}
			
		//logger.debug("order2ndErrorApproximateGradientMMFF94SumEA : " + order2ndErrorApproximateGradientMMFF94SumEA);
	}


	/**
	 *  Get the angle bending 2nd order error approximate gradient.
	 *
	 *
	 *@return           Angle bending 2nd order error approximate gradient value.
	 */
	public double[][] getAngleBending2ndOrderErrorApproximateGradient() {
		return angleBendingOrder2ndErrorApproximateGradient;
	}


	/**
	 *  Evaluate the gradient of the angle bending term for a given atoms
	 *  coordinates
	 *
	 *@param  coord3d  Current molecule coordinates.
	 */
	public void setGradientMMFF94SumEA(GVector coord3d) {
		gradientMMFF94SumEA = new GVector(coord3d.getSize());

/*		boolean[] sigmaChange = new boolean[coord3d.getSize()];
		for (int i=1; i < sigmaChange.length; i++) {sigmaChange[i] = false;}
*/
		setDeltav(coord3d);
		setAngleBendingFirstDerivative(coord3d);
		
		double sumGradientEA;
		for (int m = 0; m < gradientMMFF94SumEA.getSize(); m++) {

			sumGradientEA = 0;
			for (int l = 0; l < angleNumber; l++) {

				sumGradientEA = sumGradientEA + (k2[l] * 2 * deltav[l] + k3[l] * 3 * Math.pow(deltav[l],2)) * dDeltav[m][l];
			}
			
			gradientMMFF94SumEA.setElement(m, sumGradientEA);
		}

		//logger.debug("gradientMMFF94SumEA : " + gradientMMFF94SumEA);
	}


	/**
	 *  Get the gradient of the angle bending term.
	 *
	 *
	 *@return           Angle bending gradient value
	 */
	public GVector getGradientMMFF94SumEA() {
		return gradientMMFF94SumEA;
	}


	/**
	 *  Evaluate a 2nd order approximation of the gradient, of the angle bending term, for a given atoms
	 *  coordinates
	 *
	 *@param  coord3d  Current molecule coordinates.
	 */
	public void set2ndOrderErrorApproximateGradientMMFF94SumEA(GVector coord3d) {
		order2ndErrorApproximateGradientMMFF94SumEA = new GVector(coord3d.getSize());
		double sigma = Math.pow(0.000000000000001,0.33);
		GVector xplusSigma = new GVector(coord3d.getSize());
		GVector xminusSigma = new GVector(coord3d.getSize());
		boolean[] sigmaChange = new boolean[coord3d.getSize()];
		for (int i=1; i < sigmaChange.length; i++) {sigmaChange[i] = false;}
		
		for (int m = 0; m < order2ndErrorApproximateGradientMMFF94SumEA.getSize(); m++) {
			xplusSigma.set(coord3d);
			xplusSigma.setElement(m,coord3d.getElement(m) + sigma);
			xminusSigma.set(coord3d);
			xminusSigma.setElement(m,coord3d.getElement(m) - sigma);
			order2ndErrorApproximateGradientMMFF94SumEA.setElement(m,(functionMMFF94SumEA(xplusSigma) - functionMMFF94SumEA(xminusSigma)) / (2 * sigma));
		}
			
		//logger.debug("order2ndErrorApproximateGradientMMFF94SumEA : " + order2ndErrorApproximateGradientMMFF94SumEA);
	}


	/**
	 *  Get the 2nd order error approximate gradient of the angle bending term.
	 *
	 *
	 *@return           Angle bending 2nd order error approximate gradient value.
	 */
	public GVector get2ndOrderErrorApproximateGradientMMFF94SumEA() {
		return order2ndErrorApproximateGradientMMFF94SumEA;
	}


	/**
	 *  Evaluate a 5th order error approximation of the gradient, of the angle bending term, for a given atoms
	 *  coordinates
	 *
	 *@param  coord3d  Current molecule coordinates.
	 */
	public void set5thOrderErrorApproximateGradientMMFF94SumEA(GVector coord3d) {
		order5thErrorApproximateGradientMMFF94SumEA = 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());
		boolean[] sigmaChange = new boolean[coord3d.getSize()];
		for (int i=1; i < sigmaChange.length; i++) {sigmaChange[i] = false;}
		
		for (int m=0; m < order5thErrorApproximateGradientMMFF94SumEA.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);
			order5thErrorApproximateGradientMMFF94SumEA.setElement(m, (8 * (functionMMFF94SumEA(xplusSigma) - functionMMFF94SumEA(xminusSigma)) 
					- (functionMMFF94SumEA(xplus2Sigma) - functionMMFF94SumEA(xminus2Sigma))) / (12 * sigma));
		}
			
		//logger.debug("order5thErrorApproximateGradientMMFF94SumEA : " + order5thErrorApproximateGradientMMFF94SumEA);
	}


	/**
	 *  Get the 5 order approximate gradient of the angle bending term.
	 *
	 *@return        Angle bending 5 order approximate gradient value.
	 */
	public GVector get5thOrderErrorApproximateGradientMMFF94SumEA() {
		return order5thErrorApproximateGradientMMFF94SumEA;
	}


	/**
	 *  Calculate the angle bending second derivative respect to the cartesian coordinates of the atoms.
	 *
	 *@param  coord3d  Current molecule coordinates.
	 */
	public void setAngleBendingSecondDerivative(GVector coord3d) {
		
		ddDeltav = new double[coord3d.getSize()][][];
		
		Double forAtomNumber = null;
		int atomNumbern;
		int atomNumberm;
		int coordinaten;
		int coordinatem;
		
		double ddDeltav1;
		double ddDeltav2;
		
		double ddDeltav2a;
		double ddDeltav2b;
		
		double ddDeltav2a1;
		double ddDeltav2a2;
		double ddDeltav2a3;
		double ddDeltav2a4;

		double ddDeltav2a1a;
		double ddDeltav2a1b;
		
		double ddDeltav2a4a1;
		double ddDeltav2a4a1a;
		double ddDeltav2a4a1b;
		double ddDeltav2a4a1c;
		double ddDeltav2a4a2;

		double ddDeltav2a4b;

		double ddDeltav2a4c1;
		double ddDeltav2a4c1a;
		double ddDeltav2a4c1b;
		double ddDeltav2a4c1c;
		double ddDeltav2a4c2;

		double ddDeltav2a4d;

		setAngleBendingFirstDerivative(coord3d);
		
		for (int n=0; n<coord3d.getSize(); n++) {
			ddDeltav[n] = new double[coord3d.getSize()][];
			
			forAtomNumber = new Double(n/3);
			coordinaten = n % 3;
			//logger.debug("coordinaten = " + coordinaten);
				
			atomNumbern = forAtomNumber.intValue();
			//logger.debug("atomNumbern = " + atomNumbern);
				
			for (int m = 0; m < coord3d.getSize(); m++) {
			
				ddDeltav[n][m] = new double[angleNumber];
			
				forAtomNumber = new Double(m/3);
				coordinatem = m % 3;
				//logger.debug("coordinatem = " + coordinatem);

				atomNumberm = forAtomNumber.intValue();
				//logger.debug("atomNumberm = " + atomNumberm);

				for (int l = 0; l < angleNumber; l++) {
				
					if ((angleAtomPosition[l][0] == atomNumberm) | (angleAtomPosition[l][1] == atomNumberm) | (angleAtomPosition[l][2] == atomNumberm)) {
						if ((angleAtomPosition[l][0] == atomNumbern) | (angleAtomPosition[l][1] == atomNumbern) | (angleAtomPosition[l][2] == atomNumbern)) {

							Point3d xi = new Point3d(coord3d.getElement(3 * angleAtomPosition[l][0]), coord3d.getElement(3 * angleAtomPosition[l][0] + 1),coord3d.getElement( 3 * angleAtomPosition[l][0] + 2));
							//logger.debug("xi = " + xi);
							Point3d xj = new Point3d(coord3d.getElement(3 * angleAtomPosition[l][1]), coord3d.getElement(3 * angleAtomPosition[l][1] + 1),coord3d.getElement( 3 * angleAtomPosition[l][1] + 2));
							//logger.debug("xj = " + xj);
							Point3d xk = new Point3d(coord3d.getElement(3 * angleAtomPosition[l][2]), coord3d.getElement(3 * angleAtomPosition[l][2] + 1),coord3d.getElement( 3 * angleAtomPosition[l][2] + 2));
							//logger.debug("xk = " + xk);
				
							Vector3d xij = new Vector3d();
							xij.sub(xi,xj);
							//logger.debug("xij = " + xij);
							Vector3d xkj = new Vector3d();
							xkj.sub(xk,xj);
							//logger.debug("xkj = " + xkj);
					
							double rij = xi.distance(xj);
							//logger.debug("rij = " + rij);
							double rkj = xk.distance(xj);
							//logger.debug("rkj = " + rkj);

							ddDeltav1 = (-1/Math.sqrt(Math.pow(1-Math.pow(xij.dot(xkj)/(rij * rkj),2),3)))
									* (xij.dot(xkj)/(rij * rkj))
									* (1/(Math.pow(rij,2) * Math.pow(rkj,2)));
							
							ddDeltav2 = (-1/Math.sqrt(1-Math.pow(xij.dot(xkj)/(rij * rkj),2))) 
									* (1/(Math.pow(rij,4) * Math.pow(rkj,4)));
							
							ddDeltav2a = Math.pow(rij,2) * Math.pow(rkj,2);
							
							ddDeltav2a1 = rij * rkj;
							ddDeltav2a1a = 0;
							ddDeltav2a1b = 0;
							ddDeltav2a2 = 0;
							ddDeltav2a3 = 0;
							ddDeltav2a4 = xij.dot(xkj);
							ddDeltav2a4a1a = 0;
							ddDeltav2a4a1b = 0;
							ddDeltav2a4a1c = 0;
							ddDeltav2a4a2 = 0;
							ddDeltav2b = 0;

							ddDeltav2a4a1 = 0;

							ddDeltav2a4b = 0;

							ddDeltav2a4c1 = 0;
							ddDeltav2a4c1a = 0;
							ddDeltav2a4c1b = 0;
							ddDeltav2a4c1c = 0;
							ddDeltav2a4c2 = 0;

							ddDeltav2a4d = 0;

							//logger.debug("OK: had d1 and have the atomNumbern");
						
							if (angleAtomPosition[l][0] == atomNumberm) {

								switch (coordinatem) {
									case 0: ddDeltav1 = ddDeltav1 * ((xk.x-xj.x) * rij * rkj - (xij.dot(xkj)) * rkj * ((xi.x-xj.x)/rij));
										ddDeltav2b = (xk.x-xj.x) * rij * rkj - (xij.dot(xkj)) * rkj * ((xi.x-xj.x)/rij);
										ddDeltav2a1a = 1;
										ddDeltav2a1b = 0;
										ddDeltav2a2 = xk.x-xj.x;
										ddDeltav2a3 = rkj * ((xi.x-xj.x)/rij);
										ddDeltav2a4a1a = 1;
										ddDeltav2a4a2 = xi.x-xj.x;
										ddDeltav2a4c1a = 0;
										ddDeltav2a4c2 = 0;
										ddDeltav2a4b = (xi.x-xj.x)/rij;
										ddDeltav2a4d = 0;
										break;
									case 1:	ddDeltav1 = ddDeltav1 * ((xk.y-xj.y) * rij * rkj - (xij.dot(xkj)) * rkj * ((xi.y-xj.y)/rij));
										ddDeltav2b = (xk.y-xj.y) * rij * rkj - (xij.dot(xkj)) * rkj * ((xi.y-xj.y)/rij);
										ddDeltav2a1a = 1;
										ddDeltav2a1b = 0;
										ddDeltav2a2 = xk.y-xj.y;
										ddDeltav2a3 = rkj * ((xi.y-xj.y)/rij);
										ddDeltav2a4a1b = 1;
										ddDeltav2a4a2 = xi.y-xj.y;
										ddDeltav2a4c1b = 0;
										ddDeltav2a4c2 = 0;
										ddDeltav2a4b = (xi.y-xj.y)/rij;
										ddDeltav2a4d = 0;
										break;
									case 2: ddDeltav1 = ddDeltav1 * ((xk.z-xj.z) * rij * rkj - (xij.dot(xkj)) * rkj * ((xi.z-xj.z)/rij));
										ddDeltav2b = (xk.z-xj.z) * rij * rkj - (xij.dot(xkj)) * rkj * ((xi.z-xj.z)/rij);
										ddDeltav2a1a = 1;
										ddDeltav2a1b = 0;
										ddDeltav2a2 = xk.z-xj.z;
										ddDeltav2a3 = rkj * ((xi.z-xj.z)/rij);
										ddDeltav2a4a1c = 1;
										ddDeltav2a4a2 = xi.z-xj.z;
										ddDeltav2a4c1c = 0;
										ddDeltav2a4c2 = 0;
										ddDeltav2a4b = (xi.z-xj.z)/rij;
										ddDeltav2a4d = 0;
										break;
								}
							}
							if (angleAtomPosition[l][1] == atomNumberm) {

								switch (coordinatem) {
									case 0: ddDeltav1 = ddDeltav1 * ((2 * xj.x - xk.x - xi.x) * rij * rkj - (xij.dot(xkj)) * ((-(xi.x-xj.x)/rij) * rkj + (-(xk.x-xj.x)/rkj) * rij));
										ddDeltav2b = (2 * xj.x - xk.x - xi.x) * rij * rkj - (xij.dot(xkj)) * ((-(xi.x-xj.x)/rij) * rkj + (-(xk.x-xj.x)/rkj) * rij);
										ddDeltav2a1a = -1;
										ddDeltav2a1b = -1;
										ddDeltav2a2 = 2 * xj.x - xk.x - xi.x;
										ddDeltav2a3 = (-(xi.x-xj.x)/rij) * rkj + (-(xk.x-xj.x)/rkj) * rij;
										ddDeltav2a4a1a = -1;
										ddDeltav2a4a2 = -(xi.x-xj.x);
										ddDeltav2a4c1a = -1;
										ddDeltav2a4c2 = -(xk.x-xj.x);
										ddDeltav2a4b = -(xi.x-xj.x)/rij;
										ddDeltav2a4d = -(xk.x-xj.x)/rkj;
										break;
									case 1:	ddDeltav1 = ddDeltav1 * ((2 * xj.y - xk.y - xi.y) * rij * rkj - (xij.dot(xkj)) * ((-(xi.y-xj.y)/rij) * rkj + (-(xk.y-xj.y)/rkj) * rij));
										ddDeltav2b = (2 * xj.y - xk.y - xi.y) * rij * rkj - (xij.dot(xkj)) * ((-(xi.y-xj.y)/rij) * rkj + (-(xk.y-xj.y)/rkj) * rij);
										ddDeltav2a1a = -1;
										ddDeltav2a1b = -1;
										ddDeltav2a2 = 2 * xj.y - xk.y - xi.y;
										ddDeltav2a3 = (-(xi.y-xj.y)/rij) * rkj + (-(xk.y-xj.y)/rkj) * rij;
										ddDeltav2a4a1b = -1;
										ddDeltav2a4a2 = -(xi.y-xj.y);
										ddDeltav2a4c1b = -1;
										ddDeltav2a4c2 = -(xk.y-xj.y);
										ddDeltav2a4b = -(xi.y-xj.y)/rij;
										ddDeltav2a4d = -(xk.y-xj.y)/rkj;
										break;
									case 2: ddDeltav1 = ddDeltav1 * ((2 * xj.z - xk.z - xi.z) * rij * rkj - (xij.dot(xkj)) * ((-(xi.z-xj.z)/rij) * rkj + (-(xk.z-xj.z)/rkj) * rij));
										ddDeltav2b = (2 * xj.z - xk.z - xi.z) * rij * rkj - (xij.dot(xkj)) * ((-(xi.z-xj.z)/rij) * rkj + (-(xk.z-xj.z)/rkj) * rij);
										ddDeltav2a1a = -1;