package org.openscience.cdk.modeling.forcefield;import java.util.Hashtable;import java.util.Vector;import javax.vecmath.GMatrix;import javax.vecmath.GVector;import org.openscience.cdk.interfaces.IAtomContainer;import org.openscience.cdk.qsar.IAtomicDescriptor;import org.openscience.cdk.qsar.descriptors.atomic.BondsToAtomDescriptor;import org.openscience.cdk.qsar.result.IntegerResult;//import org.openscience.cdk.tools.LoggingTool;/**			 *  Van Der Waals Interactions calculator for the potential energy function. Include function and derivatives. * *@author     vlabarta *@cdk.created    February 17, 2005 *@cdk.module     forcefield */public class VanDerWaalsInteractions {	String functionShape = " Van Der Waals Interactions ";	GVector currentCoordinates = null;	double mmff94SumEvdW = 0;	double[] vdWE1 = null;	double[] vdWE2 = null;	double ccgSumEvdWSlaterKirkwood = 0;	double ccgSumEvdWAverage = 0;	GVector gradientMMFF94SumEvdW = null;	double[] vdWEG1 = null;	double[] vdWEG2 = null;	GVector order2ndErrorApproximateGradientMMFF94SumEvdW = null;	GVector order5thErrorApproximateGradientMMFF94SumEvdW = null;	GVector gradientCCGSumEvdWSlaterKirkwood = null;	GVector gradientCCGSumEvdWAverage = null;	GMatrix hessianMMFF94SumEvdW = null;	double[] forHessian = null;		double[][] dR = null;	// Atom distance first order derivative respect to atoms coordinates	double[][][] ddR = null;	GVector dterm1 = null;	GVector dterm2 = null;	GVector ds = null;	GVector dt = null;	GVector dIvdw = null;		//int[][] distances = null;	//Better check common atom connected	IAtomicDescriptor shortestPathBetweenToAtoms = new BondsToAtomDescriptor();	Object[] params = {new Integer(0)};		int vdwInteractionNumber;	int[][] vdWiAtomPosition = null;	double[] eSK = null; 	// vdW well depths (mmff94: Slater-Kirkwood-based formula).	double[] asteriskR = null;	// minimum-energy separation in angstroms (mmff94).	double[] r = null;	// interatomic distance	double[]  atomRadiu0;	double[] eAv = null; 	// vdW well depths (Average).	double[] capitalR = null;	// minimum-energy separation in angstroms (Average).		double bb = 0.2;	double microB = 12;	// Parameters for the ccg vdWaals function	double a = 0.07;	// buffering constants (a,b). Prevent division by 0.	double b = 0.12;	double n = 7;	double m = 14;	double nij = n;	double mij = m - n;	double[] t = null;	double[] ivdw = null;	double vdwScale14 = 1;	// Scale factor for 1-4 interactions. To take in the future from mmff94.prm files.		//private LoggingTool logger;		/**	 *  Constructor for the VanDerWaalsInteractions object	 */	public VanDerWaalsInteractions() {        		//logger = new LoggingTool(this);	}	/**	 *  Set CCG Van Der Waals parameters for the molecule.	 *	 *	 *@param  molecule       The molecule like an AtomContainer object.	 *@param  parameterSet   MMFF94 parameters set	 *@exception  Exception  Description of the Exception	 */	public void setMMFF94VanDerWaalsParameters(IAtomContainer molecule, Hashtable parameterSet) throws Exception {		//distances = wnd.getShortestPathLengthBetweenAtoms((AtomContainer) molecule);				vdwInteractionNumber = 0;		for (int i=0; i<molecule.getAtomCount(); i++) {			for (int j=i+1; j<molecule.getAtomCount(); j++) {				params[0] = new Integer(j);				shortestPathBetweenToAtoms.setParameters(params);				//if (distances[molecule.getAtomNumber(molecule.getAtomAt(i))][molecule.getAtomNumber(molecule.getAtomAt(j))]>2) {				if (((IntegerResult)shortestPathBetweenToAtoms.calculate(molecule.getAtom(i),molecule).getValue()).intValue()>2){					vdwInteractionNumber += 1;				}			}		}		//logger.debug("vdwInteractionNumber : " + vdwInteractionNumber);				Vector vdwInteractionData = null;		eSK = new double[vdwInteractionNumber];		asteriskR = new double[vdwInteractionNumber];		eAv = new double[vdwInteractionNumber];		capitalR = new double[vdwInteractionNumber];		r = new double[vdwInteractionNumber];		atomRadiu0 = new double[molecule.getAtomCount()];		t = new double[vdwInteractionNumber];		ivdw = new double[vdwInteractionNumber];		double gI;	// To eSK calculation		double gJ;	// To eSK calculation		double alphaI;	// To eSK calculation and asteriskR[l]		double alphaJ;	// To eSK calculation and asteriskR[l]		double nI;		double nJ;		double aaI;	// To calculate asteriskRI		double aaJ;	// To calculate asteriskRJ		double asteriskRI;		double asteriskRJ;		double gamma;	// To calculate asteriskR[l]		double dI;		double dJ;		double eI;		double eJ;				vdWiAtomPosition = new int[vdwInteractionNumber][];		int l = -1;		for (int i=0; i<molecule.getAtomCount(); i++) {			for (int j=i+1; j<molecule.getAtomCount(); j++) {				params[0] = new Integer(j);				shortestPathBetweenToAtoms.setParameters(params);				//if (distances[molecule.getAtomNumber(molecule.getAtomAt(i))][molecule.getAtomNumber(molecule.getAtomAt(j))]>2) {				if (((IntegerResult)shortestPathBetweenToAtoms.calculate(molecule.getAtom(i),molecule).getValue()).intValue()>2){					l += 1;					vdwInteractionData = (Vector) parameterSet.get("data" + molecule.getAtom(i).getAtomTypeName());					//logger.debug("vdwInteractionData " + l + " : " + vdwInteractionData);					aaI = ((Double) vdwInteractionData.get(6)).doubleValue();					gI = ((Double) vdwInteractionData.get(7)).doubleValue();					alphaI = ((Double) vdwInteractionData.get(1)).doubleValue();					nI = ((Double) vdwInteractionData.get(2)).doubleValue();					eI = ((Double) vdwInteractionData.get(0)).doubleValue();					asteriskRI = aaI * Math.pow(alphaI,0.25);										vdwInteractionData = (Vector) parameterSet.get("data" + molecule.getAtom(j).getAtomTypeName());					//logger.debug("vdwInteractionData : " + vdwInteractionData);					aaJ = ((Double) vdwInteractionData.get(6)).doubleValue();					gJ = ((Double) vdwInteractionData.get(7)).doubleValue();					alphaJ = ((Double) vdwInteractionData.get(1)).doubleValue();					nJ = ((Double) vdwInteractionData.get(2)).doubleValue();					eJ = ((Double) vdwInteractionData.get(0)).doubleValue();					asteriskRJ = aaJ * Math.pow(alphaJ,0.25);										if (molecule.getAtom(i).getAtomTypeName() == molecule.getAtom(j).getAtomTypeName()) {						asteriskR[l] = asteriskRI;					} else {						gamma = (asteriskRI - asteriskRJ) / (asteriskRI + asteriskRJ);						asteriskR[l] = 0.5 * (asteriskRI + asteriskRJ) * (1 + bb * (1 - Math.exp((-1) * microB * Math.pow(gamma,2))));					}					eSK[l] = ((181.16 * gI * gJ * alphaI * alphaJ) / (Math.sqrt(alphaI/nI) + Math.sqrt(alphaJ/nJ))) * 1 / Math.pow(asteriskR[l], 6);					//logger.debug("eSK = " + eSK[l]);										vdwInteractionData = (Vector) parameterSet.get("vdw" + molecule.getAtom(i).getAtomTypeName());					//logger.debug("vdwInteractionData " + l + " : " + vdwInteractionData);					atomRadiu0[i] = ((Double) vdwInteractionData.get(0)).doubleValue();					vdwInteractionData = (Vector) parameterSet.get("vdw" + molecule.getAtom(j).getAtomTypeName());					atomRadiu0[j] = ((Double) vdwInteractionData.get(0)).doubleValue();										dI = 2 * atomRadiu0[i];					dJ = 2 * atomRadiu0[j];					capitalR[l] = (dI + dJ)/2;										eAv[l] = Math.sqrt(eI * eJ);										t[l] = 1;										params[0] = new Integer(j);					shortestPathBetweenToAtoms.setParameters(params);					if (((IntegerResult)shortestPathBetweenToAtoms.calculate(molecule.getAtom(i),molecule).getValue()).intValue()==3){					//if (distances[molecule.getAtomNumber(molecule.getAtomAt(i))][molecule.getAtomNumber(molecule.getAtomAt(j))] == 3) {						ivdw[l] = vdwScale14;					}else {						ivdw[l] = 1;					}										vdWiAtomPosition[l] = new int[2];					vdWiAtomPosition[l][0] = i;					vdWiAtomPosition[l][1] = j;				}			}		}				currentCoordinates = new GVector(molecule.getAtomCount());		vdWE1 = new double[vdwInteractionNumber];		vdWE2 = new double[vdwInteractionNumber];		vdWEG1 = new double[vdwInteractionNumber];		vdWEG2 = new double[vdwInteractionNumber];	}	/**	 *  Calculate the actual Rij	 *	 *@param  coords3d  Current molecule coordinates.	 */	public void setAtomDistance(GVector coords3d) {		for (int l = 0; l < vdwInteractionNumber; l++) {			r[l] = ForceFieldTools.distanceBetweenTwoAtomsFrom3xNCoordinates(coords3d, vdWiAtomPosition[l][0], vdWiAtomPosition[l][1]);			//logger.debug("r[" + l + "]= " + r[l]);		}	}	/**	 *  Get the atom distances values (Rij).	 *	 *@return        Atom distance values.	 */	public double[] getAtomDistance() {		return r;	}	/**	 *  Evaluate the MMFF94 Van Der Waals interaction energy given the atoms cartesian coordinates.	 *	 *@param  coords3d  Current molecule coordinates.	 */	public void setFunctionMMFF94SumEvdW(GVector coords3d) {		if (currentCoordinates.equals(coords3d)) {} 		else {			currentCoordinates.set(coords3d);			setAtomDistance(coords3d);			mmff94SumEvdW = 0;			for (int l = 0; l < vdwInteractionNumber; l++) {				//logger.debug(" ");				//logger.debug("eSK[" + l + "] = " + eSK[l]);				//logger.debug("asteriskR[" + l + "] = " + asteriskR[l]);				//logger.debug("r[" + l + "] = " + r[l]);				//logger.debug("term rest with minus 2 : " + (1.12 * Math.pow(asteriskR[l],7) / (Math.pow(r[l],7) + 0.12 * Math.pow(asteriskR[l],7))));				//logger.debug("vdwInteraction energy = " + (eSK[l] *					//(Math.pow(1.07 * asteriskR[l] / (r[l] + 0.07 * asteriskR[l]) ,7)) *					//((1.12 * Math.pow(asteriskR[l],7) / (Math.pow(r[l],7) + 0.12 * Math.pow(asteriskR[l],7))) - 2)));				vdWE1[l] = eSK[l] * (Math.pow(1.07 * asteriskR[l] / (r[l] + 0.07 * asteriskR[l]) ,7));				vdWE2[l] = (1.12 * Math.pow(asteriskR[l],7) / (Math.pow(r[l],7) + 0.12 * Math.pow(asteriskR[l],7))) - 2;				mmff94SumEvdW = mmff94SumEvdW + vdWE1[l] * vdWE2[l];				//logger.debug("mmff94SumEvdW = " + mmff94SumEvdW);			}			//mmff94SumEvdW = Math.abs(mmff94SumEvdW);		}		//logger.debug("mmff94SumEvdW = " + mmff94SumEvdW);	}	/**	 *  Get the MMFF94 Van Der Waals interaction energy for the current atoms coordinates.	 *	 *@return        MMFF94 Van Der Waals interaction energy value.	 */	public double getFunctionMMFF94SumEvdW() {		return mmff94SumEvdW;	}	/**	 *  Evaluate the CCG Van Der Waals interaction term.	 *	 *@param  coords3d  Current molecule coordinates.	 *@return        CCG Van Der Waals interaction term value.	 */	public double functionCCGSumEvdWSK(GVector coords3d, double[] s) {		if (currentCoordinates.equals(coords3d)) {} 		else {			setAtomDistance(coords3d);			ccgSumEvdWSlaterKirkwood = 0;			double c;			for (int l = 0; l < vdwInteractionNumber; l++) {				c = ((1+a) * asteriskR[l]) / (r[l] + a * asteriskR[l]);				ccgSumEvdWSlaterKirkwood = ccgSumEvdWSlaterKirkwood +								(eSK[l] *								(Math.pow(c,nij)) *								((nij/mij) * ((1+b) * Math.pow(asteriskR[l],mij) / (Math.pow(r[l],mij) + b * Math.pow(asteriskR[l],mij)))   - (mij + nij)/mij) *								s[l] * t[l] * ivdw[l]);			}		}		//logger.debug("ccgSumEvdWSlaterKirkwood = " + ccgSumEvdWSlaterKirkwood);		return ccgSumEvdWSlaterKirkwood;	}	/**	 *  Evaluate the CCG Van Der Waals interaction term.	 *	 *@param  coords3d  Current molecule coordinates.	 *@return        CCG Van Der Waals interaction term value.	 */	public double functionCCGSumEvdWAv(GVector coords3d, double[] s) {		if (currentCoordinates.equals(coords3d)) {} 		else {			setAtomDistance(coords3d);			ccgSumEvdWAverage = 0;			double c;			for (int l = 0; l < vdwInteractionNumber; l++) {				c = ((1+a) * capitalR[l]) / (r[l] + a * capitalR[l]);				ccgSumEvdWAverage = ccgSumEvdWAverage +								(eAv[l] *								(Math.pow(c,nij)) *								((nij/mij) * ((1+b) * Math.pow(capitalR[l],mij) / (Math.pow(r[l],mij) + b * Math.pow(capitalR[l],mij))) - (mij + nij)/mij) *								s[l] * t[l] * ivdw[l]);											}		}		//logger.debug("ccgSumEvdWAverage = " + ccgSumEvdWAverage);		return ccgSumEvdWAverage;	}	/**	 *  Calculate the atoms distances (Rij) first derivative respect to the cartesian coordinates of the atoms.	 *	 *@param  coord3d  Current molecule coordinates.	 */	public void setAtomsDistancesFirstOrderDerivative(GVector coord3d) {				dR = new double[coord3d.getSize()][];				Double forAtomNumber = null;		int atomNumber = 0;		int coordinate;		for (int i = 0; i < dR.length; i++) {						dR[i] = new double[vdwInteractionNumber];						forAtomNumber = new Double(i/3);			coordinate = i % 3;			//logger.debug("coordinate = " + coordinate);			atomNumber = forAtomNumber.intValue();			//logger.debug("atomNumber = " + atomNumber);			for (int j = 0; j < vdwInteractionNumber; j++) {				if ((vdWiAtomPosition[j][0] == atomNumber) | (vdWiAtomPosition[j][1] == atomNumber)) {					switch (coordinate) {						//x-coordinate						case 0: dR[i][j] = (coord3d.getElement(3 * vdWiAtomPosition[j][0]) - coord3d.getElement(3 * vdWiAtomPosition[j][1]))								/ Math.sqrt(Math.pow(coord3d.getElement(3 * vdWiAtomPosition[j][0]) - coord3d.getElement(3 * vdWiAtomPosition[j][1]),2) + Math.pow(coord3d.getElement(3 * vdWiAtomPosition[j][0] + 1) - coord3d.getElement(3 * vdWiAtomPosition[j][1] + 1),2) + Math.pow(coord3d.getElement(3 * vdWiAtomPosition[j][0] + 2) - coord3d.getElement(3 * vdWiAtomPosition[j][1] + 2),2)); 							break;						//y-coordinate						case 1:	dR[i][j] = (coord3d.getElement(3 * vdWiAtomPosition[j][0] + 1) - coord3d.getElement(3 * vdWiAtomPosition[j][1] + 1))								/ Math.sqrt(Math.pow(coord3d.getElement(3 * vdWiAtomPosition[j][0]) - coord3d.getElement(3 * vdWiAtomPosition[j][1]),2) + Math.pow(coord3d.getElement(3 * vdWiAtomPosition[j][0] + 1) - coord3d.getElement(3 * vdWiAtomPosition[j][1] + 1),2) + Math.pow(coord3d.getElement(3 * vdWiAtomPosition[j][0] + 2) - coord3d.getElement(3 * vdWiAtomPosition[j][1] + 2),2)); 							break;						//z-coordinate						case 2: dR[i][j] = (coord3d.getElement(3 * vdWiAtomPosition[j][0] + 2) - coord3d.getElement(3 * vdWiAtomPosition[j][1] + 2))								/ Math.sqrt(Math.pow(coord3d.getElement(3 * vdWiAtomPosition[j][0]) - coord3d.getElement(3 * vdWiAtomPosition[j][1]),2) + Math.pow(coord3d.getElement(3 * vdWiAtomPosition[j][0] + 1) - coord3d.getElement(3 * vdWiAtomPosition[j][1] + 1),2) + Math.pow(coord3d.getElement(3 * vdWiAtomPosition[j][0] + 2) - coord3d.getElement(3 * vdWiAtomPosition[j][1] + 2),2)); 							break;					}					if (vdWiAtomPosition[j][1] == atomNumber) {						dR[i][j] = (-1) * dR[i][j];					}				} else {					dR[i][j] = 0;				}				//logger.debug("vdW Interaction " + j + " : " + "dR[" + i + "][" + j + "] = " + dR[i][j]);			}		}	}	/**	 *  Get the atoms distances first order derivative respect to the cartesian coordinates of the atoms.	 *	 *@return        Atoms distances first order derivative value [dimension(3xN)] [vdW interaction number]	 */	public double[][] getAtomsDistancesFirstDerivative() {		return dR;	}