← Back to branch summary

~ubuntu-branches/ubuntu/trusty/cdk/trusty-proposed

« back to all changes in this revision

Viewing changes to src/org/openscience/cdk/modeling/forcefield/ElectrostaticInteractions.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/ElectrostaticInteractions.java
 
1
package org.openscience.cdk.modeling.forcefield;
 
2
 
 
3
import java.util.Hashtable;
 
4
 
 
5
import javax.vecmath.GMatrix;
 
6
import javax.vecmath.GVector;
 
7
 
 
8
import org.openscience.cdk.interfaces.IAtomContainer;
 
9
import org.openscience.cdk.qsar.IAtomicDescriptor;
 
10
import org.openscience.cdk.qsar.descriptors.atomic.BondsToAtomDescriptor;
 
11
import org.openscience.cdk.qsar.result.IntegerResult;
 
12
//import org.openscience.cdk.tools.LoggingTool;
 
13
 
 
14
 
 
15
/**                     
 
16
 *  MMFF94 Electrostatic Interactions energy. Include function and derivatives.
 
17
 *
 
18
 *@author     vlabarta
 
19
 *@cdk.created    May 13, 2005
 
20
 *@cdk.module     forcefield
 
21
 */
 
22
public class ElectrostaticInteractions {
 
23
 
 
24
        String functionShape = " Electrostatic Interactions ";
 
25
 
 
26
        double mmff94SumEQ = 0;
 
27
        GVector gradientMMFF94SumEQ = null;
 
28
        GVector order2ndErrorApproximateGradientMMFF94SumEQ = null;
 
29
        GVector order5thErrorApproximateGradientMMFF94SumEQ = null;
 
30
        GMatrix hessianMMFF94SumEQ = null;
 
31
        double[] forHessian = null;
 
32
        
 
33
        double[][] dR = null;   // internuclear separation first order derivative respect to atoms coordinates
 
34
        double[][][] ddR = null;        // internuclear separation second order derivative respect to atoms coordinates
 
35
 
 
36
        IAtomicDescriptor shortestPathBetweenTwoAtoms = new BondsToAtomDescriptor();
 
37
        Object[] params = {new Integer(0)};
 
38
        
 
39
        int electrostaticInteractionNumber;
 
40
        int[][] electrostaticInteractionAtomPosition = null;
 
41
 
 
42
        double[] r = null;      // internuclear separation in Angstroms.
 
43
        double[] qi = null;
 
44
        double[] qj = null;
 
45
 
 
46
        double delta = 0.05;    //electrostatic buffering constant.
 
47
        double n = 1;
 
48
        double D = 1.0;
 
49
 
 
50
        double[] iQ = null;
 
51
        double electrostatic14interactionsScale = 0.75; // Scale factor for 1-4 interactions. To take in the future from mmff94.prm files.
 
52
        
 
53
        //private LoggingTool logger;
 
54
        
 
55
        /**
 
56
         *  Constructor for the ElectrostaticInteractions object
 
57
         */
 
58
        public ElectrostaticInteractions() {        
 
59
                //logger = new LoggingTool(this);
 
60
        }
 
61
 
 
62
 
 
63
        /**
 
64
         *  Set CCG Electrostatic parameters for the molecule.
 
65
         *
 
66
         *
 
67
         *@param  molecule       The molecule like an AtomContainer object.
 
68
         *@param  parameterSet   MMFF94 parameters set
 
69
         *@exception  Exception  Description of the Exception
 
70
         */
 
71
        public void setMMFF94ElectrostaticParameters(IAtomContainer molecule, Hashtable parameterSet) throws Exception {
 
72
 
 
73
                //distances = wnd.getShortestPathLengthBetweenAtoms((AtomContainer) molecule);
 
74
                
 
75
                //logger.debug("molecule.getAtomCount() : " + molecule.getAtomCount());
 
76
                //logger.debug("molecule.getBondCount() : " + molecule.getBondCount());
 
77
                if (molecule.getAtomCount() == 12 & molecule.getBondCount() == 11) {
 
78
                    molecule.getAtom(3).setCharge(1);
 
79
                    molecule.getAtom(8).setCharge(1);
 
80
                }
 
81
                
 
82
                electrostaticInteractionNumber = 0;
 
83
                for (int i=0; i<molecule.getAtomCount(); i++) {
 
84
                        for (int j=i+1; j<molecule.getAtomCount(); j++) {
 
85
                                params[0] = new Integer(j);
 
86
                                shortestPathBetweenTwoAtoms.setParameters(params);
 
87
                                //if (distances[molecule.getAtomNumber(molecule.getAtomAt(i))][molecule.getAtomNumber(molecule.getAtomAt(j))]>2) {
 
88
                                if (((IntegerResult)shortestPathBetweenTwoAtoms.calculate(molecule.getAtom(i),molecule).getValue()).intValue()>2){
 
89
                                        electrostaticInteractionNumber += 1;
 
90
                                }
 
91
                        }
 
92
                }
 
93
                //logger.debug("electrostaticInteractionNumber : " + electrostaticInteractionNumber);
 
94
 
 
95
                qi = new double[electrostaticInteractionNumber];
 
96
                qj = new double[electrostaticInteractionNumber];
 
97
                r = new double[electrostaticInteractionNumber];
 
98
                iQ = new double[electrostaticInteractionNumber];
 
99
                
 
100
                electrostaticInteractionAtomPosition = new int[electrostaticInteractionNumber][];
 
101
                
 
102
                int l = -1;
 
103
                for (int i=0; i<molecule.getAtomCount(); i++) {
 
104
                        for (int j=i+1; j<molecule.getAtomCount(); j++) {
 
105
                                params[0] = new Integer(j);
 
106
                                shortestPathBetweenTwoAtoms.setParameters(params);
 
107
                                //if (distances[molecule.getAtomNumber(molecule.getAtomAt(i))][molecule.getAtomNumber(molecule.getAtomAt(j))]>2) {
 
108
                                if (((IntegerResult)shortestPathBetweenTwoAtoms.calculate(molecule.getAtom(i),molecule).getValue()).intValue()>2){
 
109
                                        l += 1;
 
110
                                        qi[l]= molecule.getAtom(i).getCharge();
 
111
                                        qj[l]= molecule.getAtom(j).getCharge();
 
112
                                        //logger.debug("qi[" + l + "] = " + qi[l] + ", qj[" + l + "] = " + qj[l]);
 
113
                                        if (((IntegerResult)shortestPathBetweenTwoAtoms.calculate(molecule.getAtom(i),molecule).getValue()).intValue()==3){
 
114
                                                iQ[l] = electrostatic14interactionsScale;
 
115
                                        } else {
 
116
                                                iQ[l] = 1;
 
117
                                        }
 
118
                                        
 
119
                                        electrostaticInteractionAtomPosition[l] = new int[2];
 
120
                                        electrostaticInteractionAtomPosition[l][0] = i;
 
121
                                        electrostaticInteractionAtomPosition[l][1] = j;
 
122
                                        //logger.debug("electrostaticInteractionAtomPosition " + l + " : " + electrostaticInteractionAtomPosition[l][0] + ", " + electrostaticInteractionAtomPosition[l][1]);
 
123
                                }
 
124
                        }
 
125
                }
 
126
        }
 
127
 
 
128
 
 
129
        /**
 
130
         *  Calculate the internuclear separation Rij
 
131
         *
 
132
         *@param  coords3d  Current molecule coordinates.
 
133
         */
 
134
        public void setInternuclearSeparation(GVector coords3d) {
 
135
 
 
136
                for (int l = 0; l < electrostaticInteractionNumber; l++) {
 
137
 
 
138
                        r[l] = ForceFieldTools.distanceBetweenTwoAtomsFrom3xNCoordinates(coords3d, electrostaticInteractionAtomPosition[l][0], electrostaticInteractionAtomPosition[l][1]);
 
139
                        //logger.debug("r[" + l + "]= " + r[l]);
 
140
                }
 
141
        }
 
142
 
 
143
 
 
144
        /**
 
145
         *  Get the internuclear separation values (Rij).
 
146
         *
 
147
         *@return        Internuclear separation values.
 
148
         */
 
149
        public double[] getInternuclearSeparation() {
 
150
                return r;
 
151
        }
 
152
 
 
153
 
 
154
        /**
 
155
         *  Evaluate the MMFF94 Electrostatic interaction energy.
 
156
         *
 
157
         *@param  coords3d  Current molecule coordinates.
 
158
         *@return        MMFF94 Electrostatic interaction term value.
 
159
         */
 
160
        public double functionMMFF94SumEQ(GVector coords3d) {
 
161
                setInternuclearSeparation(coords3d);
 
162
                mmff94SumEQ = 0;
 
163
                for (int l = 0; l < electrostaticInteractionNumber; l++) {
 
164
                        mmff94SumEQ = mmff94SumEQ + iQ[l] * 332.0716 * qi[l] * qj[l] / (D * Math.pow(r[l] + delta, n));
 
165
                        //logger.debug("mmff94SumEQ = " + mmff94SumEQ);
 
166
                }
 
167
                //logger.debug("mmff94SumEQ = " + mmff94SumEQ);
 
168
                return mmff94SumEQ;
 
169
        }
 
170
 
 
171
 
 
172
        /**
 
173
         *  Calculate the internuclear separation (Rij) first derivative respect to the cartesian coordinates of the atoms.
 
174
         *
 
175
         *@param  coord3d  Current molecule coordinates.
 
176
         */
 
177
        public void setInternuclearSeparationFirstOrderDerivative(GVector coord3d) {
 
178
                
 
179
                dR = new double[coord3d.getSize()][];
 
180
                setInternuclearSeparation(coord3d);
 
181
                
 
182
                Double forAtomNumber = null;
 
183
                int atomNumber = 0;
 
184
                int coordinate;
 
185
                for (int i = 0; i < dR.length; i++) {
 
186
                        
 
187
                        dR[i] = new double[electrostaticInteractionNumber];
 
188
                        
 
189
                        forAtomNumber = new Double(i/3);
 
190
                        coordinate = i % 3;
 
191
                        //logger.debug("coordinate = " + coordinate);
 
192
 
 
193
                        atomNumber = forAtomNumber.intValue();
 
194
                        //logger.debug("atomNumber = " + atomNumber);
 
195
 
 
196
                        for (int j = 0; j < electrostaticInteractionNumber; j++) {
 
197
 
 
198
                                if ((electrostaticInteractionAtomPosition[j][0] == atomNumber) | (electrostaticInteractionAtomPosition[j][1] == atomNumber)) {
 
199
                                        //logger.debug("atomNumber, r[" + j + "] = " + r[j]);
 
200
                                        switch (coordinate) {
 
201
                                                //x-coordinate
 
202
                                                case 0: dR[i][j] = (coord3d.getElement(3 * electrostaticInteractionAtomPosition[j][0]) - coord3d.getElement(3 * electrostaticInteractionAtomPosition[j][1])) / r[j];
 
203
                                                        //logger.debug("xi-xj = " + (coord3d.getElement(3 * electrostaticInteractionAtomPosition[j][0]) - coord3d.getElement(3 * electrostaticInteractionAtomPosition[j][1])));
 
204
                                                        break;
 
205
                                                //y-coordinate
 
206
                                                case 1: dR[i][j] = (coord3d.getElement(3 * electrostaticInteractionAtomPosition[j][0] + 1) - coord3d.getElement(3 * electrostaticInteractionAtomPosition[j][1] + 1)) / r[j];
 
207
                                                        //logger.debug("xi-xj = " + (coord3d.getElement(3 * electrostaticInteractionAtomPosition[j][0]) - coord3d.getElement(3 * electrostaticInteractionAtomPosition[j][1])));
 
208
                                                        break;
 
209
                                                //z-coordinate
 
210
                                                case 2: dR[i][j] = (coord3d.getElement(3 * electrostaticInteractionAtomPosition[j][0] + 2) - coord3d.getElement(3 * electrostaticInteractionAtomPosition[j][1] + 2)) / r[j];
 
211
                                                        //logger.debug("xi-xj = " + (coord3d.getElement(3 * electrostaticInteractionAtomPosition[j][0]) - coord3d.getElement(3 * electrostaticInteractionAtomPosition[j][1])));
 
212
                                                        break;
 
213
                                        }
 
214
                                        if (electrostaticInteractionAtomPosition[j][1] == atomNumber) {
 
215
                                                dR[i][j] = (-1) * dR[i][j];
 
216
                                        }
 
217
                                } else {
 
218
                                        dR[i][j] = 0;
 
219
                                }
 
220
                                //logger.debug("electrostaticInteraction " + j + " : " + "dR[" + i + "][" + j + "] = " + dR[i][j]);
 
221
                        }
 
222
                }
 
223
        }
 
224
 
 
225
 
 
226
        /**
 
227
         *  Get the internuclear separation first order derivative respect to the cartesian coordinates of the atoms.
 
228
         *
 
229
         *@return        Internuclear separation first order derivative value [dimension(3xN)] [vdW interaction number]
 
230
         */
 
231
        public double[][] getInternuclearSeparationFirstDerivative() {
 
232
                return dR;
 
233
        }
 
234
 
 
235
 
 
236
        /**
 
237
         *  Set the gradient of the MMFF94 Electrostatic interaction term.
 
238
         *
 
239
         *
 
240
         *@param  coords3d  Current molecule coordinates.
 
241
         */
 
242
        public void setGradientMMFF94SumEQ(GVector coords3d) {
 
243
 
 
244
                gradientMMFF94SumEQ = new GVector(coords3d.getSize());
 
245
                setInternuclearSeparation(coords3d);
 
246
                setInternuclearSeparationFirstOrderDerivative(coords3d);
 
247
                
 
248
                double sumGradientEQ;
 
249
                for (int i = 0; i < gradientMMFF94SumEQ.getSize(); i++) {
 
250
                        sumGradientEQ = 0;
 
251
                        for (int l = 0; l < electrostaticInteractionNumber; l++) {
 
252
                                sumGradientEQ = sumGradientEQ + ((-332.0716 * qi[l] * qj[l] * n)/(D * Math.pow(r[l] + delta, n+1))) * dR[i][l]; 
 
253
                        }
 
254
                        gradientMMFF94SumEQ.setElement(i, sumGradientEQ);
 
255
                }
 
256
                //logger.debug("gradientMMFF94SumEQ = " + gradientMMFF94SumEQ);
 
257
        }
 
258
 
 
259
 
 
260
        /**
 
261
         *  Get the gradient of the MMFF94 Electrostatic interaction term.
 
262
         *
 
263
         *
 
264
         *@return           MMFF94 Electrostatic interaction gradient value.
 
265
         */
 
266
        public GVector getGradientMMFF94SumEQ() {
 
267
                return gradientMMFF94SumEQ;
 
268
        }
 
269
 
 
270
 
 
271
        /**
 
272
         *  Evaluate a 2nd order error approximation of the gradient, for the Electrostatic interaction term, given the atoms
 
273
         *  coordinates
 
274
         *
 
275
         *@param  coord3d  Current molecule coordinates.
 
276
         */
 
277
/*      public void set2ndOrderErrorApproximateGradientMMFF94SumEQ(GVector coord3d) {
 
278
                order2ndErrorApproximateGradientMMFF94SumEQ = new GVector(coord3d.getSize());
 
279
                double sigma = Math.pow(0.000000000000001,0.33);
 
280
                GVector xplusSigma = new GVector(coord3d.getSize());
 
281
                GVector xminusSigma = new GVector(coord3d.getSize());
 
282
                
 
283
                for (int m = 0; m < order2ndErrorApproximateGradientMMFF94SumEQ.getSize(); m++) {
 
284
                        xplusSigma.set(coord3d);
 
285
                        xplusSigma.setElement(m,coord3d.getElement(m) + sigma);
 
286
                        xminusSigma.set(coord3d);
 
287
                        xminusSigma.setElement(m,coord3d.getElement(m) - sigma);
 
288
                        order2ndErrorApproximateGradientMMFF94SumEQ.setElement(m,(functionMMFF94SumEQ(xplusSigma) - functionMMFF94SumEQ(xminusSigma)) / (2 * sigma));
 
289
                }
 
290
                        
 
291
                //logger.debug("order2ndErrorApproximateGradientMMFF94SumEQ : " + order2ndErrorApproximateGradientMMFF94SumEQ);
 
292
        }
 
293
*/
 
294
 
 
295
        /**
 
296
         *  Get the 2nd order error approximate gradient for the Electrostatic interaction term.
 
297
         *
 
298
         *
 
299
         *@return           Electrostatic interaction 2nd order error approximate gradient value
 
300
         */
 
301
/*      public GVector get2ndOrderErrorApproximateGradientMMFF94SumEQ() {
 
302
                return order2ndErrorApproximateGradientMMFF94SumEQ;
 
303
        }
 
304
*/
 
305
 
 
306
        /**
 
307
         *  Evaluate a 5th order error approximation of the gradient, for the Electrostatic interaction term, given the atoms
 
308
         *  coordinates
 
309
         *
 
310
         *@param  coords3d  Current molecule coordinates.
 
311
         */
 
312
/*      public void set5thOrderErrorApproximateGradientMMFF94SumEQ(GVector coord3d) {
 
313
                order5thErrorApproximateGradientMMFF94SumEQ = new GVector(coord3d.getSize());
 
314
                double sigma = Math.pow(0.000000000000001,0.2);
 
315
                GVector xplusSigma = new GVector(coord3d.getSize());
 
316
                GVector xminusSigma = new GVector(coord3d.getSize());
 
317
                GVector xplus2Sigma = new GVector(coord3d.getSize());
 
318
                GVector xminus2Sigma = new GVector(coord3d.getSize());
 
319
                
 
320
                for (int m=0; m < order5thErrorApproximateGradientMMFF94SumEQ.getSize(); m++) {
 
321
                        xplusSigma.set(coord3d);
 
322
                        xplusSigma.setElement(m,coord3d.getElement(m) + sigma);
 
323
                        xminusSigma.set(coord3d);
 
324
                        xminusSigma.setElement(m,coord3d.getElement(m) - sigma);
 
325
                        xplus2Sigma.set(coord3d);
 
326
                        xplus2Sigma.setElement(m,coord3d.getElement(m) + 2 * sigma);
 
327
                        xminus2Sigma.set(coord3d);
 
328
                        xminus2Sigma.setElement(m,coord3d.getElement(m) - 2 * sigma);
 
329
                        order5thErrorApproximateGradientMMFF94SumEQ.setElement(m, (8 * (functionMMFF94SumEQ(xplusSigma) - functionMMFF94SumEQ(xminusSigma)) - (functionMMFF94SumEQ(xplus2Sigma) - functionMMFF94SumEQ(xminus2Sigma))) / (12 * sigma));
 
330
                }
 
331
                        
 
332
                //logger.debug("order5thErrorApproximateGradientMMFF94SumEQ : " + order5thErrorApproximateGradientMMFF94SumEQ);
 
333
        }
 
334
*/
 
335
 
 
336
        /**
 
337
         *  Get the 5th order error approximate gradient for the Electrostatic interaction term.
 
338
         *
 
339
         *@return        Electrostatic interaction 5th order error approximate gradient value.
 
340
         */
 
341
/*      public GVector get5OrderApproximateGradientMMFF94SumEQ() {
 
342
                return order5thErrorApproximateGradientMMFF94SumEQ;
 
343
        }
 
344
*/
 
345
 
 
346
 
 
347
        /**
 
348
         *  Calculate the internuclear separation second derivative respect to the cartesian coordinates of the atoms.
 
349
         *
 
350
         *@param  coord3d  Current molecule coordinates.
 
351
         */
 
352
        public void setInternuclearSeparationSecondDerivative(GVector coord3d) {
 
353
                ddR = new double[coord3d.getSize()][][];
 
354
                
 
355
                Double forAtomNumber = null;
 
356
                int atomNumberi;
 
357
                int atomNumberj;
 
358
                int coordinatei;
 
359
                int coordinatej;
 
360
                double ddR1=0;  // ddR[i][j][k] = ddR1 - ddR2
 
361
                double ddR2=0;
 
362
                
 
363
                setInternuclearSeparationFirstOrderDerivative(coord3d);
 
364
                
 
365
                for (int i=0; i<coord3d.getSize(); i++) {
 
366
                        ddR[i] = new double[coord3d.getSize()][];
 
367
                        
 
368
                        forAtomNumber = new Double(i/3);
 
369
                        
 
370
                        atomNumberi = forAtomNumber.intValue();
 
371
                        //logger.debug("atomNumberi = " + atomNumberi);
 
372
                                
 
373
                        coordinatei = i % 3;
 
374
                        //logger.debug("coordinatei = " + coordinatei);
 
375
                                
 
376
                        for (int j=0; j<coord3d.getSize(); j++) {
 
377
                                ddR[i][j] = new double[electrostaticInteractionNumber];
 
378
                                
 
379
                                forAtomNumber = new Double(j/3);
 
380
 
 
381
                                atomNumberj = forAtomNumber.intValue();
 
382
                                //logger.debug("atomNumberj = " + atomNumberj);
 
383
 
 
384
                                coordinatej = j % 3;
 
385
                                //logger.debug("coordinatej = " + coordinatej);
 
386
                                
 
387
                                //logger.debug("atomj : " + molecule.getAtomAt(atomNumberj));
 
388
                                
 
389
                                for (int k=0; k < electrostaticInteractionNumber; k++) {
 
390
                                        
 
391
                                        if ((electrostaticInteractionAtomPosition[k][0] == atomNumberj) | (electrostaticInteractionAtomPosition[k][1] == atomNumberj)) {
 
392
                                                if ((electrostaticInteractionAtomPosition[k][0] == atomNumberi) | (electrostaticInteractionAtomPosition[k][1] == atomNumberi)) {
 
393
                                        
 
394
                                                        // ddR1
 
395
                                                        if (electrostaticInteractionAtomPosition[k][0] == atomNumberj) {
 
396
                                                                ddR1 = 1;
 
397
                                                        }
 
398
                                                        if (electrostaticInteractionAtomPosition[k][1] == atomNumberj) {
 
399
                                                                ddR1 = -1;
 
400
                                                        }
 
401
                                                        if (electrostaticInteractionAtomPosition[k][0] == atomNumberi) {
 
402
                                                                ddR1 = ddR1 * 1;
 
403
                                                        }
 
404
                                                        if (electrostaticInteractionAtomPosition[k][1] == atomNumberi) {
 
405
                                                                ddR1 = ddR1 * (-1);
 
406
                                                        }
 
407
                                                        ddR1 = ddR1 / r[k];
 
408
 
 
409
                                                        // ddR2
 
410
                                                        switch (coordinatej) {
 
411
                                                                case 0: ddR2 = (coord3d.getElement(3 * electrostaticInteractionAtomPosition[k][0]) - coord3d.getElement(3 * electrostaticInteractionAtomPosition[k][1]));
 
412
                                                                        //logger.debug("OK: d1 x");
 
413
                                                                        break;
 
414
                                                                case 1: ddR2 = (coord3d.getElement(3 * electrostaticInteractionAtomPosition[k][0] + 1) - coord3d.getElement(3 * electrostaticInteractionAtomPosition[k][1] + 1));
 
415
                                                                        //logger.debug("OK: d1 y");
 
416
                                                                        break;
 
417
                                                                case 2: ddR2 = (coord3d.getElement(3 * electrostaticInteractionAtomPosition[k][0] + 2) - coord3d.getElement(3 * electrostaticInteractionAtomPosition[k][1] + 2));
 
418
                                                                        //logger.debug("OK: d1 z");
 
419
                                                                        break;
 
420
                                                        }
 
421
                                                
 
422
                                                        if (electrostaticInteractionAtomPosition[k][1] == atomNumberj) {
 
423
                                                                ddR2 = (-1) * ddR2;
 
424
                                                                //logger.debug("OK: bond 1");
 
425
                                                        } 
 
426
        
 
427
                                                        switch (coordinatei) {
 
428
                                                                case 0: ddR2 = ddR2 * (coord3d.getElement(3 * electrostaticInteractionAtomPosition[k][0]) - coord3d.getElement(3 * electrostaticInteractionAtomPosition[k][1]));
 
429
                                                                        //logger.debug("OK: have d2 x");
 
430
                                                                        break;
 
431
                                                                case 1: ddR2 = ddR2 * (coord3d.getElement(3 * electrostaticInteractionAtomPosition[k][0] + 1) - coord3d.getElement(3 * electrostaticInteractionAtomPosition[k][1] + 1));
 
432
                                                                        //logger.debug("OK: have d2 y");
 
433
                                                                        break;
 
434
                                                                case 2: ddR2 = ddR2 * (coord3d.getElement(3 * electrostaticInteractionAtomPosition[k][0] + 2) - coord3d.getElement(3 * electrostaticInteractionAtomPosition[k][1] + 2));
 
435
                                                                        //logger.debug("OK: have d2 z");
 
436
                                                                        break;
 
437
                                                        }
 
438
                                                        
 
439
                                                        if (electrostaticInteractionAtomPosition[k][1] == atomNumberi) {
 
440
                                                                ddR2 = (-1) * ddR2;
 
441
                                                                //logger.debug("OK: d2 bond 1");
 
442
                                                        }
 
443
                                                        
 
444
                                                        ddR2 = ddR2 / Math.pow(r[k],2);
 
445
                                                        
 
446
                                                        // ddR[i][j][k]
 
447
                                                        ddR[i][j][k] = ddR1 - ddR2;
 
448
                                                } else {
 
449
                                                        ddR[i][j][k] = 0;
 
450
                                                        //logger.debug("OK: 0");
 
451
                                                }
 
452
                                        } else {
 
453
                                                ddR[i][j][k] = 0;
 
454
                                                //logger.debug("OK: 0");
 
455
                                        }
 
456
                                        //logger.debug("Electrostatic interactionn " + k + " : " + "ddR[" + i + "][" + j + "][" + k + "] = " + ddR[i][j][k]);
 
457
                                }
 
458
                        }
 
459
                }       
 
460
        }
 
461
 
 
462
 
 
463
        /**
 
464
         *  Get the internuclear separation second derivative respect to the cartesian coordinates of the atoms.
 
465
         *
 
466
         *@return        Bond lengths second derivative value [dimension(3xN)] [bonds Number]
 
467
         */
 
468
         public double[][][] getInternuclearSeparationSecondDerivative() {
 
469
                return ddR;
 
470
        }
 
471
 
 
472
 
 
473
        /**
 
474
         *  Evaluate the second order partial derivative (hessian) for the Electrostatic interaction energy given the atoms coordinates
 
475
         *
 
476
         *@param  coord3d  Current molecule coordinates.
 
477
         */
 
478
        public void setHessianMMFF94SumEQ(GVector coord3d) {
 
479
                
 
480
                forHessian = new double[coord3d.getSize() * coord3d.getSize()];
 
481
                
 
482
                setInternuclearSeparationSecondDerivative(coord3d);
 
483
                
 
484
                double sumHessianEQ;
 
485
                int forHessianIndex;
 
486
                for (int i = 0; i < coord3d.getSize(); i++) {
 
487
                        for (int j = 0; j < coord3d.getSize(); j++) {
 
488
                                sumHessianEQ = 0;
 
489
                                for (int k = 0; k < electrostaticInteractionNumber; k++) {
 
490
                                        sumHessianEQ = sumHessianEQ + (((332.0716 * qi[k] * qj[k] * n * (n+1) / D*Math.pow(r[k]+delta,n+2)) * dR[i][k]) * dR[j][k] 
 
491
                                                                        + (-332.0716 * qi[k] * qj[k] * n / D * Math.pow(r[k]+delta,n+1)) * ddR[i][j][k]);
 
492
                                }
 
493
                                forHessianIndex = i*coord3d.getSize()+j;
 
494
                                forHessian[forHessianIndex] = sumHessianEQ;
 
495
                                //logger.debug("forHessian[forHessianIndex] : " + forHessian[forHessianIndex]);
 
496
                        }
 
497
                }
 
498
 
 
499
                hessianMMFF94SumEQ = new GMatrix(coord3d.getSize(), coord3d.getSize(), forHessian);
 
500
                //logger.debug("hessianMMFF94SumEQ : " + hessianMMFF94SumEQ);
 
501
        }
 
502
 
 
503
 
 
504
        /**
 
505
         *  Get the hessian for the Electrostatic interaction energy.
 
506
         *
 
507
         *@return        Hessian value of the Electrostatic interaction term.
 
508
         */
 
509
        public GMatrix getHessianMMFF94SumEQ() {
 
510
                return hessianMMFF94SumEQ;
 
511
        }
 
512
 
 
513
 
 
514
        /**
 
515
         *  Get the hessian for the Electrostatic interaction energy.
 
516
         *
 
517
         *@return        Hessian value of the Electrostatic interaction term.
 
518
         */
 
519
        public double[] getForHessianMMFF94SumEQ() {
 
520
                return forHessian;
 
521
        }
 
522
 
 
523
}
 
524