anglebending.java

化学图形处理软件

										ddDeltav2a1b = -1;
										ddDeltav2a2 = 2 * xj.z - xk.z - xi.z;
										ddDeltav2a3 = (-(xi.z-xj.z)/rij) * rkj + (-(xk.z-xj.z)/rkj) * rij;
										ddDeltav2a4a1c = -1;
										ddDeltav2a4a2 = -(xi.z-xj.z);
										ddDeltav2a4c1c = -1;
										ddDeltav2a4c2 = -(xk.z-xj.z);
										ddDeltav2a4b = -(xi.z-xj.z)/rij;
										ddDeltav2a4d = -(xk.z-xj.z)/rkj;
										break;
								}
							}
							if (angleAtomPosition[l][2] == atomNumberm) {

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

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

								switch (coordinaten) {
									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 = ddDeltav2b * (2 * rij * (-(xi.x-xj.x)/rij) * Math.pow(rkj,2) + Math.pow(rij,2) * 2 * rkj * (-(xk.x-xj.x)/rkj));
										ddDeltav2a1a = ddDeltav2a1a * -1;
										ddDeltav2a1b = ddDeltav2a1b * -1;
										ddDeltav2a2 = ddDeltav2a2 * ((-(xi.x-xj.x)/rij) * rkj + (-(xk.x-xj.x)/rkj) * rij);
										ddDeltav2a3 = ddDeltav2a3 * (2 * xj.x - xk.x - xi.x);
										ddDeltav2a4a1a = ddDeltav2a4a1a * (-1);
										ddDeltav2a4a2 = ddDeltav2a4a2 * (-(xi.x-xj.x)/rij);
										ddDeltav2a4c1a = ddDeltav2a4c1a * (-1);
										ddDeltav2a4c2 = ddDeltav2a4c2 * (-(xk.x-xj.x)/rkj);
										ddDeltav2a4b = ddDeltav2a4b * (-(xk.x-xj.x)/rkj);
										ddDeltav2a4d = ddDeltav2a4d * (-(xi.x-xj.x)/rij);
										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 = ddDeltav2b * (2 * rij * (-(xi.y-xj.y)/rij) * Math.pow(rkj,2) + Math.pow(rij,2) * 2 * rkj * (-(xk.y-xj.y)/rkj));
										ddDeltav2a1a = ddDeltav2a1a * -1;
										ddDeltav2a1b = ddDeltav2a1b * -1;
										ddDeltav2a2 = ddDeltav2a2 * ((-(xi.y-xj.y)/rij) * rkj + (-(xk.y-xj.y)/rkj) * rij);
										ddDeltav2a3 = ddDeltav2a3 * (2 * xj.y - xk.y - xi.y);
										ddDeltav2a4a1b = ddDeltav2a4a1b * (-1);
										ddDeltav2a4a2 = ddDeltav2a4a2 * (-(xi.y-xj.y)/rij);
										ddDeltav2a4c1b = ddDeltav2a4c1b * (-1);
										ddDeltav2a4c2 = ddDeltav2a4c2 * (-(xk.y-xj.y)/rkj);
										ddDeltav2a4b = ddDeltav2a4b * (-(xk.y-xj.y)/rkj);
										ddDeltav2a4d = ddDeltav2a4d * (-(xi.y-xj.y)/rij);
										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 = ddDeltav2b * (2 * rij * (-(xi.z-xj.z)/rij) * Math.pow(rkj,2) + Math.pow(rij,2) * 2 * rkj * (-(xk.z-xj.z)/rkj));
										ddDeltav2a1a = ddDeltav2a1a * -1;
										ddDeltav2a1b = ddDeltav2a1b * -1;
										ddDeltav2a2 = ddDeltav2a2 * ((-(xi.z-xj.z)/rij) * rkj + (-(xk.z-xj.z)/rkj) * rij);
										ddDeltav2a3 = ddDeltav2a3 * (2 * xj.z - xk.z - xi.z);
										ddDeltav2a4a1c = ddDeltav2a4a1c * (-1);
										ddDeltav2a4a2 = ddDeltav2a4a2 * (-(xi.z-xj.z)/rij);
										ddDeltav2a4c1c = ddDeltav2a4c1c * (-1);
										ddDeltav2a4c2 = ddDeltav2a4c2 * (-(xk.z-xj.z)/rkj);
										ddDeltav2a4b = ddDeltav2a4b * (-(xk.z-xj.z)/rkj);
										ddDeltav2a4d = ddDeltav2a4d * (-(xi.z-xj.z)/rij);
										break;
								}
							}
							if (angleAtomPosition[l][2] == atomNumbern) {

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

							ddDeltav2a4a1 = (ddDeltav2a4a1a + ddDeltav2a4a1b + ddDeltav2a4a1c) * rij;
							ddDeltav2a4c1 = (ddDeltav2a4c1a + ddDeltav2a4c1b + ddDeltav2a4c1c) * rkj;
							ddDeltav2a4 = ddDeltav2a4 * (((ddDeltav2a4a1 - ddDeltav2a4a2) / Math.pow(rij,2)) * rkj + ddDeltav2a4b + ((ddDeltav2a4c1 + ddDeltav2a4c2) / Math.pow(rkj,2)) * rij + ddDeltav2a4d);
							ddDeltav2 = ddDeltav2 * (((ddDeltav2a1a + ddDeltav2a1b) * ddDeltav2a1 + ddDeltav2a2 - ddDeltav2a3 - ddDeltav2a4) * ddDeltav2a - ddDeltav2b);
							ddDeltav[n][m][l] = ddDeltav1 + ddDeltav2;
						}
					}
					else {
						ddDeltav[n][m][l] = 0;
					}
					//logger.debug("ddDeltav[" + n + "][" + m + "][" + l + "] = " + ddDeltav[n][m][l]);
				}
			}
		}
		
	}


	/**
	 *  Get the angle bending second derivative respect to the cartesian coordinates of the atoms.
	 *
	 *@return        Delta angle bending second derivative value [dimension(3xN)] [angles Number]
	 */
	public double[][][] getAngleBendingSecondDerivative() {
		return ddDeltav;
	}


	/**
	 *  Evaluate a 2nd order approximation of the Hessian, for the angle bending,
	 *  given the atoms coordinates.
	 *
	 *@param  coord3d  Current molecule coordinates.
	 */
	public void setAngleBending2ndOrderErrorApproximateHessian(GVector coord3d) {
		angleBendingOrder2ndErrorApproximateHessian = new double[coord3d.getSize()][][];
		
		double sigma = Math.pow(0.000000000000001,0.33);
		GVector xminusSigma = new GVector(coord3d.getSize());
		GVector xplusSigma = new GVector(coord3d.getSize());
		double[][] gradientAtXminusSigma = null;
		double[][] gradientAtXplusSigma = null;
		
		for (int i = 0; i < coord3d.getSize(); i++) {
			xminusSigma.set(coord3d);
			xminusSigma.setElement(i,coord3d.getElement(i) - sigma);
			setAngleBending2ndOrderErrorApproximateGradient(xminusSigma);
			gradientAtXminusSigma = this.getAngleBending2ndOrderErrorApproximateGradient();
			xplusSigma.set(coord3d);
			xplusSigma.setElement(i,coord3d.getElement(i) + sigma);
			setAngleBending2ndOrderErrorApproximateGradient(xplusSigma);
			gradientAtXplusSigma = this.getAngleBending2ndOrderErrorApproximateGradient();
			angleBendingOrder2ndErrorApproximateHessian[i] = new double[coord3d.getSize()][];
			for (int j = 0; j < coord3d.getSize(); j++) {
				angleBendingOrder2ndErrorApproximateHessian[i][j] = new double[angleNumber];
				for (int k=0; k < angleNumber; k++) {
					angleBendingOrder2ndErrorApproximateHessian[i][j][k] = (gradientAtXplusSigma[j][k] - gradientAtXminusSigma[j][k]) / (2 * sigma);
				}
			}
		}
	}


	/**
	 *  Get the 2nd order error approximate Hessian for the angle bending.
	 *
	 *
	 *@return           Angle bending 2nd order error approximate Hessian values.
	 */
	public double[][][] getAngleBending2ndOrderErrorApproximateHessian() {
		return angleBendingOrder2ndErrorApproximateHessian;
	}


	/**
	 *  Evaluate the hessian for the angle bending.
	 *
	 *@param  coord3d  Current molecule coordinates.
	 */
	public void setHessianMMFF94SumEA(GVector coord3d) {

		double[] forHessian = new double[coord3d.getSize() * coord3d.getSize()];
		
		/*boolean[] sigmaChange = new boolean[coord3d.getSize()];
		for (int i=1; i < sigmaChange.length; i++) {sigmaChange[i] = false;}
		*/
		
		setDeltav(coord3d);
		setAngleBendingSecondDerivative(coord3d);
		
		double sumHessianEA = 0;

		int forHessianIndex;

		for (int n = 0; n < coord3d.getSize(); n++) {
			for (int m= 0; m < coord3d.getSize(); m++) {
				for (int l = 0; l < angleNumber; l++) {
					sumHessianEA = sumHessianEA + (2 * k2[l] + 6 * k3[l] * deltav[l]) * dDeltav[n][l] * dDeltav[m][l]
							+ (k2[l] * 2 * deltav[l] + k3[l] * 3 * Math.pow(deltav[l],2)) * ddDeltav[n][m][l];
				}

				forHessianIndex = n*coord3d.getSize()+m;
				forHessian[forHessianIndex] = sumHessianEA;
			}
		}
		
		/*for (int n = 0; n < forHessian.length; n++) {
			logger.debug(forHessian[n] + ", ");
			if (n % 6 == 5) {
				logger.debug("");
			}
		}*/
		
		hessianMMFF94SumEA = new GMatrix(coord3d.getSize(), coord3d.getSize(), forHessian);
		//logger.debug("hessianMMFF94SumEA : " + hessianMMFF94SumEA);
		
		NewtonRaphsonMethod nrm = new NewtonRaphsonMethod();
		nrm.hessianEigenValues(forHessian, coord3d.getSize());
		
	}


	/**
	 *  Get the hessian for the angle bending.
	 *
	 *@return        Hessian value of the angle bending term.
	 */
	public GMatrix getHessianMMFF94SumEA() {
		return hessianMMFF94SumEA;
	}


	/**
	 *  Evaluate a 2nd order approximation of the Hessian, for the angle bending energy term,
	 *  given the atoms coordinates.
	 *
	 *@param  coord3d  Current molecule coordinates.
	 */
	public void set2ndOrderErrorApproximateHessianMMFF94SumEA(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);
			setGradientMMFF94SumEA(xminusSigma);
			gradientAtXminusSigma.set(gradientMMFF94SumEA);
			xplusSigma.set(coord3d);
			xplusSigma.setElement(i,coord3d.getElement(i) + sigma);
			setGradientMMFF94SumEA(xplusSigma);
			gradientAtXplusSigma.set(gradientMMFF94SumEA);
			for (int j = 0; j < coord3d.getSize(); j++) {
				forHessianIndex = i*coord3d.getSize()+j;
				forOrder2ndErrorApproximateHessian[forHessianIndex] = (gradientAtXplusSigma.getElement(j) - gradientAtXminusSigma.getElement(j)) / (2 * sigma);
				//(functionMMFF94SumEA(xplusSigma) - 2 * fx + functionMMFF94SumEA(xminusSigma)) / Math.pow(sigma,2);
			}
		}
		
		order2ndErrorApproximateHessianMMFF94SumEA = new GMatrix(coord3d.getSize(), coord3d.getSize());
		order2ndErrorApproximateHessianMMFF94SumEA.set(forOrder2ndErrorApproximateHessian);
		//logger.debug("order2ndErrorApproximateHessianMMFF94SumEA : " + order2ndErrorApproximateHessianMMFF94SumEA);
	}


	/**
	 *  Get the 2nd order error approximate Hessian for the angle bending energy term.
	 *
	 *
	 *@return           Angle bending energy 2nd order error approximate Hessian value.
	 */
	public GMatrix get2ndOrderErrorApproximateHessianMMFF94SumEA() {
		return order2ndErrorApproximateHessianMMFF94SumEA;
	}

}