← Back to branch summary

~ubuntu-branches/ubuntu/maverick/cdk/maverick

« back to all changes in this revision

Viewing changes to src/org/openscience/cdk/modeling/forcefield/NewtonRaphsonMethod.java

  • Committer: Bazaar Package Importer
  • Author(s): Paul Cager
  • Date: 2008-04-09 21:17:53 UTC
  • Revision ID: james.westby@ubuntu.com-20080409211753-46lmjw5z8mx5pd8d
Tags: upstream-1.0.2
Import upstream version 1.0.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/org/openscience/cdk/modeling/forcefield/NewtonRaphsonMethod.java
 
1
/* $Revision: 8201 $ $Author: egonw $ $Date: 2007-04-16 10:40:19 +0200 (Mon, 16 Apr 2007) $
 
2
 *
 
3
 * Copyright (C) 2005-2007  Violeta Labarta <vlabarta@users.sf.net>
 
4
 *
 
5
 * Contact: cdk-devel@lists.sourceforge.net
 
6
 *
 
7
 * This program is free software; you can redistribute it and/or
 
8
 * modify it under the terms of the GNU Lesser General Public License
 
9
 * as published by the Free Software Foundation; either version 2.1
 
10
 * of the License, or (at your option) any later version.
 
11
 * All we ask is that proper credit is given for our work, which includes
 
12
 * - but is not limited to - adding the above copyright notice to the beginning
 
13
 * of your source code files, and to any copyright notice that you may distribute
 
14
 * with programs based on this work.
 
15
 *
 
16
 * This program is distributed in the hope that it will be useful,
 
17
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 
18
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
19
 * GNU Lesser General Public License for more details.
 
20
 *
 
21
 * You should have received a copy of the GNU Lesser General Public License
 
22
 * along with this program; if not, write to the Free Software
 
23
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
 
24
 */
 
25
package org.openscience.cdk.modeling.forcefield;
 
26
 
 
27
import javax.vecmath.GMatrix;
 
28
import javax.vecmath.GVector;
 
29
 
 
30
import org.openscience.cdk.tools.LoggingTool;
 
31
 
 
32
import Jama.Matrix;
 
33
 
 
34
/**
 
35
 * Methods of Newton-Raphson approach.
 
36
 *
 
37
 * @author        vlabarta
 
38
 * @cdk.created   2005-06-01
 
39
 * @cdk.module    forcefield
 
40
 * 
 
41
 * @cdk.keyword   Newton-Raphson
 
42
 */
 
43
public class NewtonRaphsonMethod {
 
44
        GVector gradientPerInverseHessianVector = null;
 
45
        Matrix matrixForDeterminatCalculation = null;
 
46
        private LoggingTool logger;
 
47
 
 
48
 
 
49
        /**
 
50
         *  Constructor for the NR object
 
51
         */
 
52
        public NewtonRaphsonMethod() {
 
53
                logger = new LoggingTool(this);
 
54
        }
 
55
 
 
56
 
 
57
        /**
 
58
         *  Calculate the eigen values for the hessian matrix.
 
59
         *
 
60
         *@param  forMatrix  Hessian matrix
 
61
         *@param  size       coordinates dimension
 
62
         */
 
63
        public void hessianEigenValues(double[] forMatrix, int size) {
 
64
                Matrix A = new Matrix(forMatrix, size);
 
65
                Matrix As = A.plus(A.transpose());      // Simetric matrix: As = 1/2 * (A + AT);
 
66
                As.timesEquals(0.5); 
 
67
                //logger.debug("Simetric matrix Hs = 1/2 * (H + HT) = ");
 
68
                //As.print(As.getRowDimension(), As.getColumnDimension());
 
69
                double[] realEigenvalues = As.eig().getRealEigenvalues();
 
70
                double[] imagEigenvalues = As.eig().getImagEigenvalues();
 
71
                System.out.println(" ");
 
72
                System.out.println("Hs EigenValues :");
 
73
                for (int i=0; i < As.getColumnDimension(); i++) {
 
74
                        System.out.println("Eigen value " + i + ": real part = " + realEigenvalues[i]);
 
75
                        System.out.println(", imaginary part = " + imagEigenvalues[i]);
 
76
                 }
 
77
 
 
78
        }
 
79
 
 
80
    public void determinat(GVector gradientk, GMatrix hessiank) {
 
81
                //logger.debug(" ");
 
82
                //logger.debug("calculate hessian determinat: ");
 
83
                double[][] forDeterminatCalculation = new double[gradientk.getSize()][];
 
84
                for (int i = 0; i < gradientk.getSize(); i++) {
 
85
                        forDeterminatCalculation[i] = new double[gradientk.getSize()];
 
86
                        for (int j = 0; j < forDeterminatCalculation[i].length; j++) {
 
87
                                forDeterminatCalculation[i][j] = hessiank.getElement(i, j);
 
88
                        }
 
89
                }
 
90
 
 
91
                //logger.debug("gradientk.getSize() = " + gradientk.getSize());
 
92
                /*
 
93
                 *  if (gradientk.getSize() == 36) {
 
94
                 *  logger.debug();
 
95
                 *  for (int i = 0; i < forDeterminatCalculation.length; i++) {
 
96
                 *  for (int j = 0; j < forDeterminatCalculation[i].length; j++) {
 
97
                 *  logger.debug(forDeterminatCalculation[i][j] + " ");
 
98
                 *  }
 
99
                 *  logger.debug();
 
100
                 *  }
 
101
                 *  }
 
102
                 */
 
103
                matrixForDeterminatCalculation = new Matrix(forDeterminatCalculation);
 
104
                //matrixForDeterminatCalculation.print(gradientk.getSize(), gradientk.getSize());
 
105
 
 
106
                //logger.debug("matrixForDeterminatCalculation.det() = " + matrixForDeterminatCalculation.det());
 
107
 
 
108
                return;
 
109
        }
 
110
 
 
111
        public void gradientPerInverseHessian(GVector gradientk, GMatrix hessiank) {
 
112
                this.determinat(gradientk, hessiank);
 
113
                if (matrixForDeterminatCalculation.det() != 0) {
 
114
                        hessiank.invert();
 
115
                        //logger.debug("hessiank.invert() = " + hessiank);
 
116
                        gradientPerInverseHessianVector = new GVector(gradientk.getSize());
 
117
                        gradientPerInverseHessianVector.mul(gradientk, hessiank);
 
118
                } else {
 
119
                        logger.debug("The Newton-Raphson method can't be execute because the hessian can't be inverted");
 
120
                }
 
121
                return;
 
122
        }
 
123
 
 
124
 
 
125
        /**
 
126
         *  Gets the gradientPerInverseHessian attribute of the NewtonRaphsonMethod
 
127
         *  object
 
128
         *
 
129
         *@return    The gradientPerInverseHessian value
 
130
         */
 
131
        public GVector getGradientPerInverseHessian() {
 
132
                return gradientPerInverseHessianVector;
 
133
        }
 
134
 
 
135
}
 
136