bondstretching.java

化学图形处理软件

package org.openscience.cdk.modeling.forcefield;

import org.openscience.cdk.interfaces.IAtomContainer;
import org.openscience.cdk.interfaces.IBond;
import org.openscience.cdk.modeling.builder3d.MMFF94ParametersCall;

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;


/**
 *  Bond Stretching calculator for the potential energy function. Include function and derivatives.
 *
 *@author     vlabarta
 *@cdk.created    January 27, 2005
 *@cdk.module     forcefield
 */
public class BondStretching {

	String functionShape = " Bond Stretching ";
	
		
	double mmff94SumEB = 0;
	GVector gradientMMFF94SumEB = null;
	GMatrix hessianMMFF94SumEB = null;
	double[] forHessian = null;
	GMatrix order2ndErrorApproximateHessianMMFF94SumEB = null;
	double[] forOrder2ndErrorApproximateHessian = null;

	int bondsNumber;
	int[][] bondAtomPosition = null;

	double[] r0 = null;	// Force field parameters
	double[] k2 = null;
	double[] k3 = null;
	double[] k4 = null;
	double cs = -2;
	
	double[] r = null;	// The actual bond lengths
	double[] deltar = null;	// The difference between actual and reference bond lengths
	
	double[][] dDeltar = null;
	double[][][] ddDeltar = null;

	//private LoggingTool logger;


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


	/**
	 *  Set MMFF94 reference bond lengths r0IJ and the constants k2, k3, k4 for
	 *  each i-j bond in the molecule.
	 *
	 *@param  molecule       The molecule like an AtomContainer object.
	 *@param  parameterSet   MMFF94 parameters set
	 *@exception  Exception  Description of the Exception
	 */
	public void setMMFF94BondStretchingParameters(IAtomContainer molecule, Hashtable parameterSet) throws Exception {

		//logger.debug("setMMFF94BondStretchingParameters");
		
		bondsNumber = molecule.getBondCount();
		//logger.debug("bondsNumber = " + bondsNumber);
		bondAtomPosition = new int[molecule.getBondCount()][];

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

		r0 = new double[molecule.getBondCount()];
		k2 = new double[molecule.getBondCount()];
		k3 = new double[molecule.getBondCount()];
		k4 = new double[molecule.getBondCount()];

        String bondType;
        Iterator bonds = molecule.bonds();
        int i = 0;
        while (bonds.hasNext()) {
            IBond bond = (IBond) bonds.next();

            //atomsInBond = bonds[i].getatoms();

            bondType = bond.getProperty("MMFF94 bond type").toString();
            //logger.debug("bondType " + i + " = " + bondType);

            bondAtomPosition[i] = new int[bond.getAtomCount()];

            for (int j = 0; j < bond.getAtomCount(); j++) {
                bondAtomPosition[i][j] = molecule.getAtomNumber(bond.getAtom(j));
            }

            /*System.out.println("bond " + i + " : " + bondType + " " + bonds[i].getAtom(0).getAtomTypeName() + "(" + bondAtomPosition[i][0] + "), " +
                       bonds[i].getAtom(1).getAtomTypeName() + "(" + bondAtomPosition[i][1] + ")");
               */
            bondData = pc.getBondData(bondType, bond.getAtom(0).getAtomTypeName(), bond.getAtom(1).getAtomTypeName());
            //logger.debug("bondData : " + bondData);
            r0[i] = ((Double) bondData.get(0)).doubleValue();
            k2[i] = ((Double) bondData.get(1)).doubleValue();
            k3[i] = ((Double) bondData.get(2)).doubleValue();
            k4[i] = ((Double) bondData.get(3)).doubleValue();

            i++;
        }

        r = new double[molecule.getBondCount()];
		deltar = new double[molecule.getBondCount()];
		

	}


	/**
	 *  Calculate the actual bond distance rij and the difference with the reference bond distances.
	 *
	 *@param  coord3d  Current molecule coordinates.
	 */
	public void calculateDeltar(GVector coord3d) {

		//logger.debug("deltar.length = " + deltar.length);
		for (int i = 0; i < bondAtomPosition.length; i++) {
			r[i] = ForceFieldTools.distanceBetweenTwoAtomsFrom3xNCoordinates(coord3d, bondAtomPosition[i][0], bondAtomPosition[i][1]);
			deltar[i] = r[i] - r0[i];
			//if (deltar[i] > 0) {
			//	deltar[i] = (-1) * deltar[i];
			//}
		}
	}


	/**
	 *  Evaluate the MMFF94 bond stretching term for the given atoms cartesian coordinates.
	 *
	 *@param  coord3d  Current molecule coordinates.
	 *@return        bond stretching value
	 */
	public double functionMMFF94SumEB(GVector coord3d) {
		/*for (int i=0; i < bondsNumber; i++) {
			logger.debug("before:	deltar[" + i + "] = " + deltar[i]);
		}*/
		
		calculateDeltar(coord3d);

		/*for (int i=0; i < bondsNumber; i++) {
			logger.debug("after:	deltar[" + i + "] = " + deltar[i]);
		}*/
		
		mmff94SumEB = 0;

		for (int i = 0; i < bondsNumber; i++) {
			mmff94SumEB = mmff94SumEB + k2[i] * Math.pow(deltar[i],2) 
							+ k3[i] * Math.pow(deltar[i],3) + k4[i] * Math.pow(deltar[i],4);
		}

		//logger.debug("mmff94SumEB = " + mmff94SumEB);
		
		return mmff94SumEB;
	}


	/**
	 *  Set the bond lengths first derivative respect to the cartesian
	 *  coordinates of the atoms.
	 *
	 *@param  coord3d           Current molecule coordinates.
	 *@param  deltar           Difference between the current bonds and the reference bonds.
	 *@param  bondAtomPosition  Position of the bending atoms in the atoms coordinates (0:N, N: atoms number).
	 */
	public void setBondLengthsFirstDerivative(GVector coord3d, double[] deltar, int[][] bondAtomPosition) {

		dDeltar = new double[coord3d.getSize()][];

		Double forAtomNumber = null;
		int atomNumber = 0;
		int coordinate;
		for (int i = 0; i < dDeltar.length; i++) {

			dDeltar[i] = new double[deltar.length];

			forAtomNumber = new Double(i / 3);
			coordinate = i % 3;
			//logger.debug("coordinate = " + coordinate);

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

			for (int j = 0; j < deltar.length; j++) {

				if ((bondAtomPosition[j][0] == atomNumber) | (bondAtomPosition[j][1] == atomNumber)) {
					switch (coordinate) {
									//x-coordinate
									case 0:
										dDeltar[i][j] = (coord3d.getElement(3 * bondAtomPosition[j][0]) - coord3d.getElement(3 * bondAtomPosition[j][1]))
												 / Math.sqrt(Math.pow(coord3d.getElement(3 * bondAtomPosition[j][0]) - coord3d.getElement(3 * bondAtomPosition[j][1]), 2) + Math.pow(coord3d.getElement(3 * bondAtomPosition[j][0] + 1) - coord3d.getElement(3 * bondAtomPosition[j][1] + 1), 2) + Math.pow(coord3d.getElement(3 * bondAtomPosition[j][0] + 2) - coord3d.getElement(3 * bondAtomPosition[j][1] + 2), 2));
										break;
									//y-coordinate
									case 1:
										dDeltar[i][j] = (coord3d.getElement(3 * bondAtomPosition[j][0] + 1) - coord3d.getElement(3 * bondAtomPosition[j][1] + 1))
												 / Math.sqrt(Math.pow(coord3d.getElement(3 * bondAtomPosition[j][0]) - coord3d.getElement(3 * bondAtomPosition[j][1]), 2) + Math.pow(coord3d.getElement(3 * bondAtomPosition[j][0] + 1) - coord3d.getElement(3 * bondAtomPosition[j][1] + 1), 2) + Math.pow(coord3d.getElement(3 * bondAtomPosition[j][0] + 2) - coord3d.getElement(3 * bondAtomPosition[j][1] + 2), 2));
										break;
									//z-coordinate
									case 2:
										dDeltar[i][j] = (coord3d.getElement(3 * bondAtomPosition[j][0] + 2) - coord3d.getElement(3 * bondAtomPosition[j][1] + 2))
												 / Math.sqrt(Math.pow(coord3d.getElement(3 * bondAtomPosition[j][0]) - coord3d.getElement(3 * bondAtomPosition[j][1]), 2) + Math.pow(coord3d.getElement(3 * bondAtomPosition[j][0] + 1) - coord3d.getElement(3 * bondAtomPosition[j][1] + 1), 2) + Math.pow(coord3d.getElement(3 * bondAtomPosition[j][0] + 2) - coord3d.getElement(3 * bondAtomPosition[j][1] + 2), 2));
										break;
					}
					if (bondAtomPosition[j][1] == atomNumber) {
						dDeltar[i][j] = (-1) * dDeltar[i][j];
					}
				} else {
					dDeltar[i][j] = 0;
				}
				//logger.debug("bond " + j + " : " + "dDeltar[" + i + "][" + j + "] = " + dDeltar[i][j]);
			}
		}
	}


	/**
	 *  Get the bond lengths first derivative respect to the cartesian coordinates of the
	 *  atoms.
	 *
	 *@return    Delta bond lengths derivative value [dimension(3xN)] [bonds Number]
	 *      
	 */
	public double[][] getBondLengthsFirstDerivative() {
		return dDeltar;
	}


	/**
	 *  Evaluate the first order partial derivative for the bond stretching given the atoms coordinates
	 *
	 *@param  coord3d  Current molecule coordinates.
	 */
	public void setGradientMMFF94SumEB(GVector coord3d) {
		
		gradientMMFF94SumEB = new GVector(coord3d.getSize());
		
		calculateDeltar(coord3d);
		setBondLengthsFirstDerivative(coord3d, deltar, bondAtomPosition);
		
		double sumGradientEB;
		for (int i = 0; i < gradientMMFF94SumEB.getSize(); i++) {
			
			sumGradientEB = 0;
			for (int j = 0; j < bondsNumber; j++) {