stretchbendinteractions.java

化学图形处理软件

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.interfaces.IBond;
import org.openscience.cdk.interfaces.IAtom;
import org.openscience.cdk.tools.manipulator.AtomContainerManipulator;
import org.openscience.cdk.modeling.builder3d.MMFF94ParametersCall;
import org.openscience.cdk.qsar.IAtomicDescriptor;
import org.openscience.cdk.qsar.descriptors.atomic.PeriodicTablePositionDescriptor;
import org.openscience.cdk.qsar.result.IntegerResult;
import org.openscience.cdk.tools.LoggingTool;

/**
 *  Stretch-Bend Interaction calculator for the potential energy function.
 *  Include function and derivatives.
 *
 *@author         vlabarta
 *@cdk.created    2005-02-15
 *@cdk.module     forcefield
 */
public class StretchBendInteractions {

	String functionShape = " Stretch-Bend Interactions ";

	double mmff94SumEBA = 0;
	GVector gradientMMFF94SumEBA = null;
	GVector order2ndErrorApproximateGradientMMFF94SumEBA = null;
	GVector order5thErrorApproximateGradientMMFF94SumEBA = null;
	GMatrix hessianMMFF94SumEBA = null;
	GVector currentCoordinates = null;
	GVector gradientCurrentCoordinates = null;

	double[][] dDeltarij = null;
	double[][] dDeltarkj = null;
	double[][] dDeltav = null;

	int[][] bondijAtomPosition = null;
	int[][] bondkjAtomPosition = null;
	double[] r0IJ = null;
	double[] r0KJ = null;
	double[] kbaIJK = null;
	double[] kbaKJI = null;
	double[] rij = null;
	double[] rkj = null;
	double[] deltarij = null;
	double[] deltarkj = null;

	BondStretching bs = new BondStretching();
	AngleBending ab = new AngleBending();
	private LoggingTool logger;

	GVector moleculeCurrentCoordinates = null;
	boolean[] changeAtomCoordinates = null;
	int changedCoordinates;


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


	/**
	 *  Set MMFF94 reference bond lengths r0IJ and r0JK and stretch-bend
	 *  interaction constants kbaIJK and kbaKJI 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 setMMFF94StretchBendParameters(IAtomContainer molecule, Hashtable parameterSet, boolean angleBendingFlag) throws Exception {

		//logger.debug("setMMFF94StretchBendParameters");
		
		ab.setMMFF94AngleBendingParameters(molecule, parameterSet, angleBendingFlag);

		IAtom[] atomConnected = null;

		Vector stretchBendInteractionsData = null;
		Vector bondData = null;
		MMFF94ParametersCall pc = new MMFF94ParametersCall();
		pc.initialize(parameterSet);

		bondijAtomPosition = new int[ab.angleNumber][];
		bondkjAtomPosition = new int[ab.angleNumber][];
		r0IJ = new double[ab.angleNumber];
		r0KJ = new double[ab.angleNumber];
		kbaIJK = new double[ab.angleNumber];
		kbaKJI = new double[ab.angleNumber];

		String strbndType;
		String angleType;
		String bondIJType;
		String bondKJType;

		IBond bondIJ = null;
		IBond bondKJ = null;
		
		IAtomicDescriptor descriptor  = new PeriodicTablePositionDescriptor();
		int iR = 0;
		int jR = 0;
		int kR = 0;

		
		int l = -1;
		for (int j = 0; j < molecule.getAtomCount(); j++) {
			
			atomConnected = AtomContainerManipulator.getAtomArray(molecule.getConnectedAtomsList(molecule.getAtom(j)));
			
			if (atomConnected.length > 1) {
				
				for (int i = 0; i < atomConnected.length; i++) {
					
					for (int k = i + 1; k < atomConnected.length; k++) {
						
						l += 1;
						
						bondIJ = molecule.getBond(atomConnected[i], molecule.getAtom(j));
						bondIJType = bondIJ.getProperty("MMFF94 bond type").toString();
						
						bondKJ = molecule.getBond(atomConnected[k], molecule.getAtom(j));
						bondKJType = bondKJ.getProperty("MMFF94 bond type").toString();
						
						angleType = "0";
						if ((bondIJType == "1") | (bondKJType == "1")) {
							angleType = "1";
						}  
						if ((bondIJType == "1") & (bondKJType == "1")) {
							angleType = "2";
						}  
						
						//logger.debug("bondIJType = " + bondIJType + ", bondKJType = " + bondKJType + ", angleType = " + angleType);
						
						strbndType = "0";
						if ((angleType == "0") & (bondIJType == "0") & (bondKJType == "0")) {strbndType = "0";}
						else if ((angleType == "1") & (bondIJType == "1") & (bondKJType == "0")) {strbndType = "1";}
						else if ((angleType == "1") & (bondIJType == "0") & (bondKJType == "1")) {strbndType = "2";}
						else if ((angleType == "2") & (bondIJType == "1") & (bondKJType == "1")) {strbndType = "3";}
						else if ((angleType == "4") & (bondIJType == "0") & (bondKJType == "0")) {strbndType = "4";}
						else if ((angleType == "3") & (bondIJType == "0") & (bondKJType == "0")) {strbndType = "5";}
						else if ((angleType == "5") & (bondIJType == "1") & (bondKJType == "0")) {strbndType = "6";}
						else if ((angleType == "5") & (bondIJType == "0") & (bondKJType == "1")) {strbndType = "7";}
						else if ((angleType == "6") & (bondIJType == "1") & (bondKJType == "1")) {strbndType = "8";}
						else if ((angleType == "7") & (bondIJType == "1") & (bondKJType == "0")) {strbndType = "9";}
						else if ((angleType == "7") & (bondIJType == "0") & (bondKJType == "1")) {strbndType = "10";}
						else if ((angleType == "8") & (bondIJType == "1") & (bondKJType == "1")) {strbndType = "11";}
						
						//logger.debug("strbnd: " + strbndType + ", " + atomConnected[i].getAtomTypeName() + "(" + molecule.getAtomNumber(atomConnected[i]) + "), " + molecule.getAtom(j).getAtomTypeName() + "(" + molecule.getAtomNumber(molecule.getAtom(j)) + "), " + ((IAtom)atomConnected.get(k)).getAtomTypeName() + "(" + molecule.getAtomNumber((IAtom)atomConnected.get(k)) + ")");
						stretchBendInteractionsData = pc.getBondAngleInteractionData(strbndType, atomConnected[i].getAtomTypeName(), molecule.getAtom(j).getAtomTypeName(), atomConnected[k].getAtomTypeName());
						
						if (stretchBendInteractionsData == null) {
							if (angleType == "1") {
								if (strbndType == "1") {strbndType = "2";}
								else {strbndType = "1";}
								//logger.debug("strbnd: " + strbndType + ", " + ((IAtom)atomConnected.get(i)).getAtomTypeName() + "(" + molecule.getAtomNumber((IAtom)atomConnected.get(i)) + "), " + molecule.getAtom(j).getAtomTypeName() + "(" + molecule.getAtomNumber(molecule.getAtom(j)) + "), " + ((IAtom)atomConnected.get(k)).getAtomTypeName() + "(" + molecule.getAtomNumber((IAtom)atomConnected.get(k)) + ")");
								stretchBendInteractionsData = pc.getBondAngleInteractionData(strbndType, atomConnected[i].getAtomTypeName(), molecule.getAtom(j).getAtomTypeName(), atomConnected[k].getAtomTypeName());
							}
						}
						
						if (stretchBendInteractionsData == null) {
							iR = ((IntegerResult)descriptor.calculate(atomConnected[i],molecule).getValue()).intValue();
							jR = ((IntegerResult)descriptor.calculate(molecule.getAtom(j),molecule).getValue()).intValue();
							kR = ((IntegerResult)descriptor.calculate(atomConnected[k],molecule).getValue()).intValue();
							stretchBendInteractionsData = pc.getDefaultStretchBendData(iR, jR, kR);
						} 

						//logger.debug("stretchBendInteractionsData : " + stretchBendInteractionsData);
						kbaIJK[l] = ((Double) stretchBendInteractionsData.get(0)).doubleValue();
						kbaKJI[l] = ((Double) stretchBendInteractionsData.get(1)).doubleValue();

						//logger.debug("kbaIJK[" + l + "] = " + kbaIJK[l]);
						//logger.debug("kbaKJI[" + l + "] = " + kbaKJI[l]);

						
						bondData = pc.getBondData(bondIJType, atomConnected[i].getAtomTypeName(), molecule.getAtom(j).getAtomTypeName());
						r0IJ[l] = ((Double) bondData.get(0)).doubleValue();
						bondData = pc.getBondData(bondKJType, atomConnected[k].getAtomTypeName(), molecule.getAtom(j).getAtomTypeName());
						r0KJ[l] = ((Double) bondData.get(0)).doubleValue();
						
						bondijAtomPosition[l] = new int[2];
						bondijAtomPosition[l][0] = molecule.getAtomNumber(atomConnected[i]);
						bondijAtomPosition[l][1] = j;
						
						bondkjAtomPosition[l] = new int[2];
						bondkjAtomPosition[l][0] = molecule.getAtomNumber(atomConnected[k]);
						bondkjAtomPosition[l][1] = j;
					}
				}
			}
		}
		rij = new double[ab.angleNumber];
		rkj = new double[ab.angleNumber];
		deltarij = new double[ab.angleNumber];
		deltarkj = new double[ab.angleNumber];
		currentCoordinates = new GVector(3 * molecule.getAtomCount());
		gradientCurrentCoordinates = new GVector(3 * molecule.getAtomCount());
		gradientMMFF94SumEBA = new GVector(3 * molecule.getAtomCount());
		dDeltarij = new double[3 * molecule.getAtomCount()][];
		dDeltarkj = new double[3 * molecule.getAtomCount()][];
		dDeltav = new double[3 * molecule.getAtomCount()][];
		hessianMMFF94SumEBA = new GMatrix(3 * molecule.getAtomCount(), 3 * molecule.getAtomCount());
		for (int i = 0; i < 3 * molecule.getAtomCount(); i++) {
			dDeltarij[i] = new double[ab.angleNumber];
			dDeltarkj[i] = new double[ab.angleNumber];
			dDeltav[i] = new double[ab.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()];

	}


	/**
	 *  Calculate the current bond distances rij and rkj for each angle j, and the
	 *  difference with the reference bonds.
	 *
	 *@param  coords3d  Current molecule coordinates.
	 */
	public void setDeltarijAndDeltarkj(GVector coords3d) {

		changedCoordinates = 0;
		//logger.debug("Setting Deltarij and Deltarkj");
		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 i = 0; i < ab.angleNumber; i++) {
			if ((changeAtomCoordinates[ab.angleAtomPosition[i][0]] == true) | 
					(changeAtomCoordinates[ab.angleAtomPosition[i][1]] == true))		{
			
				rij[i] = ForceFieldTools.distanceBetweenTwoAtomsFrom3xNCoordinates(coords3d, ab.angleAtomPosition[i][1], ab.angleAtomPosition[i][0]);
				deltarij[i] = rij[i] - r0IJ[i];
				//logger.debug("deltarij[" + i + "] = " + deltarij[i]);
			}
			//else {System.out.println("deltarij[" + i + "] was no recalculated");}
			if ((changeAtomCoordinates[ab.angleAtomPosition[i][1]] == true) | 
					(changeAtomCoordinates[ab.angleAtomPosition[i][2]] == true))		{
			
				rkj[i] = ForceFieldTools.distanceBetweenTwoAtomsFrom3xNCoordinates(coords3d, ab.angleAtomPosition[i][1], ab.angleAtomPosition[i][2]);
				deltarkj[i] = rkj[i] - r0KJ[i];
				//logger.debug("deltarkj[" + i + "] = " + deltarkj[i]);