← 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/StretchBendInteractions.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/StretchBendInteractions.java
 
1
package org.openscience.cdk.modeling.forcefield;
 
2
 
 
3
import java.util.Hashtable;
 
4
import java.util.Vector;
 
5
 
 
6
import javax.vecmath.GMatrix;
 
7
import javax.vecmath.GVector;
 
8
 
 
9
import org.openscience.cdk.interfaces.IAtomContainer;
 
10
import org.openscience.cdk.interfaces.IBond;
 
11
import org.openscience.cdk.interfaces.IAtom;
 
12
import org.openscience.cdk.tools.manipulator.AtomContainerManipulator;
 
13
import org.openscience.cdk.modeling.builder3d.MMFF94ParametersCall;
 
14
import org.openscience.cdk.qsar.IAtomicDescriptor;
 
15
import org.openscience.cdk.qsar.descriptors.atomic.PeriodicTablePositionDescriptor;
 
16
import org.openscience.cdk.qsar.result.IntegerResult;
 
17
import org.openscience.cdk.tools.LoggingTool;
 
18
 
 
19
/**
 
20
 *  Stretch-Bend Interaction calculator for the potential energy function.
 
21
 *  Include function and derivatives.
 
22
 *
 
23
 *@author         vlabarta
 
24
 *@cdk.created    2005-02-15
 
25
 *@cdk.module     forcefield
 
26
 */
 
27
public class StretchBendInteractions {
 
28
 
 
29
        String functionShape = " Stretch-Bend Interactions ";
 
30
 
 
31
        double mmff94SumEBA = 0;
 
32
        GVector gradientMMFF94SumEBA = null;
 
33
        GVector order2ndErrorApproximateGradientMMFF94SumEBA = null;
 
34
        GVector order5thErrorApproximateGradientMMFF94SumEBA = null;
 
35
        GMatrix hessianMMFF94SumEBA = null;
 
36
        GVector currentCoordinates = null;
 
37
        GVector gradientCurrentCoordinates = null;
 
38
 
 
39
        double[][] dDeltarij = null;
 
40
        double[][] dDeltarkj = null;
 
41
        double[][] dDeltav = null;
 
42
 
 
43
        int[][] bondijAtomPosition = null;
 
44
        int[][] bondkjAtomPosition = null;
 
45
        double[] r0IJ = null;
 
46
        double[] r0KJ = null;
 
47
        double[] kbaIJK = null;
 
48
        double[] kbaKJI = null;
 
49
        double[] rij = null;
 
50
        double[] rkj = null;
 
51
        double[] deltarij = null;
 
52
        double[] deltarkj = null;
 
53
 
 
54
        BondStretching bs = new BondStretching();
 
55
        AngleBending ab = new AngleBending();
 
56
        private LoggingTool logger;
 
57
 
 
58
        GVector moleculeCurrentCoordinates = null;
 
59
        boolean[] changeAtomCoordinates = null;
 
60
        int changedCoordinates;
 
61
 
 
62
 
 
63
        /**
 
64
         *  Constructor for the StretchBendInteractions object
 
65
         */
 
66
        public StretchBendInteractions() {
 
67
                logger = new LoggingTool(this);
 
68
        }
 
69
 
 
70
 
 
71
        /**
 
72
         *  Set MMFF94 reference bond lengths r0IJ and r0JK and stretch-bend
 
73
         *  interaction constants kbaIJK and kbaKJI for each i-j-k angle in the
 
74
         *  molecule.
 
75
         *
 
76
         *@param  molecule       The molecule like an AtomContainer object.
 
77
         *@param  parameterSet   MMFF94 parameters set
 
78
         *@exception  Exception  Description of the Exception
 
79
         */
 
80
        public void setMMFF94StretchBendParameters(IAtomContainer molecule, Hashtable parameterSet, boolean angleBendingFlag) throws Exception {
 
81
 
 
82
                //logger.debug("setMMFF94StretchBendParameters");
 
83
                
 
84
                ab.setMMFF94AngleBendingParameters(molecule, parameterSet, angleBendingFlag);
 
85
 
 
86
                IAtom[] atomConnected = null;
 
87
 
 
88
                Vector stretchBendInteractionsData = null;
 
89
                Vector bondData = null;
 
90
                MMFF94ParametersCall pc = new MMFF94ParametersCall();
 
91
                pc.initialize(parameterSet);
 
92
 
 
93
                bondijAtomPosition = new int[ab.angleNumber][];
 
94
                bondkjAtomPosition = new int[ab.angleNumber][];
 
95
                r0IJ = new double[ab.angleNumber];
 
96
                r0KJ = new double[ab.angleNumber];
 
97
                kbaIJK = new double[ab.angleNumber];
 
98
                kbaKJI = new double[ab.angleNumber];
 
99
 
 
100
                String strbndType;
 
101
                String angleType;
 
102
                String bondIJType;
 
103
                String bondKJType;
 
104
 
 
105
                IBond bondIJ = null;
 
106
                IBond bondKJ = null;
 
107
                
 
108
                IAtomicDescriptor descriptor  = new PeriodicTablePositionDescriptor();
 
109
                int iR = 0;
 
110
                int jR = 0;
 
111
                int kR = 0;
 
112
 
 
113
                
 
114
                int l = -1;
 
115
                for (int j = 0; j < molecule.getAtomCount(); j++) {
 
116
                        
 
117
                        atomConnected = AtomContainerManipulator.getAtomArray(molecule.getConnectedAtomsList(molecule.getAtom(j)));
 
118
                        
 
119
                        if (atomConnected.length > 1) {
 
120
                                
 
121
                                for (int i = 0; i < atomConnected.length; i++) {
 
122
                                        
 
123
                                        for (int k = i + 1; k < atomConnected.length; k++) {
 
124
                                                
 
125
                                                l += 1;
 
126
                                                
 
127
                                                bondIJ = molecule.getBond(atomConnected[i], molecule.getAtom(j));
 
128
                                                bondIJType = bondIJ.getProperty("MMFF94 bond type").toString();
 
129
                                                
 
130
                                                bondKJ = molecule.getBond(atomConnected[k], molecule.getAtom(j));
 
131
                                                bondKJType = bondKJ.getProperty("MMFF94 bond type").toString();
 
132
                                                
 
133
                                                angleType = "0";
 
134
                                                if ((bondIJType == "1") | (bondKJType == "1")) {
 
135
                                                        angleType = "1";
 
136
                                                }  
 
137
                                                if ((bondIJType == "1") & (bondKJType == "1")) {
 
138
                                                        angleType = "2";
 
139
                                                }  
 
140
                                                
 
141
                                                //logger.debug("bondIJType = " + bondIJType + ", bondKJType = " + bondKJType + ", angleType = " + angleType);
 
142
                                                
 
143
                                                strbndType = "0";
 
144
                                                if ((angleType == "0") & (bondIJType == "0") & (bondKJType == "0")) {strbndType = "0";}
 
145
                                                else if ((angleType == "1") & (bondIJType == "1") & (bondKJType == "0")) {strbndType = "1";}
 
146
                                                else if ((angleType == "1") & (bondIJType == "0") & (bondKJType == "1")) {strbndType = "2";}
 
147
                                                else if ((angleType == "2") & (bondIJType == "1") & (bondKJType == "1")) {strbndType = "3";}
 
148
                                                else if ((angleType == "4") & (bondIJType == "0") & (bondKJType == "0")) {strbndType = "4";}
 
149
                                                else if ((angleType == "3") & (bondIJType == "0") & (bondKJType == "0")) {strbndType = "5";}
 
150
                                                else if ((angleType == "5") & (bondIJType == "1") & (bondKJType == "0")) {strbndType = "6";}
 
151
                                                else if ((angleType == "5") & (bondIJType == "0") & (bondKJType == "1")) {strbndType = "7";}
 
152
                                                else if ((angleType == "6") & (bondIJType == "1") & (bondKJType == "1")) {strbndType = "8";}
 
153
                                                else if ((angleType == "7") & (bondIJType == "1") & (bondKJType == "0")) {strbndType = "9";}
 
154
                                                else if ((angleType == "7") & (bondIJType == "0") & (bondKJType == "1")) {strbndType = "10";}
 
155
                                                else if ((angleType == "8") & (bondIJType == "1") & (bondKJType == "1")) {strbndType = "11";}
 
156
                                                
 
157
                                                //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)) + ")");
 
158
                                                stretchBendInteractionsData = pc.getBondAngleInteractionData(strbndType, atomConnected[i].getAtomTypeName(), molecule.getAtom(j).getAtomTypeName(), atomConnected[k].getAtomTypeName());
 
159
                                                
 
160
                                                if (stretchBendInteractionsData == null) {
 
161
                                                        if (angleType == "1") {
 
162
                                                                if (strbndType == "1") {strbndType = "2";}
 
163
                                                                else {strbndType = "1";}
 
164
                                                                //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)) + ")");
 
165
                                                                stretchBendInteractionsData = pc.getBondAngleInteractionData(strbndType, atomConnected[i].getAtomTypeName(), molecule.getAtom(j).getAtomTypeName(), atomConnected[k].getAtomTypeName());
 
166
                                                        }
 
167
                                                }
 
168
                                                
 
169
                                                if (stretchBendInteractionsData == null) {
 
170
                                                        iR = ((IntegerResult)descriptor.calculate(atomConnected[i],molecule).getValue()).intValue();
 
171
                                                        jR = ((IntegerResult)descriptor.calculate(molecule.getAtom(j),molecule).getValue()).intValue();
 
172
                                                        kR = ((IntegerResult)descriptor.calculate(atomConnected[k],molecule).getValue()).intValue();
 
173
                                                        stretchBendInteractionsData = pc.getDefaultStretchBendData(iR, jR, kR);
 
174
                                                } 
 
175
 
 
176
                                                //logger.debug("stretchBendInteractionsData : " + stretchBendInteractionsData);
 
177
                                                kbaIJK[l] = ((Double) stretchBendInteractionsData.get(0)).doubleValue();
 
178
                                                kbaKJI[l] = ((Double) stretchBendInteractionsData.get(1)).doubleValue();
 
179
 
 
180
                                                //logger.debug("kbaIJK[" + l + "] = " + kbaIJK[l]);
 
181
                                                //logger.debug("kbaKJI[" + l + "] = " + kbaKJI[l]);
 
182
 
 
183
                                                
 
184
                                                bondData = pc.getBondData(bondIJType, atomConnected[i].getAtomTypeName(), molecule.getAtom(j).getAtomTypeName());
 
185
                                                r0IJ[l] = ((Double) bondData.get(0)).doubleValue();
 
186
                                                bondData = pc.getBondData(bondKJType, atomConnected[k].getAtomTypeName(), molecule.getAtom(j).getAtomTypeName());
 
187
                                                r0KJ[l] = ((Double) bondData.get(0)).doubleValue();
 
188
                                                
 
189
                                                bondijAtomPosition[l] = new int[2];
 
190
                                                bondijAtomPosition[l][0] = molecule.getAtomNumber(atomConnected[i]);
 
191
                                                bondijAtomPosition[l][1] = j;
 
192
                                                
 
193
                                                bondkjAtomPosition[l] = new int[2];
 
194
                                                bondkjAtomPosition[l][0] = molecule.getAtomNumber(atomConnected[k]);
 
195
                                                bondkjAtomPosition[l][1] = j;
 
196
                                        }
 
197
                                }
 
198
                        }
 
199
                }
 
200
                rij = new double[ab.angleNumber];
 
201
                rkj = new double[ab.angleNumber];
 
202
                deltarij = new double[ab.angleNumber];
 
203
                deltarkj = new double[ab.angleNumber];
 
204
                currentCoordinates = new GVector(3 * molecule.getAtomCount());
 
205
                gradientCurrentCoordinates = new GVector(3 * molecule.getAtomCount());
 
206
                gradientMMFF94SumEBA = new GVector(3 * molecule.getAtomCount());
 
207
                dDeltarij = new double[3 * molecule.getAtomCount()][];
 
208
                dDeltarkj = new double[3 * molecule.getAtomCount()][];
 
209
                dDeltav = new double[3 * molecule.getAtomCount()][];
 
210
                hessianMMFF94SumEBA = new GMatrix(3 * molecule.getAtomCount(), 3 * molecule.getAtomCount());
 
211
                for (int i = 0; i < 3 * molecule.getAtomCount(); i++) {
 
212
                        dDeltarij[i] = new double[ab.angleNumber];
 
213
                        dDeltarkj[i] = new double[ab.angleNumber];
 
214
                        dDeltav[i] = new double[ab.angleNumber];
 
215
                }
 
216
                
 
217
                this.moleculeCurrentCoordinates = new GVector(3 * molecule.getAtomCount());
 
218
                for (int i=0; i<moleculeCurrentCoordinates.getSize(); i++) {
 
219
                        this.moleculeCurrentCoordinates.setElement(i,1E10);
 
220
                } 
 
221
 
 
222
                this.changeAtomCoordinates = new boolean[molecule.getAtomCount()];
 
223
 
 
224
        }
 
225
 
 
226
 
 
227
        /**
 
228
         *  Calculate the current bond distances rij and rkj for each angle j, and the
 
229
         *  difference with the reference bonds.
 
230
         *
 
231
         *@param  coords3d  Current molecule coordinates.
 
232
         */
 
233
        public void setDeltarijAndDeltarkj(GVector coords3d) {
 
234
 
 
235
                changedCoordinates = 0;
 
236
                //logger.debug("Setting Deltarij and Deltarkj");
 
237
                for (int i=0; i < changeAtomCoordinates.length; i++) {
 
238
                        this.changeAtomCoordinates[i] = false;
 
239
                }
 
240
                this.moleculeCurrentCoordinates.sub(coords3d);
 
241
                
 
242
                for (int i = 0; i < this.moleculeCurrentCoordinates.getSize(); i++) {
 
243
                        //logger.debug("moleculeCurrentCoordinates " + i + " = " + this.moleculeCurrentCoordinates.getElement(i));
 
244
                        if (Math.abs(this.moleculeCurrentCoordinates.getElement(i)) > 0) {
 
245
                                changeAtomCoordinates[i/3] = true;
 
246
                                changedCoordinates = changedCoordinates + 1;
 
247
                                //logger.debug("changeAtomCoordinates[" + i/3 + "] = " + changeAtomCoordinates[i/3]);
 
248
                                i = i + (2 - i % 3);
 
249
                        }
 
250
                }
 
251
 
 
252
                for (int i = 0; i < ab.angleNumber; i++) {
 
253
                        if ((changeAtomCoordinates[ab.angleAtomPosition[i][0]] == true) | 
 
254
                                        (changeAtomCoordinates[ab.angleAtomPosition[i][1]] == true))            {
 
255
                        
 
256
                                rij[i] = ForceFieldTools.distanceBetweenTwoAtomsFrom3xNCoordinates(coords3d, ab.angleAtomPosition[i][1], ab.angleAtomPosition[i][0]);
 
257
                                deltarij[i] = rij[i] - r0IJ[i];
 
258
                                //logger.debug("deltarij[" + i + "] = " + deltarij[i]);
 
259
                        }
 
260
                        //else {System.out.println("deltarij[" + i + "] was no recalculated");}
 
261
                        if ((changeAtomCoordinates[ab.angleAtomPosition[i][1]] == true) | 
 
262
                                        (changeAtomCoordinates[ab.angleAtomPosition[i][2]] == true))            {
 
263
                        
 
264
                                rkj[i] = ForceFieldTools.distanceBetweenTwoAtomsFrom3xNCoordinates(coords3d, ab.angleAtomPosition[i][1], ab.angleAtomPosition[i][2]);
 
265
                                deltarkj[i] = rkj[i] - r0KJ[i];
 
266
                                //logger.debug("deltarkj[" + i + "] = " + deltarkj[i]);
 
267
                        }
 
268
                        //else {System.out.println("deltarkj[" + i + "] was no recalculated");}
 
269
                }
 
270
                /*if    (changedCoordinates == changeAtomCoordinates.length) {
 
271
                        for (int m = 0; m < ab.angleNumber; m++) {
 
272
                                System.out.println("phi[" + m + "] = " + Math.toDegrees(phi[m]));
 
273
                        }
 
274
                }
 
275
                */
 
276
                moleculeCurrentCoordinates.set(coords3d);
 
277
        }
 
278
 
 
279
 
 
280
        
 
281
        /**
 
282
         *  Set the MMFF94 stretch-bend interaction term given the atoms cartesian
 
283
         *  coordinates.
 
284
         *
 
285
         *@param  coords3d  Current molecule coordinates.
 
286
         */
 
287
        public void setFunctionMMFF94SumEBA(GVector coords3d) {
 
288
                //ab.setAngleBendingFlag(false);
 
289
                if (currentCoordinates.equals(coords3d)) {
 
290
                } else {
 
291
                        setDeltarijAndDeltarkj(coords3d);
 
292
                        ab.setDeltav(coords3d);
 
293
                        mmff94SumEBA = 0;
 
294
                        for (int j = 0; j < ab.angleNumber; j++) {
 
295
                                //logger.debug("kbaIJK[" + j + "] = " + kbaIJK[j]);
 
296
                                //logger.debug("kbaKJI[" + j + "] = " + kbaKJI[j]);
 
297
                                //logger.debug("deltarij[" + j + "] = " + deltarij[j]);
 
298
                                //logger.debug("deltarkj[" + j + "] = " + deltarkj[j]);
 
299
                                //logger.debug("ab.deltav[" + j + "] = " + ab.deltav[j]);
 
300
                                mmff94SumEBA = mmff94SumEBA + 2.51210 * (kbaIJK[j] * deltarij[j] + kbaKJI[j] * deltarkj[j]) * ab.deltav[j];
 
301
                                //logger.debug("mmff94SumEBA = " + mmff94SumEBA);
 
302
                        }
 
303
                        //mmff94SumEBA = Math.abs(mmff94SumEBA);
 
304
                        //logger.debug("mmff94SumEBA = " + mmff94SumEBA);
 
305
                        currentCoordinates.set(coords3d);
 
306
                }
 
307
        }
 
308
 
 
309
 
 
310
        /**
 
311
         *  Get the MMFF94 stretch-bend interaction term.
 
312
         *
 
313
         *@return    MMFF94 stretch-bend interaction term value.
 
314
         */
 
315
        public double getFunctionMMFF94SumEBA() {
 
316
                return mmff94SumEBA;
 
317
        }
 
318
 
 
319
 
 
320
        /**
 
321
         *  Evaluate the gradient of the stretch-bend interaction term.
 
322
         *
 
323
         *@param  coords3d  Current molecule coordinates.
 
324
         */
 
325
        public void setGradientMMFF94SumEBA(GVector coords3d) {
 
326
 
 
327
                if (currentCoordinates.equals(coords3d)) {
 
328
                } else {
 
329
                        setFunctionMMFF94SumEBA(coords3d);
 
330
                }
 
331
 
 
332
                bs.setBondLengthsFirstDerivative(coords3d, deltarij, bondijAtomPosition);
 
333
                dDeltarij = bs.getBondLengthsFirstDerivative();
 
334
                bs.setBondLengthsFirstDerivative(coords3d, deltarkj, bondkjAtomPosition);
 
335
                dDeltarkj = bs.getBondLengthsFirstDerivative();
 
336
                ab.setAngleBending2ndOrderErrorApproximateGradient(coords3d);
 
337
                dDeltav = ab.getAngleBending2ndOrderErrorApproximateGradient();
 
338
 
 
339
                if (dDeltav == null) {logger.debug("setGradient: dDeltav null");} 
 
340
                double sumGradientEBA;
 
341
                for (int i = 0; i < gradientMMFF94SumEBA.getSize(); i++) {
 
342
                        sumGradientEBA = 0;
 
343
                        for (int j = 0; j < ab.angleNumber; j++) {
 
344
                                sumGradientEBA = sumGradientEBA + (kbaIJK[j] * dDeltarij[i][j] + kbaKJI[j] * dDeltarkj[i][j]) * ab.deltav[j]
 
345
                                                 + (kbaIJK[j] * deltarij[j] + kbaKJI[j] * deltarkj[j]) * ab.angleBendingOrder2ndErrorApproximateGradient[i][j];
 
346
                        }
 
347
                        sumGradientEBA = sumGradientEBA * 2.51210;
 
348
 
 
349
                        gradientMMFF94SumEBA.setElement(i, sumGradientEBA);
 
350
                        gradientCurrentCoordinates.set(coords3d);
 
351
                }
 
352
                //logger.debug("gradientMMFF94SumEBA = " + gradientMMFF94SumEBA);
 
353
        }
 
354
 
 
355
 
 
356
        /**
 
357
         *  Get the gradient of the stretch-bend interaction term.
 
358
         *
 
359
         *@return    stretch-bend interaction gradient value.
 
360
         */
 
361
        public GVector getGradientMMFF94SumEBA() {
 
362
                return gradientMMFF94SumEBA;
 
363
        }
 
364
 
 
365
 
 
366
        /**
 
367
         *  Evaluate a 2nd order approximation of the gradient for the stretch-bend interaction term, 
 
368
         *  given the atoms coordinates.
 
369
         *
 
370
         *@param  coord3d  Current molecule coordinates.
 
371
         */
 
372
        public void set2ndOrderErrorApproximateGradientMMFF94SumEBA(GVector coord3d) {
 
373
                order2ndErrorApproximateGradientMMFF94SumEBA = new GVector(coord3d.getSize());
 
374
                double sigma = Math.pow(0.000000000000001,0.33);
 
375
                GVector xplusSigma = new GVector(coord3d.getSize());
 
376
                GVector xminusSigma = new GVector(coord3d.getSize());
 
377
                double fInXplusSigma = 0;
 
378
                double fInXminusSigma = 0;
 
379
 
 
380
                for (int m = 0; m < order2ndErrorApproximateGradientMMFF94SumEBA.getSize(); m++) {
 
381
                        xplusSigma.set(coord3d);
 
382
                        xplusSigma.setElement(m,coord3d.getElement(m) + sigma);
 
383
                        this.setFunctionMMFF94SumEBA(xplusSigma);
 
384
                        fInXplusSigma = this.getFunctionMMFF94SumEBA();
 
385
                        xminusSigma.set(coord3d);
 
386
                        xminusSigma.setElement(m,coord3d.getElement(m) - sigma);
 
387
                        this.setFunctionMMFF94SumEBA(xminusSigma);
 
388
                        fInXminusSigma = this.getFunctionMMFF94SumEBA();
 
389
                        order2ndErrorApproximateGradientMMFF94SumEBA.setElement(m,(fInXplusSigma - fInXminusSigma) / (2 * sigma));
 
390
                }
 
391
                        
 
392
                //logger.debug("order2ndErrorApproximateGradientMMFF94SumEBA : " + order2ndErrorApproximateGradientMMFF94SumEBA);
 
393
        }
 
394
 
 
395
 
 
396
        /**
 
397
         *  Get the 2nd order error approximate gradient for the stretch-bend term.
 
398
         *
 
399
         *
 
400
         *@return           Stretch-bend interaction 2nd order error approximate gradient value.
 
401
         */
 
402
        public GVector get2ndOrderErrorApproximateGradientMMFF94SumEBA() {
 
403
                return order2ndErrorApproximateGradientMMFF94SumEBA;
 
404
        }
 
405
 
 
406
 
 
407
        /**
 
408
         *  Evaluate a 5th order error approximation of the gradient, of the stretch-bend interaction term, for a given atoms
 
409
         *  coordinates
 
410
         *
 
411
         *@param  coord3d  Current molecule coordinates.
 
412
         */
 
413
        public void set5thOrderErrorApproximateGradientMMFF94SumEBA(GVector coord3d) {
 
414
                order5thErrorApproximateGradientMMFF94SumEBA = new GVector(coord3d.getSize());
 
415
                double sigma = Math.pow(0.000000000000001,0.2);
 
416
                GVector xplusSigma = new GVector(coord3d.getSize());
 
417
                GVector xminusSigma = new GVector(coord3d.getSize());
 
418
                GVector xplus2Sigma = new GVector(coord3d.getSize());
 
419
                GVector xminus2Sigma = new GVector(coord3d.getSize());
 
420
                double fInXplusSigma = 0;
 
421
                double fInXminusSigma = 0;
 
422
                double fInXplus2Sigma = 0;
 
423
                double fInXminus2Sigma = 0;
 
424
                
 
425
                for (int m=0; m < order5thErrorApproximateGradientMMFF94SumEBA.getSize(); m++) {
 
426
                        xplusSigma.set(coord3d);
 
427
                        xplusSigma.setElement(m,coord3d.getElement(m) + sigma);
 
428
                        this.setFunctionMMFF94SumEBA(xplusSigma);
 
429
                        fInXplusSigma = this.getFunctionMMFF94SumEBA();
 
430
                        xminusSigma.set(coord3d);
 
431
                        xminusSigma.setElement(m,coord3d.getElement(m) - sigma);
 
432
                        this.setFunctionMMFF94SumEBA(xminusSigma);
 
433
                        fInXminusSigma = this.getFunctionMMFF94SumEBA();
 
434
                        xplus2Sigma.set(coord3d);
 
435
                        xplus2Sigma.setElement(m,coord3d.getElement(m) + 2 * sigma);
 
436
                        this.setFunctionMMFF94SumEBA(xplus2Sigma);
 
437
                        fInXplus2Sigma = this.getFunctionMMFF94SumEBA();
 
438
                        xminus2Sigma.set(coord3d);
 
439
                        xminus2Sigma.setElement(m,coord3d.getElement(m) - 2 * sigma);
 
440
                        this.setFunctionMMFF94SumEBA(xminus2Sigma);
 
441
                        fInXminus2Sigma = this.getFunctionMMFF94SumEBA();
 
442
                        order5thErrorApproximateGradientMMFF94SumEBA.setElement(m, (8 * (fInXplusSigma - fInXminusSigma) - (fInXplus2Sigma - fInXminus2Sigma)) / (12 * sigma));
 
443
                }
 
444
                        
 
445
                //logger.debug("order5thErrorApproximateGradientMMFF94SumEBA : " + order5thErrorApproximateGradientMMFF94SumEBA);
 
446
        }
 
447
 
 
448
 
 
449
        /**
 
450
         *  Get the 5 order approximate gradient of the stretch-bend interaction term.
 
451
         *
 
452
         *@return        stretch-bend interaction 5 order approximate gradient value.
 
453
         */
 
454
        public GVector get5thOrderErrorApproximateGradientMMFF94SumEBA() {
 
455
                return order5thErrorApproximateGradientMMFF94SumEBA;
 
456
        }
 
457
 
 
458
 
 
459
        /**
 
460
         *  Evaluate the hessian of the stretch-bend interaction.
 
461
         *
 
462
         *@param  coords3d  Current molecule coordinates.
 
463
         */
 
464
        public void setHessianMMFF94SumEBA(GVector coords3d) {
 
465
 
 
466
                double[] forHessian = new double[coords3d.getSize() * coords3d.getSize()];
 
467
 
 
468
                if (currentCoordinates.equals(coords3d)) {
 
469
                } else {
 
470
                        setFunctionMMFF94SumEBA(coords3d);
 
471
                }
 
472
 
 
473
                if (dDeltarij == null) {logger.debug("dDeltarij null");} 
 
474
                if (gradientCurrentCoordinates.equals(coords3d) == false) {
 
475
                        bs.setBondLengthsFirstDerivative(coords3d, deltarij, bondijAtomPosition);
 
476
                        dDeltarij = bs.getBondLengthsFirstDerivative();
 
477
                }
 
478
                if (dDeltarkj == null) {logger.debug("dDeltarkj null");} 
 
479
                if (gradientCurrentCoordinates.equals(coords3d) == false) {
 
480
                        bs.setBondLengthsFirstDerivative(coords3d, deltarkj, bondkjAtomPosition);
 
481
                        dDeltarkj = bs.getBondLengthsFirstDerivative();
 
482
                }
 
483
                if (dDeltav == null) {logger.debug("setHessian: dDeltav null");} 
 
484
                if (gradientCurrentCoordinates.equals(coords3d) == false) {
 
485
                        //logger.debug("ab.setAngleBending2ndOrderErrorApproximateGradient()");
 
486
                        ab.setAngleBending2ndOrderErrorApproximateGradient(coords3d);
 
487
                        dDeltav = ab.getAngleBending2ndOrderErrorApproximateGradient();
 
488
                }
 
489
                
 
490
                ab.setAngleBending2ndOrderErrorApproximateHessian(coords3d);
 
491
                double[][][] ddDeltav = ab.getAngleBending2ndOrderErrorApproximateHessian();
 
492
                bs.setBondLengthsSecondDerivative(coords3d, deltarij, bondijAtomPosition);
 
493
                double[][][] ddDeltarij = bs.getBondLengthsSecondDerivative();
 
494
                bs.setBondLengthsSecondDerivative(coords3d, deltarkj, bondkjAtomPosition);
 
495
                double[][][] ddDeltarkj = bs.getBondLengthsSecondDerivative();
 
496
                
 
497
                if (dDeltav == null) {logger.debug("setHessian: dDeltav null");} 
 
498
                //logger.debug("ab.angleNumber = " + ab.angleNumber);
 
499
                double sumHessianEBA;
 
500
                int forHessianIndex;
 
501
                for (int i = 0; i < coords3d.getSize(); i++) {
 
502
                        for (int j = 0; j < coords3d.getSize(); j++) {
 
503
                                forHessianIndex = i*coords3d.getSize()+j;
 
504
                                sumHessianEBA = 0;
 
505
                                for (int k = 0; k < ab.angleNumber; k++) {
 
506
                                        sumHessianEBA = sumHessianEBA + (kbaIJK[k] * ddDeltarij[i][j][k] + kbaKJI[k] * ddDeltarkj[i][j][k]) * ab.deltav[k]
 
507
                                                 + (kbaIJK[k] * dDeltarij[j][k] + kbaKJI[k] * dDeltarkj[j][k]) * dDeltav[j][k]
 
508
                                                 + (kbaIJK[k] * dDeltarij[j][k] + kbaKJI[k] * dDeltarkj[j][k]) * dDeltav[j][k]
 
509
                                                 + (kbaIJK[k] * deltarij[k] + kbaKJI[k] * deltarkj[k]) * ddDeltav[i][j][k];
 
510
                                }
 
511
                                forHessian[forHessianIndex] = sumHessianEBA;
 
512
                        }
 
513
                }
 
514
 
 
515
                hessianMMFF94SumEBA.setSize(coords3d.getSize(), coords3d.getSize());
 
516
                hessianMMFF94SumEBA.set(forHessian);
 
517
                //logger.debug("hessianMMFF94SumEBA : " + hessianMMFF94SumEBA);
 
518
        }
 
519
 
 
520
 
 
521
        /**
 
522
         *  Get the hessian of the stretch-bend interaction.
 
523
         *
 
524
         *@return    Hessian value of the stretch-bend interaction term.
 
525
         */
 
526
        public GMatrix getHessianMMFF94SumEBA() {
 
527
                return hessianMMFF94SumEBA;
 
528
        }
 
529
 
 
530
}
 
531