Menu
▾
▴
[Cdk-commits] CVS: cdk/src/org/openscience/cdk/modeling/forcefield AngleBending.java,1.12,1.13 BondStretching.java,1.16,1.17 ConjugateGradientMethod.java,1.12,1.13 GeometricMinimizer.java,1.21,1.22 LineSearch.java,1.9,1.10 MMFF94EnergyFunction.java,1.6,1.7 NewtonRaphsonMethod.java,1.5,1.6 PotentialFunction.java,1.9,1.10 SteepestDescentsMethod.java,1.9,1.10 StretchBendInteractions.java,1.11,1.12 Torsions.java,1.7,1.8 VanDerWaalsInteractions.java,1.11,1.12
Update of /cvsroot/cdk/cdk/src/org/openscience/cdk/modeling/forcefield In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv13478/src/org/openscience/cdk/modeling/forcefield Modified Files: AngleBending.java BondStretching.java ConjugateGradientMethod.java GeometricMinimizer.java LineSearch.java MMFF94EnergyFunction.java NewtonRaphsonMethod.java PotentialFunction.java SteepestDescentsMethod.java StretchBendInteractions.java Torsions.java VanDerWaalsInteractions.java Log Message: forcefield: second derivative for all energy term completed, stretch-bend interactions energy term completed Index: AngleBending.java =================================================================== RCS file: /cvsroot/cdk/cdk/src/org/openscience/cdk/modeling/forcefield/AngleBending.java,v retrieving revision 1.12 retrieving revision 1.13 diff -u -r1.12 -r1.13 --- AngleBending.java 4 May 2005 14:18:04 -0000 1.12 +++ AngleBending.java 13 Jun 2005 13:17:16 -0000 1.13 @@ -25,13 +25,18 @@ String functionShape = " Angle bending "; double mmff94SumEA = 0; - GVector gradientMMFF94SumEA = new GVector(3); - GVector order2ndApproximateGradientMMFF94SumEA = new GVector(3); - GVector order5ApproximateGradientMMFF94SumEA = new GVector(3); - GMatrix hessianMMFF94SumEA = new GMatrix(3,3); + 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; @@ -122,7 +127,6 @@ v = new double[angleNumber]; deltav = new double[angleNumber]; - } @@ -131,7 +135,7 @@ * *@param coord3d Current molecule coordinates. */ - public void calculateDeltav(GVector coord3d) { + public void setDeltav(GVector coord3d) { //logger.debug("deltav.length = " + deltav.length); for (int i = 0; i < angleNumber; i++) { @@ -150,13 +154,23 @@ /** + * 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) { - calculateDeltav(coord3d); + setDeltav(coord3d); mmff94SumEA = 0; for (int i = 0; i < angleNumber; i++) { //logger.debug("k2[" + i + "]= " + k2[i]); @@ -283,14 +297,58 @@ /** + * Set a 2nd order approximation of the gradient, for the angle bending, given the atoms + * coordinates + * + *@param coord3d Current molecule coordinates. + */ + public void setAngleBending2ndOrderErrorApproximateGradient(GVector coord3d) { + angleBendingOrder2ndErrorApproximateGradient = new double[coord3d.getSize()][]; + double sigma = Math.pow(0.000000000000001,0.33); + GVector xplusSigma = new GVector(coord3d.getSize()); + GVector xminusSigma = new GVector(coord3d.getSize()); + + double[] deltavInXplusSigma = null; + double[] deltavInXminusSigma = null; + for (int m = 0; m < angleBendingOrder2ndErrorApproximateGradient.length; m++) { + angleBendingOrder2ndErrorApproximateGradient[m] = new double[angleNumber]; + xplusSigma.set(coord3d); + xplusSigma.setElement(m,coord3d.getElement(m) + sigma); + setDeltav(xplusSigma); + deltavInXplusSigma = getDeltav(); + xminusSigma.set(coord3d); + xminusSigma.setElement(m,coord3d.getElement(m) - sigma); + setDeltav(xminusSigma); + deltavInXminusSigma = getDeltav(); + for (int l=0; l < angleNumber; l++) { + angleBendingOrder2ndErrorApproximateGradient[m][l] = (deltavInXplusSigma[l] - deltavInXminusSigma[l]) / (2 * sigma); + } + } + + //logger.debug("order2ndErrorApproximateGradientMMFF94SumEA : " + order2ndErrorApproximateGradientMMFF94SumEA); + } + + + /** + * Get the angle bending 2nd order error approximate gradient. + * + * + *@return Angle bending 2nd order error approximate gradient value. + */ + public double[][] getAngleBending2ndOrderErrorApproximateGradient() { + return angleBendingOrder2ndErrorApproximateGradient; + } + + + /** * Evaluate the gradient of the angle bending term for a given atoms * coordinates * *@param coord3d Current molecule coordinates. */ public void setGradientMMFF94SumEA(GVector coord3d) { - gradientMMFF94SumEA.setSize(coord3d.getSize()); - calculateDeltav(coord3d); + gradientMMFF94SumEA = new GVector(coord3d.getSize()); + setDeltav(coord3d); setAngleBendingFirstDerivative(coord3d); double sumGradientEA; @@ -326,50 +384,50 @@ * *@param coord3d Current molecule coordinates. */ - public void set2ndOrderApproximateGradientMMFF94SumEA(GVector coord3d) { - order2ndApproximateGradientMMFF94SumEA.setSize(coord3d.getSize()); + public void set2ndOrderErrorApproximateGradientMMFF94SumEA(GVector coord3d) { + order2ndErrorApproximateGradientMMFF94SumEA = new GVector(coord3d.getSize()); double sigma = Math.pow(0.000000000000001,0.33); GVector xplusSigma = new GVector(coord3d.getSize()); GVector xminusSigma = new GVector(coord3d.getSize()); - for (int m = 0; m < order2ndApproximateGradientMMFF94SumEA.getSize(); m++) { + for (int m = 0; m < order2ndErrorApproximateGradientMMFF94SumEA.getSize(); m++) { xplusSigma.set(coord3d); xplusSigma.setElement(m,coord3d.getElement(m) + sigma); xminusSigma.set(coord3d); xminusSigma.setElement(m,coord3d.getElement(m) - sigma); - order2ndApproximateGradientMMFF94SumEA.setElement(m,(functionMMFF94SumEA(xplusSigma) - functionMMFF94SumEA(xminusSigma)) / (2 * sigma)); + order2ndErrorApproximateGradientMMFF94SumEA.setElement(m,(functionMMFF94SumEA(xplusSigma) - functionMMFF94SumEA(xminusSigma)) / (2 * sigma)); } - //logger.debug("order2ndApproximateGradientMMFF94SumEA : " + order2ndApproximateGradientMMFF94SumEA); + //logger.debug("order2ndErrorApproximateGradientMMFF94SumEA : " + order2ndErrorApproximateGradientMMFF94SumEA); } /** - * Get the 2nd order approximate gradient of the angle bending term. + * Get the 2nd order error approximate gradient of the angle bending term. * * - *@return Angle bending 2nd order approximate gradient value. + *@return Angle bending 2nd order error approximate gradient value. */ - public GVector get2ndOrderApproximateGradientMMFF94SumEA() { - return order2ndApproximateGradientMMFF94SumEA; + public GVector get2ndOrderErrorApproximateGradientMMFF94SumEA() { + return order2ndErrorApproximateGradientMMFF94SumEA; } /** - * Evaluate an 5 order approximation of the gradient, of the angle bending term, for a given atoms + * Evaluate a 5th order error approximation of the gradient, of the angle bending term, for a given atoms * coordinates * *@param coords3d Current molecule coordinates. */ - public void set5OrderApproximateGradientMMFF94SumEA(GVector coord3d) { - order5ApproximateGradientMMFF94SumEA.setSize(coord3d.getSize()); + public void set5thOrderErrorApproximateGradientMMFF94SumEA(GVector coord3d) { + order5thErrorApproximateGradientMMFF94SumEA = new GVector(coord3d.getSize()); double sigma = Math.pow(0.000000000000001,0.2); GVector xplusSigma = new GVector(coord3d.getSize()); GVector xminusSigma = new GVector(coord3d.getSize()); GVector xplus2Sigma = new GVector(coord3d.getSize()); GVector xminus2Sigma = new GVector(coord3d.getSize()); - for (int m=0; m < order5ApproximateGradientMMFF94SumEA.getSize(); m++) { + for (int m=0; m < order5thErrorApproximateGradientMMFF94SumEA.getSize(); m++) { xplusSigma.set(coord3d); xplusSigma.setElement(m,coord3d.getElement(m) + sigma); xminusSigma.set(coord3d); @@ -378,10 +436,10 @@ xplus2Sigma.setElement(m,coord3d.getElement(m) + 2 * sigma); xminus2Sigma.set(coord3d); xminus2Sigma.setElement(m,coord3d.getElement(m) - 2 * sigma); - order5ApproximateGradientMMFF94SumEA.setElement(m, (8 * (functionMMFF94SumEA(xplusSigma) - functionMMFF94SumEA(xminusSigma)) - (functionMMFF94SumEA(xplus2Sigma) - functionMMFF94SumEA(xminus2Sigma))) / (12 * sigma)); + order5thErrorApproximateGradientMMFF94SumEA.setElement(m, (8 * (functionMMFF94SumEA(xplusSigma) - functionMMFF94SumEA(xminusSigma)) - (functionMMFF94SumEA(xplus2Sigma) - functionMMFF94SumEA(xminus2Sigma))) / (12 * sigma)); } - //logger.debug("order5ApproximateGradientMMFF94SumEA : " + order5ApproximateGradientMMFF94SumEA); + //logger.debug("order5thErrorApproximateGradientMMFF94SumEA : " + order5thErrorApproximateGradientMMFF94SumEA); } @@ -390,8 +448,8 @@ * *@return Angle bending 5 order approximate gradient value. */ - public GVector get5OrderApproximateGradientMMFF94SumEA() { - return order5ApproximateGradientMMFF94SumEA; + public GVector get5thOrderErrorApproximateGradientMMFF94SumEA() { + return order5thErrorApproximateGradientMMFF94SumEA; } @@ -807,7 +865,7 @@ /** - * Get the bond lengths second derivative respect to the cartesian coordinates of the atoms. + * Get the angle bending second derivative respect to the cartesian coordinates of the atoms. * *@return Delta angle bending second derivative value [dimension(3xN)] [angles Number] */ @@ -817,6 +875,52 @@ /** + * Evaluate a 2nd order approximation of the Hessian, for the angle bending, + * given the atoms coordinates. + * + *@param coord3d Current molecule coordinates. + */ + public void setAngleBending2ndOrderErrorApproximateHessian(GVector coord3d) { + angleBendingOrder2ndErrorApproximateHessian = new double[coord3d.getSize()][][]; + + double sigma = Math.pow(0.000000000000001,0.33); + GVector xminusSigma = new GVector(coord3d.getSize()); + GVector xplusSigma = new GVector(coord3d.getSize()); + double[][] gradientAtXminusSigma = null; + double[][] gradientAtXplusSigma = null; + + for (int i = 0; i < coord3d.getSize(); i++) { + xminusSigma.set(coord3d); + xminusSigma.setElement(i,coord3d.getElement(i) - sigma); + setAngleBending2ndOrderErrorApproximateGradient(xminusSigma); + gradientAtXminusSigma = this.getAngleBending2ndOrderErrorApproximateGradient(); + xplusSigma.set(coord3d); + xplusSigma.setElement(i,coord3d.getElement(i) + sigma); + setAngleBending2ndOrderErrorApproximateGradient(xplusSigma); + gradientAtXplusSigma = this.getAngleBending2ndOrderErrorApproximateGradient(); + angleBendingOrder2ndErrorApproximateHessian[i] = new double[coord3d.getSize()][]; + for (int j = 0; j < coord3d.getSize(); j++) { + angleBendingOrder2ndErrorApproximateHessian[i][j] = new double[angleNumber]; + for (int k=0; k < angleNumber; k++) { + angleBendingOrder2ndErrorApproximateHessian[i][j][k] = (gradientAtXplusSigma[j][k] - gradientAtXminusSigma[j][k]) / (2 * sigma); + } + } + } + } + + + /** + * Get the 2nd order error approximate Hessian for the angle bending. + * + * + *@return Angle bending 2nd order error approximate Hessian values. + */ + public double[][][] getAngleBending2ndOrderErrorApproximateHessian() { + return angleBendingOrder2ndErrorApproximateHessian; + } + + + /** * Evaluate the hessian for the angle bending. * *@param coord3d Current molecule coordinates. @@ -825,7 +929,7 @@ double[] forHessian = new double[coord3d.getSize() * coord3d.getSize()]; - calculateDeltav(coord3d); + setDeltav(coord3d); setAngleBendingSecondDerivative(coord3d); double sumHessianEA = 0; @@ -868,5 +972,54 @@ return hessianMMFF94SumEA; } + + /** + * Evaluate a 2nd order approximation of the Hessian, for the angle bending energy term, + * given the atoms coordinates. + * + *@param coord3d Current molecule coordinates. + */ + public void set2ndOrderErrorApproximateHessianMMFF94SumEA(GVector coord3d) { + forOrder2ndErrorApproximateHessian = new double[coord3d.getSize() * coord3d.getSize()]; + + double sigma = Math.pow(0.000000000000001,0.33); + GVector xminusSigma = new GVector(coord3d.getSize()); + GVector xplusSigma = new GVector(coord3d.getSize()); + GVector gradientAtXminusSigma = new GVector(coord3d.getSize()); + GVector gradientAtXplusSigma = new GVector(coord3d.getSize()); + + int forHessianIndex; + for (int i = 0; i < coord3d.getSize(); i++) { + xminusSigma.set(coord3d); + xminusSigma.setElement(i,coord3d.getElement(i) - sigma); + setGradientMMFF94SumEA(xminusSigma); + gradientAtXminusSigma.set(gradientMMFF94SumEA); + xplusSigma.set(coord3d); + xplusSigma.setElement(i,coord3d.getElement(i) + sigma); + setGradientMMFF94SumEA(xplusSigma); + gradientAtXplusSigma.set(gradientMMFF94SumEA); + for (int j = 0; j < coord3d.getSize(); j++) { + forHessianIndex = i*coord3d.getSize()+j; + forOrder2ndErrorApproximateHessian[forHessianIndex] = (gradientAtXplusSigma.getElement(j) - gradientAtXminusSigma.getElement(j)) / (2 * sigma); + //(functionMMFF94SumEA(xplusSigma) - 2 * fx + functionMMFF94SumEA(xminusSigma)) / Math.pow(sigma,2); + } + } + + order2ndErrorApproximateHessianMMFF94SumEA = new GMatrix(coord3d.getSize(), coord3d.getSize()); + order2ndErrorApproximateHessianMMFF94SumEA.set(forOrder2ndErrorApproximateHessian); + //logger.debug("order2ndErrorApproximateHessianMMFF94SumEA : " + order2ndErrorApproximateHessianMMFF94SumEA); + } + + + /** + * Get the 2nd order error approximate Hessian for the angle bending energy term. + * + * + *@return Angle bending energy 2nd order error approximate Hessian value. + */ + public GMatrix get2ndOrderErrorApproximateHessianMMFF94SumEA() { + return order2ndErrorApproximateHessianMMFF94SumEA; + } + } Index: BondStretching.java =================================================================== RCS file: /cvsroot/cdk/cdk/src/org/openscience/cdk/modeling/forcefield/BondStretching.java,v retrieving revision 1.16 retrieving revision 1.17 diff -u -r1.16 -r1.17 --- BondStretching.java 4 May 2005 14:18:06 -0000 1.16 +++ BondStretching.java 13 Jun 2005 13:17:16 -0000 1.17 @@ -26,8 +26,11 @@ ForceFieldTools ffTools = new ForceFieldTools(); double mmff94SumEB = 0; - GVector gradientMMFF94SumEB = new GVector(3); - GMatrix hessianMMFF94SumEB = new GMatrix(3,3); + GVector gradientMMFF94SumEB = null; + GMatrix hessianMMFF94SumEB = null; + double[] forHessian = null; + GMatrix order2ndErrorApproximateHessianMMFF94SumEB = null; + double[] forOrder2ndErrorApproximateHessian = null; int bondsNumber; int[][] bondAtomPosition = null; @@ -124,7 +127,7 @@ /** - * Evaluate the MMFF94 bond stretching term for the given atoms coordinates. + * Evaluate the MMFF94 bond stretching term for the given atoms cartesian coordinates. * *@param coord3d Current molecule coordinates. *@return bond stretching value @@ -152,44 +155,50 @@ /** - * Calculate the bond lengths first derivative respect to the cartesian coordinates of the atoms. + * Set the bond lengths first derivative respect to the cartesian + * coordinates of the atoms. * - *@param coord3d Current molecule coordinates. + *@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) { - + 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[bondsNumber]; - - forAtomNumber = new Double(i/3); + + 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 < bondsNumber; j++) { + 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; + //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]; @@ -204,9 +213,11 @@ /** - * Get the bond lengths derivative respect to the cartesian coordinates of the atoms. + * Get the bond lengths first derivative respect to the cartesian coordinates of the + * atoms. * - *@return Delta bond lengths derivative value [dimension(3xN)] [bonds Number] + *@return Delta bond lengths derivative value [dimension(3xN)] [bonds Number] + * */ public double[][] getBondLengthsFirstDerivative() { return dDeltar; @@ -214,16 +225,16 @@ /** - * Evaluate the gradient for the bond stretching in a given atoms coordinates + * 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.setSize(coord3d.getSize()); + gradientMMFF94SumEB = new GVector(coord3d.getSize()); calculateDeltar(coord3d); - setBondLengthsFirstDerivative(coord3d); + setBondLengthsFirstDerivative(coord3d, deltar, bondAtomPosition); double sumGradientEB; for (int i = 0; i < gradientMMFF94SumEB.getSize(); i++) { @@ -257,7 +268,7 @@ * *@param coord3d Current molecule coordinates. */ - public void setBondLengthsSecondDerivative(GVector coord3d) { + public void setBondLengthsSecondDerivative(GVector coord3d, double[] deltar, int[][] bondAtomPosition) { ddDeltar = new double[coord3d.getSize()][][]; Double forAtomNumber = null; @@ -268,7 +279,12 @@ double ddDeltar1=0; // ddDeltar[i][j][k] = ddDeltar1 - ddDeltar2 double ddDeltar2=0; - setBondLengthsFirstDerivative(coord3d); + setBondLengthsFirstDerivative(coord3d, deltar, bondAtomPosition); + //logger.debug("bondAtomPosition.length = " + bondAtomPosition.length); + double[] rTemp = new double[bondAtomPosition.length]; + for (int i = 0; i < bondAtomPosition.length; i++) { + rTemp[i] = ffTools.distanceBetweenTwoAtomsFrom3xNCoordinates(coord3d, bondAtomPosition[i][0], bondAtomPosition[i][1]); + } for (int i=0; i<coord3d.getSize(); i++) { ddDeltar[i] = new double[coord3d.getSize()][]; @@ -282,7 +298,7 @@ //logger.debug("coordinatei = " + coordinatei); for (int j=0; j<coord3d.getSize(); j++) { - ddDeltar[i][j] = new double[bondsNumber]; + ddDeltar[i][j] = new double[deltar.length]; forAtomNumber = new Double(j/3); @@ -294,7 +310,7 @@ //logger.debug("atomj : " + molecule.getAtomAt(atomNumberj)); - for (int k=0; k < bondsNumber; k++) { + for (int k=0; k < deltar.length; k++) { if ((bondAtomPosition[k][0] == atomNumberj) | (bondAtomPosition[k][1] == atomNumberj)) { if ((bondAtomPosition[k][0] == atomNumberi) | (bondAtomPosition[k][1] == atomNumberi)) { @@ -312,7 +328,7 @@ if (bondAtomPosition[k][1] == atomNumberi) { ddDeltar1 = ddDeltar1 * (-1); } - ddDeltar1 = ddDeltar1 / r[k]; + ddDeltar1 = ddDeltar1 / rTemp[k]; // ddDeltar2 switch (coordinatej) { @@ -349,7 +365,7 @@ //logger.debug("OK: d2 bond 1"); } - ddDeltar2 = ddDeltar2 / Math.pow(r[k],2); + ddDeltar2 = ddDeltar2 / Math.pow(rTemp[k],2); // ddDeltar[i][j][k] ddDeltar[i][j][k] = ddDeltar1 - ddDeltar2; @@ -379,16 +395,16 @@ /** - * Evaluate the hessian for the bond stretching. + * Evaluate the second order partial derivative (hessian) for the bond stretching given the atoms coordinates * *@param coord3d Current molecule coordinates. */ public void setHessianMMFF94SumEB(GVector coord3d) { - double[] forHessian = new double[coord3d.getSize() * coord3d.getSize()]; + forHessian = new double[coord3d.getSize() * coord3d.getSize()]; calculateDeltar(coord3d); - setBondLengthsSecondDerivative(coord3d); + setBondLengthsSecondDerivative(coord3d, deltar, bondAtomPosition); double sumHessianEB; int forHessianIndex; @@ -404,8 +420,7 @@ } } - hessianMMFF94SumEB.setSize(coord3d.getSize(), coord3d.getSize()); - hessianMMFF94SumEB.set(forHessian); + hessianMMFF94SumEB = new GMatrix(coord3d.getSize(), coord3d.getSize(),forHessian); //logger.debug("hessianMMFF94SumEB : " + hessianMMFF94SumEB); } @@ -420,5 +435,63 @@ } + /** + * Get the hessian for the bond stretching. + * + *@return Hessian value of the bond stretching term. + */ + public double[] getForHessianMMFF94SumEB() { + return forHessian; + } + + + /** + * Evaluate a 2nd order approximation of the Hessian, for the bond stretching energy term, + * given the atoms coordinates. + * + *@param coord3d Current molecule coordinates. + */ + public void set2ndOrderErrorApproximateHessianMMFF94SumEB(GVector coord3d) { + forOrder2ndErrorApproximateHessian = new double[coord3d.getSize() * coord3d.getSize()]; + + double sigma = Math.pow(0.000000000000001,0.33); + GVector xminusSigma = new GVector(coord3d.getSize()); + GVector xplusSigma = new GVector(coord3d.getSize()); + GVector gradientAtXminusSigma = new GVector(coord3d.getSize()); + GVector gradientAtXplusSigma = new GVector(coord3d.getSize()); + + int forHessianIndex; + for (int i = 0; i < coord3d.getSize(); i++) { + xminusSigma.set(coord3d); + xminusSigma.setElement(i,coord3d.getElement(i) - sigma); + setGradientMMFF94SumEB(xminusSigma); + gradientAtXminusSigma.set(gradientMMFF94SumEB); + xplusSigma.set(coord3d); + xplusSigma.setElement(i,coord3d.getElement(i) + sigma); + setGradientMMFF94SumEB(xplusSigma); + gradientAtXplusSigma.set(gradientMMFF94SumEB); + for (int j = 0; j < coord3d.getSize(); j++) { + forHessianIndex = i*coord3d.getSize()+j; + forOrder2ndErrorApproximateHessian[forHessianIndex] = (gradientAtXplusSigma.getElement(j) - gradientAtXminusSigma.getElement(j)) / (2 * sigma); + //(functionMMFF94SumEB(xplusSigma) - 2 * fx + functionMMFF94SumEB(xminusSigma)) / Math.pow(sigma,2); + } + } + + order2ndErrorApproximateHessianMMFF94SumEB = new GMatrix(coord3d.getSize(), coord3d.getSize()); + order2ndErrorApproximateHessianMMFF94SumEB.set(forOrder2ndErrorApproximateHessian); + //logger.debug("order2ndErrorApproximateHessianMMFF94SumEB : " + order2ndErrorApproximateHessianMMFF94SumEB); + } + + + /** + * Get the 2nd order error approximate Hessian for the bond stretching term. + * + * + *@return Bond stretching 2nd order error approximate Hessian value. + */ + public GMatrix get2ndOrderErrorApproximateHessianMMFF94SumEB() { + return order2ndErrorApproximateHessianMMFF94SumEB; + } + } Index: ConjugateGradientMethod.java =================================================================== RCS file: /cvsroot/cdk/cdk/src/org/openscience/cdk/modeling/forcefield/ConjugateGradientMethod.java,v retrieving revision 1.12 retrieving revision 1.13 diff -u -r1.12 -r1.13 --- ConjugateGradientMethod.java 4 May 2005 14:18:07 -0000 1.12 +++ ConjugateGradientMethod.java 13 Jun 2005 13:17:16 -0000 1.13 @@ -53,7 +53,9 @@ public void setvk(GVector gk, int iterNumber) { if (iterNumber != 1) { - if (gk.angle(gkminus1) > 1) { + //logger.debug("gk.angle(gkminus1) = " + gk.angle(gkminus1)); + //if (gk.angle(gkminus1) > 1) { + if (iterNumber % 25 != 0) { vkminus1.set(vk); setuk(gkminus1,gk); vkminus1.scale(uk); Index: GeometricMinimizer.java =================================================================== RCS file: /cvsroot/cdk/cdk/src/org/openscience/cdk/modeling/forcefield/GeometricMinimizer.java,v retrieving revision 1.21 retrieving revision 1.22 diff -u -r1.21 -r1.22 --- GeometricMinimizer.java 4 May 2005 14:18:07 -0000 1.21 +++ GeometricMinimizer.java 13 Jun 2005 13:17:16 -0000 1.22 @@ -29,6 +29,7 @@ GVector kCoordinates = null; GVector kplus1Coordinates = null; + GVector g0 = null; int atomNumbers = 1; double fxk = 0; @@ -49,6 +50,7 @@ double d = 0; ForceFieldTools ffTools = new ForceFieldTools(); + NewtonRaphsonMethod nrm = new NewtonRaphsonMethod(); private LoggingTool logger; Molecule molecule; @@ -109,7 +111,6 @@ } public void initializeMinimizationParameters(GVector initialCoord) { - //kCoordinates.setSize(dimension); //kCoordinates.set(initialCoord); kCoordinates=initialCoord; @@ -123,8 +124,9 @@ iterationNumberk = -1; atomNumbers = initialCoord.getSize()/3; - - + + g0 = new GVector(kCoordinates.getSize()); + g0.zero(); } @@ -133,6 +135,7 @@ */ public void checkConvergence(double convergenceCriterion) { + //logger.debug("Checking convergence : "); RMS = 0; RMS = gradientkplus1.dot(gradientkplus1); RMS = RMS / kCoordinates.getSize(); @@ -202,12 +205,16 @@ GVector sk = new GVector(gradientk.getSize()); GVector skplus1 = new GVector(gradientk.getSize()); - LineSearch ls = new LineSearch(); + if (gradientk.equals(g0)) { + convergence = true; + kplus1Coordinates.set(kCoordinates); + } //logger.debug(""); //logger.debug("FORCEFIELDTESTS steepestDescentTest"); SteepestDescentsMethod sdm = new SteepestDescentsMethod(kCoordinates); + LineSearch ls = new LineSearch(); while ((iterationNumberkplus1 < SDMaximumIteration) & (convergence == false)) { @@ -226,7 +233,7 @@ sk.set(skplus1); } //logger.debug("Search direction: "); - sdm.setSk(gradientk, iterationNumberkplus1); + sdm.setSk(gradientk); skplus1.set(sdm.sk); ls.setLineSearchLambda(kCoordinates, sdm.sk, forceField); @@ -240,14 +247,14 @@ checkConvergence(SDconvergenceCriterion); - if (iterationNumberkplus1 % 50 == 0) { +/* if (iterationNumberkplus1 % 50 == 0) { logger.debug(""); - //logger.debug("f(x" + iterationNumberk + ") = " + fxk); - logger.debug("f(x" + iterationNumberkplus1 + ") = " + fxkplus1); - logger.debug("fxkplus1 - fxk = " + (fxkplus1 - fxk)); + logger.debug("f(x" + iterationNumberk + ") = " + fxk); +*/ logger.debug("f(x" + iterationNumberkplus1 + ") = " + fxkplus1); +/* logger.debug("fxkplus1 - fxk = " + (fxkplus1 - fxk)); logger.debug("gradientkplus1, gradientk angle = " + gradientkplus1.angle(gradientk)); } - +*/ /*if (iterationNumberkplus1 != 1) { logger.debug("sk+1.sk = " + skplus1.dot(sk)); logger.debug("gk+1.gk = " + gradientkplus1.dot(gradientk)); @@ -255,7 +262,7 @@ if (molecule != null){ //logger.debug("STEEPESTDM: kplus1Coordinates:"+kplus1Coordinates); - ffTools.assignCoordinatesToMolecule(kplus1Coordinates, molecule); + ffTools.assignCoordinatesToMolecule(kplus1Coordinates, molecule); } } steepestDescentsMinimum.set(kplus1Coordinates); @@ -309,22 +316,27 @@ initializeMinimizationParameters(initialCoordinates); fxk = forceField.energyFunction(initialCoordinates); - logger.debug("f(x0) = " + fxk); + //logger.debug("f(x0) = " + fxk); forceField.setEnergyGradient(kCoordinates); gradientk.set(forceField.getEnergyGradient()); - //logger.debug("gradient at iteration 1 : g1 = " + gradientk); + logger.debug("gradient at iteration 1 : g1 = " + gradientk); conjugateGradientMinimum = new GVector(kCoordinates); - LineSearch ls = new LineSearch(); + + if (gradientk.equals(g0)) { + convergence = true; + kplus1Coordinates.set(kCoordinates); + } ConjugateGradientMethod cgm = new ConjugateGradientMethod(); + LineSearch ls = new LineSearch(); while ((iterationNumberkplus1 < CGMaximumIteration) & (convergence == false)) { iterationNumberkplus1 += 1; iterationNumberk += 1; - //logger.debug(""); + logger.debug(""); //logger.debug(""); //logger.debug("CG Iteration number: " + iterationNumberkplus1); @@ -340,33 +352,38 @@ setkplus1Coordinates(cgm.vk, ls.lineSearchLambda); fxkplus1 = forceField.energyFunction(kplus1Coordinates); //logger.debug("x" + iterationNumberkplus1 + " = " + kplus1Coordinates + " "); - //logger.debug("f(x" + iterationNumberkplus1 + ") = " + fxkplus1); + logger.debug("f(x" + iterationNumberkplus1 + ") = " + fxkplus1); forceField.setEnergyGradient(kplus1Coordinates); gradientkplus1.set(forceField.getEnergyGradient()); checkConvergence(CGconvergenceCriterion); + //logger.debug(" "); + //logger.debug("convergence = " + convergence); - if (iterationNumberkplus1 % 50 == 0 | iterationNumberkplus1 % 50 == 49 | iterationNumberkplus1 % 50 == 1) { - logger.debug(""); +// if (iterationNumberkplus1 % 50 == 0 | iterationNumberkplus1 % 50 == 49 | iterationNumberkplus1 % 50 == 1) { + //logger.debug(""); //logger.debug("f(x" + iterationNumberk + ") = " + fxk); - logger.debug("f(x" + iterationNumberkplus1 + ") = " + fxkplus1); - logger.debug("fxkplus1 - fxk = " + (fxkplus1 - fxk)); - logger.debug("gradientkplus1, gradientk angle = " + gradientkplus1.angle(gradientk)); - } + //logger.debug("f(x" + iterationNumberkplus1 + ") = " + fxkplus1); +// logger.debug("fxkplus1 - fxk = " + (fxkplus1 - fxk)); +// logger.debug("gradientkplus1, gradientk angle = " + gradientkplus1.angle(gradientk)); +// } /*if (iterationNumberkplus1 != 1) { logger.debug("gk+1.gk = " + gradientkplus1.dot(gradientk)); }*/ - if (molecule !=null){ //logger.debug("CGM: kplus1Coordinates:"+kplus1Coordinates); - ffTools.assignCoordinatesToMolecule(kplus1Coordinates, molecule); + ffTools.assignCoordinatesToMolecule(kplus1Coordinates, molecule); } - + + //forceField.setEnergyHessian(kplus1Coordinates); + //nrm.determinat(gradientkplus1, forceField.getEnergyHessian()); + //nrm.hessianEigenValues(forceField.getForEnergyHessian(), kplus1Coordinates.getSize()); + } conjugateGradientMinimum.set(kplus1Coordinates); - forceField.setEnergyGradient(conjugateGradientMinimum); + //logger.debug("conjugateGradientMinimum, forceField.getEnergyGradient().norm() = " + forceField.getEnergyGradient().norm()); //logger.debug("The CG minimum energy is at: " + conjugateGradientMinimum); //logger.debug("f(minimum) = " + fxkplus1); @@ -415,7 +432,6 @@ //logger.debug(""); //logger.debug("FORCEFIELDTESTS NewtonRaphsonTest"); - NewtonRaphsonMethod nrm = new NewtonRaphsonMethod(); GMatrix hessian = new GMatrix(initialCoordinates.getSize(),initialCoordinates.getSize()); while ((iterationNumberkplus1 < NRMaximumIteration) & (convergence == false)) { Index: LineSearch.java =================================================================== RCS file: /cvsroot/cdk/cdk/src/org/openscience/cdk/modeling/forcefield/LineSearch.java,v retrieving revision 1.9 retrieving revision 1.10 diff -u -r1.9 -r1.10 --- LineSearch.java 4 May 2005 14:18:08 -0000 1.9 +++ LineSearch.java 13 Jun 2005 13:17:16 -0000 1.10 @@ -20,13 +20,14 @@ GVector direction = new GVector(3); double lambdaa = 0; - double lambdab = 1; + double lambdab = 0.5; double lambdac = 0; double fxa = 0; double fxb = 0; double fxc = 0; double parabolMinimumLambda = 0; double lineSearchLambda = 0; + double lambdabOld = 0; private LoggingTool logger; @@ -39,58 +40,98 @@ /** - * Bracket the minimum using the Golden Section Search. - * The bracketing phase determines the range of points on the search line to be consider for the interpolation. - * Look for 3 points along the line where the energy of the middle point is lower than the energy of the two outer points. - * The bracket corresponds to an interval specifying the range of values of Lambda. + * Initially Bracketing a Minimum: Look for 3 points along the line search + * where the energy of the middle point is lower than the energy of the two outer points. + * (Numerical recipes in C++. The Art of Scientific Computing. Second Edition.) * *@param kPoint Current point, xk *@param searchDirection Search direction *@param forceFieldFunction Potential energy function */ - public void goldenSectionSearch(GVector kPoint, GVector searchDirection, PotentialFunction forceFieldFunction) { + public void initiallyBracketingAMinimum(GVector kPoint, GVector searchDirection, PotentialFunction forceFieldFunction) { - //logger.debug("Bracketing the minimum using the Golden Section Search:"); + //logger.debug("Initially Bracketing a Minimum:"); - double lambda1 = 0; GVector x1 = new GVector(kPoint); double fx1 = forceFieldFunction.energyFunction(x1); - double lambda4 = 100 * lambdab; - double lambda2 = 0.3819660112 * (lambda4 - lambda1); - double lambda3 = 0.6180339887498948482 * (lambda4 - lambda1); - - direction.set(searchDirection); direction.scale(lambda2); direction.add(kPoint); - GVector x2 = new GVector(direction); double fx2 = forceFieldFunction.energyFunction(x2); - direction.set(searchDirection); direction.scale(lambda3); direction.add(kPoint); - GVector x3 = new GVector(direction); double fx3 = forceFieldFunction.energyFunction(x3); - direction.set(searchDirection); direction.scale(lambda4); direction.add(kPoint); - GVector x4 = new GVector(direction); double fx4 = forceFieldFunction.energyFunction(x4); + direction.set(searchDirection); + lambdaa = 0; + if (lambdab < 0.5) {} + else { + lambdab = 0.5; + } + + double gold = 1.618034; //Default ratio by which successive intervals are magnified. + double glimit = 100; //Maximum magnification allowed for a parabolic-fit step. + double tiny = 1.0E-20; //It is used to prevent any possible division by zero. + + double lambdau = 0; + double r = 0; + double q = 0; + double fxu = 0; + double ulim = 0; + + GVector xa = new GVector(kPoint); fxa = forceFieldFunction.energyFunction(xa); + GVector xb = new GVector(kPoint.getSize()); xb.scaleAdd(lambdab, direction, kPoint); //Sets the value of this vector to the scalar multiplication by s of vector v1 plus vector v2 (this = s*v1 + v2) + fxb = forceFieldFunction.energyFunction(xb); + + if (fxb >fxa) { + double lambdat = lambdaa; GVector xt = new GVector(xa); double fxt = fxa; + lambdaa = lambdab; xa.set(xb); fxa = fxb; + lambdab = lambdat; xb.set(xt); fxb = fxt; + } - boolean finish = false; - while (finish != true) { - if (fx2 < fx3) {//we can bracket the minimum by the interval [x1, x3] - lambda4 = lambda3; x4.set(x3); fx4 = fx3; - lambda3 = lambda2; x3.set(x2); fx3 = fx2; - lambda2 = lambda1 + 0.3819660112 * (lambda4 - lambda1); - direction.set(searchDirection); direction.scale(lambda2); direction.add(kPoint); - x2.set(direction); fx2 = forceFieldFunction.energyFunction(x2); - } else {//we can bracket the minimum by the interval [x2, x4] - lambda1 = lambda2; x1.set(x2); fx1 = fx2; - lambda2 = lambda3; x2.set(x3); fx2 = fx3; - lambda3 = lambda1 + 0.6180339887498948482 * (lambda4 - lambda1); - direction.set(searchDirection); direction.scale(lambda3); direction.add(kPoint); - x3.set(direction); fx3 = forceFieldFunction.energyFunction(x3); - } - if (fx4-fx1<0.0000001) { - finish = true; - if (fx2 < fx3) {lineSearchLambda = lambda2;} - else {lineSearchLambda = lambda3;} + lambdac = lambdab + gold * (lambdab - lambdaa); + GVector xc = new GVector(kPoint.getSize()); xc.scaleAdd(lambdac, direction, kPoint); + fxc = forceFieldFunction.energyFunction(xc); + + double sign = 0; + GVector xu = new GVector(kPoint.getSize()); + + while (fxb > fxc) { + r = (lambdab-lambdaa) * (fxb-fxc); + q = (lambdab-lambdac) * (fxb-fxa); + if (q-r > 0) {sign = 1;} + else {sign = -1;} + lambdau = lambdab - ((lambdab-lambdac) * q - (lambdab-lambdaa) * r) / (2.0 * sign * Math.max(Math.abs(q-r),tiny)); + ulim = lambdab + glimit * (lambdac - lambdab); //We won't go farther than this. + + //Test various possibilities: + if ((lambdab - lambdau) * (lambdau - lambdac) > 0.0) { //Parabolic u is between b and c + xu.scaleAdd(lambdau, direction, kPoint); fxu = forceFieldFunction.energyFunction(xu); + if (fxu < fxc) { // Got a minimum between b and c + lambdaa = lambdab; fxa = fxb; + lambdab = lambdau; fxb = fxu; + } else { + if (fxu > fxb) { // Got a minimum between a and u + lambdac = lambdau; + fxc = fxu; + } + } + lambdau = lambdac + gold * (lambdac-lambdab); // Parabolic fit was not use. Use default magnification. + xu.scaleAdd(lambdau, direction, kPoint); fxu = forceFieldFunction.energyFunction(xu); + } else if ((lambdac-lambdau) * (lambdau-ulim) > 0.0) { // Parabolic fit is between c and its allowed limit. + xu.scaleAdd(lambdau, direction, kPoint); fxu = forceFieldFunction.energyFunction(xu); + if (fxu < fxc) { + //shft3(lambdab,lambdac,lambdau,lambdau+gold*(lambdau-lambdac)); //? + //shft3(fxb,fxc,fxu,func(u)); + } + } else if (lambdau > ulim) { // Limit parabolic u to maximum allowed value + lambdau = ulim; + xu.scaleAdd(lambdau, direction, kPoint); fxu = forceFieldFunction.energyFunction(xu); + } else { // Reject parabolic u, use default magnification. + lambdau = lambdac + gold * (lambdac-lambdab); + xu.scaleAdd(lambdau, direction, kPoint); fxu = forceFieldFunction.energyFunction(xu); } - System.out.println(); - System.out.println("lambda1= " + lambda1 + " ; fx1 = " + fx1); - System.out.println("lambda2= " + lambda2 + " ; fx2 = " + fx2); - System.out.println("lambda3= " + lambda3 + " ; fx3 = " + fx3); - System.out.println("lambda4= " + lambda4 + " ; fx4 = " + fx4); - System.out.println("finish = " + finish); + //shft3(lambdaa,lambdab,lambdac,lambdau); // ? Eliminate oldest point and continue + //shft3(fxa,fxb,fxc,fxu); } + + //logger.debug("lambdaa = " + lambdaa); + //logger.debug("fxa = " + fxa); + //logger.debug("lambdab = " + lambdab); + //logger.debug("fxb = " + fxb); + //logger.debug("lambdac = " + lambdac); + //logger.debug("fxc = " + fxc); + return; } @@ -112,7 +153,10 @@ GVector xa = new GVector(kPoint); fxa = forceFieldFunction.energyFunction(xa); - lambdab = 2 * lambdab; + if (lambdab < 0.5) {} + else { + lambdab = 0.5; + } GVector xb = new GVector(kPoint); direction.set(searchDirection); direction.scale(lambdab); @@ -130,23 +174,25 @@ if (fxb < fxa) { + //logger.debug("Brent's exponential search"); + lambdac = 1.2 * lambdab; xc.set(kPoint); direction.set(searchDirection); direction.scale(lambdac); xc.add(direction); fxc = forceFieldFunction.energyFunction(xc); - + //logger.debug("lambdaa = " + lambdaa + " "); //logger.debug("fxa = " + fxa); //logger.debug("lambdab = " + lambdab + " "); //logger.debug("fxb = " + fxb); //logger.debug("lambdac = " + lambdac + " "); //logger.debug("fxc = " + fxc); - + while (finish == false) { - - if (fxc > fxb) { + + if (fxc >= fxb) { finish = true; } else { @@ -155,10 +201,11 @@ fxa = fxb; xb.set(xc); + lambdabOld = lambdab; lambdab = lambdac; fxb = fxc; - lambdac = 1.2 * lambdac; + lambdac = 1.618 * (lambdac-lambdabOld) + lambdac; // Brent's exponential search xc.set(kPoint); direction.set(searchDirection); direction.scale(lambdac); @@ -171,37 +218,107 @@ //logger.debug("fxb = " + fxb); //logger.debug("lambdac = " + lambdac + " "); //logger.debug("fxc = " + fxc); - + + if (fxc < 0.0000000000000001) { + finish = true; + lineSearchLambda = lambdab; + //logger.debug("fxb < 0.0000000000000001"); + //logger.debug("lambdaa = " + lambdaa + " "); + //logger.debug("fxa = " + fxa); + //logger.debug("lambdab = " + lambdab + " "); + //logger.debug("fxb = " + fxb); + //logger.debug("lambdac = " + lambdac + " "); + //logger.debug("fxc = " + fxc); + } } } } else { - while (finish == false) { - - xc.set(xb); - lambdac = lambdab; - fxc = fxb; - - lambdab = lambdab / 2; - xb.set(kPoint); - direction.set(searchDirection); - direction.scale(lambdab); - xb.add(direction); - fxb = forceFieldFunction.energyFunction(xb); - - //logger.debug("lambdaa = " + lambdaa + " "); - //logger.debug("fxa = " + fxa); - //logger.debug("lambdab = " + lambdab + " "); - //logger.debug("fxb = " + fxb); - //logger.debug("lambdac = " + lambdac + " "); - //logger.debug("fxc = " + fxc); - - if (fxb < fxa) { + //logger.debug("Golden Section Method"); + xc.set(xb); + lambdac = lambdab; + fxc = fxb; + + double lambdab1 = 0.3819660112 * (lambdac - lambdaa); + double lambdab2 = 0.6180339887498948482 * (lambdac - lambdaa); + + GVector xb1 = new GVector(kPoint); + direction.set(searchDirection); + direction.scale(lambdab1); + xb1.add(direction); + double fxb1 = forceFieldFunction.energyFunction(xb1); + + GVector xb2 = new GVector(kPoint); + direction.set(searchDirection); + direction.scale(lambdab2); + xb2.add(direction); + double fxb2 = forceFieldFunction.energyFunction(xb2); + + //logger.debug("lambdaa = " + lambdaa + " "); + //logger.debug("fxa = " + fxa); + //logger.debug("lambdab1 = " + lambdab1 + " "); + //logger.debug("fxb1 = " + fxb1); + //logger.debug("lambdab2 = " + lambdab2 + " "); + //logger.debug("fxb2 = " + fxb2); + //logger.debug("lambdac = " + lambdac + " "); + //logger.debug("fxc = " + fxc); + + while (finish != true) { + if (fxb1 <= fxb2) {//we can bracket the minimum by the interval [lambdaa, lambdab2] + lambdac = lambdab2; xc.set(xb2); fxc = fxb2; + if (fxa > fxb1) { + finish = true; + lambdab = lambdab1; xb.set(xb1); fxb = fxb1; + } else { + lambdab2 = lambdab1; xb2.set(xb1); fxb2 = fxb1; + lambdab1 = lambdaa + 0.3819660112 * (lambdac - lambdaa); + xb1.set(kPoint); direction.set(searchDirection); direction.scale(lambdab1); + xb1.add(direction); fxb1 = forceFieldFunction.energyFunction(xb1); + } + } else {//we can bracket the minimum by the interval [lambdab1, lambdac] + if (fxa > fxb1) { + lambdaa = lambdab1; xa.set(xb1); fxa = fxb1; + if (fxb2 < fxc) { + finish = true; + lambdab = lambdab2; xb.set(xb2); fxb = fxb2; + } else { + lambdab1 = lambdab2; xb1.set(xb2); fxb1 = fxb2; + lambdab2 = lambdaa + 0.6180339887498948482 * (lambdac - lambdaa); + xb2.set(kPoint); direction.set(searchDirection); direction.scale(lambdab2); + xb2.add(direction); fxb2 = forceFieldFunction.energyFunction(xb2); + } + } else {//we can bracket the minimum by the interval [lambdaa, lambdab2] + lambdac = lambdab2; xc.set(xb2); fxc = fxb2; + lambdab2 = lambdab1; xb2.set(xb1); fxb2 = fxb1; + lambdab1 = lambdaa + 0.3819660112 * (lambdac - lambdaa); + xb1.set(kPoint); direction.set(searchDirection); direction.scale(lambdab1); + xb1.add(direction); fxb1 = forceFieldFunction.energyFunction(xb1); + } + } + if (Math.abs(fxc-fxa) < 0.000000000000001) { finish = true; + logger.debug("fxc-fxa < 0.00000000001"); + if (fxb1 < fxb2) {lineSearchLambda = lambdab1;} + else {lineSearchLambda = lambdab2;} } - if (lambdab < 0.0000000000000001) { + if (lambdab1 < 0.0000000000000001) { finish = true; + logger.debug("lambdab < 0.0000000000000001"); + lineSearchLambda = lambdaa; + //logger.debug("lambdaa = " + lambdaa + " "); + //logger.debug("fxa = " + fxa); + //logger.debug("lambdab = " + lambdab1 + " "); + //logger.debug("fxb = " + fxb1); + //logger.debug("lambdac = " + lambdac + " "); + //logger.debug("fxc = " + fxc); } + + //logger.debug(" "); + //logger.debug("lambdaa= " + lambdaa + " ; fxa = " + fxa); + //logger.debug("lambdab1= " + lambdab1 + " ; fxb1 = " + fxb1); + //logger.debug("lambdab2= " + lambdab2 + " ; fxb2 = " + fxb2); + //logger.debug("lambdac= " + lambdac + " ; fxc = " + fxc); + //logger.debug("finish = " + finish); } } return; @@ -224,27 +341,7 @@ *@return The stepSize value */ public double getStepSize() { - return lambdab; - } - - - /** - * xk+1= Xk + Lambdak Sk - * - *@param oldCoordinates Coordinates of the previous step, k - *@param currentStepSize Step size estimated - *@return New coordinates of the atoms, k+1 step - */ - public GVector setCoordinates(GVector oldCoordinates, double currentStepSize) { - GVector coordinates = new GVector(oldCoordinates.getSize()); - return coordinates; - } - - - /** - */ - public void setStepSize() { - return; + return lineSearchLambda; } @@ -261,9 +358,13 @@ direction.setSize(searchDirection.getSize()); -// if (lambdab > 0.00000001) { - bracketingTheMinimum(kPoint, searchDirection, forceFieldFunction); - if (fxb < fxa & fxb < fxc) { + bracketingTheMinimum(kPoint, searchDirection, forceFieldFunction); + if (fxb < fxa & fxb < fxc) { + if (Math.abs(fxc-fxa)/Math.abs(lambdac-lambdaa) < 1) { + goldenSectionMethod(kPoint, searchDirection, forceFieldFunction); + } else { + //logger.debug("Quadratic interpolation"); + GVector xI = new GVector(kPoint.getSize()); double fxI = 0; double fMinimumMinus1 = fxa; @@ -311,16 +412,67 @@ } }*/ lineSearchLambda = lambdab; - }else {lineSearchLambda = 0;} - //}else { - // goldenSectionSearch(kPoint, searchDirection, forceFieldFunction); - //} + } + } + } + + + /** + * Bracket the minimum using the Golden Section Search. + * The bracketing phase determines the range of points on the search line to be consider for the interpolation. + * Look for 3 points along the line where the energy of the middle point is lower than the energy of the two outer points. + * The bracket corresponds to an interval specifying the range of values of Lambda. + * + *@param kPoint Current point, xk + *@param searchDirection Search direction + *@param forceFieldFunction Potential energy function + */ + public void goldenSectionMethod(GVector kPoint, GVector searchDirection, PotentialFunction forceFieldFunction) { + + logger.debug("Golden Section Search"); - //logger.debug("xb = " + xb + " "); - //logger.debug("fxb = " + fxb); - - //logger.debug("Successive parabolic interpolations : "); + double lambda1 = 0; GVector x1 = new GVector(kPoint); double fx1 = forceFieldFunction.energyFunction(x1); + double lambda4 = lambdac; + double lambda2 = 0.3819660112 * (lambda4 - lambda1); + double lambda3 = 0.6180339887498948482 * (lambda4 - lambda1); + direction.set(searchDirection); direction.scale(lambda2); direction.add(kPoint); + GVector x2 = new GVector(direction); double fx2 = forceFieldFunction.energyFunction(x2); + direction.set(searchDirection); direction.scale(lambda3); direction.add(kPoint); + GVector x3 = new GVector(direction); double fx3 = forceFieldFunction.energyFunction(x3); + direction.set(searchDirection); direction.scale(lambda4); direction.add(kPoint); + GVector x4 = new GVector(direction); double fx4 = forceFieldFunction.energyFunction(x4); + + boolean finish = false; + while (finish != true) { + if (fx2 < fx3) {//we can bracket the minimum by the interval [x1, x3] + lambda4 = lambda3; x4.set(x3); fx4 = fx3; + lambda3 = lambda2; x3.set(x2); fx3 = fx2; + lambda2 = lambda1 + 0.3819660112 * (lambda4 - lambda1); + direction.set(searchDirection); direction.scale(lambda2); direction.add(kPoint); + x2.set(direction); fx2 = forceFieldFunction.energyFunction(x2); + } else {//we can bracket the minimum by the interval [x2, x4] + lambda1 = lambda2; x1.set(x2); fx1 = fx2; + lambda2 = lambda3; x2.set(x3); fx2 = fx3; + lambda3 = lambda1 + 0.6180339887498948482 * (lambda4 - lambda1); + direction.set(searchDirection); direction.scale(lambda3); direction.add(kPoint); + x3.set(direction); fx3 = forceFieldFunction.energyFunction(x3); + } + if (fx4-fx1 < 0.000000001) { + finish = true; + if (fx2 < fx3) {lineSearchLambda = lambda2;} + else {lineSearchLambda = lambda3;} + } + //logger.debug(); + //logger.debug("lambda1= " + lambda1 + " ; fx1 = " + fx1); + //logger.debug("lambda2= " + lambda2 + " ; fx2 = " + fx2); + //logger.debug("lambda3= " + lambda3 + " ; fx3 = " + fx3); + //logger.debug("lambda4= " + lambda4 + " ; fx4 = " + fx4); + //logger.debug("finish = " + finish); + } + return; } + + } Index: MMFF94EnergyFunction.java =================================================================== RCS file: /cvsroot/cdk/cdk/src/org/openscience/cdk/modeling/forcefield/MMFF94EnergyFunction.java,v retrieving revision 1.6 retrieving revision 1.7 diff -u -r1.6 -r1.7 --- MMFF94EnergyFunction.java 4 May 2005 14:18:08 -0000 1.6 +++ MMFF94EnergyFunction.java 13 Jun 2005 13:17:16 -0000 1.7 @@ -20,20 +20,22 @@ public class MMFF94EnergyFunction implements PotentialFunction { String energyFunctionShape = " MMFF94 energy "; double energy = 0; - GVector energyGradient = new GVector(3); - GVector order2ndApproximateEnergyGradient = new GVector(3); - GVector order5ApproximateEnergyGradient = new GVector(3); + GVector energyGradient = null; + GVector order2ndErrorApproximateEnergyGradient = new GVector(3); + GVector order5thErrorApproximateEnergyGradient = new GVector(3); GMatrix energyHessian = null; - + double[] forHessian = null; ForceFieldTools fft = new ForceFieldTools(); private LoggingTool logger; BondStretching bs = new BondStretching(); AngleBending ab = new AngleBending(); - //StretchBendInteractions sbi = new StretchBendInteractions(); + StretchBendInteractions sbi = new StretchBendInteractions(); Torsions t =new Torsions(); - //VanDerWaalsInteractions vdwi = new VanDerWaalsInteractions(); + VanDerWaalsInteractions vdwi = new VanDerWaalsInteractions(); + ElectrostaticInteractions ei = new ElectrostaticInteractions(); + int count = 0; /** @@ -44,9 +46,10 @@ //logger.debug(molecule.getAtomCount() + " "+mmff94Tables.size()); bs.setMMFF94BondStretchingParameters(molecule, mmff94Tables); ab.setMMFF94AngleBendingParameters(molecule, mmff94Tables); - //sbi.setMMFF94StretchBendParameters(molecule, mmff94Tables); + sbi.setMMFF94StretchBendParameters(molecule, mmff94Tables); t.setMMFF94TorsionsParameters(molecule, mmff94Tables); - //vdwi.setMMFF94VanDerWaalsParameters(molecule, mmff94Tables); + vdwi.setMMFF94VanDerWaalsParameters(molecule, mmff94Tables); + ei.setMMFF94ElectrostaticParameters(molecule, mmff94Tables); logger = new LoggingTool(this); } @@ -58,17 +61,28 @@ *@return MMFF94 energy function value. */ public double energyFunction(GVector coords3d) { - //logger.debug("bs.functionMMFF94SumEB(coords3d) = " + bs.functionMMFF94SumEB(coords3d)); - //logger.debug("ab.functionMMFF94SumEA(coords3d) = " + ab.functionMMFF94SumEA(coords3d)); - //logger.debug("sbi.functionMMFF94SumEBA(coords3d) = " + sbi.functionMMFF94SumEBA(coords3d)); - //logger.debug("t.functionMMFF94SumET(coords3d) = " + t.functionMMFF94SumET(coords3d)); - //logger.debug("vdwi.functionMMFF94SumEvdW(coords3d) = " + vdwi.functionMMFF94SumEvdW(coords3d)); - - energy = bs.functionMMFF94SumEB(coords3d) - + ab.functionMMFF94SumEA(coords3d) // + sbi.functionMMFF94SumEBA(coords3d) - + t.functionMMFF94SumET(coords3d) ; //+ vdwi.functionMMFF94SumEvdW(coords3d); + vdwi.setFunctionMMFF94SumEvdW(coords3d); + sbi.setFunctionMMFF94SumEBA(coords3d); - //logger.debug("energy = " + energy); + energy = bs.functionMMFF94SumEB(coords3d) + + ab.functionMMFF94SumEA(coords3d) + + sbi.getFunctionMMFF94SumEBA() + + t.functionMMFF94SumET(coords3d) + + vdwi.getFunctionMMFF94SumEvdW() + + ei.functionMMFF94SumEQ(coords3d) + ; + + count += 1; + if (count == 0 | count % 20 == 0) { + logger.debug("count = " + count); + logger.debug("bs.functionMMFF94SumEB(coords3d) = " + bs.functionMMFF94SumEB(coords3d)); + logger.debug("ab.functionMMFF94SumEA(coords3d) = " + ab.functionMMFF94SumEA(coords3d)); + logger.debug("sbi.functionMMFF94SumEBA(coords3d) = " + sbi.getFunctionMMFF94SumEBA()); + logger.debug("t.functionMMFF94SumET(coords3d) = " + t.functionMMFF94SumET(coords3d)); + logger.debug("vdwi.functionMMFF94SumEvdW(coords3d) = " + vdwi.getFunctionMMFF94SumEvdW()); + logger.debug("ei.functionMMFF94SumEQ(coords3d) = " + ei.functionMMFF94SumEQ(coords3d)); + logger.debug("energy = " + energy); + } return energy; } @@ -79,18 +93,19 @@ *@param coords3d Current molecule coordinates. */ public void setEnergyGradient(GVector coords3d) { - //setOrder2ndApproximateEnergyGradient(coords3d); + //setOrder2ndErrorApproximateEnergyGradient(coords3d); //logger.debug("coords3d : " + coords3d); - energyGradient.setSize(coords3d.getSize()); + energyGradient = new GVector(coords3d.getSize()); bs.setGradientMMFF94SumEB(coords3d); - ab.set2ndOrderApproximateGradientMMFF94SumEA(coords3d); - //ab.set5thOrderApproximateGradientMMFF94SumEA(coords3d); - //sbi.setGradientMMFF94SumEBA(coords3d); - t.set2ndOrderApproximateGradientMMFF94SumET(coords3d); - //t.set5thOrderApproximateGradientMMFF94SumET(coords3d); - //vdwi.setGradientMMFF94SumEvdW(coords3d); + ab.set2ndOrderErrorApproximateGradientMMFF94SumEA(coords3d); + //ab.set5thOrderErrorApproximateGradientMMFF94SumEA(coords3d); + sbi.setGradientMMFF94SumEBA(coords3d); + t.set2ndOrderErrorApproximateGradientMMFF94SumET(coords3d); + //t.set5thOrderErrorApproximateGradientMMFF94SumET(coords3d); + vdwi.setGradientMMFF94SumEvdW(coords3d); + ei.setGradientMMFF94SumEQ(coords3d); //logger.debug("bs.getGradientMMFF94SumEB() = " + bs.getGradientMMFF94SumEB()); //logger.debug("ab.getGradientMMFF94SumEA() = " + ab.getGradientMMFF94SumEA()); @@ -98,8 +113,12 @@ for (int i=0; i < energyGradient.getSize(); i++) { energyGradient.setElement(i, bs.getGradientMMFF94SumEB().getElement(i) - + ab.get2ndOrderApproximateGradientMMFF94SumEA().getElement(i) // + sbi.getGradientMMFF94SumEBA().getElement(i) - + t.get2ndOrderApproximateGradientMMFF94SumET().getElement(i) ); // + vdwi.getGradientMMFF94SumEvdW().getElement(i)); + + ab.get2ndOrderErrorApproximateGradientMMFF94SumEA().getElement(i) + + sbi.getGradientMMFF94SumEBA().getElement(i) + + t.get2ndOrderErrorApproximateGradientMMFF94SumET().getElement(i) + + vdwi.getGradientMMFF94SumEvdW().getElement(i) + + ei.getGradientMMFF94SumEQ().getElement(i) + ); } } @@ -111,58 +130,58 @@ */ public GVector getEnergyGradient() { return energyGradient; - //return order2ndApproximateEnergyGradient; + //return order2ndErrorApproximateEnergyGradient; } /** - * Evaluate the order 2 approximate gradient for the MMFF94 energy function in a given 3xN point. + * Evaluate the 2nd order error approximate gradient for the MMFF94 energy function in a given 3xN point. * *@param coords3d Current molecule coordinates. */ - public void setOrder2ndApproximateEnergyGradient(GVector coords3d) { + public void setOrder2ndErrorApproximateEnergyGradient(GVector coords3d) { //logger.debug("coords3d : " + coords3d); - order2ndApproximateEnergyGradient.setSize(coords3d.getSize()); + order2ndErrorApproximateEnergyGradient.setSize(coords3d.getSize()); double sigma = Math.pow(0.0000000000000001,0.3333); GVector xplusSigma = new GVector(coords3d.getSize()); GVector xminusSigma = new GVector(coords3d.getSize()); - for (int i=0; i < order2ndApproximateEnergyGradient.getSize(); i++) { + for (int i=0; i < order2ndErrorApproximateEnergyGradient.getSize(); i++) { xplusSigma.set(coords3d); xplusSigma.setElement(i,coords3d.getElement(i) + sigma); xminusSigma.set(coords3d); xminusSigma.setElement(i,coords3d.getElement(i) - sigma); - order2ndApproximateEnergyGradient.setElement(i, (energyFunction(xplusSigma) - energyFunction(xminusSigma)) / (2 * sigma)); + order2ndErrorApproximateEnergyGradient.setElement(i, (energyFunction(xplusSigma) - energyFunction(xminusSigma)) / (2 * sigma)); } - //logger.debug("order2ndApproximateEnergyGradient : " + order2ndApproximateEnergyGradient); + //logger.debug("order2ndErrorApproximateEnergyGradient : " + order2ndErrorApproximateEnergyGradient); } /** - * Get the order 2 approximate gradient for the MMFF94 energy function in a given 3xN point. + * Get the 2nd o... [truncated message content] |