anglebending.java

化学图形处理软件

package org.openscience.cdk.modeling.forcefield;

import java.util.Hashtable;
import java.util.Vector;

import javax.vecmath.GMatrix;
import javax.vecmath.GVector;
import javax.vecmath.Point3d;
import javax.vecmath.Vector3d;

import org.openscience.cdk.interfaces.IAtomContainer;
import org.openscience.cdk.interfaces.IAtom;
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.LoggingTool;


/**
 *  Angle bending calculator for the potential energy function. Include function and derivatives.
 *
 *@author     vlabarta
 *@cdk.created    February 8, 2005
 *@cdk.module     forcefield
 */
public class AngleBending {

	String functionShape = " Angle bending ";

	double mmff94SumEA = 0;
	GVector gradientMMFF94SumEA = null;
	GVector order2ndErrorApproximateGradientMMFF94SumEA = null;
	GVector order5thErrorApproximateGradientMMFF94SumEA = null;
	GMatrix hessianMMFF94SumEA = null;
	double[] forHessian = null;
	GMatrix order2ndErrorApproximateHessianMMFF94SumEA = null;
	double[] forOrder2ndErrorApproximateHessian = null;

	double[][] dDeltav = null;
	double[][] angleBendingOrder2ndErrorApproximateGradient = null;
	double[][][] ddDeltav = null;
	double[][][] angleBendingOrder2ndErrorApproximateHessian = null;

	int angleNumber = 0;
	int[][] angleAtomPosition = null;

	double[] v0 = null;
	double[] k2 = null;
	double[] k3 = null;
	double[] k4 = null;
	double cb = -0.007;
	double[] currentCoordinates_v = null;
	double[] currentCoordinates_deltav = null;
	double[] v = null;
	double[] deltav = null;
	
	boolean angleBending;
	
	GVector moleculeCurrentCoordinates = null;
	boolean[] changeAtomCoordinates = null;
	int changedCoordinates;



	/**
	 *  Constructor for the AngleBending object
	 */
	public AngleBending() {
	        //logger = new LoggingTool(this);
	}

	public void setAngleBendingFlag(boolean flag){
		angleBending=flag;
	}
	
	/**
	 *  Set MMFF94 reference angle v0IJK and the constants k2, k3, k4 for each
	 *  i-j-k angle in the molecule.
	 *
	 *@param  molecule       The molecule like an AtomContainer object.
	 *@param  parameterSet   MMFF94 parameters set
	 *@exception  Exception  Description of the Exception
	 */
	public void setMMFF94AngleBendingParameters(IAtomContainer molecule, Hashtable parameterSet, boolean angleBendingFlag ) throws Exception {

		//logger.debug("setMMFF94AngleBendingParameters");		
		
		IAtom[] atomConnected = null;
		angleBending=angleBendingFlag;
		for (int i = 0; i < molecule.getAtomCount(); i++) {
			atomConnected = AtomContainerManipulator.getAtomArray(molecule.getConnectedAtomsList(molecule.getAtom(i)));
			if (atomConnected.length > 1) {
				for (int j = 0; j < atomConnected.length; j++) {
					for (int k = j+1; k < atomConnected.length; k++) {
						angleNumber += 1;
					}
				}
			}
		}
		//logger.debug("angleNumber = " + angleNumber);

		Vector angleData = null;
		MMFF94ParametersCall pc = new MMFF94ParametersCall();
		pc.initialize(parameterSet);
		
		v0 = new double[angleNumber];
		k2 = new double[angleNumber];
		k3 = new double[angleNumber];
		k4 = new double[angleNumber];
		
		angleAtomPosition = new int[angleNumber][];

		String angleType;
		IBond bondIJ = null;
		IBond bondKJ = null;
		String bondIJType;
		String bondKJType;
		int l = -1;
		for (int i = 0; i < molecule.getAtomCount(); i++) {
			atomConnected = AtomContainerManipulator.getAtomArray(molecule.getConnectedAtomsList(molecule.getAtom(i)));
			if (atomConnected.length > 1) {
				for (int j = 0; j < atomConnected.length; j++) {
					for (int k = j+1; k < atomConnected.length; k++) {
						l += 1;
						bondIJ = molecule.getBond(atomConnected[j], molecule.getAtom(i));
						bondIJType = bondIJ.getProperty("MMFF94 bond type").toString();
						//logger.debug("bondIJType = " + bondIJType);

						bondKJ = molecule.getBond(atomConnected[k], molecule.getAtom(i));
						bondKJType = bondKJ.getProperty("MMFF94 bond type").toString();
						//logger.debug("bondKJType = " + bondKJType);
						
						angleType = "0";
						if ((bondIJType == "1") | (bondKJType == "1")) {
							angleType = "1";
						}  
						if ((bondIJType == "1") & (bondKJType == "1")) {
							angleType = "2";
						}  
						
						//logger.debug(angleType + ", " + atomConnected[j].getAtomTypeName() + ", " + molecule.getAtom(i).getAtomTypeName() + ", " + atomConnected[k].getAtomTypeName());
						angleData = pc.getAngleData(angleType, atomConnected[j].getAtomTypeName(), molecule.getAtom(i).getAtomTypeName(), atomConnected[k].getAtomTypeName());
						//logger.debug("angleData : " + angleData);
						v0[l] = ((Double) angleData.get(0)).doubleValue();
						k2[l] = ((Double) angleData.get(1)).doubleValue();
						k3[l] = ((Double) angleData.get(2)).doubleValue();
						//k4[l] = ((Double) angleData.get(3)).doubleValue();

						//logger.debug("v0[" + l + "] = " + v0[l]);
						//logger.debug("k2[" + l + "] = " + k2[l]);
						//logger.debug("k3[" + l + "] = " + k3[l]);
						//logger.debug("k4[" + l + "] = " + k4[l]);
						
						angleAtomPosition[l] = new int[3];
						angleAtomPosition[l][0] = molecule.getAtomNumber(atomConnected[j]);
						angleAtomPosition[l][1] = i;
						angleAtomPosition[l][2] = molecule.getAtomNumber(atomConnected[k]);

					}
				}
			}
		}

		v = new double[angleNumber];
		deltav = new double[angleNumber];

		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()];
		for (int i=0; i < molecule.getAtomCount(); i++) {
			this.changeAtomCoordinates[i] = false;
		}

	}


	/**
	 *  Calculate the actual bond angles vijk and the difference with the reference angles.
	 *
	 *@param  coord3d  Current molecule coordinates.
	 */
	public void setDeltav(GVector coord3d) {
		changedCoordinates = 0;
		//logger.debug("Setting Deltav");
		for (int i=0; i < changeAtomCoordinates.length; i++) {
			this.changeAtomCoordinates[i] = false;
		}
		this.moleculeCurrentCoordinates.sub(coord3d);
		for (int i = 0; i < this.moleculeCurrentCoordinates.getSize(); i++) {
			//logger.debug("this.moleculeCurrentCoordinates.getElement(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);
			}
		}
		//logger.debug("currentCoordinates_deltav.length = " + currentCoordinates_deltav.length);

		for (int i = 0; i < angleNumber; i++) {
			if ((changeAtomCoordinates[angleAtomPosition[i][0]] == true) | 
				(changeAtomCoordinates[angleAtomPosition[i][1]] == true) | 
				(changeAtomCoordinates[angleAtomPosition[i][2]] == true))		{
				
				v[i] = ForceFieldTools.angleBetweenTwoBondsFrom3xNCoordinates(coord3d,angleAtomPosition[i][0],angleAtomPosition[i][1],angleAtomPosition[i][2]);
				//logger.debug("currentCoordinates_v[" + i + "] = " + currentCoordinates_v[i]);
				//logger.debug("v0[" + i + "] = " + v0[i]);
				deltav[i] = v[i] - v0[i];
				if (deltav[i] > 0 & angleBending) {
					deltav[i]= (-1) * deltav[i]; 
				}else if (deltav[i] < 0 & !angleBending){
					deltav[i]= (-1) * deltav[i]; 
				}
				/*if (Math.abs(currentCoordinates_deltav[i]) < 0.05) {
				 logger.debug("currentCoordinates_deltav[" + i + "]= " + currentCoordinates_deltav[i]);
				 }*/
			}
			//else {System.out.println("v[" + i + "] remain the same");}
		/*if 	(changedCoordinates == changeAtomCoordinates.length) {
		 		for (int m = 0; m < torsionNumber; m++) {
			 		System.out.println("phi[" + m + "] = " + Math.toDegrees(phi[m]));
			 	}
		}
		*/
		}
		moleculeCurrentCoordinates.set(coord3d);
	}



	/**
	 *  Get the current difference between the bond angles vijk and the reference angles.
	 *
	 *@return  Difference between the current bond angles vijk and the reference angles.
	 */
	public double[] getDeltav() {
		return deltav;	
	}


	/**
	 *  Evaluate the MMFF94 angle bending term for the given atoms coordinates
	 *
	 *@param  coord3d  Current molecule coordinates.
	 *@return        MMFF94 angle bending term value
	 */
	public double functionMMFF94SumEA(GVector coord3d) {
		//this.angleBending=true;
		this.setDeltav(coord3d);
		mmff94SumEA = 0;
		for (int i = 0; i < angleNumber; i++) {
			//logger.debug("k2[" + i + "]= " + k2[i]);
			//logger.debug("k3[" + i + "]= " + k3[i]);
			//logger.debug("currentCoordinates_deltav[" + i + "]= " + currentCoordinates_deltav[i]);
			//logger.debug("For Angle " + i + " : " + k2[i] * Math.pow(currentCoordinates_deltav[i],2) + k3[i] * Math.pow(currentCoordinates_deltav[i],3));

			mmff94SumEA = mmff94SumEA + k2[i] * Math.pow(deltav[i],2) 
						+ k3[i] * Math.pow(deltav[i],3);
						
			//mmff94SumEA = Math.abs(mmff94SumEA);
			
			//logger.debug("mmff94SumEA = " + mmff94SumEA);
		}
		//mmff94SumEA = Math.abs(mmff94SumEA);
		//logger.debug("mmff94SumEA = " + mmff94SumEA);
		return mmff94SumEA;
	}


	/**
	 *  Calculate the angle bending first derivative respect to the cartesian coordinates of the atoms.
	 *
	 *@param  coord3d  Current molecule coordinates.
	 */
	public void setAngleBendingFirstDerivative(GVector coord3d) {
		
		dDeltav = new double[coord3d.getSize()][];
		
		Double forAtomNumber = null;
		int atomNumber = 0;
		int coordinate;

		for (int m = 0; m < dDeltav.length; m++) {
			
			dDeltav[m] = new double[angleNumber];
			
			forAtomNumber = new Double(m/3);
			coordinate = m % 3;
			//logger.debug("coordinate = " + coordinate);

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

			for (int l = 0; l < angleNumber; l++) {
				
				if ((angleAtomPosition[l][0] == atomNumber) | (angleAtomPosition[l][1] == atomNumber) | (angleAtomPosition[l][2] == atomNumber)) {

					Point3d xi = new Point3d(coord3d.getElement(3 * angleAtomPosition[l][0]), coord3d.getElement(3 * angleAtomPosition[l][0] + 1),coord3d.getElement( 3 * angleAtomPosition[l][0] + 2));
					//logger.debug("xi = " + xi);
					Point3d xj = new Point3d(coord3d.getElement(3 * angleAtomPosition[l][1]), coord3d.getElement(3 * angleAtomPosition[l][1] + 1),coord3d.getElement( 3 * angleAtomPosition[l][1] + 2));
					//logger.debug("xj = " + xj);
					Point3d xk = new Point3d(coord3d.getElement(3 * angleAtomPosition[l][2]), coord3d.getElement(3 * angleAtomPosition[l][2] + 1),coord3d.getElement( 3 * angleAtomPosition[l][2] + 2));
					//logger.debug("xk = " + xk);
				
					Vector3d xij = new Vector3d();
					xij.sub(xi,xj);
					//logger.debug("xij = " + xij);
					Vector3d xkj = new Vector3d();
					xkj.sub(xk,xj);
					//logger.debug("xkj = " + xkj);
					
					double rji = xj.distance(xi);
					//logger.debug("rji = " + rji);
					double rjk = xj.distance(xk);
					//logger.debug("rji = " + rjk);

					dDeltav[m][l] = (-1/Math.sqrt(1-Math.pow(xij.dot(xkj)/(rji * rjk),2))) * (1/(Math.pow(rji,2) * Math.pow(rjk,2))); 

					if (angleAtomPosition[l][0] == atomNumber) {

						switch (coordinate) {
							case 0: dDeltav[m][l] = dDeltav[m][l] * ((xk.x-xj.x) * rji * rjk - (xij.dot(xkj)) * rjk * (-(xj.x-xi.x)/rji));
								break;
							case 1:	dDeltav[m][l] = dDeltav[m][l] * ((xk.y-xj.y) * rji * rjk - (xij.dot(xkj)) * rjk * (-(xj.y-xi.y)/rji));
								break;
							case 2: dDeltav[m][l] = dDeltav[m][l] * ((xk.z-xj.z) * rji * rjk - (xij.dot(xkj)) * rjk * (-(xj.z-xi.z)/rji));
								break;
						}
					}
					if (angleAtomPosition[l][1] == atomNumber) {

						switch (coordinate) {
							case 0: dDeltav[m][l] = dDeltav[m][l] * ((2 * xj.x - xk.x - xi.x) * rji * rjk - (xij.dot(xkj)) * (((xj.x-xi.x)/rji) * rjk + ((xj.x-xk.x)/rjk) * rji));
								break;
							case 1:	dDeltav[m][l] = dDeltav[m][l] * ((2 * xj.y - xk.y - xi.y) * rji * rjk - (xij.dot(xkj)) * (((xj.y-xi.y)/rji) * rjk + ((xj.y-xk.y)/rjk) * rji));
								break;
							case 2: dDeltav[m][l] = dDeltav[m][l] * ((2 * xj.z - xk.z - xi.z) * rji * rjk - (xij.dot(xkj)) * (((xj.z-xi.z)/rji) * rjk + ((xj.z-xk.z)/rjk) * rji));
								break;
						}
					}

					if (angleAtomPosition[l][2] == atomNumber) {

						switch (coordinate) {
							case 0: dDeltav[m][l] = dDeltav[m][l] * ((xi.x-xj.x) * rji * rjk - (xij.dot(xkj)) * rji * (-(xj.x-xk.x)/rjk));
								break;
							case 1:	dDeltav[m][l] = dDeltav[m][l] * ((xi.y-xj.y) * rji * rjk - (xij.dot(xkj)) * rji * (-(xj.y-xk.y)/rjk));
								break;
							case 2: dDeltav[m][l] = dDeltav[m][l] * ((xi.z-xj.z) * rji * rjk - (xij.dot(xkj)) * rji * (-(xj.z-xk.z)/rjk));
								break;
						}
					}

				}
				else {
					dDeltav[m][l] = 0;
				}
				//logger.debug("dDeltav[" + m + "][" + l + "] = " + dDeltav[m][l]);
			}
		}
		
	}


	/**
	 *  Get the bond lengths first derivative respect to the cartesian coordinates of the atoms.
	 *
	 *@return        Delta angle bending first derivative value [dimension(3xN)] [angles Number]
	 */
	public double[][] getAngleBendingFirstDerivative() {
		return dDeltav;