3
* $Date: 2007-01-04 18:46:10 +0100 (Thu, 04 Jan 2007) $
6
* Copyright (C) 2001-2007 The Chemistry Development Kit (CDK) project
8
* Contact: cdk-devel@lists.sf.net
10
* This program is free software; you can redistribute it and/or
11
* modify it under the terms of the GNU Lesser General Public License
12
* as published by the Free Software Foundation; either version 2.1
13
* of the License, or (at your option) any later version.
14
* All we ask is that proper credit is given for our work, which includes
15
* - but is not limited to - adding the above copyright notice to the beginning
16
* of your source code files, and to any copyright notice that you may distribute
17
* with programs based on this work.
19
* This program is distributed in the hope that it will be useful,
20
* but WITHOUT ANY WARRANTY; without even the implied warranty of
21
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22
* GNU Lesser General Public License for more details.
24
* You should have received a copy of the GNU Lesser General Public License
25
* along with this program; if not, write to the Free Software
26
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
29
package org.openscience.cdk.math.qm;
31
import org.openscience.cdk.math.Matrix;
32
import org.openscience.cdk.math.Vector;
35
* This class contains the information to use gauss function as a base
36
* for calculation of quantum mechanics. The function is defined as:
38
* f(x,y,z) = (x-rx)^nx * (y-ry)^ny * (z-rz)^nz * exp(-alpha*(r-ri)^2)
42
* S = <phi_i|phi_j><br>
43
* J = <d/dr phi_i | d/dr phi_j><br>
44
* V = <phi_i | 1/r | phi_j><br>
46
* @author Stephan Michels <stephan@vern.chem.tu-berlin.de>
47
* @cdk.created 2001-06-14
49
* @cdk.keyword Gaussian basis set
51
public class GaussiansBasis implements IBasis
53
private int count; // number of the basis functions
54
private int[] nx; // [base]
55
private int[] ny; // [base]
56
private int[] nz; // [base]
57
private double[] alpha; // [base]
58
private double[] norm; // Normalize the bases
59
private Vector[] r; // [basis] Positions of base functions
61
private int count_atoms; // number of the atoms
62
private Vector[] rN; // [atom] Positions of the atoms
63
private int[] oz; // [atom] Atomic numbers of the atoms
66
private double minx = 0; private double maxx = 0;
67
private double miny = 0; private double maxy = 0;
68
private double minz = 0; private double maxz = 0;
70
public GaussiansBasis()
75
* Set up basis with gauss funktions
76
* f(x,y,z) = (x-rx)^nx * (y-ry)^ny * (z-rz)^nz * exp(-alpha*(r-ri)^2).
78
* @param atoms The atom will need to calculate the core potential
80
public GaussiansBasis(int[] nx, int[] ny, int[] nz, double[] alpha, Vector[] r, org.openscience.cdk.interfaces.IAtom[] atoms)
82
setBasis(nx, ny, nz, alpha, r, atoms);
86
* Set up basis with gauss funktions
87
* f(x,y,z) = (x-rx)^nx * (y-ry)^ny * (z-rz)^nz * exp(-alpha*(r-ri)^2).
89
* @param atoms The atom will need to calculate the core potential
91
protected void setBasis(int[] nx, int[] ny, int[] nz, double[] alpha, Vector[] r, org.openscience.cdk.interfaces.IAtom[] atoms)
95
//count_atoms = molecule.getSize();
96
count_atoms = atoms.length;
98
// logger.debug("Count of atoms: "+count_atoms);
100
this.rN = new Vector[count_atoms];
101
this.oz = new int[count_atoms];
102
for(i=0; i<count_atoms; i++)
104
this.rN[i] = (new Vector(atoms[i].getPoint3d())).mul(1.8897);
105
this.oz[i] = atoms[i].getAtomicNumber();
106
// logger.debug((i+1)+".Atom Z="+this.oz[i]+" r="+(new Vector(atoms[i].getPoint3d()))+"[angstrom]");
110
count = Math.min(nx.length,
112
Math.min(nz.length,alpha.length)));
114
// logger.debug("Count of bases: "+count);
116
this.nx = new int[count];
117
this.ny = new int[count];
118
this.nz = new int[count];
119
this.alpha = new double[count];
120
//this.atoms = new int[count];
121
this.norm = new double[count];
122
this.r = new Vector[count];
123
this.oz = new int[count];
125
for(i=0; i<count; i++)
130
this.alpha[i] = alpha[i];
131
//this.atoms[i] = atoms[i];
132
//this.r[i] = new Vector(atoms[i].getPoint3d()).mul(1.8897);
133
this.r[i] = r[i].mul(1.8897);
135
norm[i] = Math.sqrt(calcS(i,i));
141
// logger.debug((i+1)+".Base nx="+nx[i]+" ny="+ny[i]+" nz="+nz[i]+" alpha="+
142
// alpha[i]+" r="+r[i]+" norm="+norm[i]);
146
minx = Math.min(minx,this.r[i].vector[0]); maxx = Math.max(maxx,this.r[i].vector[0]);
147
miny = Math.min(miny,this.r[i].vector[1]); maxy = Math.max(maxy,this.r[i].vector[1]);
148
minz = Math.min(minz,this.r[i].vector[2]); maxz = Math.max(maxz,this.r[i].vector[2]);
152
minx = r[0].vector[0]; maxx = r[0].vector[0];
153
miny = r[0].vector[1]; maxy = r[0].vector[1];
154
minz = r[0].vector[2]; maxz = r[0].vector[2];
158
minx -= 2; maxx += 2;
159
miny -= 2; maxy += 2;
160
minz -= 2; maxz += 2;
166
* Gets the number of base vectors.
174
* Gets the dimension of the volume, which describes the base.
176
public double getMinX()
182
* Gets the dimension of the volume, which describes the base.
184
public double getMaxX()
190
* Gets the dimension of the volume, which describes the base.
192
public double getMinY()
198
* Gets the dimension of the volume, which describes the base.
200
public double getMaxY()
206
* Gets the dimension of the volume, which describes the base.
208
public double getMinZ()
214
* Gets the dimension of the volume, which describes the base.
216
public double getMaxZ()
222
* Calculates the function value an (x,y,z).
223
* @param index The number of the base
225
public double getValue(int index, double x, double y, double z)
227
double dx = (x*1.8897)-r[index].vector[0];
228
double dy = (y*1.8897)-r[index].vector[1];
229
double dz = (z*1.8897)-r[index].vector[2];
232
for(i=0; i<nx[index]; i++)
234
for(i=0; i<ny[index]; i++)
236
for(i=0; i<nz[index]; i++)
238
return result*Math.exp(-alpha[index]*(dx*dx+dy*dy+dz*dz));
242
* Calculates the function values.
243
* @param index The number of the base
245
public Vector getValues(int index, Matrix m)
250
double x,y,z,dx,dy,dz,mx,my,mz;
252
x = m.matrix[0][0]; y = m.matrix[1][0]; z = m.matrix[2][0];
254
dx = (x*1.8897)-r[index].vector[0];
255
dy = (y*1.8897)-r[index].vector[1];
256
dz = (z*1.8897)-r[index].vector[2];
258
Vector result = new Vector(m.columns);
262
for(i=0; i<nx[index]; i++)
266
for(i=0; i<ny[index]; i++)
270
for(i=0; i<nz[index]; i++)
273
dx *= dx; dy *= dy; dz *= dz;
275
result.vector[0] = mx*my*mz*Math.exp(-alpha[index]*(dx+dy+dz));
277
for(j=1; j<m.columns; j++)
279
if (x!=m.matrix[0][j])
282
dx = (x*1.8897)-r[index].vector[0];
284
for(i=0; i<nx[index]; i++)
289
if (y!=m.matrix[1][j])
292
dy = (y*1.8897)-r[index].vector[1];
294
for(i=0; i<ny[index]; i++)
299
if (z!=m.matrix[2][j])
302
dz = (z*1.8897)-r[index].vector[2];
304
for(i=0; i<nz[index]; i++)
309
result.vector[j] = mx*my*mz*Math.exp(-alpha[index]*(dx+dy+dz));
316
* Gets the position of a base.
318
public Vector getPosition(int index)
320
return r[index].duplicate().mul(0.52918);
323
public double calcD(double normi, double normj, double alphai, double alphaj, Vector ri, Vector rj)
325
double dx = ri.vector[0]-rj.vector[0];
326
double dy = ri.vector[1]-rj.vector[1];
327
double dz = ri.vector[2]-rj.vector[2];
328
return Math.exp(-((alphai*alphaj)/(alphai+alphaj))*(dx*dx+dy*dy+dz*dz))*normi*normj;
332
* Transfer equation for the calculation of the overlap integral..
334
private double calcI(int ni, int nj, double alphai, double alphaj, double xi, double xj)
336
if ((ni<0) || (nj<0))
338
System.err.println("Error [Basis.calcI()]: nj="+nj);
339
return Double.NaN; // Falls fehlerhafte Parameter
341
double[][] I = new double[nj+1][];
343
double alphaij = alphai+alphaj;
344
double xij = (alphai*xi+alphaj*xj)/alphaij;
346
I[0] = new double[(nj+ni+1)];
347
I[0][0] = Math.sqrt(Math.PI)/Math.sqrt(alphaij); // I(0,0)=G(0)
351
I[0][1] = -(xi-xij)*I[0][0]; // I(0,1)=G(1)
354
for(int i=2; i<=(nj+ni); i++)
355
I[0][i] = ((i-1)/(2*alphaij))*I[0][i-2]-(xi-xij)*I[0][i-1];
357
for(int j=1; j<=nj; j++)
359
I[j] = new double[(nj+ni+1)-j];
361
for(int i=0; i<=(nj+ni-j); i++)
362
I[j][i] = I[nj-1][ni+1]+(xi-xj)*I[nj-1][ni];
369
public double calcS(int i, int j)
371
//logger.debug("S: i="+i+" j="+j+" r[i]="+r[i]+" r[j]="+r[j]);
372
return calcD(norm[i], norm[j], alpha[i],alpha[j],r[i],r[j]) *
373
calcI(nx[i],nx[j],alpha[i],alpha[j],r[i].vector[0],r[j].vector[0]) *
374
calcI(ny[i],ny[j],alpha[i],alpha[j],r[i].vector[1],r[j].vector[1]) *
375
calcI(nz[i],nz[j],alpha[i],alpha[j],r[i].vector[2],r[j].vector[2]);
379
* Transfer equation the the calculation of the impulse
381
public double calcJ(int ni, int nj, double alphai, double alphaj, double xi, double xj)
384
return -4d*alphai*alphai* calcI(ni+2,nj,alphai,alphaj,xi,xj)
385
+2d*alphai*(2*ni+1)* calcI(ni ,nj,alphai,alphaj,xi,xj)
386
-ni*(ni-1)* calcI(ni-2,nj,alphai,alphaj,xi,xj);
388
return -4d*alphai*alphai* calcI(3 ,nj,alphai,alphaj,xi,xj)
389
+6d*alphai* calcI(1 ,nj,alphai,alphaj,xi,xj);
391
return -4d*alphai*alphai* calcI(2 ,nj,alphai,alphaj,xi,xj)
392
+2d*alphai* calcI(0 ,nj,alphai,alphaj,xi,xj);
394
System.err.println("Error [Basis.calcJ]: ni="+ni);
395
return Double.NaN; // Falls fehlerhafte Parameter
398
public double calcJ(int i, int j)
400
return calcD(norm[i], norm[j], alpha[i],alpha[j],r[i],r[j])*
401
(calcJ(nx[i],nx[j],alpha[i],alpha[j],r[i].vector[0],r[j].vector[0])*
402
calcI(ny[i],ny[j],alpha[i],alpha[j],r[i].vector[1],r[j].vector[1])*
403
calcI(nz[i],nz[j],alpha[i],alpha[j],r[i].vector[2],r[j].vector[2])+
405
calcI(nx[i],nx[j],alpha[i],alpha[j],r[i].vector[0],r[j].vector[0])*
406
calcJ(ny[i],ny[j],alpha[i],alpha[j],r[i].vector[1],r[j].vector[1])*
407
calcI(nz[i],nz[j],alpha[i],alpha[j],r[i].vector[2],r[j].vector[2])+
409
calcI(nx[i],nx[j],alpha[i],alpha[j],r[i].vector[0],r[j].vector[0])*
410
calcI(ny[i],ny[j],alpha[i],alpha[j],r[i].vector[1],r[j].vector[1])*
411
calcJ(nz[i],nz[j],alpha[i],alpha[j],r[i].vector[2],r[j].vector[2]));
415
* Transfer equation for the calculation of core potentials
417
public double calcG(int n, double t, double alphai, double alphaj, double xi, double xj, double xN)
420
return ((n-1)/(2*(alphai+alphaj)))* calcG(n-2, t, alphai, alphaj, xi, xj, xN)
421
-(n-1)*t*t* calcG(n-2, t, alphai, alphaj, xi, xj, xN)
422
+(((alphai*xi+alphaj*xj)/(alphai+alphaj))-xi)* calcG(n-1, t, alphai, alphaj, xi, xj, xN)
423
+(((alphai*xi+alphaj*xj)/(alphai+alphaj))-xN)*t*t* calcG(n-1, t, alphai, alphaj, xi, xj, xN);
426
return (((alphai*xi+alphaj*xj)/(alphai+alphaj))-xi)* calcG(0, t, alphai, alphaj, xi, xj, xN)
427
+(((alphai*xi+alphaj*xj)/(alphai+alphaj))-xN)*t*t* calcG(0, t, alphai, alphaj, xi, xj, xN);
430
return Math.sqrt(Math.PI)/Math.sqrt(alphai+alphaj);
432
System.err.println("Error [Basis.calcG]: n="+n);
433
return Double.NaN; // Falls fehlerhafte Parameter
437
* Transfer equation for the calculation of core potentials.
439
private double calcI(int ni, int nj, double t, double alphai, double alphaj,
440
double xi, double xj, double xN)
443
return calcI(ni+1, nj-1, t, alphai, alphaj, xi, xj, xN)+
444
(xi-xj)*calcI(ni, nj-1, t, alphai, alphaj, xi, xj, xN);
447
return calcG(ni, t, alphai, alphaj, xi, xj, xN);
449
System.err.println("Error [Basis.calcI()]: nj="+nj);
450
return Double.NaN; // Falls fehlerhafte Parameter
454
* Calculates the core potential.
455
* It use a 10 point Simpson formula.
457
* @param i Index of the first base
458
* @param j Index of the second base
459
* @param rN Position the core potential
460
* @param Z Atomic number of the nucleous
462
public double calcV(int i, int j, Vector rN, double Z)
472
double alphaij = alpha[i]+alpha[j];
474
double rxij = (alpha[i]*r[i].vector[0]+alpha[j]*r[j].vector[0])/alphaij;
475
double ryij = (alpha[i]*r[i].vector[1]+alpha[j]*r[j].vector[1])/alphaij;
476
double rzij = (alpha[i]*r[i].vector[2]+alpha[j]*r[j].vector[2])/alphaij;
478
double X = alphaij*((rxij-rN.vector[0])*(rxij-rN.vector[0]) +
479
(ryij-rN.vector[1])*(ryij-rN.vector[1]) +
480
(rzij-rN.vector[2])*(rzij-rN.vector[2]));
482
double C = 2*calcD(norm[i], norm[j],
483
alpha[i], alpha[j], r[i], r[j])*Math.sqrt(alphaij)/Math.sqrt(Math.PI);
486
for(f = 1; f<steps; f=f+2)
489
sum1 += Math.exp(-X*t*t)*calcI(nx[i], nx[j], t, alpha[i], alpha[j],
490
r[i].vector[0], r[j].vector[0], rN.vector[0]) *
491
calcI(ny[i], ny[j], t, alpha[i], alpha[j],
492
r[i].vector[1], r[j].vector[1], rN.vector[1]) *
493
calcI(nz[i], nz[j], t, alpha[i], alpha[j],
494
r[i].vector[2], r[j].vector[2], rN.vector[2]);
498
for(f = 2; f<steps; f=f+2)
501
sum2 += Math.exp(-X*t*t)*calcI(nx[i], nx[j], t, alpha[i], alpha[j],
502
r[i].vector[0], r[j].vector[0], rN.vector[0]) *
503
calcI(ny[i], ny[j], t, alpha[i], alpha[j],
504
r[i].vector[1], r[j].vector[1], rN.vector[1]) *
505
calcI(nz[i], nz[j], t, alpha[i], alpha[j],
506
r[i].vector[2], r[j].vector[2], rN.vector[2]);
510
f1 = Math.exp(-X*t*t)*calcI(nx[i], nx[j], t, alpha[i], alpha[j],
511
r[i].vector[0], r[j].vector[0], rN.vector[0]) *
512
calcI(ny[i], ny[j], t, alpha[i], alpha[j],
513
r[i].vector[1], r[j].vector[1], rN.vector[1]) *
514
calcI(nz[i], nz[j], t, alpha[i], alpha[j],
515
r[i].vector[2], r[j].vector[2], rN.vector[2]);
518
f2 = Math.exp(-X*t*t)*calcI(nx[i], nx[j], t, alpha[i], alpha[j],
519
r[i].vector[0], r[j].vector[0], rN.vector[0]) *
520
calcI(ny[i], ny[j], t, alpha[i], alpha[j],
521
r[i].vector[1], r[j].vector[1], rN.vector[1]) *
522
calcI(nz[i], nz[j], t, alpha[i], alpha[j],
523
r[i].vector[2], r[j].vector[2], rN.vector[2]);
525
return (h/3)*(f1 + 4*sum1 + 2*sum2 + f2)*Z*C;
529
* Calculates the core potential.
530
* It use a 10 point Simpson formula.
532
* @param i Index of the first base
533
* @param j Index of the second base
535
public double calcV(int i, int j)
538
for(int k=0; k<count_atoms; k++)
540
//logger.debug("k="+k+" r="+r[k]);
541
result += calcV(i,j, rN[k], oz[k]);
543
return -result; // Vorsicht negatives Vorzeichen
547
* Transfer equation for a four center integral.
549
public double calcG(int n, int m, double u,
550
double alphai, double alphaj, double alphak, double alphal, double xi, double xj, double xk, double xl)
554
// logger.debug("Error(CalcG):Bad parameter n="+n+" m="+m);
558
double alphaij = alphai+alphaj;
559
double alphakl = alphak+alphal;
561
double xij = (alphai*xi+alphaj*xj)/alphaij;
562
double xkl = (alphak*xk+alphal*xl)/alphakl;
564
double C00 = (xij-xi)-((u*u*alphakl*(xij-xkl))/(u*u*(alphaij+alphakl)+alphaij*alphakl));
565
double Cs00 = (xkl-xk)+((u*u*alphaij*(xij-xkl))/(u*u*(alphaij+alphakl)+alphaij*alphakl));
566
double B00 = u*u/(2*(u*u*(alphaij+alphakl)+alphaij*alphakl));
567
double B10 = (u*u+alphakl)/(2*(u*u*(alphaij+alphakl)+alphaij*alphakl));
568
double Bs01 = (u*u+alphaij)/(2*(u*u*(alphaij+alphakl)+alphaij*alphakl));
570
double[][] G = new double[n+1][m+1];
582
G[i][0] = (i-1)*B10 *G[i-2][0]+
591
G[0][i] = (i-1)*Bs01 *G[0][i-2]+
597
G[1][i] = i*B00 *G[0][i-1]+
603
G[i][j] = (i-1)*B10 *G[i-2][j]+
611
* Transfer equation for a four center integral.
613
public double calcI(int ni, int nj, int nk, int nl, double u,
614
double alphai, double alphaj, double alphak, double alphal,
615
double xi, double xj, double xk, double xl)
618
return calcI(ni+1,nj-1,nk ,nl ,u,alphai,alphaj,alphak,alphal,xi,xj,xk,xl)+
619
(xj-xi)*calcI(ni ,nj-1,nk ,nl ,u,alphai,alphaj,alphak,alphal,xi,xj,xk,xl);
622
return calcI(ni ,nj ,nk+1,nl-1,u,alphai,alphaj,alphak,alphal,xi,xj,xk,xl)+
623
(xl-xk)*calcI(ni ,nj ,nk ,nl-1,u,alphai,alphaj,alphak,alphal,xi,xj,xk,xl);
625
if ((ni==0) && (nj==0) && (nk==0) && (nl==0))
628
if ((nj==0) && (nl==0))
629
return calcG(ni,nk,u,alphai,alphaj,alphak,alphal,xi,xj,xk,xl);
631
// logger.debug("Error(CalcI):Bad parameter ni="+ni+" nj="+nj+" nk="+nk+" nl="+nl);
635
public double calcI(int i, int j, int k, int l)
642
// Berechnen der Integration nach Simson
647
double alphaij = alpha[i]+alpha[j];
648
double alphakl = alpha[k]+alpha[l];
650
double rxij = (alpha[i]*r[i].vector[0]+alpha[j]*r[j].vector[0])/alphaij;
651
double ryij = (alpha[i]*r[i].vector[1]+alpha[j]*r[j].vector[1])/alphaij;
652
double rzij = (alpha[i]*r[i].vector[2]+alpha[j]*r[j].vector[2])/alphaij;
654
double rxkl = (alpha[k]*r[k].vector[0]+alpha[l]*r[l].vector[0])/alphakl;
655
double rykl = (alpha[k]*r[k].vector[1]+alpha[l]*r[l].vector[1])/alphakl;
656
double rzkl = (alpha[k]*r[k].vector[2]+alpha[l]*r[l].vector[2])/alphakl;
658
double alpha0 = alphaij*alphakl/(alphaij+alphakl);
660
double X = alpha0*((rxij-rxkl)*(rxij-rxkl) +
661
(ryij-rykl)*(ryij-rzkl) +
662
(rzij-rzkl)*(rzij-rzkl));
664
double C = (Math.PI*Math.PI*Math.PI/Math.pow((alpha[i]+alpha[j])*(alpha[k]+alpha[l]),1.5))*
666
calcD(norm[i], norm[j], alpha[i], alpha[j], r[i], r[j])*
667
calcD(norm[k], norm[l], alpha[k], alpha[l], r[k], r[l])*
668
(2d/Math.sqrt(Math.PI));
671
for(f = 1; f<steps; f=f+2)
674
sum1 += Math.exp(-X*t*t)*calcI(nx[i], nx[j], nx[k], nx[l], t, alpha[i], alpha[j], alpha[k], alpha[l],
675
r[i].vector[0], r[j].vector[0], r[k].vector[0], r[l].vector[0]) *
676
calcI(ny[i], ny[j], ny[k], ny[l], t, alpha[i], alpha[j], alpha[k], alpha[l],
677
r[i].vector[1], r[j].vector[1], r[k].vector[1], r[l].vector[1]) *
678
calcI(nz[i], nz[j], nz[k], nz[l], t, alpha[i], alpha[j], alpha[k], alpha[l],
679
r[i].vector[2], r[j].vector[2], r[k].vector[2], r[l].vector[2]);
684
for(f = 2; f<steps; f=f+2)
687
sum2 += Math.exp(-X*t*t)*calcI(nx[i], nx[j], nx[k], nx[l], t, alpha[i], alpha[j], alpha[k], alpha[l],
688
r[i].vector[0], r[j].vector[0], r[k].vector[0], r[l].vector[0]) *
689
calcI(ny[i], ny[j], ny[k], ny[l], t, alpha[i], alpha[j], alpha[k], alpha[l],
690
r[i].vector[1], r[j].vector[1], r[k].vector[1], r[l].vector[1]) *
691
calcI(nz[i], nz[j], nz[k], nz[l], t, alpha[i], alpha[j], alpha[k], alpha[l],
692
r[i].vector[2], r[j].vector[2], r[k].vector[2], r[l].vector[2]);
696
f1 = Math.exp(-X*t*t)*calcI(nx[i], nx[j], nx[k], nx[l], t, alpha[i], alpha[j], alpha[k], alpha[l],
697
r[i].vector[0], r[j].vector[0], r[k].vector[0], r[l].vector[0]) *
698
calcI(ny[i], ny[j], ny[k], ny[l], t, alpha[i], alpha[j], alpha[k], alpha[l],
699
r[i].vector[1], r[j].vector[1], r[k].vector[1], r[l].vector[1]) *
700
calcI(nz[i], nz[j], nz[k], nz[l], t, alpha[i], alpha[j], alpha[k], alpha[l],
701
r[i].vector[2], r[j].vector[2], r[k].vector[2], r[l].vector[2]);
704
f2 = Math.exp(-X*t*t)*calcI(nx[i], nx[j], nx[k], nx[l], t, alpha[i], alpha[j], alpha[k], alpha[l],
705
r[i].vector[0], r[j].vector[0], r[k].vector[0], r[l].vector[0]) *
706
calcI(ny[i], ny[j], ny[k], ny[l], t, alpha[i], alpha[j], alpha[k], alpha[l],
707
r[i].vector[1], r[j].vector[1], r[k].vector[1], r[l].vector[1]) *
708
calcI(nz[i], nz[j], nz[k], nz[l], t, alpha[i], alpha[j], alpha[k], alpha[l],
709
r[i].vector[2], r[j].vector[2], r[k].vector[2], r[l].vector[2]);
711
return C * (h/3)*(f1 + 4*sum1 + 2*sum2 + f2);
added
removed
src/org/openscience/cdk/math/qm/GaussiansBasis.java
