~ubuntu-branches/ubuntu/vivid/psicode/vivid

« back to all changes in this revision

Viewing changes to src/bin/cints/GIAO_Deriv/giao_oe_deriv.cc

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2008-06-07 16:49:57 UTC
  • mfrom: (2.1.2 hardy)
  • Revision ID: james.westby@ubuntu.com-20080607164957-8pifvb133yjlkagn
Tags: 3.3.0-3
* debian/rules (DEB_MAKE_CHECK_TARGET): Do not abort test suite on
  failures.
* debian/rules (DEB_CONFIGURE_EXTRA_FLAGS): Set ${bindir} to /usr/lib/psi.
* debian/rules (install/psi3): Move psi3 file to /usr/bin.
* debian/patches/07_464867_move_executables.dpatch: New patch, add
  /usr/lib/psi to the $PATH, so that the moved executables are found.
  (closes: #464867)
* debian/patches/00list: Adjusted.

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
extern "C" {
 
2
#include<math.h>
 
3
#include<stdio.h>
 
4
#include<stdlib.h>
 
5
#include<string.h>
 
6
#include<libipv1/ip_lib.h>
 
7
#include<libiwl/iwl.h>
 
8
#include<libciomr/libciomr.h>
 
9
#include<libint/libint.h>
 
10
#include<psifiles.h>
 
11
 
 
12
#include"defines.h"
 
13
#define EXTERN
 
14
#include"global.h"
 
15
#include"taylor_fm_eval.h"
 
16
#include"oe_osrr.h"
 
17
#include"small_fns.h"
 
18
}
 
19
 
 
20
/*-----------------------------------------------------------------------------------------
 
21
  This function computes some derivatives of integrals over GIAO Gaussians with respect to
 
22
  E and B fields (at E=0, B=0, i.e. in absence of external fields) and writes them out.
 
23
  
 
24
                                   t
 
25
  d2h/dBdE = -0.5 i <mu| (AB x r) r  |nu>
 
26
  
 
27
  d h/dB   =  0.5   <mu| L   + i (AB x r) h |nu>
 
28
                          B
 
29
 
 
30
  d h/dE   =  - <mu| r |nu>   -- simply electric dipole integral, already computed, hence
 
31
                                 not computed here
 
32
 -----------------------------------------------------------------------------------------*/
 
33
extern "C" void giao_oe_deriv()
 
34
{
 
35
 
 
36
#ifdef USE_TAYLOR_FM
 
37
  init_Taylor_Fm_Eval(BasisSet.max_am*4+1,UserOptions.cutoff);
 
38
#endif  
 
39
 
 
40
  /*--- allocate room for the one-e matrices ---*/
 
41
  double** D2HDBXDEX = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
 
42
  double** D2HDBXDEY = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
 
43
  double** D2HDBXDEZ = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
 
44
  double** D2HDBYDEX = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
 
45
  double** D2HDBYDEY = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
 
46
  double** D2HDBYDEZ = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
 
47
  double** D2HDBZDEX = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
 
48
  double** D2HDBZDEY = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
 
49
  double** D2HDBZDEZ = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
 
50
  double** DHDBX = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
 
51
  double** DHDBY = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
 
52
  double** DHDBZ = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
 
53
  double** DSDBX = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
 
54
  double** DSDBY = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
 
55
  double** DSDBZ = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
 
56
 
 
57
  /*--- allocate storage for shell blocks of one electron integrals ---*/
 
58
  int dimension = ioff[BasisSet.max_am];
 
59
  double** d2HdBxdExtemp = block_matrix(dimension,dimension);
 
60
  double** d2HdBxdEytemp = block_matrix(dimension,dimension);
 
61
  double** d2HdBxdEztemp = block_matrix(dimension,dimension);
 
62
  double** d2HdBydExtemp = block_matrix(dimension,dimension);
 
63
  double** d2HdBydEytemp = block_matrix(dimension,dimension);
 
64
  double** d2HdBydEztemp = block_matrix(dimension,dimension);
 
65
  double** d2HdBzdExtemp = block_matrix(dimension,dimension);
 
66
  double** d2HdBzdEytemp = block_matrix(dimension,dimension);
 
67
  double** d2HdBzdEztemp = block_matrix(dimension,dimension);
 
68
  double** dHdBxtemp = block_matrix(dimension,dimension);
 
69
  double** dHdBytemp = block_matrix(dimension,dimension);
 
70
  double** dHdBztemp = block_matrix(dimension,dimension);
 
71
  double** dSdBxtemp = block_matrix(dimension,dimension);
 
72
  double** dSdBytemp = block_matrix(dimension,dimension);
 
73
  double** dSdBztemp = block_matrix(dimension,dimension);
 
74
 
 
75
  /* Note the "+1" on the row dimension -- we are computing kinetic energy
 
76
     integrals with one extra quantum in the bra function here */
 
77
  double** OIX = block_matrix(BasisSet.max_am+1,BasisSet.max_am+2);
 
78
  double** OIY = block_matrix(BasisSet.max_am+1,BasisSet.max_am+2);
 
79
  double** OIZ = block_matrix(BasisSet.max_am+1,BasisSet.max_am+2);
 
80
  /* same here -- need to compute integrals with extra quantum in the bra */
 
81
  int indmax_0 = (BasisSet.max_am-1)*BasisSet.max_am*BasisSet.max_am+1;
 
82
  int indmax_p1 = (BasisSet.max_am)*(BasisSet.max_am+1)*(BasisSet.max_am+1)+1;
 
83
  double*** AI0 = init_box(indmax_p1,indmax_0,2*BasisSet.max_am+3);
 
84
  
 
85
  for (int si=0; si<BasisSet.num_shells; si++){
 
86
    int am_i = BasisSet.shells[si].am-1;
 
87
    int ni = ioff[BasisSet.shells[si].am];
 
88
    struct coordinates A = Molecule.centers[BasisSet.shells[si].center-1];
 
89
    int ioffset = BasisSet.shells[si].fao - 1;
 
90
    int izm = 1;
 
91
    int iym = am_i+2;     // Note +2 instead of +1 here -- we'll be computing nuclear attraction
 
92
                          // integrals with extra quantum in bra!
 
93
    int ixm = iym*iym;
 
94
 
 
95
    for (int sj=0; sj<BasisSet.num_shells; sj++){
 
96
      int nj = ioff[BasisSet.shells[sj].am];
 
97
      int am_j = BasisSet.shells[sj].am-1;
 
98
      struct coordinates B = Molecule.centers[BasisSet.shells[sj].center-1];
 
99
      int joffset = BasisSet.shells[sj].fao - 1;
 
100
      int jzm = 1;
 
101
      int jym = am_j+1;
 
102
      int jxm = jym*jym;
 
103
 
 
104
      struct shell_pair* sp = &(BasisSet.shell_pairs[si][sj]);
 
105
      struct coordinates AB;
 
106
      AB.x = sp->AB[0];
 
107
      AB.y = sp->AB[1];
 
108
      AB.z = sp->AB[2];
 
109
      double ab2 = AB.x * AB.x;
 
110
      ab2 += AB.y * AB.y;
 
111
      ab2 += AB.z * AB.z;
 
112
 
 
113
      /*--- zero the temporary storage for accumulating contractions ---*/
 
114
      memset(d2HdBxdExtemp[0],0,sizeof(double)*dimension*dimension);
 
115
      memset(d2HdBxdEytemp[0],0,sizeof(double)*dimension*dimension);
 
116
      memset(d2HdBxdEztemp[0],0,sizeof(double)*dimension*dimension);
 
117
      memset(d2HdBydExtemp[0],0,sizeof(double)*dimension*dimension);
 
118
      memset(d2HdBydEytemp[0],0,sizeof(double)*dimension*dimension);
 
119
      memset(d2HdBydEztemp[0],0,sizeof(double)*dimension*dimension);
 
120
      memset(d2HdBzdExtemp[0],0,sizeof(double)*dimension*dimension);
 
121
      memset(d2HdBzdEytemp[0],0,sizeof(double)*dimension*dimension);
 
122
      memset(d2HdBzdEztemp[0],0,sizeof(double)*dimension*dimension);
 
123
      memset(dHdBxtemp[0],0,sizeof(double)*dimension*dimension);
 
124
      memset(dHdBytemp[0],0,sizeof(double)*dimension*dimension);
 
125
      memset(dHdBztemp[0],0,sizeof(double)*dimension*dimension);
 
126
      memset(dSdBxtemp[0],0,sizeof(double)*dimension*dimension);
 
127
      memset(dSdBytemp[0],0,sizeof(double)*dimension*dimension);
 
128
      memset(dSdBztemp[0],0,sizeof(double)*dimension*dimension);
 
129
      
 
130
      /*--- contract by primitives here ---*/
 
131
      for (int i = 0; i < BasisSet.shells[si].n_prims; i++) {
 
132
        double a1 = sp->a1[i];
 
133
        double inorm = sp->inorm[i];
 
134
        for (int j = 0; j < BasisSet.shells[sj].n_prims; j++) {
 
135
          double a2 = sp->a2[j];
 
136
          double gam = sp->gamma[i][j];
 
137
          double jnorm = sp->jnorm[j];
 
138
          struct coordinates PA, PB;
 
139
          PA.x = sp->PA[i][j][0];
 
140
          PA.y = sp->PA[i][j][1];
 
141
          PA.z = sp->PA[i][j][2];
 
142
          PB.x = sp->PB[i][j][0];
 
143
          PB.y = sp->PB[i][j][1];
 
144
          PB.z = sp->PB[i][j][2];
 
145
          double oog = 1.0/gam;
 
146
          double over_pf = exp(-a1*a2*ab2*oog)*sqrt(M_PI*oog)*M_PI*oog*inorm*jnorm;
 
147
 
 
148
          OI_OSrecurs(OIX,OIY,OIZ,PA,PB,gam,am_i+1,am_j+2);
 
149
 
 
150
          /*--- create all am components of si ---*/
 
151
          int ai = 0;
 
152
          for(int ii = 0; ii <= am_i; ii++){
 
153
            int l1 = am_i - ii;
 
154
            for(int jj = 0; jj <= ii; jj++){
 
155
              int m1 = ii - jj;
 
156
              int n1 = jj;
 
157
              /*--- create all am components of sj ---*/
 
158
              int aj = 0;
 
159
              for(int kk = 0; kk <= am_j; kk++){
 
160
                int l2 = am_j - kk;
 
161
                for(int ll = 0; ll <= kk; ll++){
 
162
                  int m2 = kk - ll;
 
163
                  int n2 = ll;
 
164
 
 
165
                  double x00 = OIX[l1][l2];
 
166
                  double y00 = OIY[m1][m2];
 
167
                  double z00 = OIZ[n1][n2];
 
168
                  
 
169
                  double x01 = OIX[l1][l2+1];
 
170
                  double y01 = OIY[m1][m2+1];
 
171
                  double z01 = OIZ[n1][n2+1];
 
172
                  
 
173
                  double x10 = OIX[l1+1][l2];
 
174
                  double y10 = OIY[m1+1][m2];
 
175
                  double z10 = OIZ[n1+1][n2];
 
176
                  
 
177
                  double x11 = OIX[l1+1][l2+1];
 
178
                  double y11 = OIY[m1+1][m2+1];
 
179
                  double z11 = OIZ[n1+1][n2+1];
 
180
 
 
181
                  double x0 = x00;
 
182
                  double y0 = y00;
 
183
                  double z0 = z00;
 
184
                  double x1 = x01 + x00*B.x;
 
185
                  double y1 = y01 + y00*B.y;
 
186
                  double z1 = z01 + z00*B.z;
 
187
                  double x2 = x11 + x10*B.x + x01*A.x + x00*A.x*B.x;
 
188
                  double y2 = y11 + y10*B.y + y01*A.y + y00*A.y*B.y;
 
189
                  double z2 = z11 + z10*B.z + z01*A.z + z00*A.z*B.z;
 
190
 
 
191
                  double dSdB_pfac = 0.5*over_pf;
 
192
                  double x = x1*y0*z0;
 
193
                  double y = x0*y1*z0;
 
194
                  double z = x0*y0*z1;
 
195
                  dSdBxtemp[ai][aj] += dSdB_pfac * (AB.y * z - AB.z * y);
 
196
                  dSdBytemp[ai][aj] += dSdB_pfac * (AB.z * x - AB.x * z);
 
197
                  dSdBztemp[ai][aj] += dSdB_pfac * (AB.x * y - AB.y * x);
 
198
 
 
199
/*                  dSdBxtemp[ai][aj] += 2.0*dSdB_pfac * x;
 
200
                  dSdBytemp[ai][aj] += 2.0*dSdB_pfac * y;
 
201
                  dSdBztemp[ai][aj] += 2.0*dSdB_pfac * z;*/
 
202
                  
 
203
                  double d2BE_pfac = -0.5*over_pf;
 
204
                  double qxx = x2*y0*z0;
 
205
                  double qxy = x1*y1*z0;
 
206
                  double qxz = x1*y0*z1;
 
207
                  double qyx = qxy;
 
208
                  double qyy = x0*y2*z0;
 
209
                  double qyz = x0*y1*z1;
 
210
                  double qzx = qxz;
 
211
                  double qzy = qyz;
 
212
                  double qzz = x0*y0*z2;
 
213
                  d2HdBxdExtemp[ai][aj] += d2BE_pfac * ( AB.y * qzx - AB.z * qyx);
 
214
                  d2HdBxdEytemp[ai][aj] += d2BE_pfac * ( AB.y * qzy - AB.z * qyy);
 
215
                  d2HdBxdEztemp[ai][aj] += d2BE_pfac * ( AB.y * qzz - AB.z * qyz);
 
216
                  d2HdBydExtemp[ai][aj] += d2BE_pfac * ( AB.z * qxx - AB.x * qzx);
 
217
                  d2HdBydEytemp[ai][aj] += d2BE_pfac * ( AB.z * qxy - AB.x * qzy);
 
218
                  d2HdBydEztemp[ai][aj] += d2BE_pfac * ( AB.z * qxz - AB.x * qzz);
 
219
                  d2HdBzdExtemp[ai][aj] += d2BE_pfac * ( AB.x * qyx - AB.y * qxx);
 
220
                  d2HdBzdEytemp[ai][aj] += d2BE_pfac * ( AB.x * qyy - AB.y * qxy);
 
221
                  d2HdBzdEztemp[ai][aj] += d2BE_pfac * ( AB.x * qyz - AB.y * qxz);
 
222
                  
 
223
                  double nx = -2.0*a2*x01;
 
224
                  if (l2 >= 1)
 
225
                    nx += l2*OIX[l1][l2-1];
 
226
                  double ny = -2.0*a2*y01;
 
227
                  if (m2 >= 1)
 
228
                    ny += m2*OIY[m1][m2-1];
 
229
                  double nz = -2.0*a2*z01;
 
230
                  if (n2 >= 1)
 
231
                    nz += n2*OIZ[n1][n2-1];
 
232
 
 
233
                  double dB_pfac = -0.5*over_pf;
 
234
                  dHdBxtemp[ai][aj] += dB_pfac * x00 * ( y01 * nz - z01 * ny);
 
235
                  dHdBytemp[ai][aj] += dB_pfac * y00 * ( z01 * nx - x01 * nz);
 
236
                  dHdBztemp[ai][aj] += dB_pfac * z00 * ( x01 * ny - y01 * nx);
 
237
 
 
238
                  double tx00 = a2*(2*l2+1)*OIX[l1][l2] - 2*a2*a2*OIX[l1][l2+2];
 
239
                  if (l2 >= 2)
 
240
                    tx00 -= 0.5*l2*(l2-1)*OIX[l1][l2-2];
 
241
                  double ty00 = a2*(2*m2+1)*OIY[m1][m2] - 2*a2*a2*OIY[m1][m2+2];
 
242
                  if (m2 >= 2)
 
243
                    ty00 -= 0.5*m2*(m2-1)*OIY[m1][m2-2];
 
244
                  double tz00 = a2*(2*n2+1)*OIZ[n1][n2] - 2*a2*a2*OIZ[n1][n2+2];
 
245
                  if (n2 >= 2)
 
246
                    tz00 -= 0.5*n2*(n2-1)*OIZ[n1][n2-2];
 
247
 
 
248
                  double tx10 = a2*(2*l2+1)*(OIX[l1+1][l2] + A.x*OIX[l1][l2]) - 2*a2*a2*(OIX[l1+1][l2+2] + A.x*OIX[l1][l2+2]);
 
249
                  if (l2 >= 2)
 
250
                    tx10 -= 0.5*l2*(l2-1)*(OIX[l1+1][l2-2] + A.x*OIX[l1][l2-2]);
 
251
                  double ty10 = a2*(2*m2+1)*(OIY[m1+1][m2] + A.y*OIY[m1][m2]) - 2*a2*a2*(OIY[m1+1][m2+2] + A.y*OIY[m1][m2+2]);
 
252
                  if (m2 >= 2)
 
253
                    ty10 -= 0.5*m2*(m2-1)*(OIY[m1+1][m2-2] + A.y*OIY[m1][m2-2]);
 
254
                  double tz10 = a2*(2*n2+1)*(OIZ[n1+1][n2] + A.z*OIZ[n1][n2]) - 2*a2*a2*(OIZ[n1+1][n2+2] + A.z*OIZ[n1][n2+2]);
 
255
                  if (n2 >= 2)
 
256
                    tz10 -= 0.5*n2*(n2-1)*(OIZ[n1+1][n2-2] + A.z*OIZ[n1][n2-2]);
 
257
 
 
258
                  double xT = tx10*y0*z0 + x1*ty00*z0 + x1*y0*tz00;
 
259
                  double yT = tx00*y1*z0 + x0*ty10*z0 + x0*y1*tz00;
 
260
                  double zT = tx00*y0*z1 + x0*ty00*z1 + x0*y0*tz10;
 
261
 
 
262
                  dB_pfac = 0.5*over_pf;
 
263
                  dHdBxtemp[ai][aj] += dB_pfac * ( AB.y * zT - AB.z * yT);
 
264
                  dHdBytemp[ai][aj] += dB_pfac * ( AB.z * xT - AB.x * zT);
 
265
                  dHdBztemp[ai][aj] += dB_pfac * ( AB.x * yT - AB.y * xT);
 
266
 
 
267
/*                  dHdBxtemp[ai][aj] += 2.0*dB_pfac * (tx00*y0*z0 + x0*ty00*z0+x0*y0*tz00);
 
268
                  dHdBztemp[ai][aj] += 2.0*dB_pfac * (tx00*y0*z0 + x0*ty00*z0+x0*y0*tz00);*/
 
269
                  
 
270
                  /*double T = tx00*y0*z0 + x0*ty00*z0 + x0*y0*tz00;
 
271
                  dHdBytemp[ai][aj] += 2.0*dB_pfac*T;*/
 
272
                  
 
273
                  aj++;
 
274
                }
 
275
              }
 
276
              ai++;
 
277
            }
 
278
          } /*--- end cartesian components for (si,sj) with primitives (i,j) ---*/
 
279
          
 
280
          for(int atom=0;atom<Molecule.num_atoms;atom++) {
 
281
            struct coordinates PC;
 
282
            PC.x = sp->P[i][j][0] - Molecule.centers[atom].x;
 
283
            PC.y = sp->P[i][j][1] - Molecule.centers[atom].y;
 
284
            PC.z = sp->P[i][j][2] - Molecule.centers[atom].z;
 
285
            AI_OSrecurs(AI0,PA,PB,PC,gam,am_i+1,am_j);
 
286
 
 
287
            double dB_pfac = 0.5 * over_pf * (-Molecule.centers[atom].Z_nuc);
 
288
 
 
289
            int ai = 0;
 
290
            for(int ii = 0; ii <= am_i; ii++){
 
291
              int l1 = am_i - ii;
 
292
              for(int jj = 0; jj <= ii; jj++){
 
293
                int m1 = ii - jj;
 
294
                int n1 = jj ;
 
295
                int iind_000 = n1*izm + m1*iym + l1*ixm;
 
296
                int iind_100 = n1*izm + m1*iym + (l1+1)*ixm;
 
297
                int iind_010 = n1*izm + (m1+1)*iym + l1*ixm;
 
298
                int iind_001 = (n1+1)*izm + m1*iym + l1*ixm;
 
299
                /*--- create all am components of sj ---*/
 
300
                int aj = 0;
 
301
                for(int kk = 0; kk <= am_j; kk++){
 
302
                  int l2 = am_j - kk;
 
303
                  for(int ll = 0; ll <= kk; ll++){
 
304
                    int m2 = kk - ll;
 
305
                    int n2 = ll ;
 
306
 
 
307
                    int jind = n2*jzm + m2*jym + l2*jxm;
 
308
 
 
309
                    double V = AI0[iind_000][jind][0];
 
310
                    double xV = AI0[iind_100][jind][0] + A.x * V;
 
311
                    double yV = AI0[iind_010][jind][0] + A.y * V;
 
312
                    double zV = AI0[iind_001][jind][0] + A.z * V;
 
313
 
 
314
                    dHdBxtemp[ai][aj] += dB_pfac * ( AB.y * zV - AB.z * yV);
 
315
                    dHdBytemp[ai][aj] += dB_pfac * ( AB.z * xV - AB.x * zV);
 
316
                    dHdBztemp[ai][aj] += dB_pfac * ( AB.x * yV - AB.y * xV);
 
317
 
 
318
/*                    dHdBxtemp[ai][aj] += dB_pfac * xV;
 
319
                    dHdBytemp[ai][aj] += dB_pfac * yV;
 
320
                    dHdBztemp[ai][aj] += dB_pfac * zV;*/
 
321
                  
 
322
                    aj++;
 
323
                  }
 
324
                }  
 
325
                ai++;
 
326
              }
 
327
            } /*--- end cartesian components for (si,sj) with primitives (i,j) ---*/
 
328
          }
 
329
          
 
330
        }
 
331
      } /*--- end primitive contraction ---*/
 
332
    
 
333
      /*--- Normalize the contracted integrals ---*/
 
334
      double* ptr1 = GTOs.bf_norm[am_i];
 
335
      double* ptr2 = GTOs.bf_norm[am_j];
 
336
      for(int i=0; i<ni; i++) {
 
337
        double norm1 = ptr1[i];
 
338
        for(int j=0; j<nj; j++) {
 
339
          double norm12 = norm1*ptr2[j];
 
340
          d2HdBxdExtemp[i][j] *= norm12;
 
341
          d2HdBxdEytemp[i][j] *= norm12;
 
342
          d2HdBxdEztemp[i][j] *= norm12;
 
343
          d2HdBydExtemp[i][j] *= norm12;
 
344
          d2HdBydEytemp[i][j] *= norm12;
 
345
          d2HdBydEztemp[i][j] *= norm12;
 
346
          d2HdBzdExtemp[i][j] *= norm12;
 
347
          d2HdBzdEytemp[i][j] *= norm12;
 
348
          d2HdBzdEztemp[i][j] *= norm12;
 
349
          dHdBxtemp[i][j] *= norm12;
 
350
          dHdBytemp[i][j] *= norm12;
 
351
          dHdBztemp[i][j] *= norm12;
 
352
          dSdBxtemp[i][j] *= norm12;
 
353
          dSdBytemp[i][j] *= norm12;
 
354
          dSdBztemp[i][j] *= norm12;
 
355
        }
 
356
      }
 
357
 
 
358
      for(int i=0;i<ni;i++)
 
359
        for(int j=0;j<nj;j++) {
 
360
          int ii = ioffset + i;
 
361
          int jj = joffset + j;
 
362
          D2HDBXDEX[ii][jj] = d2HdBxdExtemp[i][j];
 
363
          D2HDBXDEY[ii][jj] = d2HdBxdEytemp[i][j];
 
364
          D2HDBXDEZ[ii][jj] = d2HdBxdEztemp[i][j];
 
365
          D2HDBYDEX[ii][jj] = d2HdBydExtemp[i][j];
 
366
          D2HDBYDEY[ii][jj] = d2HdBydEytemp[i][j];
 
367
          D2HDBYDEZ[ii][jj] = d2HdBydEztemp[i][j];
 
368
          D2HDBZDEX[ii][jj] = d2HdBzdExtemp[i][j];
 
369
          D2HDBZDEY[ii][jj] = d2HdBzdEytemp[i][j];
 
370
          D2HDBZDEZ[ii][jj] = d2HdBzdEztemp[i][j];
 
371
          DHDBX[ii][jj] = dHdBxtemp[i][j];
 
372
          DHDBY[ii][jj] = dHdBytemp[i][j];
 
373
          DHDBZ[ii][jj] = dHdBztemp[i][j];
 
374
          DSDBX[ii][jj] = dSdBxtemp[i][j];
 
375
          DSDBY[ii][jj] = dSdBytemp[i][j];
 
376
          DSDBZ[ii][jj] = dSdBztemp[i][j];
 
377
        }
 
378
    }
 
379
  }  /*--- This shell pair is done ---*/
 
380
 
 
381
  /*--- flush it all away ---*/
 
382
  fflush(outfile);
 
383
  free_block(OIX);
 
384
  free_block(OIY);
 
385
  free_block(OIZ);
 
386
  free_box(AI0,indmax_p1,indmax_0);
 
387
  dimension=BasisSet.num_ao*BasisSet.num_ao;
 
388
  iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_D2HDBDE_XX, dimension,D2HDBXDEX[0]);
 
389
  iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_D2HDBDE_XY, dimension,D2HDBXDEY[0]);
 
390
  iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_D2HDBDE_XZ, dimension,D2HDBXDEZ[0]);
 
391
  iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_D2HDBDE_YX, dimension,D2HDBYDEX[0]);
 
392
  iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_D2HDBDE_YY, dimension,D2HDBYDEY[0]);
 
393
  iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_D2HDBDE_YZ, dimension,D2HDBYDEZ[0]);
 
394
  iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_D2HDBDE_ZX, dimension,D2HDBZDEX[0]);
 
395
  iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_D2HDBDE_ZY, dimension,D2HDBZDEY[0]);
 
396
  iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_D2HDBDE_ZZ, dimension,D2HDBZDEZ[0]);
 
397
  iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_DHDB_X, dimension,DHDBX[0]);
 
398
  iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_DHDB_Y, dimension,DHDBY[0]);
 
399
  iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_DHDB_Z, dimension,DHDBZ[0]);
 
400
  iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_DSDB_X, dimension,DSDBX[0]);
 
401
  iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_DSDB_Y, dimension,DSDBY[0]);
 
402
  iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_DSDB_Z, dimension,DSDBZ[0]);
 
403
  if (UserOptions.print_lvl >= PRINT_OEI) {
 
404
    fprintf(outfile,"  -dS/dB_x AO integrals:\n\n");
 
405
    print_mat(DSDBX,BasisSet.num_ao,BasisSet.num_ao,outfile);
 
406
    fprintf(outfile,"  -dS/dB_y AO integrals:\n\n");
 
407
    print_mat(DSDBY,BasisSet.num_ao,BasisSet.num_ao,outfile);
 
408
    fprintf(outfile,"  -dS/dB_z AO integrals:\n\n");
 
409
    print_mat(DSDBZ,BasisSet.num_ao,BasisSet.num_ao,outfile);
 
410
    fprintf(outfile,"  -dh/dB_x AO integrals:\n\n");
 
411
    print_mat(DHDBX,BasisSet.num_ao,BasisSet.num_ao,outfile);
 
412
    fprintf(outfile,"  -dh/dB_y AO integrals:\n\n");
 
413
    print_mat(DHDBY,BasisSet.num_ao,BasisSet.num_ao,outfile);
 
414
    fprintf(outfile,"  -dh/dB_z AO integrals:\n\n");
 
415
    print_mat(DHDBZ,BasisSet.num_ao,BasisSet.num_ao,outfile);
 
416
/*    fprintf(outfile,"  -d2h/dB_x dE_x AO integrals:\n\n");
 
417
    print_mat(D2HDBXDEX,BasisSet.num_ao,BasisSet.num_ao,outfile);
 
418
    fprintf(outfile,"  -d2h/dB_x dE_y AO integrals:\n\n");
 
419
    print_mat(D2HDBXDEY,BasisSet.num_ao,BasisSet.num_ao,outfile);
 
420
    fprintf(outfile,"  -d2h/dB_x dE_z AO integrals:\n\n");
 
421
    print_mat(D2HDBXDEZ,BasisSet.num_ao,BasisSet.num_ao,outfile);
 
422
    fprintf(outfile,"  -d2h/dB_y dE_x AO integrals:\n\n");
 
423
    print_mat(D2HDBYDEX,BasisSet.num_ao,BasisSet.num_ao,outfile);
 
424
    fprintf(outfile,"  -d2h/dB_y dE_y AO integrals:\n\n");
 
425
    print_mat(D2HDBYDEY,BasisSet.num_ao,BasisSet.num_ao,outfile);
 
426
    fprintf(outfile,"  -d2h/dB_y dE_z AO integrals:\n\n");
 
427
    print_mat(D2HDBYDEZ,BasisSet.num_ao,BasisSet.num_ao,outfile);
 
428
    fprintf(outfile,"  -d2h/dB_z dE_x AO integrals:\n\n");
 
429
    print_mat(D2HDBZDEX,BasisSet.num_ao,BasisSet.num_ao,outfile);
 
430
    fprintf(outfile,"  -d2h/dB_z dE_y AO integrals:\n\n");
 
431
    print_mat(D2HDBZDEY,BasisSet.num_ao,BasisSet.num_ao,outfile);
 
432
    fprintf(outfile,"  -d2h/dB_z dE_z AO integrals:\n\n");
 
433
    print_mat(D2HDBZDEZ,BasisSet.num_ao,BasisSet.num_ao,outfile);*/
 
434
    fprintf(outfile,"\n");
 
435
  }
 
436
 
 
437
  free_block(d2HdBxdExtemp);
 
438
  free_block(d2HdBxdEytemp);
 
439
  free_block(d2HdBxdEztemp);
 
440
  free_block(d2HdBydExtemp);
 
441
  free_block(d2HdBydEytemp);
 
442
  free_block(d2HdBydEztemp);
 
443
  free_block(d2HdBzdExtemp);
 
444
  free_block(d2HdBzdEytemp);
 
445
  free_block(d2HdBzdEztemp);
 
446
  free_block(dHdBxtemp);
 
447
  free_block(dHdBytemp);
 
448
  free_block(dHdBztemp);
 
449
  free_block(dSdBxtemp);
 
450
  free_block(dSdBytemp);
 
451
  free_block(dSdBztemp);
 
452
 
 
453
  free_block(D2HDBXDEX);
 
454
  free_block(D2HDBXDEY);
 
455
  free_block(D2HDBXDEZ);
 
456
  free_block(D2HDBYDEX);
 
457
  free_block(D2HDBYDEY);
 
458
  free_block(D2HDBYDEZ);
 
459
  free_block(D2HDBZDEX);
 
460
  free_block(D2HDBZDEY);
 
461
  free_block(D2HDBZDEZ);
 
462
  free_block(DHDBX);
 
463
  free_block(DHDBY);
 
464
  free_block(DHDBZ);
 
465
  free_block(DSDBX);
 
466
  free_block(DSDBY);
 
467
  free_block(DSDBZ);
 
468
 
 
469
#ifdef USE_TAYLOR_FM
 
470
  free_Taylor_Fm_Eval();
 
471
#endif
 
472
 
 
473
  return;
 
474
}