/**	 *  Set the gradient of the MMFF94 Van Der Waals interaction term.	 *	 *	 *@param  coords3d  Current molecule coordinates.	 */	public void setGradientMMFF94SumEvdW(GVector coords3d) {		gradientMMFF94SumEvdW = new GVector(coords3d.getSize());		if (currentCoordinates.equals(coords3d)) {} 		else {setFunctionMMFF94SumEvdW(coords3d);}		setAtomsDistancesFirstOrderDerivative(coords3d);		for (int l = 0; l < vdwInteractionNumber; l++) {			vdWEG1[l] = eSK[l] * 7 * Math.pow((1.07 * asteriskR[l]) / (r[l] + 0.07 * asteriskR[l]),6) * 					1.07 * asteriskR[l] * (-1) * (1/Math.pow((r[l] + 0.07 * asteriskR[l]),2));			vdWEG2[l] = 1.12 * Math.pow(asteriskR[l],7) * 					(-1) * (1/Math.pow((Math.pow(r[l],7) + 0.12 * Math.pow(asteriskR[l],7)),2)) * 7 * Math.pow(r[l],6);		}				double sumGradientEvdW;		for (int i = 0; i < gradientMMFF94SumEvdW.getSize(); i++) {			sumGradientEvdW = 0;			for (int l = 0; l < vdwInteractionNumber; l++) {				sumGradientEvdW = sumGradientEvdW + vdWEG1[l] * dR[i][l] * vdWE2[l] + vdWE1[l] * vdWEG2[l]  * dR[i][l];			}			gradientMMFF94SumEvdW.setElement(i, sumGradientEvdW);		}		//logger.debug("gradientMMFF94SumEvdW = " + gradientMMFF94SumEvdW);	}	/**	 *  Get the gradient of the MMFF94 Van Der Waals interaction term.	 *	 *	 *@return           MMFF94 Van Der Waals interaction gradient value.	 */	public GVector getGradientMMFF94SumEvdW() {		return gradientMMFF94SumEvdW;	}	/**	 *  Evaluate a 2nd order error approximation of the gradient, for the Van Der Waals interaction term, given the atoms	 *  coordinates	 *	 *@param  coord3d  Current molecule coordinates.	 */	public void set2ndOrderErrorApproximateGradientMMFF94SumEvdW(GVector coord3d) {		order2ndErrorApproximateGradientMMFF94SumEvdW = 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 fxplusSigma = 0;		double fxminusSigma = 0;				for (int m = 0; m < order2ndErrorApproximateGradientMMFF94SumEvdW.getSize(); m++) {			xplusSigma.set(coord3d);			xplusSigma.setElement(m,coord3d.getElement(m) + sigma);			setFunctionMMFF94SumEvdW(xplusSigma);			fxplusSigma = getFunctionMMFF94SumEvdW();			xminusSigma.set(coord3d);			xminusSigma.setElement(m,coord3d.getElement(m) - sigma);			setFunctionMMFF94SumEvdW(xminusSigma);			fxminusSigma = getFunctionMMFF94SumEvdW();			order2ndErrorApproximateGradientMMFF94SumEvdW.setElement(m,(fxplusSigma - fxminusSigma) / (2 * sigma));		}					//logger.debug("order2ndErrorApproximateGradientMMFF94SumEvdW : " + order2ndErrorApproximateGradientMMFF94SumEvdW);	}	/**	 *  Get the 2nd order error approximate gradient for the Van Der Waals interaction term.	 *	 *	 *@return           Van Der Waals interaction 2nd order error approximate gradient value	 */	public GVector get2ndOrderErrorApproximateGradientMMFF94SumEvdW() {		return order2ndErrorApproximateGradientMMFF94SumEvdW;	}	/**	 *  Evaluate a 5th order error approximation of the gradient, for the Van Der Waals interaction term, given the atoms	 *  coordinates	 *	 *@param  coord3d  Current molecule coordinates.	 */	public void set5thOrderErrorApproximateGradientMMFF94SumEvdW(GVector coord3d) {		order5thErrorApproximateGradientMMFF94SumEvdW = 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 fxplusSigma = 0;		double fxminusSigma = 0;		double fxplus2Sigma = 0;		double fxminus2Sigma = 0;				for (int m=0; m < order5thErrorApproximateGradientMMFF94SumEvdW.getSize(); m++) {			xplusSigma.set(coord3d);			xplusSigma.setElement(m,coord3d.getElement(m) + sigma);			setFunctionMMFF94SumEvdW(xplusSigma);			fxplusSigma = getFunctionMMFF94SumEvdW();			xminusSigma.set(coord3d);			xminusSigma.setElement(m,coord3d.getElement(m) - sigma);			setFunctionMMFF94SumEvdW(xminusSigma);			fxminusSigma = getFunctionMMFF94SumEvdW();			xplus2Sigma.set(coord3d);			xplus2Sigma.setElement(m,coord3d.getElement(m) + 2 * sigma);			setFunctionMMFF94SumEvdW(xplus2Sigma);			fxplus2Sigma = getFunctionMMFF94SumEvdW();			xminus2Sigma.set(coord3d);			xminus2Sigma.setElement(m,coord3d.getElement(m) - 2 * sigma);			setFunctionMMFF94SumEvdW(xminus2Sigma);			fxminus2Sigma = getFunctionMMFF94SumEvdW();			order5thErrorApproximateGradientMMFF94SumEvdW.setElement(m, (8 * (fxplusSigma - fxminusSigma) - (fxplus2Sigma - fxminus2Sigma)) / (12 * sigma));		}					//logger.debug("order5thErrorApproximateGradientMMFF94SumEvdW : " + order5thErrorApproximateGradientMMFF94SumEvdW);	}	/**	 *  Get the 5th order error approximate gradient for the Van Der Waals interaction term.	 *	 *@return        Torsion 5th order error approximate gradient value.	 */	public GVector get5OrderApproximateGradientMMFF94SumEvdW() {		return order5thErrorApproximateGradientMMFF94SumEvdW;	}	/**	 *  Evaluate the gradient of the CCG Van Der Waals interaction term.	 *  	 *	 *@param  coords3d  Current molecule coordinates.	 *@return           CCG Van Der Waals interaction gradient value.	 *//*	public GVector gradientMMFF94SumEvdW(GVector coords3d, double[] s) {		gradientCCGSumEvdWSlaterKirkwood.setSize(molecule.getAtomCount() * 3);		setAtomDistance(coords3d);		dR.setSize(molecule.getAtomCount() * 3);		dterm1.setSize(molecule.getAtomCount() * 3);		dterm2.setSize(molecule.getAtomCount() * 3);				double c;		double[] term1 = new double[vdwInteractionNumber];		double[] term2 = new double[vdwInteractionNumber];		double sumGradientEvdW;		for (int i = 0; i < gradientCCGSumEvdWSlaterKirkwood.getSize(); i++) {			dterm1.setElement(i,1);                 // dterm1 : partial derivative of term1. To change in the future			dterm2.setElement(i,1);                 // dterm2 : partial derivative of term2. To change in the future			ds.setElement(i,1);                 // ds : partial derivative of s. To change in the future			dt.setElement(i,1);                 // dt : partial derivative of t. To change in the future			dIvdw.setElement(i,1);                 // dIvdw : partial derivative of Ivdw. To change in the future			sumGradientEvdW = 0;			for (int l = 0; l < vdwInteractionNumber; l++) {								c = ((1+a) * asteriskR[l]) / (r[l] + a * asteriskR[l]);				term1[l] = Math.pow(c,nij);				term2[l] = (nij/mij) * ((1+b) * Math.pow(asteriskR[l],mij) / (Math.pow(r[l],mij) + b * Math.pow(asteriskR[l],mij)))   - (mij + nij)/mij;								sumGradientEvdW = sumGradientEvdW + (deSK.getElement(i) * term1[l] * term2[l] + eSK[l] * (dterm1.getElement(i) * term2[l] + term1[l] * dterm2.getElement(i))) * s[l] * t[l] * ivdw[l] + 								(eSK[l] * term1[l] * term2[l]) * (ds.getElement(i) * t[l] * ivdw[l] + s[l] * (dt.getElement(i) * ivdw[l] + t[l] * dIvdw.getElement(i)));			}			sumGradientEvdW = sumGradientEvdW * 2.51210;						gradientCCGSumEvdWSlaterKirkwood.setElement(i, sumGradientEvdW);		}		//logger.debug("gradientCCGSumEvdWSlaterKirkwood = " + gradientCCGSumEvdWSlaterKirkwood);		return gradientCCGSumEvdWSlaterKirkwood;	}*/	/**	 *  Evaluate the hessian of the CCG Van Der Waals interaction term.	 *	 *@param  coords3d  Current molecule coordinates.	 *@return        Hessian value of the CCG Van Der Waals interaction term.	 *//*	public GMatrix hessian(GVector coords3d) {		double[] forHessian = new double[coords3d.getSize() * coords3d.getSize()];		setAtomDistance(coords3d);		double sumHessianEvdW = 0;		GMatrix ddR = new GMatrix(coords3d.getSize(),coords3d.getSize());		ddR.setZero();		for (int i = 0; i < forHessian.length; i++) {			for (int j = 0; j < vdwInteractionNumber; j++) {				sumHessianEvdW = sumHessianEvdW + 1;			}			forHessian[i] = sumHessianEvdW;		}		hessianMMFF94SumEvdW.setSize(coords3d.getSize(), coords3d.getSize());		hessianMMFF94SumEvdW.set(forHessian); 		//logger.debug("hessianMMFF94SumEvdW : " + hessianMMFF94SumEvdW);		return hessianMMFF94SumEvdW;	}*/	/**	 *  Calculate the bond lengths second derivative respect to the cartesian coordinates of the atoms.	 *	 *@param  coord3d  Current molecule coordinates.	 */	public void setBondLengthsSecondDerivative(GVector coord3d) {		ddR = new double[coord3d.getSize()][][];				Double forAtomNumber = null;		int atomNumberi;		int atomNumberj;		int coordinatei;		int coordinatej;		double ddR1=0;	// ddR[i][j][k] = ddR1 - ddR2		double ddR2=0;				setAtomsDistancesFirstOrderDerivative(coord3d);				for (int i=0; i<coord3d.getSize(); i++) {			ddR[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++) {				ddR[i][j] = new double[vdwInteractionNumber];								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 < vdwInteractionNumber; k++) {										if ((vdWiAtomPosition[k][0] == atomNumberj) | (vdWiAtomPosition[k][1] == atomNumberj)) {						if ((vdWiAtomPosition[k][0] == atomNumberi) | (vdWiAtomPosition[k][1] == atomNumberi)) {												// ddR1							if (vdWiAtomPosition[k][0] == atomNumberj) {								ddR1 = 1;							}							if (vdWiAtomPosition[k][1] == atomNumberj) {								ddR1 = -1;							}							if (vdWiAtomPosition[k][0] == atomNumberi) {								ddR1 = ddR1 * 1;							}							if (vdWiAtomPosition[k][1] == atomNumberi) {								ddR1 = ddR1 * (-1);							}							ddR1 = ddR1 / r[k];							// ddR2							switch (coordinatej) {								case 0: ddR2 = (coord3d.getElement(3 * vdWiAtomPosition[k][0]) - coord3d.getElement(3 * vdWiAtomPosition[k][1]));									//logger.debug("OK: d1 x");									break;								case 1:	ddR2 = (coord3d.getElement(3 * vdWiAtomPosition[k][0] + 1) - coord3d.getElement(3 * vdWiAtomPosition[k][1] + 1));									//logger.debug("OK: d1 y");									break;								case 2:	ddR2 = (coord3d.getElement(3 * vdWiAtomPosition[k][0] + 2) - coord3d.getElement(3 * vdWiAtomPosition[k][1] + 2));									//logger.debug("OK: d1 z");									break;							}													if (vdWiAtomPosition[k][1] == atomNumberj) {								ddR2 = (-1) * ddR2;								//logger.debug("OK: bond 1");							} 								switch (coordinatei) {								case 0: ddR2 = ddR2 * (coord3d.getElement(3 * vdWiAtomPosition[k][0]) - coord3d.getElement(3 * vdWiAtomPosition[k][1]));									//logger.debug("OK: have d2 x");									break;								case 1:	ddR2 = ddR2 * (coord3d.getElement(3 * vdWiAtomPosition[k][0] + 1) - coord3d.getElement(3 * vdWiAtomPosition[k][1] + 1));									//logger.debug("OK: have d2 y");									break;								case 2: ddR2 = ddR2 * (coord3d.getElement(3 * vdWiAtomPosition[k][0] + 2) - coord3d.getElement(3 * vdWiAtomPosition[k][1] + 2));									//logger.debug("OK: have d2 z");									break;							}														if (vdWiAtomPosition[k][1] == atomNumberi) {								ddR2 = (-1) * ddR2;								//logger.debug("OK: d2 bond 1");							}														ddR2 = ddR2 / Math.pow(r[k],2);														// ddR[i][j][k]							ddR[i][j][k] = ddR1 - ddR2;						} else {							ddR[i][j][k] = 0;							//logger.debug("OK: 0");						}					} else {						ddR[i][j][k] = 0;						//logger.debug("OK: 0");					}					//logger.debug("bond " + k + " : " + "ddR[" + i + "][" + j + "][" + k + "] = " + ddR[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 ddR;	}	/**	 *  Evaluate the second order partial derivative (hessian) for the van der Waals interactions given the atoms coordinates	 *	 *@param  coord3d  Current molecule coordinates.	 */	public void setHessianMMFF94SumEvdW(GVector coord3d) {				forHessian = new double[coord3d.getSize() * coord3d.getSize()];				if (currentCoordinates.equals(coord3d)) {} 		else {setFunctionMMFF94SumEvdW(coord3d);}				setBondLengthsSecondDerivative(coord3d);				double sumHessianEvdW;		int forHessianIndex;		double vdWEHessian1 = 0;		double vdWEHessian2 = 0;		double vdWESD1 = 0;		double vdWESD2 = 0;		for (int i = 0; i < coord3d.getSize(); i++) {			for (int j = 0; j < coord3d.getSize(); j++) {				sumHessianEvdW = 0;				for (int k = 0; k < vdwInteractionNumber; k++) {					vdWESD1 = 89.47 * eSK[k] * Math.pow(asteriskR[k], 7) * (1/Math.pow(r[k] + 0.07 * asteriskR[k], 9)) * dR[i][k];					vdWESD2 = -7.84 * Math.pow(asteriskR[k],7) * (6 * Math.pow(r[k],5) * dR[i][k] * Math.pow(Math.pow(r[k],7) + 0.12 * Math.pow(asteriskR[k],7),2) - Math.pow(r[k],12) * 14 * (Math.pow(r[k],7) + 0.12 * Math.pow(asteriskR[k],7)) * dR[i][k])/Math.pow(Math.pow(r[k],7) + 0.12 * Math.pow(asteriskR[k],7),4);					vdWEHessian1 = (vdWESD1 * dR[j][k] + vdWEG1[k] * ddR[i][j][k]) * vdWE2[k] + vdWEG1[k] * (vdWEG2[k] * dR[i][k]);					vdWEHessian2 = (vdWEG1[k] * dR[i][k]) * (vdWEG2[k]) + (vdWE1[k]) * (vdWESD2 * dR[j][k] + vdWEG2[k] * ddR[i][j][k]);					sumHessianEvdW = sumHessianEvdW + (vdWEHessian1 + vdWEHessian2);				}				forHessianIndex = i*coord3d.getSize()+j;				forHessian[forHessianIndex] = sumHessianEvdW;			}		}				hessianMMFF94SumEvdW = new GMatrix(coord3d.getSize(), coord3d.getSize(), forHessian);		//logger.debug("hessianMMFF94SumEvdW : " + hessianMMFF94SumEvdW);	}	/**	 *  Get the hessian for the van der Waals interactions.	 *	 *@return        Hessian value of the van der Waals interactions term.	 */	public GMatrix getHessianMMFF94SumEvdW() {		return hessianMMFF94SumEvdW;	}	/**	 *  Get the hessian for the van der Waals interactions.	 *	 *@return        Hessian value of the van der Waals interactions term.	 */	public double[] getForHessianMMFF94SumEvdW() {		return forHessian;	}}