← Back to branch summary

~ubuntu-branches/ubuntu/karmic/psicode/karmic

« back to all changes in this revision

Viewing changes to src/lib/libmints/kinetic.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/kinetic.cc
 
1
#include <libciomr/libciomr.h>
 
2
 
 
3
#include <libmints/basisset.h>
 
4
#include <libmints/gshell.h>
 
5
#include <libmints/overlap.h>
 
6
#include <libmints/kinetic.h>
 
7
#include <physconst.h>
 
8
 
 
9
#define MAX(a, b) ((a) > (b) ? (a) : (b))
 
10
 
 
11
using namespace psi;
 
12
 
 
13
// Initialize overlap_recur_ to +1 basis set angular momentum
 
14
KineticInt::KineticInt(IntegralFactory* integral, Ref<BasisSet>& bs1, Ref<BasisSet>& bs2, int deriv) :
 
15
    OneBodyInt(integral, bs1, bs2, deriv), overlap_recur_(bs1->max_am()+1+deriv, bs2->max_am()+1+deriv)
 
16
{
 
17
    int maxam1 = bs1_->max_am();
 
18
    int maxam2 = bs2_->max_am();
 
19
    
 
20
    int maxnao1 = (maxam1+1)*(maxam1+2)/2;
 
21
    int maxnao2 = (maxam2+1)*(maxam2+2)/2;
 
22
    if (deriv == 1) {
 
23
        maxnao1 *= 3 * natom_;
 
24
        maxnao2 *= 3 * natom_;
 
25
    }
 
26
    
 
27
    buffer_ = new double[maxnao1*maxnao2];
 
28
}
 
29
 
 
30
KineticInt::~KineticInt()
 
31
{
 
32
    delete[] buffer_;
 
33
}
 
34
 
 
35
void KineticInt::compute_shell(int sh1, int sh2)
 
36
{
 
37
    compute_pair(bs1_->shell(sh1), bs2_->shell(sh2));
 
38
}
 
39
 
 
40
void KineticInt::compute_shell_deriv1(int sh1, int sh2)
 
41
{
 
42
    compute_pair_deriv1(bs1_->shell(sh1), bs2_->shell(sh2));
 
43
}
 
44
 
 
45
// The engine only supports segmented basis sets
 
46
void KineticInt::compute_pair(Ref<GaussianShell>& s1, Ref<GaussianShell>& s2)
 
47
{
 
48
    int ao12;
 
49
    int am1 = s1->am(0);
 
50
    int am2 = s2->am(0);
 
51
    int nprim1 = s1->nprimitive();
 
52
    int nprim2 = s2->nprimitive();
 
53
    double A[3], B[3];
 
54
    A[0] = s1->center()[0];
 
55
    A[1] = s1->center()[1];
 
56
    A[2] = s1->center()[2];
 
57
    B[0] = s2->center()[0];
 
58
    B[1] = s2->center()[1];
 
59
    B[2] = s2->center()[2];
 
60
    
 
61
    // compute intermediates
 
62
    double AB2 = 0.0;
 
63
    AB2 += (A[0] - B[0]) * (A[0] - B[0]);
 
64
    AB2 += (A[1] - B[1]) * (A[1] - B[1]);
 
65
    AB2 += (A[2] - B[2]) * (A[2] - B[2]);
 
66
    
 
67
    memset(buffer_, 0, s1->ncartesian() * s2->ncartesian() * sizeof(double));
 
68
    
 
69
    double **x = overlap_recur_.x();
 
70
    double **y = overlap_recur_.y();
 
71
    double **z = overlap_recur_.z();
 
72
    
 
73
    for (int p1=0; p1<nprim1; ++p1) {
 
74
        double a1 = s1->exp(p1);
 
75
        double c1 = s1->coef(0, p1);
 
76
        for (int p2=0; p2<nprim2; ++p2) {
 
77
            double a2 = s2->exp(p2);
 
78
            double c2 = s2->coef(0, p2);
 
79
            double gamma = a1 + a2;
 
80
            double oog = 1.0/gamma;
 
81
            
 
82
            double PA[3], PB[3];
 
83
            double P[3];
 
84
            
 
85
            P[0] = (a1*A[0] + a2*B[0])*oog;
 
86
            P[1] = (a1*A[1] + a2*B[1])*oog;
 
87
            P[2] = (a1*A[2] + a2*B[2])*oog;
 
88
            PA[0] = P[0] - A[0];
 
89
            PA[1] = P[1] - A[1];
 
90
            PA[2] = P[2] - A[2];
 
91
            PB[0] = P[0] - B[0];
 
92
            PB[1] = P[1] - B[1];
 
93
            PB[2] = P[2] - B[2];
 
94
            
 
95
            double over_pf = exp(-a1*a2*AB2*oog) * sqrt(M_PI*oog) * M_PI * oog * c1 * c2;
 
96
            
 
97
            // Do recursion
 
98
            overlap_recur_.compute(PA, PB, gamma, am1+1, am2+1);
 
99
            
 
100
            ao12 = 0;
 
101
            for(int ii = 0; ii <= am1; ii++) {
 
102
                int l1 = am1 - ii;
 
103
                for(int jj = 0; jj <= ii; jj++) {
 
104
                    int m1 = ii - jj;
 
105
                    int n1 = jj;
 
106
                    /*--- create all am components of sj ---*/
 
107
                    for(int kk = 0; kk <= am2; kk++) {
 
108
                        int l2 = am2 - kk;
 
109
                        for(int ll = 0; ll <= kk; ll++) {
 
110
                            int m2 = kk - ll;
 
111
                            int n2 = ll;
 
112
 
 
113
                            double I1, I2, I3, I4;
 
114
                            
 
115
                            I1 = (l1 == 0 || l2 == 0) ? 0.0 : x[l1-1][l2-1] * y[m1][m2] * z[n1][n2] * over_pf;
 
116
                            I2 = x[l1+1][l2+1] * y[m1][m2] * z[n1][n2] * over_pf;
 
117
                            I3 = (l2 == 0) ? 0.0 : x[l1+1][l2-1] * y[m1][m2] * z[n1][n2] * over_pf;
 
118
                            I4 = (l1 == 0) ? 0.0 : x[l1-1][l2+1] * y[m1][m2] * z[n1][n2] * over_pf;
 
119
                            // fprintf(outfile, "I1=%f, I2=%f, I3=%f, I4=%f\n", I1, I2, I3, I4);
 
120
                            double Ix = 0.5 * l1 * l2 * I1 + 2.0 * a1 * a2 * I2 - a1 * l2 * I3 - l1 * a2 * I4;
 
121
 
 
122
                            I1 = (m1 == 0 || m2 == 0) ? 0.0 : x[l1][l2] * y[m1-1][m2-1] * z[n1][n2] * over_pf;
 
123
                            I2 = x[l1][l2] * y[m1+1][m2+1] * z[n1][n2] * over_pf;
 
124
                            I3 = (m2 == 0) ? 0.0 : x[l1][l2] * y[m1+1][m2-1] * z[n1][n2] * over_pf;
 
125
                            I4 = (m1 == 0) ? 0.0 : x[l1][l2] * y[m1-1][m2+1] * z[n1][n2] * over_pf;
 
126
                            // fprintf(outfile, "I1=%f, I2=%f, I3=%f, I4=%f\n", I1, I2, I3, I4);
 
127
                            double Iy = 0.5 * m1 * m2 * I1 + 2.0 * a1 * a2 * I2 - a1 * m2 * I3 - m1 * a2 * I4;
 
128
                            
 
129
                            I1 = (n1 == 0 || n2 == 0) ? 0.0 : x[l1][l2] * y[m1][m2] * z[n1-1][n2-1] * over_pf;
 
130
                            I2 = x[l1][l2] * y[m1][m2] * z[n1+1][n2+1] * over_pf;
 
131
                            I3 = (n2 == 0) ? 0.0 : x[l1][l2] * y[m1][m2] * z[n1+1][n2-1] * over_pf;
 
132
                            I4 = (n1 == 0) ? 0.0 : x[l1][l2] * y[m1][m2] * z[n1-1][n2+1] * over_pf;
 
133
                            // fprintf(outfile, "I1=%f, I2=%f, I3=%f, I4=%f\n", I1, I2, I3, I4);
 
134
                            double Iz = 0.5 * n1 * n2 * I1 + 2.0 * a1 * a2 * I2 - a1 * n2 * I3 - n1 * a2 * I4;
 
135
                            
 
136
                            // buffer_[ao12++] += over_pf * (Ix + Iy + Iz);
 
137
                            buffer_[ao12++] += (Ix + Iy + Iz);
 
138
                            // fprintf(outfile, "Ix=%f, Iy=%f Iz=%f\n", Ix, Iy, Iz);
 
139
                        }
 
140
                    }
 
141
                }
 
142
            }
 
143
        }
 
144
    }
 
145
    
 
146
    // Integrals are done. Normalize for angular momentum
 
147
    normalize_am(s1, s2);
 
148
    
 
149
    // for (int i=0; i<s1->ncartesian() * s2->ncartesian(); ++i) {
 
150
    //     fprintf(outfile, "integral (am=%d|am=%d) = %f\n", am1, am2, buffer_[i]);
 
151
    // }
 
152
    // Spherical harmonic transformation
 
153
    // Wrapped up in the AO to SO transformation (I think)
 
154
    // spherical_transform_1e(s1, s2);
 
155
}
 
156
 
 
157
static double ke_int(double **x, double **y, double **z, double a1, int l1, int m1, int n1, 
 
158
              double a2, int l2, int m2, int n2)
 
159
{
 
160
    double I1, I2, I3, I4;
 
161
    
 
162
    I1 = (l1 == 0 || l2 == 0) ? 0.0 : x[l1-1][l2-1] * y[m1][m2] * z[n1][n2];
 
163
    I2 = x[l1+1][l2+1] * y[m1][m2] * z[n1][n2];
 
164
    I3 = (l2 == 0) ? 0.0 : x[l1+1][l2-1] * y[m1][m2] * z[n1][n2];
 
165
    I4 = (l1 == 0) ? 0.0 : x[l1-1][l2+1] * y[m1][m2] * z[n1][n2];
 
166
    double Ix = 0.5 * l1 * l2 * I1 + 2.0 * a1 * a2 * I2 - a1 * l2 * I3 - l1 * a2 * I4;
 
167
 
 
168
    I1 = (m1 == 0 || m2 == 0) ? 0.0 : x[l1][l2] * y[m1-1][m2-1] * z[n1][n2];
 
169
    I2 = x[l1][l2] * y[m1+1][m2+1] * z[n1][n2];
 
170
    I3 = (m2 == 0) ? 0.0 : x[l1][l2] * y[m1+1][m2-1] * z[n1][n2];
 
171
    I4 = (m1 == 0) ? 0.0 : x[l1][l2] * y[m1-1][m2+1] * z[n1][n2];
 
172
    double Iy = 0.5 * m1 * m2 * I1 + 2.0 * a1 * a2 * I2 - a1 * m2 * I3 - m1 * a2 * I4;
 
173
    
 
174
    I1 = (n1 == 0 || n2 == 0) ? 0.0 : x[l1][l2] * y[m1][m2] * z[n1-1][n2-1];
 
175
    I2 = x[l1][l2] * y[m1][m2] * z[n1+1][n2+1];
 
176
    I3 = (n2 == 0) ? 0.0 : x[l1][l2] * y[m1][m2] * z[n1+1][n2-1];
 
177
    I4 = (n1 == 0) ? 0.0 : x[l1][l2] * y[m1][m2] * z[n1-1][n2+1];
 
178
    double Iz = 0.5 * n1 * n2 * I1 + 2.0 * a1 * a2 * I2 - a1 * n2 * I3 - n1 * a2 * I4;
 
179
    
 
180
    return (Ix + Iy + Iz);
 
181
}
 
182
 
 
183
// The engine only supports segmented basis sets
 
184
void KineticInt::compute_pair_deriv1(Ref<GaussianShell>& s1, Ref<GaussianShell>& s2)
 
185
{
 
186
    int ao12;
 
187
    int am1 = s1->am(0);
 
188
    int am2 = s2->am(0);
 
189
    int nprim1 = s1->nprimitive();
 
190
    int nprim2 = s2->nprimitive();
 
191
    double A[3], B[3];
 
192
    A[0] = s1->center()[0];
 
193
    A[1] = s1->center()[1];
 
194
    A[2] = s1->center()[2];
 
195
    B[0] = s2->center()[0];
 
196
    B[1] = s2->center()[1];
 
197
    B[2] = s2->center()[2];
 
198
    
 
199
    size_t size = s1->ncartesian() * s2->ncartesian();
 
200
    int center_i = s1->ncenter()*3*size;
 
201
    int center_j = s2->ncenter()*3*size;
 
202
    
 
203
    // compute intermediates
 
204
    double AB2 = 0.0;
 
205
    AB2 += (A[0] - B[0]) * (A[0] - B[0]);
 
206
    AB2 += (A[1] - B[1]) * (A[1] - B[1]);
 
207
    AB2 += (A[2] - B[2]) * (A[2] - B[2]);
 
208
    
 
209
    memset(buffer_, 0, 3 * natom_ * s1->ncartesian() * s2->ncartesian() * sizeof(double));
 
210
    
 
211
    double **x = overlap_recur_.x();
 
212
    double **y = overlap_recur_.y();
 
213
    double **z = overlap_recur_.z();
 
214
    
 
215
    for (int p1=0; p1<nprim1; ++p1) {
 
216
        double a1 = s1->exp(p1);
 
217
        double c1 = s1->coef(0, p1);
 
218
        for (int p2=0; p2<nprim2; ++p2) {
 
219
            double a2 = s2->exp(p2);
 
220
            double c2 = s2->coef(0, p2);
 
221
            double gamma = a1 + a2;
 
222
            double oog = 1.0/gamma;
 
223
            
 
224
            double PA[3], PB[3];
 
225
            double P[3];
 
226
            
 
227
            P[0] = (a1*A[0] + a2*B[0])*oog;
 
228
            P[1] = (a1*A[1] + a2*B[1])*oog;
 
229
            P[2] = (a1*A[2] + a2*B[2])*oog;
 
230
            PA[0] = P[0] - A[0];
 
231
            PA[1] = P[1] - A[1];
 
232
            PA[2] = P[2] - A[2];
 
233
            PB[0] = P[0] - B[0];
 
234
            PB[1] = P[1] - B[1];
 
235
            PB[2] = P[2] - B[2];
 
236
            
 
237
            double over_pf = exp(-a1*a2*AB2*oog) * sqrt(M_PI*oog) * M_PI * oog * c1 * c2;
 
238
            
 
239
            // Do recursion
 
240
            overlap_recur_.compute(PA, PB, gamma, am1+2, am2+2);
 
241
            
 
242
            ao12 = 0;
 
243
            for(int ii = 0; ii <= am1; ii++) {
 
244
                int l1 = am1 - ii;
 
245
                for(int jj = 0; jj <= ii; jj++) {
 
246
                    int m1 = ii - jj;
 
247
                    int n1 = jj;
 
248
                    /*--- create all am components of sj ---*/
 
249
                    for(int kk = 0; kk <= am2; kk++) {
 
250
                        int l2 = am2 - kk;
 
251
                        for(int ll = 0; ll <= kk; ll++) {
 
252
                            int m2 = kk - ll;
 
253
                            int n2 = ll;
 
254
 
 
255
                            // x on center i
 
256
                            buffer_[center_i+(0*size)+ao12] += 2.0 * a1 * ke_int(x, y, z, a1, l1+1, m1, n1, a2, l2, m2, n2) * over_pf;
 
257
                            if (l1)
 
258
                                buffer_[center_i+(0*size)+ao12] -= ke_int(x, y, z, a1, l1-1, m1, n1, a2, l2, m2, n2) * over_pf;
 
259
                            // y on center i
 
260
                            buffer_[center_i+(1*size)+ao12] += 2.0 * a1 * ke_int(x, y, z, a1, l1, m1+1, n1, a2, l2, m2, n2) * over_pf;
 
261
                            if (m1)
 
262
                                buffer_[center_i+(1*size)+ao12] -= ke_int(x, y, z, a1, l1, m1-1, n1, a2, l2, m2, n2) * over_pf;
 
263
                            // z on center i
 
264
                            buffer_[center_i+(2*size)+ao12] += 2.0 * a1 * ke_int(x, y, z, a1, l1, m1, n1+1, a2, l2, m2, n2) * over_pf;
 
265
                            if (n1)
 
266
                                buffer_[center_i+(2*size)+ao12] -= ke_int(x, y, z, a1, l1, m1, n1-1, a2, l2, m2, n2) * over_pf;
 
267
 
 
268
                            // x on center j
 
269
                            buffer_[center_j+(0*size)+ao12] += 2.0 * a2 * ke_int(x, y, z, a1, l1, m1, n1, a2, l2+1, m2, n2) * over_pf;
 
270
                            if (l2)
 
271
                                buffer_[center_j+(0*size)+ao12] -= ke_int(x, y, z, a1, l1, m1, n1, a2, l2-1, m2, n2) * over_pf;
 
272
                            // y on center j
 
273
                            buffer_[center_j+(1*size)+ao12] += 2.0 * a2 * ke_int(x, y, z, a1, l1, m1, n1, a2, l2, m2+1, n2) * over_pf;
 
274
                            if (m2)
 
275
                                buffer_[center_j+(1*size)+ao12] -= ke_int(x, y, z, a1, l1, m1, n1, a2, l2, m2-1, n2) * over_pf;
 
276
                            // z on center j
 
277
                            buffer_[center_j+(2*size)+ao12] += 2.0 * a2 * ke_int(x, y, z, a1, l1, m1, n1, a2, l2, m2, n2+1) * over_pf;
 
278
                            if (n2)
 
279
                                buffer_[center_j+(2*size)+ao12] -= ke_int(x, y, z, a1, l1, m1, n1, a2, l2, m2, n2-1) * over_pf;
 
280
                                
 
281
                            ao12++;
 
282
                        }
 
283
                    }
 
284
                }
 
285
            }
 
286
        }
 
287
    }
 
288
    
 
289
    // Integrals are done. Normalize for angular momentum
 
290
    normalize_am(s1, s2);
 
291
    
 
292
    // Spherical harmonic transformation
 
293
    // Wrapped up in the AO to SO transformation (I think)
 
294
}