package org.openscience.cdk.modeling.forcefield;import org.openscience.cdk.interfaces.IAtom;import org.openscience.cdk.interfaces.IAtomContainer;import org.openscience.cdk.interfaces.IBond;import org.openscience.cdk.modeling.builder3d.MMFF94ParametersCall;import org.openscience.cdk.tools.manipulator.AtomContainerManipulator;import org.openscience.cdk.tools.manipulator.BondManipulator;import javax.vecmath.GMatrix;import javax.vecmath.GVector;import java.util.Hashtable;import java.util.Iterator;import java.util.Vector;//import org.openscience.cdk.tools.LoggingTool;/** *  Torsions calculator for the potential energy function. Include function and derivatives. * *@author     vlabarta *@cdk.created    March 2, 2005 *@cdk.module     forcefield */public class Torsions {	String functionShape = " Torsions ";	double mmff94SumET = 0;	GVector gradientMMFF94SumET = new GVector(3);	GVector dPhi = new GVector(3);		GVector order2ndErrorApproximateGradientMMFF94SumET = new GVector(3);	GVector order5thErrorApproximateGradientMMFF94SumET = new GVector(3);	GVector xplusSigma = null;	GVector xminusSigma = null;	double sigma = Math.pow(0.000000000000001,0.33);		GMatrix hessianMMFF94SumET = null;	double[] forHessian = null;	GMatrix order2ndErrorApproximateHessianMMFF94SumET = null;	double[] forOrder2ndErrorApproximateHessian = null;	int torsionNumber = 0;	int[][] torsionAtomPosition = null;	double[] v1 = null;	double[] v2 = null;	double[] v3 = null;	double[] phi = null;		IBond[] bond = null;	IAtom[] atomInBond = null;	IBond[] bondConnectedBefore = null;	IBond[] bondConnectedAfter = null;	//private LoggingTool logger;	GVector moleculeCurrentCoordinates = null;	boolean[] changeAtomCoordinates = null;	int changedCoordinates;	/**	 *  Constructor for the Torsions object	 */	public Torsions() {        		//logger = new LoggingTool(this);	}	/**	 *  Set MMFF94 constants V1, V2 and V3 for each i-j, j-k and k-l bonded pairs in the molecule.	 *	 *	 *@param  molecule       The molecule like an AtomContainer object.	 *@param  parameterSet   MMFF94 parameters set	 *@exception  Exception  Description of the Exception	 */	public void setMMFF94TorsionsParameters(IAtomContainer molecule, Hashtable parameterSet) throws Exception {		//logger.debug("setMMFF94TorsionsParameters");        // looks like we need the bonds in an array for the rest of the class        bond = new IBond[molecule.getBondCount()];        int counter = 0;        Iterator bonds = molecule.bonds();        while (bonds.hasNext()) {            IBond aBond = (IBond) bonds.next();            bond[counter] = aBond;            counter++;        }                		for (int b=0; b<bond.length; b++) {			atomInBond = BondManipulator.getAtomArray(bond[b]);			bondConnectedBefore = AtomContainerManipulator.getBondArray(molecule.getConnectedBondsList(atomInBond[0]));			if (bondConnectedBefore.length > 1) {				bondConnectedAfter = AtomContainerManipulator.getBondArray(molecule.getConnectedBondsList(atomInBond[1]));				if (bondConnectedAfter.length > 1) {					for (int bb=0; bb<bondConnectedBefore.length; bb++) {						if (bondConnectedBefore[bb].compare(bond[b])) {}						else {							for (int ba=0; ba<bondConnectedAfter.length; ba++) {								if (bondConnectedAfter[ba].compare(bond[b])) {}								else {									if (bondConnectedBefore[bb].isConnectedTo(bondConnectedAfter[ba])) {}									else {										torsionNumber += 1;										//logger.debug("atomi(" + torsionNumber + ") : " + bondConnectedBefore[bb].getConnectedAtom(atomInBond[0]).getAtomTypeName());										//logger.debug("atomj(" + torsionNumber + ") : " + atomInBond[0].getAtomTypeName());										//logger.debug("atomk(" + torsionNumber + ") : " + atomInBond[1].getAtomTypeName());										//logger.debug("atoml(" + torsionNumber + ") : " + bondConnectedAfter[ba].getConnectedAtom(atomInBond[1]).getAtomTypeName());									}									}							}						}					}				}			}		}		//logger.debug("torsionNumber = " + torsionNumber);		Vector torsionsData = null;		MMFF94ParametersCall pc = new MMFF94ParametersCall();		pc.initialize(parameterSet);				v1 = new double[torsionNumber];		v2 = new double[torsionNumber];		v3 = new double[torsionNumber];		torsionAtomPosition = new int[torsionNumber][];		String torsionType;		int m = -1;		for (int b=0; b<bond.length; b++) {			atomInBond = BondManipulator.getAtomArray(bond[b]);			bondConnectedBefore = AtomContainerManipulator.getBondArray(molecule.getConnectedBondsList(atomInBond[0]));			if (bondConnectedBefore.length > 1) {				bondConnectedAfter = AtomContainerManipulator.getBondArray(molecule.getConnectedBondsList(atomInBond[1]));				if (bondConnectedAfter.length > 1) {					for (int bb=0; bb<bondConnectedBefore.length; bb++) {						if (bondConnectedBefore[bb].compare(bond[b])) {}						else {							for (int ba=0; ba<bondConnectedAfter.length; ba++) {								if (bondConnectedAfter[ba].compare(bond[b])) {}								else {									if (bondConnectedBefore[bb].isConnectedTo(bondConnectedAfter[ba])) {}									else {										m += 1;										torsionAtomPosition[m] = new int[4];										torsionAtomPosition[m][0] = molecule.getAtomNumber(bondConnectedBefore[bb].getConnectedAtom(atomInBond[0]));										torsionAtomPosition[m][1] = molecule.getAtomNumber(atomInBond[0]);										torsionAtomPosition[m][2] = molecule.getAtomNumber(atomInBond[1]);										torsionAtomPosition[m][3] = molecule.getAtomNumber(bondConnectedAfter[ba].getConnectedAtom(atomInBond[1]));																			/*System.out.println("torsion " + m + " : " + bondConnectedBefore[bb].getConnectedAtom(atomInBond[0]).getFlag(CDKConstants.ISINRING) + "(" + torsionAtomPosition[m][0] + "), " + 												atomInBond[0].getFlag(CDKConstants.ISINRING) + "(" + torsionAtomPosition[m][1] + "), " + atomInBond[1].getFlag(CDKConstants.ISINRING) + "(" + torsionAtomPosition[m][2] + "), " + 												bondConnectedAfter[ba].getConnectedAtom(atomInBond[1]).getFlag(CDKConstants.ISINRING) + "(" + torsionAtomPosition[m][3] + ")");											    */										/*System.out.println("torsionAtomPosition[" + m + "]: " + bondConnectedBefore[bb].getConnectedAtom(atomInBond[0]).getSymbol() 											+ ", "+ atomInBond[0].getSymbol() + ", " + atomInBond[1].getSymbol() + ", " 											+ bondConnectedAfter[ba].getConnectedAtom(atomInBond[1]).getSymbol());									    */ 										torsionType = "0";										if (bond[b].getProperty("MMFF94 bond type").toString() == "1") {											torsionType = "1";										}										else if ((bond[b].getProperty("MMFF94 bond type").toString() == "0") & 												((bondConnectedBefore[bb].getProperty("MMFF94 bond type").toString() == "1") |												(bondConnectedAfter[ba].getProperty("MMFF94 bond type").toString() == "1"))) {											torsionType = "2";										}										/*System.out.println("torsion " + m + " : " + torsionType + " " + bondConnectedBefore[bb].getConnectedAtom(atomInBond[0]).getAtomTypeName() + "(" + torsionAtomPosition[m][0] + "), " + 										atomInBond[0].getAtomTypeName() + "(" + torsionAtomPosition[m][1] + "), " + atomInBond[1].getAtomTypeName() + "(" + torsionAtomPosition[m][2] + "), " + 										bondConnectedAfter[ba].getConnectedAtom(atomInBond[1]).getAtomTypeName() + "(" + torsionAtomPosition[m][3] + ")");										*/										torsionsData = pc.getTorsionData(torsionType, bondConnectedBefore[bb].getConnectedAtom(atomInBond[0]).getAtomTypeName(), 												atomInBond[0].getAtomTypeName(), atomInBond[1].getAtomTypeName(), bondConnectedAfter[ba].getConnectedAtom(atomInBond[1]).getAtomTypeName());																			//logger.debug("torsionsData " + m + ": " + torsionsData);										v1[m] = ((Double) torsionsData.get(0)).doubleValue();										v2[m] = /*(-1) * */((Double) torsionsData.get(1)).doubleValue();										v3[m] = ((Double) torsionsData.get(2)).doubleValue();									}								}								}						}					}				}			}		}		phi = new double[torsionNumber];		this.moleculeCurrentCoordinates = new GVector(3 * molecule.getAtomCount());		for (int i=0; i<moleculeCurrentCoordinates.getSize(); i++) {			this.moleculeCurrentCoordinates.setElement(i,1E10);		} 		this.changeAtomCoordinates = new boolean[molecule.getAtomCount()];	}	/**	 *  Calculate the actual phi	 *	 *@param  coords3d  Current molecule coordinates.	 */	public void setPhi(GVector coords3d) {		changedCoordinates = 0;		//logger.debug("Setting Phi");		for (int i=0; i < changeAtomCoordinates.length; i++) {			this.changeAtomCoordinates[i] = false;		}		this.moleculeCurrentCoordinates.sub(coords3d);		for (int i = 0; i < this.moleculeCurrentCoordinates.getSize(); i++) {			//logger.debug("moleculeCurrentCoordinates " + i + " = " + this.moleculeCurrentCoordinates.getElement(i));			if (Math.abs(this.moleculeCurrentCoordinates.getElement(i)) > 0) {				changeAtomCoordinates[i/3] = true;				changedCoordinates = changedCoordinates + 1;				//logger.debug("changeAtomCoordinates[" + i/3 + "] = " + changeAtomCoordinates[i/3]);				i = i + (2 - i % 3);			}		}		for (int m = 0; m < torsionNumber; m++) {			if ((changeAtomCoordinates[torsionAtomPosition[m][0]] == true) | 					(changeAtomCoordinates[torsionAtomPosition[m][1]] == true) | 					(changeAtomCoordinates[torsionAtomPosition[m][2]] == true) |					(changeAtomCoordinates[torsionAtomPosition[m][3]] == true))		{							phi[m] = ForceFieldTools.torsionAngleFrom3xNCoordinates(coords3d, torsionAtomPosition[m][0], torsionAtomPosition[m][1], 							torsionAtomPosition[m][2], torsionAtomPosition[m][3]);			} 			//else {System.out.println("phi was no recalculated");}		}		/*if 	(changedCoordinates == changeAtomCoordinates.length) {			for (int m = 0; m < torsionNumber; m++) {				System.out.println("phi[" + m + "] = " + Math.toDegrees(phi[m]));			}