bondstretching.java

化学图形处理软件

				//logger.debug("dDeltar = " + dDeltar[i][j]);
				//logger.debug("gradient " + i + "bond " + j + " : " + (k2[j] * 2 * deltar[j] + k3[j] * 3 * Math.pow(deltar[j],2) + k4[j] * 4 * Math.pow(deltar[j],3)) * dDeltar[i][j]);
				sumGradientEB = sumGradientEB + (k2[j] * 2 * deltar[j] + k3[j] * 3 * Math.pow(deltar[j],2) + k4[j] * 4 * Math.pow(deltar[j],3)) * dDeltar[i][j];
				//logger.debug(sumGradientEB);
			}
			//sumGradientEB = (-1) * sumGradientEB;
			gradientMMFF94SumEB.setElement(i, sumGradientEB);
		}
		//logger.debug("gradientMMFF94 = " + gradientMMFF94SumEB);
	}


	/**
	 *  Get the gradient for the bond stretching in a given atoms coordinates
	 *
	 *@return           Bond stretching gradient value
	 */
	public GVector getGradientMMFF94SumEB() {
		return gradientMMFF94SumEB;
	}


	/**
	 *  Calculate the bond lengths second derivative respect to the cartesian coordinates of the atoms.
	 *
	 *@param  coord3d  Current molecule coordinates.
	 */
	public void setBondLengthsSecondDerivative(GVector coord3d, double[] deltar, int[][] bondAtomPosition) {
		ddDeltar = new double[coord3d.getSize()][][];
		
		Double forAtomNumber = null;
		int atomNumberi;
		int atomNumberj;
		int coordinatei;
		int coordinatej;
		double ddDeltar1=0;	// ddDeltar[i][j][k] = ddDeltar1 - ddDeltar2
		double ddDeltar2=0;
		
		setBondLengthsFirstDerivative(coord3d, deltar, bondAtomPosition);
		//logger.debug("bondAtomPosition.length = " + bondAtomPosition.length);
		double[] rTemp = new double[bondAtomPosition.length];
		for (int i = 0; i < bondAtomPosition.length; i++) {
			rTemp[i] = ForceFieldTools.distanceBetweenTwoAtomsFrom3xNCoordinates(coord3d, bondAtomPosition[i][0], bondAtomPosition[i][1]);
		}
		
		for (int i=0; i<coord3d.getSize(); i++) {
			ddDeltar[i] = new double[coord3d.getSize()][];
			
			forAtomNumber = new Double(i/3);
			
			atomNumberi = forAtomNumber.intValue();
			//logger.debug("atomNumberi = " + atomNumberi);
				
			coordinatei = i % 3;
			//logger.debug("coordinatei = " + coordinatei);
				
			for (int j=0; j<coord3d.getSize(); j++) {
				ddDeltar[i][j] = new double[deltar.length];
				
				forAtomNumber = new Double(j/3);

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

				coordinatej = j % 3;
				//logger.debug("coordinatej = " + coordinatej);
				
				//logger.debug("atomj : " + molecule.getAtomAt(atomNumberj));
				
				for (int k=0; k < deltar.length; k++) {
					
					if ((bondAtomPosition[k][0] == atomNumberj) | (bondAtomPosition[k][1] == atomNumberj)) {
						if ((bondAtomPosition[k][0] == atomNumberi) | (bondAtomPosition[k][1] == atomNumberi)) {
					
							// ddDeltar1
							if (bondAtomPosition[k][0] == atomNumberj) {
								ddDeltar1 = 1;
							}
							if (bondAtomPosition[k][1] == atomNumberj) {
								ddDeltar1 = -1;
							}
							if (bondAtomPosition[k][0] == atomNumberi) {
								ddDeltar1 = ddDeltar1 * 1;
							}
							if (bondAtomPosition[k][1] == atomNumberi) {
								ddDeltar1 = ddDeltar1 * (-1);
							}
							ddDeltar1 = ddDeltar1 / rTemp[k];

							// ddDeltar2
							switch (coordinatej) {
								case 0: ddDeltar2 = (coord3d.getElement(3 * bondAtomPosition[k][0]) - coord3d.getElement(3 * bondAtomPosition[k][1]));
									//logger.debug("OK: d1 x");
									break;
								case 1:	ddDeltar2 = (coord3d.getElement(3 * bondAtomPosition[k][0] + 1) - coord3d.getElement(3 * bondAtomPosition[k][1] + 1));
									//logger.debug("OK: d1 y");
									break;
								case 2:	ddDeltar2 = (coord3d.getElement(3 * bondAtomPosition[k][0] + 2) - coord3d.getElement(3 * bondAtomPosition[k][1] + 2));
									//logger.debug("OK: d1 z");
									break;
							}
						
							if (bondAtomPosition[k][1] == atomNumberj) {
								ddDeltar2 = (-1) * ddDeltar2;
								//logger.debug("OK: bond 1");
							} 
	
							switch (coordinatei) {
								case 0: ddDeltar2 = ddDeltar2 * (coord3d.getElement(3 * bondAtomPosition[k][0]) - coord3d.getElement(3 * bondAtomPosition[k][1]));
									//logger.debug("OK: have d2 x");
									break;
								case 1:	ddDeltar2 = ddDeltar2 * (coord3d.getElement(3 * bondAtomPosition[k][0] + 1) - coord3d.getElement(3 * bondAtomPosition[k][1] + 1));
									//logger.debug("OK: have d2 y");
									break;
								case 2: ddDeltar2 = ddDeltar2 * (coord3d.getElement(3 * bondAtomPosition[k][0] + 2) - coord3d.getElement(3 * bondAtomPosition[k][1] + 2));
									//logger.debug("OK: have d2 z");
									break;
							}
							
							if (bondAtomPosition[k][1] == atomNumberi) {
								ddDeltar2 = (-1) * ddDeltar2;
								//logger.debug("OK: d2 bond 1");
							}
							
							ddDeltar2 = ddDeltar2 / Math.pow(rTemp[k],2);
							
							// ddDeltar[i][j][k]
							ddDeltar[i][j][k] = ddDeltar1 - ddDeltar2;
						} else {
							ddDeltar[i][j][k] = 0;
							//logger.debug("OK: 0");
						}
					} else {
						ddDeltar[i][j][k] = 0;
						//logger.debug("OK: 0");
					}
					//logger.debug("bond " + k + " : " + "ddDeltar[" + i + "][" + j + "][" + k + "] = " + ddDeltar[i][j][k]);
				}
			}
		}	
	}


	/**
	 *  Get the bond lengths second derivative respect to the cartesian coordinates of the atoms.
	 *
	 *@return        Bond lengths second derivative value [dimension(3xN)] [bonds Number]
	 */
	 public double[][][] getBondLengthsSecondDerivative() {
		return ddDeltar;
	}


	/**
	 *  Evaluate the second order partial derivative (hessian) for the bond stretching given the atoms coordinates
	 *
	 *@param  coord3d  Current molecule coordinates.
	 */
	public void setHessianMMFF94SumEB(GVector coord3d) {
		
		forHessian = new double[coord3d.getSize() * coord3d.getSize()];
		
		calculateDeltar(coord3d);
		setBondLengthsSecondDerivative(coord3d, deltar, bondAtomPosition);
		
		double sumHessianEB;
		int forHessianIndex;
		for (int i = 0; i < coord3d.getSize(); i++) {
			for (int j = 0; j < coord3d.getSize(); j++) {
				sumHessianEB = 0;
				for (int k = 0; k < bondsNumber; k++) {
					sumHessianEB = sumHessianEB + (2 * k2[k] + 6 * k3[k] * deltar[k] + 12 * k4[k] * Math.pow(deltar[k],2)) * dDeltar[i][k] * dDeltar[j][k]
								+ (k2[k] * 2 * deltar[k] + k3[k] * 3 * Math.pow(deltar[k],2) + k4[k] * 4 * Math.pow(deltar[k],3)) * ddDeltar[i][j][k];
				}
				forHessianIndex = i*coord3d.getSize()+j;
				forHessian[forHessianIndex] = sumHessianEB;
			}
		}

		hessianMMFF94SumEB = new GMatrix(coord3d.getSize(), coord3d.getSize(),forHessian);
		//logger.debug("hessianMMFF94SumEB : " + hessianMMFF94SumEB);
	}


	/**
	 *  Get the hessian for the bond stretching.
	 *
	 *@return        Hessian value of the bond stretching term.
	 */
	public GMatrix getHessianMMFF94SumEB() {
		return hessianMMFF94SumEB;
	}


	/**
	 *  Get the hessian for the bond stretching.
	 *
	 *@return        Hessian value of the bond stretching term.
	 */
	public double[] getForHessianMMFF94SumEB() {
		return forHessian;
	}


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


	/**
	 *  Get the 2nd order error approximate Hessian for the bond stretching term.
	 *
	 *
	 *@return           Bond stretching 2nd order error approximate Hessian value.
	 */
	public GMatrix get2ndOrderErrorApproximateHessianMMFF94SumEB() {
		return order2ndErrorApproximateHessianMMFF94SumEB;
	}

}