← Back to branch summary

~ubuntu-branches/ubuntu/quantal/psicode/quantal

« back to all changes in this revision

Viewing changes to src/lib/libmints/eri.cc

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2009-02-23 00:12:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20090223001202-rutldoy3dimfpesc
Tags: 3.4.0-1
* New upstream release.

[ Michael Banck ]
* debian/patches/01_DESTDIR.dpatch: Refreshed.
* debian/patches/02_FHS.dpatch: Removed, applied upstream.
* debian/patches/03_debian_docdir: Likewise.
* debian/patches/04_man.dpatch: Likewise.
* debian/patches/06_466828_fix_gcc_43_ftbfs.dpatch: Likewise.
* debian/patches/07_464867_move_executables: Fixed and refreshed.
* debian/patches/00list: Adjusted.
* debian/control: Improved description.
* debian/patches-held: Removed.
* debian/rules (install/psi3): Do not ship the ruby bindings for now.

[ Daniel Leidert ]
* debian/rules: Fix txtdir via DEB_MAKE_INSTALL_TARGET.
* debian/patches/01_DESTDIR.dpatch: Refreshed.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/lib/libmints/eri.cc
 
1
#include <stdexcept>
 
2
#include <string>
 
3
#include <libciomr/libciomr.h>
 
4
#include <libqt/qt.h>
 
5
 
 
6
#include <libmints/basisset.h>
 
7
#include <libmints/gshell.h>
 
8
#include <libmints/overlap.h>
 
9
#include <libmints/eri.h>
 
10
#include <libmints//wavefunction.h>
 
11
#include <physconst.h>
 
12
 
 
13
#define MAX(a, b) ((a) > (b) ? (a) : (b))
 
14
#define EPS 1.0e-17
 
15
#define TABLESIZE 121
 
16
 
 
17
using namespace psi;
 
18
 
 
19
void calc_f(double *F, int n, double t)
 
20
{
 
21
    int i, m;
 
22
    int m2;
 
23
    double t2;
 
24
    double num;
 
25
    double sum;
 
26
    double term1;
 
27
    static double K = 1.0/M_2_SQRTPI;
 
28
    double et;
 
29
 
 
30
    if (t>20.0){
 
31
        t2 = 2*t;
 
32
        et = exp(-t);
 
33
        t = sqrt(t);
 
34
        F[0] = K*erf(t)/t;
 
35
        for(m=0; m<=n-1; m++){
 
36
            F[m+1] = ((2*m + 1)*F[m] - et)/(t2);
 
37
        }
 
38
    }
 
39
    else {
 
40
        et = exp(-t);
 
41
        t2 = 2*t;
 
42
        m2 = 2*n;
 
43
        num = df[m2];
 
44
        i=0;
 
45
        sum = 1.0/(m2+1);
 
46
        do{
 
47
            i++;
 
48
            num = num*t2;
 
49
            term1 = num/df[m2+2*i+2];
 
50
            sum += term1;
 
51
        } while (fabs(term1) > EPS && i < MAX_FAC);
 
52
        F[n] = sum*et;
 
53
        for(m=n-1;m>=0;m--){
 
54
            F[m] = (t2*F[m+1] + et)/(2*m+1);
 
55
        }
 
56
    }
 
57
}
 
58
 
 
59
ERI::ERI(IntegralFactory* integral, Ref<BasisSet> &bs1, Ref<BasisSet> &bs2, Ref<BasisSet> &bs3, Ref<BasisSet> &bs4, int deriv) 
 
60
    : TwoBodyInt(integral, bs1, bs2, bs3, bs4, deriv)
 
61
{
 
62
    // Initialize libint static data
 
63
    init_libint_base();
 
64
    if (deriv_)
 
65
        init_libderiv_base();
 
66
    
 
67
    // Figure out some information to initialize libint/libderiv with
 
68
    // 1. Maximum angular momentum
 
69
    int max_am = MAX(MAX(bs1_->max_am(), bs2_->max_am()), MAX(bs3_->max_am(), bs4_->max_am()));
 
70
    // 2. Maximum number of primitive combinations
 
71
    int max_nprim = bs1_->max_nprimitive() * bs2_->max_nprimitive() * bs3_->max_nprimitive() * bs4_->max_nprimitive();
 
72
    // 3. Maximum Cartesian class size
 
73
    max_cart_ = ioff[bs1_->max_am()] * ioff[bs2_->max_am()] * ioff[bs3_->max_am()] * ioff[bs4_->max_am()] +1;
 
74
    
 
75
    // Make sure libint is compiled to handle our max AM
 
76
    if (max_am >= LIBINT_MAX_AM) {
 
77
        fprintf(stderr, "ERROR: ERI - libint cannot handle angular momentum this high.\n       Recompile libint for higher angular momentum, then recompile this program.\n");
 
78
        abort();
 
79
    }
 
80
    if (deriv_ == 1 && max_am >= LIBDERIV_MAX_AM1) {
 
81
        fprintf(stderr, "ERROR: ERI - libderiv cannot handle angular momentum this high.\n     Recompile libderiv for higher angular momentum, then recompile this program.\n");
 
82
        abort();
 
83
    }
 
84
    
 
85
    // Initialize libint
 
86
    init_libint(&libint_, max_am, max_nprim);   
 
87
    // and libderiv, if needed
 
88
    if (deriv_)
 
89
        init_libderiv1(&libderiv_, max_am, max_nprim, max_cart_-1);
 
90
        
 
91
    size_t size = INT_NCART(bs1_->max_am()) * INT_NCART(bs2_->max_am()) *
 
92
                  INT_NCART(bs3_->max_am()) * INT_NCART(bs4_->max_am());
 
93
                  
 
94
    // Used in pure_transform
 
95
    tformbuf_ = new double[size];
 
96
    memset(tformbuf_, 0, sizeof(double)*size);
 
97
    
 
98
    if (deriv_ == 1)
 
99
        size *= 3*natom_;
 
100
 
 
101
    target_ = new double[size];
 
102
    memset(target_, 0, sizeof(double)*size);
 
103
    
 
104
    source_ = new double[size];
 
105
    memset(source_, 0, sizeof(double)*size);
 
106
 
 
107
    init_fjt(4*max_am + DERIV_LVL);
 
108
}
 
109
 
 
110
ERI::~ERI()
 
111
{
 
112
    delete[] tformbuf_;
 
113
    delete[] target_;
 
114
    delete[] source_;
 
115
    delete[] denom_;
 
116
    free_block(d_);
 
117
    free_libint(&libint_);
 
118
    if (deriv_)
 
119
        free_libderiv(&libderiv_);
 
120
}
 
121
 
 
122
// Taken from CINTS
 
123
void ERI::init_fjt(int max)
 
124
{
 
125
    int i,j;
 
126
    double denom,d2jmax1,r2jmax1,wval,d2wval,sum,term,rexpw;
 
127
 
 
128
    int n1 = max+7;
 
129
    int n2 = TABLESIZE;
 
130
    d_ = block_matrix(n1, n2);
 
131
 
 
132
    /* Tabulate the gamma function for t(=wval)=0.0. */
 
133
    denom = 1.0;
 
134
    for (i=0; i<n1; i++) {
 
135
        d_[i][0] = 1.0/denom;
 
136
        denom += 2.0;
 
137
    }
 
138
 
 
139
    /* Tabulate the gamma function from t(=wval)=0.1, to 12.0. */
 
140
    d2jmax1 = 2.0*(n1-1) + 1.0;
 
141
    r2jmax1 = 1.0/d2jmax1;
 
142
    for (i=1; i<TABLESIZE; i++) {
 
143
        wval = 0.1 * i;
 
144
        d2wval = 2.0 * wval;
 
145
        term = r2jmax1;
 
146
        sum = term;
 
147
        denom = d2jmax1;
 
148
        for (j=2; j<=200; j++) {
 
149
            denom = denom + 2.0;
 
150
            term = term * d2wval / denom;
 
151
            sum = sum + term;
 
152
            if (term <= 1.0e-15) break;
 
153
        }
 
154
        rexpw = exp(-wval);
 
155
 
 
156
      /* Fill in the values for the highest j gtable entries (top row). */
 
157
        d_[n1-1][i] = rexpw * sum;
 
158
 
 
159
      /* Work down the table filling in the rest of the column. */
 
160
        denom = d2jmax1;
 
161
        for (j=n1 - 2; j>=0; j--) {
 
162
            denom = denom - 2.0;
 
163
            d_[j][i] = (d_[j+1][i]*d2wval + rexpw)/denom;
 
164
        }
 
165
    }
 
166
 
 
167
    /* Form some denominators, so divisions can be eliminated below. */
 
168
    denom_ = new double[max+1];
 
169
    denom_[0] = 0.0;
 
170
    for (i=1; i<=max; i++) {
 
171
        denom_[i] = 1.0/(2*i - 1);
 
172
    }
 
173
 
 
174
    wval_infinity_ = 2*max + 37.0;
 
175
    itable_infinity_ = (int) (10 * wval_infinity_);
 
176
}
 
177
 
 
178
void ERI::int_fjt(double *F, int J, double wval)
 
179
{
 
180
    const double sqrpih =  0.886226925452758;
 
181
    const double coef2 =  0.5000000000000000;
 
182
    const double coef3 = -0.1666666666666667;
 
183
    const double coef4 =  0.0416666666666667;
 
184
    const double coef5 = -0.0083333333333333;
 
185
    const double coef6 =  0.0013888888888889;
 
186
    const double gfac30 =  0.4999489092;
 
187
    const double gfac31 = -0.2473631686;
 
188
    const double gfac32 =  0.321180909;
 
189
    const double gfac33 = -0.3811559346;
 
190
    const double gfac20 =  0.4998436875;
 
191
    const double gfac21 = -0.24249438;
 
192
    const double gfac22 =  0.24642845;
 
193
    const double gfac10 =  0.499093162;
 
194
    const double gfac11 = -0.2152832;
 
195
    const double gfac00 = -0.490;
 
196
 
 
197
    double wdif, d2wal, rexpw, /* denom, */ gval, factor, rwval, term;
 
198
    int i, itable, irange;
 
199
 
 
200
    /* Compute an index into the table. */
 
201
    /* The test is needed to avoid floating point exceptions for
 
202
    * large values of wval. */
 
203
    if (wval > wval_infinity_) {
 
204
        itable = itable_infinity_;
 
205
    }
 
206
    else {
 
207
        itable = (int) (10.0 * wval);
 
208
    }
 
209
 
 
210
    /* If itable is small enough use the table to compute int_fjttable. */
 
211
    if (itable < TABLESIZE) {
 
212
 
 
213
        wdif = wval - 0.1 * itable;
 
214
 
 
215
      /* Compute fjt for J. */
 
216
        F[J] = (((((coef6 * d_[J+6][itable]*wdif
 
217
            + coef5 * d_[J+5][itable])*wdif
 
218
            + coef4 * d_[J+4][itable])*wdif
 
219
            + coef3 * d_[J+3][itable])*wdif
 
220
            + coef2 * d_[J+2][itable])*wdif
 
221
            -  d_[J+1][itable])*wdif
 
222
            +  d_[J][itable];
 
223
 
 
224
      /* Compute the rest of the fjt. */
 
225
        d2wal = 2.0 * wval;
 
226
        rexpw = exp(-wval);
 
227
      /* denom = 2*J + 1; */
 
228
        for (i=J; i>0; i--) {
 
229
        /* denom = denom - 2.0; */
 
230
            F[i-1] = (d2wal*F[i] + rexpw)*denom_[i];
 
231
        }
 
232
    }
 
233
    /* If wval <= 2*J + 36.0, use the following formula. */
 
234
    else if (itable <= 20*J + 360) {
 
235
        rwval = 1.0/wval;
 
236
        rexpw = exp(-wval);
 
237
 
 
238
      /* Subdivide wval into 6 ranges. */
 
239
        irange = itable/30 - 3;
 
240
        if (irange == 1) {
 
241
            gval = gfac30 + rwval*(gfac31 + rwval*(gfac32 + rwval*gfac33));
 
242
            F[0] = sqrpih*sqrt(rwval) - rexpw*gval*rwval;
 
243
        }
 
244
        else if (irange == 2) {
 
245
            gval = gfac20 + rwval*(gfac21 + rwval*gfac22);
 
246
            F[0] = sqrpih*sqrt(rwval) - rexpw*gval*rwval;
 
247
        }
 
248
        else if (irange == 3 || irange == 4) {
 
249
            gval = gfac10 + rwval*gfac11;
 
250
            F[0] = sqrpih*sqrt(rwval) - rexpw*gval*rwval;
 
251
        }
 
252
        else if (irange == 5 || irange == 6) {
 
253
            gval = gfac00;
 
254
            F[0] = sqrpih*sqrt(rwval) - rexpw*gval*rwval;
 
255
        }
 
256
        else {
 
257
            F[0] = sqrpih*sqrt(rwval);
 
258
        }
 
259
 
 
260
      /* Compute the rest of the int_fjttable from table->d[0]. */
 
261
        factor = 0.5 * rwval;
 
262
        term = factor * rexpw;
 
263
        for (i=1; i<=J; i++) {
 
264
            F[i] = factor * F[i-1] - term;
 
265
            factor = rwval + factor;
 
266
        }
 
267
    }
 
268
    /* For large values of wval use this algorithm: */
 
269
    else {
 
270
        rwval = 1.0/wval;
 
271
        F[0] = sqrpih*sqrt(rwval);
 
272
        factor = 0.5 * rwval;
 
273
        for (i=1; i<=J; i++) {
 
274
            F[i] = factor * F[i-1];
 
275
            factor = rwval + factor;
 
276
        }
 
277
    }
 
278
}
 
279
 
 
280
void ERI::compute_shell(int sh1, int sh2, int sh3, int sh4)
 
281
{
 
282
    // Need to ensure the ordering asked by the user is valid for libint
 
283
    // compute_quartet does NOT check this. SEGFAULTS should occur if order
 
284
    // is not guaranteed.
 
285
    int s1, s2, s3, s4;
 
286
    int am1, am2, am3, am4, temp;
 
287
    bool p13p24 = false, p12 = false, p34 = false;
 
288
    
 
289
    // AM used for ordering
 
290
    am1 = bs1_->shell(sh1)->am(0);
 
291
    am2 = bs2_->shell(sh2)->am(0);
 
292
    am3 = bs3_->shell(sh3)->am(0);
 
293
    am4 = bs4_->shell(sh4)->am(0);
 
294
    
 
295
    int n1 = bs1_->shell(sh1)->nfunction(0);
 
296
    int n2 = bs1_->shell(sh2)->nfunction(0);
 
297
    int n3 = bs1_->shell(sh3)->nfunction(0);
 
298
    int n4 = bs1_->shell(sh4)->nfunction(0);
 
299
    
 
300
    // l(a) >= l(b), l(c) >= l(d), and l(c) + l(d) >= l(a) + l(b).
 
301
    if (am1 >= am2) {
 
302
        s1 = sh1;
 
303
        s2 = sh2;
 
304
    } else {
 
305
        s1 = sh2;
 
306
        s2 = sh1;
 
307
        p12 = true;
 
308
    }
 
309
    
 
310
    if (am3 >= am4) {
 
311
        s3 = sh3;
 
312
        s4 = sh4;
 
313
    } else {
 
314
        s3 = sh4;
 
315
        s4 = sh3;  
 
316
        p34 = true;
 
317
    }
 
318
    
 
319
    if ((am1 + am2) > (am3 + am4)) {
 
320
        // Swap s1 and s2 with s3 and s4
 
321
        temp = s1;
 
322
        s1 = s3;
 
323
        s3 = temp;
 
324
        
 
325
        temp = s2;
 
326
        s2 = s4;
 
327
        s4 = temp;
 
328
        p13p24 = true;
 
329
    }
 
330
    
 
331
    // s1, s2, s3, s4 contain the shells to do in libint order
 
332
    compute_quartet(s1, s2, s3, s4);
 
333
        
 
334
    // Permute integrals back, if needed
 
335
    if (p12 || p34 || p13p24)
 
336
        permute_target(source_, target_, s1, s2, s3, s4, p12, p34, p13p24);
 
337
    else {
 
338
        // copy the integrals to the target_
 
339
        memcpy(target_, source_, n1 * n2 * n3 * n4 *sizeof(double));
 
340
    }
 
341
}
 
342
 
 
343
void ERI::compute_quartet(int sh1, int sh2, int sh3, int sh4)
 
344
{
 
345
    Ref<GaussianShell> s1, s2, s3, s4;
 
346
    
 
347
    s1 = bs1_->shell(sh1);
 
348
    s2 = bs2_->shell(sh2);
 
349
    s3 = bs3_->shell(sh3);
 
350
    s4 = bs4_->shell(sh4);
 
351
    
 
352
    int am1 = s1->am(0);
 
353
    int am2 = s2->am(0);
 
354
    int am3 = s3->am(0);
 
355
    int am4 = s4->am(0);
 
356
    int am = am1 + am2 + am3 + am4; // total am
 
357
    int nprim1 = s1->nprimitive();
 
358
    int nprim2 = s2->nprimitive();
 
359
    int nprim3 = s3->nprimitive();
 
360
    int nprim4 = s4->nprimitive();
 
361
    size_t nprim, nprim_combination = nprim1 * nprim2 * nprim3 * nprim4;
 
362
    double A[3], B[3], C[3], D[3];
 
363
    A[0] = s1->center()[0];
 
364
    A[1] = s1->center()[1];
 
365
    A[2] = s1->center()[2];
 
366
    B[0] = s2->center()[0];
 
367
    B[1] = s2->center()[1];
 
368
    B[2] = s2->center()[2];
 
369
    C[0] = s3->center()[0];
 
370
    C[1] = s3->center()[1];
 
371
    C[2] = s3->center()[2];
 
372
    D[0] = s4->center()[0];
 
373
    D[1] = s4->center()[1];
 
374
    D[2] = s4->center()[2];
 
375
    
 
376
    // compute intermediates
 
377
    double AB2 = 0.0;
 
378
    AB2 += (A[0] - B[0]) * (A[0] - B[0]);
 
379
    AB2 += (A[1] - B[1]) * (A[1] - B[1]);
 
380
    AB2 += (A[2] - B[2]) * (A[2] - B[2]);
 
381
    double CD2 = 0.0;
 
382
    CD2 += (C[0] - D[0]) * (C[0] - D[0]);
 
383
    CD2 += (C[1] - D[1]) * (C[1] - D[1]);
 
384
    CD2 += (C[2] - D[2]) * (C[2] - D[2]);
 
385
    
 
386
    libint_.AB[0] = A[0] - B[0];
 
387
    libint_.AB[1] = A[1] - B[1];
 
388
    libint_.AB[2] = A[2] - B[2];
 
389
    libint_.CD[0] = C[0] - D[0];
 
390
    libint_.CD[1] = C[1] - D[1];
 
391
    libint_.CD[2] = C[2] - D[2];
 
392
        
 
393
    // Prepare all the data needed by libint
 
394
    nprim = 0;
 
395
    for (int p1=0; p1<nprim1; ++p1) {
 
396
        double a1 = s1->exp(p1);
 
397
        double c1 = s1->coef(0, p1);
 
398
        for (int p2=0; p2<nprim2; ++p2) {
 
399
            double a2 = s2->exp(p2);
 
400
            double c2 = s2->coef(0, p2);
 
401
            double zeta = a1 + a2;
 
402
            double ooz = 1.0/zeta;
 
403
            double oo2z = 1.0/(2.0 * zeta);
 
404
            
 
405
            double PA[3], PB[3];
 
406
            double P[3];
 
407
            
 
408
            P[0] = (a1*A[0] + a2*B[0])*ooz;
 
409
            P[1] = (a1*A[1] + a2*B[1])*ooz;
 
410
            P[2] = (a1*A[2] + a2*B[2])*ooz;
 
411
            PA[0] = P[0] - A[0];
 
412
            PA[1] = P[1] - A[1];
 
413
            PA[2] = P[2] - A[2];
 
414
            PB[0] = P[0] - B[0];
 
415
            PB[1] = P[1] - B[1];
 
416
            PB[2] = P[2] - B[2];
 
417
            
 
418
            double Sab = pow(M_PI*ooz, 3.0/2.0) * exp(-a1*a2*ooz*AB2) * c1 * c2;
 
419
            
 
420
            for (int p3=0; p3<nprim3; ++p3) {
 
421
                double a3 = s3->exp(p3);
 
422
                double c3 = s3->coef(0, p3);
 
423
                for (int p4=0; p4<nprim4; ++p4) {
 
424
                    double a4 = s4->exp(p4);
 
425
                    double c4 = s4->coef(0, p4);
 
426
                    double nu = a3 + a4;
 
427
                    double oon = 1.0/nu;
 
428
                    double oo2n = 1.0/(2.0*nu);
 
429
                    double oo2zn = 1.0/(2.0*(zeta+nu));
 
430
                    double rho = (zeta*nu)/(zeta+nu);
 
431
                    double oo2rho = 1 / (2.0*rho);
 
432
                    
 
433
                    double QC[3], QD[3], WP[3], WQ[3], PQ[3];
 
434
                    double Q[3], W[3];
 
435
 
 
436
                    Q[0] = (a3*C[0] + a4*D[0])*oon;
 
437
                    Q[1] = (a3*C[1] + a4*D[1])*oon;
 
438
                    Q[2] = (a3*C[2] + a4*D[2])*oon;
 
439
                    QC[0] = Q[0] - C[0];
 
440
                    QC[1] = Q[1] - C[1];
 
441
                    QC[2] = Q[2] - C[2];
 
442
                    QD[0] = Q[0] - D[0];
 
443
                    QD[1] = Q[1] - D[1];
 
444
                    QD[2] = Q[2] - D[2];
 
445
                    PQ[0] = P[0] - Q[0];
 
446
                    PQ[1] = P[1] - Q[1];
 
447
                    PQ[2] = P[2] - Q[2];
 
448
                    
 
449
                    double PQ2 = 0.0;
 
450
                    PQ2 += (P[0] - Q[0]) * (P[0] - Q[0]);
 
451
                    PQ2 += (P[1] - Q[1]) * (P[1] - Q[1]);
 
452
                    PQ2 += (P[2] - Q[2]) * (P[2] - Q[2]);
 
453
                    
 
454
                    W[0] = (zeta*P[0] + nu*Q[0]) / (zeta + nu);
 
455
                    W[1] = (zeta*P[1] + nu*Q[1]) / (zeta + nu);
 
456
                    W[2] = (zeta*P[2] + nu*Q[2]) / (zeta + nu);
 
457
                    WP[0] = W[0] - P[0];
 
458
                    WP[1] = W[1] - P[1];
 
459
                    WP[2] = W[2] - P[2];
 
460
                    WQ[0] = W[0] - Q[0];
 
461
                    WQ[1] = W[1] - Q[1];
 
462
                    WQ[2] = W[2] - Q[2];
 
463
                    
 
464
                    for (int i=0; i<3; ++i) {
 
465
                        libint_.PrimQuartet[nprim].U[0][i] = PA[i];
 
466
                        libint_.PrimQuartet[nprim].U[2][i] = QC[i];
 
467
                        libint_.PrimQuartet[nprim].U[4][i] = WP[i];
 
468
                        libint_.PrimQuartet[nprim].U[5][i] = WQ[i];
 
469
                    }
 
470
                    libint_.PrimQuartet[nprim].oo2z = oo2z;
 
471
                    libint_.PrimQuartet[nprim].oo2n = oo2n;
 
472
                    libint_.PrimQuartet[nprim].oo2zn = oo2zn;
 
473
                    libint_.PrimQuartet[nprim].poz = rho * ooz;
 
474
                    libint_.PrimQuartet[nprim].pon = rho * oon;
 
475
                    libint_.PrimQuartet[nprim].oo2p = oo2rho;
 
476
                    
 
477
                    double T = rho * PQ2;
 
478
                    calc_f(libint_.PrimQuartet[nprim].F, am+1, T);
 
479
                    
 
480
                    // Modify F to include overlap of ab and cd, eqs 14, 15, 16 of libint manual
 
481
                    double Scd = pow(M_PI*oon, 3.0/2.0) * exp(-a3*a4*oon*CD2) * c3 * c4;
 
482
                    double val = 2.0 * sqrt(rho * M_1_PI) * Sab * Scd;
 
483
                    for (int i=0; i<am+1; ++i) {
 
484
                        libint_.PrimQuartet[nprim].F[i] *= val;
 
485
                    }
 
486
                    nprim++;
 
487
                }
 
488
            }
 
489
        }
 
490
    }
 
491
    
 
492
    // How many are there?
 
493
    size_t size = INT_NCART(am1) * INT_NCART(am2) * INT_NCART(am3) * INT_NCART(am4);
 
494
    
 
495
    // Compute the integral
 
496
    if (am) {
 
497
        REALTYPE *target_ints;
 
498
 
 
499
        target_ints = build_eri[am1][am2][am3][am4](&libint_, nprim);
 
500
 
 
501
        memcpy(source_, target_ints, sizeof(double)*size);
 
502
    } else {
 
503
        // Handle (ss|ss)
 
504
        double temp = 0.0;
 
505
        for (size_t i=0; i<nprim_combination; ++i)
 
506
            temp += (double)libint_.PrimQuartet[i].F[0];
 
507
        source_[0] = temp;
 
508
    }
 
509
 
 
510
    // Normalize the integrals for angular momentum
 
511
    normalize_am(s1, s2, s3, s4);
 
512
     
 
513
    // Transform the integrals to the spherical basis
 
514
    pure_transform(sh1, sh2, sh3, sh4);
 
515
    
 
516
    // Results are in source_
 
517
}
 
518
 
 
519
void ERI::compute_shell_deriv1(int sh1, int sh2, int sh3, int sh4)
 
520
{
 
521
    if (deriv_ >= 1) {
 
522
        fprintf(stderr, "ERROR - ERI: ERI object not initialized to handle derivatives.\n");
 
523
        abort();
 
524
    }
 
525
    // Need to ensure the ordering asked by the user is valid for libint
 
526
    // compute_quartet does NOT check this. SEGFAULTS should occur if order
 
527
    // is not guaranteed.
 
528
    int s1, s2, s3, s4;
 
529
    int am1, am2, am3, am4, temp;
 
530
    bool p13p24 = false, p12 = false, p34 = false;
 
531
    
 
532
    // AM used for ordering
 
533
    am1 = bs1_->shell(sh1)->am(0);
 
534
    am2 = bs2_->shell(sh2)->am(0);
 
535
    am3 = bs3_->shell(sh3)->am(0);
 
536
    am4 = bs4_->shell(sh4)->am(0);
 
537
    
 
538
    int n1 = bs1_->shell(sh1)->nfunction(0);
 
539
    int n2 = bs1_->shell(sh2)->nfunction(0);
 
540
    int n3 = bs1_->shell(sh3)->nfunction(0);
 
541
    int n4 = bs1_->shell(sh4)->nfunction(0);
 
542
    
 
543
    // l(a) >= l(b), l(c) >= l(d), and l(c) + l(d) >= l(a) + l(b).
 
544
    if (am1 >= am2) {
 
545
        s1 = sh1;
 
546
        s2 = sh2;
 
547
    } else {
 
548
        s1 = sh2;
 
549
        s2 = sh1;
 
550
        p12 = true;
 
551
    }
 
552
    
 
553
    if (am3 >= am4) {
 
554
        s3 = sh3;
 
555
        s4 = sh4;
 
556
    } else {
 
557
        s3 = sh4;
 
558
        s4 = sh3;  
 
559
        p34 = true;
 
560
    }
 
561
    
 
562
    if ((am1 + am2) > (am3 + am4)) {
 
563
        // Swap s1 and s2 with s3 and s4
 
564
        temp = s1;
 
565
        s1 = s3;
 
566
        s3 = temp;
 
567
        
 
568
        temp = s2;
 
569
        s2 = s4;
 
570
        s4 = temp;
 
571
        p13p24 = true;
 
572
    }
 
573
    
 
574
    // s1, s2, s3, s4 contain the shells to do in libderive order
 
575
    compute_quartet_deriv1(s1, s2, s3, s4);    // compute 9 sets of integral derivatives
 
576
        
 
577
    size_t size = n1 * n2 * n3 * n4;
 
578
    // Permute integrals back, if needed
 
579
    if (p12 || p34 || p13p24) {
 
580
        // 3n of them
 
581
        for (int i=0; i<3*natom_; ++i)
 
582
            permute_target(source_+(i*size), target_+(i*size), s1, s2, s3, s4, p12, p34, p13p24);
 
583
    }
 
584
    else {
 
585
        // copy the integrals to the target_, 3n of them
 
586
        memcpy(target_, source_, 3 * natom_ * size *sizeof(double));
 
587
    }
 
588
}
 
589
 
 
590
void ERI::compute_quartet_deriv1(int sh1, int sh2, int sh3, int sh4)
 
591
{
 
592
    Ref<GaussianShell> s1, s2, s3, s4;
 
593
    
 
594
    s1 = bs1_->shell(sh1);
 
595
    s2 = bs2_->shell(sh2);
 
596
    s3 = bs3_->shell(sh3);
 
597
    s4 = bs4_->shell(sh4);
 
598
    
 
599
    int am1 = s1->am(0);
 
600
    int am2 = s2->am(0);
 
601
    int am3 = s3->am(0);
 
602
    int am4 = s4->am(0);
 
603
    int am = am1 + am2 + am3 + am4; // total am
 
604
    int nprim1 = s1->nprimitive();
 
605
    int nprim2 = s2->nprimitive();
 
606
    int nprim3 = s3->nprimitive();
 
607
    int nprim4 = s4->nprimitive();
 
608
    size_t nprim;
 
609
    double A[3], B[3], C[3], D[3];
 
610
    A[0] = s1->center()[0];
 
611
    A[1] = s1->center()[1];
 
612
    A[2] = s1->center()[2];
 
613
    B[0] = s2->center()[0];
 
614
    B[1] = s2->center()[1];
 
615
    B[2] = s2->center()[2];
 
616
    C[0] = s3->center()[0];
 
617
    C[1] = s3->center()[1];
 
618
    C[2] = s3->center()[2];
 
619
    D[0] = s4->center()[0];
 
620
    D[1] = s4->center()[1];
 
621
    D[2] = s4->center()[2];
 
622
    
 
623
    // compute intermediates
 
624
    double AB2 = 0.0;
 
625
    AB2 += (A[0] - B[0]) * (A[0] - B[0]);
 
626
    AB2 += (A[1] - B[1]) * (A[1] - B[1]);
 
627
    AB2 += (A[2] - B[2]) * (A[2] - B[2]);
 
628
    double CD2 = 0.0;
 
629
    CD2 += (C[0] - D[0]) * (C[0] - D[0]);
 
630
    CD2 += (C[1] - D[1]) * (C[1] - D[1]);
 
631
    CD2 += (C[2] - D[2]) * (C[2] - D[2]);
 
632
    
 
633
    libderiv_.AB[0] = A[0] - B[0];
 
634
    libderiv_.AB[1] = A[1] - B[1];
 
635
    libderiv_.AB[2] = A[2] - B[2];
 
636
    libderiv_.CD[0] = C[0] - D[0];
 
637
    libderiv_.CD[1] = C[1] - D[1];
 
638
    libderiv_.CD[2] = C[2] - D[2];
 
639
        
 
640
    // Prepare all the data needed by libint
 
641
    nprim = 0;
 
642
    for (int p1=0; p1<nprim1; ++p1) {
 
643
        double a1 = s1->exp(p1);
 
644
        double c1 = s1->coef(0, p1);
 
645
        for (int p2=0; p2<nprim2; ++p2) {
 
646
            double a2 = s2->exp(p2);
 
647
            double c2 = s2->coef(0, p2);
 
648
            double zeta = a1 + a2;
 
649
            double ooz = 1.0/zeta;
 
650
            double oo2z = 1.0/(2.0 * zeta);
 
651
            
 
652
            double PA[3], PB[3];
 
653
            double P[3];
 
654
            
 
655
            P[0] = (a1*A[0] + a2*B[0])*ooz;
 
656
            P[1] = (a1*A[1] + a2*B[1])*ooz;
 
657
            P[2] = (a1*A[2] + a2*B[2])*ooz;
 
658
            PA[0] = P[0] - A[0];
 
659
            PA[1] = P[1] - A[1];
 
660
            PA[2] = P[2] - A[2];
 
661
            PB[0] = P[0] - B[0];
 
662
            PB[1] = P[1] - B[1];
 
663
            PB[2] = P[2] - B[2];
 
664
            
 
665
            double Sab = pow(M_PI*ooz, 3.0/2.0) * exp(-a1*a2*ooz*AB2) * c1 * c2;
 
666
            
 
667
            for (int p3=0; p3<nprim3; ++p3) {
 
668
                double a3 = s3->exp(p3);
 
669
                double c3 = s3->coef(0, p3);
 
670
                for (int p4=0; p4<nprim4; ++p4) {
 
671
                    double a4 = s4->exp(p4);
 
672
                    double c4 = s4->coef(0, p4);
 
673
                    double nu = a3 + a4;
 
674
                    double oon = 1.0/nu;
 
675
                    double oo2n = 1.0/(2.0*nu);
 
676
                    double oo2zn = 1.0/(2.0*(zeta+nu));
 
677
                    double rho = (zeta*nu)/(zeta+nu);
 
678
                    double oo2rho = 1 / (2.0*rho);
 
679
                    
 
680
                    double QC[3], QD[3], WP[3], WQ[3], PQ[3];
 
681
                    double Q[3], W[3];
 
682
 
 
683
                    Q[0] = (a3*C[0] + a4*D[0])*oon;
 
684
                    Q[1] = (a3*C[1] + a4*D[1])*oon;
 
685
                    Q[2] = (a3*C[2] + a4*D[2])*oon;
 
686
                    QC[0] = Q[0] - C[0];
 
687
                    QC[1] = Q[1] - C[1];
 
688
                    QC[2] = Q[2] - C[2];
 
689
                    QD[0] = Q[0] - D[0];
 
690
                    QD[1] = Q[1] - D[1];
 
691
                    QD[2] = Q[2] - D[2];
 
692
                    PQ[0] = P[0] - Q[0];
 
693
                    PQ[1] = P[1] - Q[1];
 
694
                    PQ[2] = P[2] - Q[2];
 
695
                    
 
696
                    double PQ2 = 0.0;
 
697
                    PQ2 += (P[0] - Q[0]) * (P[0] - Q[0]);
 
698
                    PQ2 += (P[1] - Q[1]) * (P[1] - Q[1]);
 
699
                    PQ2 += (P[2] - Q[2]) * (P[2] - Q[2]);
 
700
                    
 
701
                    W[0] = (zeta*P[0] + nu*Q[0]) / (zeta + nu);
 
702
                    W[1] = (zeta*P[1] + nu*Q[1]) / (zeta + nu);
 
703
                    W[2] = (zeta*P[2] + nu*Q[2]) / (zeta + nu);
 
704
                    WP[0] = W[0] - P[0];
 
705
                    WP[1] = W[1] - P[1];
 
706
                    WP[2] = W[2] - P[2];
 
707
                    WQ[0] = W[0] - Q[0];
 
708
                    WQ[1] = W[1] - Q[1];
 
709
                    WQ[2] = W[2] - Q[2];
 
710
                    
 
711
                    for (int i=0; i<3; ++i) {
 
712
                        libderiv_.PrimQuartet[nprim].U[0][i] = PA[i];
 
713
                        libderiv_.PrimQuartet[nprim].U[1][i] = PB[i];
 
714
                        libderiv_.PrimQuartet[nprim].U[2][i] = QC[i];
 
715
                        libderiv_.PrimQuartet[nprim].U[3][i] = QD[i];
 
716
                        libderiv_.PrimQuartet[nprim].U[4][i] = WP[i];
 
717
                        libderiv_.PrimQuartet[nprim].U[5][i] = WQ[i];
 
718
                    }
 
719
                    libderiv_.PrimQuartet[nprim].oo2z = oo2z;
 
720
                    libderiv_.PrimQuartet[nprim].oo2n = oo2n;
 
721
                    libderiv_.PrimQuartet[nprim].oo2zn = oo2zn;
 
722
                    libderiv_.PrimQuartet[nprim].poz = rho * ooz;
 
723
                    libderiv_.PrimQuartet[nprim].pon = rho * oon;
 
724
                    // libderiv_.PrimQuartet[nprim].oo2p = oo2rho;   // NOT SET IN CINTS
 
725
                    libderiv_.PrimQuartet[nprim].twozeta_a = 2.0 * a1;
 
726
                    libderiv_.PrimQuartet[nprim].twozeta_b = 2.0 * a2;
 
727
                    libderiv_.PrimQuartet[nprim].twozeta_c = 2.0 * a3;
 
728
                    libderiv_.PrimQuartet[nprim].twozeta_d = 2.0 * a4;
 
729
                    
 
730
                    double T = rho * PQ2;
 
731
                    int_fjt(libderiv_.PrimQuartet[nprim].F, am+DERIV_LVL, T);
 
732
                    
 
733
                    // Modify F to include overlap of ab and cd, eqs 14, 15, 16 of libint manual
 
734
                    double Scd = pow(M_PI*oon, 3.0/2.0) * exp(-a3*a4*oon*CD2) * c3 * c4;
 
735
                    double val = 2.0 * sqrt(rho * M_1_PI) * Sab * Scd;
 
736
                    for (int i=0; i<=am+DERIV_LVL; ++i) {
 
737
                        libderiv_.PrimQuartet[nprim].F[i] *= val;
 
738
                    }
 
739
                   
 
740
                    nprim++;
 
741
                }
 
742
            }
 
743
        }
 
744
    }
 
745
    
 
746
    // How many are there?
 
747
    size_t size = INT_NCART(am1) * INT_NCART(am2) * INT_NCART(am3) * INT_NCART(am4);
 
748
    
 
749
    // Compute the integral
 
750
    build_deriv1_eri[am1][am2][am3][am4](&libderiv_, nprim);
 
751
 
 
752
    // Sum results into derivs_
 
753
    int center_i = s1->ncenter()*3*size;
 
754
    int center_j = s2->ncenter()*3*size;
 
755
    int center_k = s3->ncenter()*3*size;
 
756
    int center_l = s4->ncenter()*3*size;
 
757
    
 
758
    memset(source_, 0, sizeof(double) * size * 3 * natom_);
 
759
    
 
760
    for (size_t i=0; i < size; ++i) {
 
761
        source_[center_i+(0*size)+i] += libderiv_.ABCD[0][i];
 
762
        source_[center_i+(1*size)+i] += libderiv_.ABCD[1][i];
 
763
        source_[center_i+(2*size)+i] += libderiv_.ABCD[2][i];
 
764
 
 
765
        // Use translational invariance to determine center_j derivatives
 
766
        source_[center_j+(0*size)+i] -= (libderiv_.ABCD[0][i] + libderiv_.ABCD[6][i] + libderiv_.ABCD[9][i]);
 
767
        source_[center_j+(1*size)+i] -= (libderiv_.ABCD[1][i] + libderiv_.ABCD[7][i] + libderiv_.ABCD[10][i]);
 
768
        source_[center_j+(2*size)+i] -= (libderiv_.ABCD[2][i] + libderiv_.ABCD[8][i] + libderiv_.ABCD[11][i]);
 
769
        
 
770
        source_[center_k+(0*size)+i] += libderiv_.ABCD[6][i];
 
771
        source_[center_k+(1*size)+i] += libderiv_.ABCD[7][i];
 
772
        source_[center_k+(2*size)+i] += libderiv_.ABCD[8][i];
 
773
        
 
774
        source_[center_l+(0*size)+i] += libderiv_.ABCD[9][i];
 
775
        source_[center_l+(1*size)+i] += libderiv_.ABCD[10][i];
 
776
        source_[center_l+(2*size)+i] += libderiv_.ABCD[11][i];        
 
777
    }
 
778
    
 
779
    // Normalize the 3n sets of integrals
 
780
    normalize_am(s1, s2, s3, s4, 3*natom_);
 
781
        
 
782
    // Transform the integrals to the spherical basis
 
783
    // pure_transform can only do 1 at a time. We need to do 3n.
 
784
    for (int i=0; i<3*natom_; ++i)
 
785
        pure_transform(sh1, sh2, sh3, sh4, i);
 
786
 
 
787
    // Results are in source_
 
788
}
 
789