← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/lib/libdpd/cc3_sigma_RHF.c

  • 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
http://bugs.debian.org/464867
* 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:
src/lib/libdpd/cc3_sigma_RHF.c
 
1
#include <stdio.h>
 
2
#include <stdlib.h>
 
3
#include <string.h>
 
4
#include <math.h>
 
5
#include <libciomr/libciomr.h>
 
6
#include <libdpd/dpd.h>
 
7
#include <libqt/qt.h>
 
8
#include <ccfiles.h>
 
9
 
 
10
/*
 
11
  This function computes contributions to singles and doubles of
 
12
  matrix elements of triples:
 
13
    SIA   <-- <S|(Dints)           <T|(Wmbij,Wabei) CMNEF |0> |T> / (w-wt)
 
14
    SIjAb <-- <D|(FME,WmAEf,WMnIe) <T|(Wmbij,Wabei) CMNEF |0> |T> / (w-wt)
 
15
  Irrep variables are:
 
16
    GS <--            GW                  GWX3    ^ GC
 
17
                                  (           GX3            )
 
18
  These are used to make X3 quantity in T3_RHF:
 
19
    CIjAb, WAbEi, WMbIj, fIJ2, fAB2, omega
 
20
*/
 
21
 
 
22
void cc3_sigma_RHF(dpdbuf4 *CIjAb, dpdbuf4 *WAbEi, dpdbuf4 *WMbIj,
 
23
    int do_singles, dpdbuf4 *Dints, dpdfile2 *SIA,
 
24
    int do_doubles, dpdfile2 *FME, dpdbuf4 *WmAEf, dpdbuf4 *WMnIe,
 
25
    dpdbuf4 *SIjAb, int *occpi, int *occ_off, int *virtpi, int *vir_off,
 
26
    double omega, FILE *outfile)
 
27
{
 
28
  int h, nirreps;
 
29
  int Gi, Gj, Gk, Gl, Ga, Gb, Gc, Gd;
 
30
  int i, j, k, l, a, b, c, d;
 
31
  int I, J, K, L, A, B, C, D;
 
32
  int kj, jk, ji, ij, ik, ki;
 
33
  int Gkj, Gjk, Gji, Gij, Gik, Gki, Gkd;
 
34
  int Gijk, GS, GC, GWX3, GW, GX3, nrows,ncols;
 
35
  int ab, ba, ac, ca, bc, cb;
 
36
  int Gab, Gba, Gac, Gca, Gbc, Gcb, Gid, Gjd;
 
37
  int id, jd, kd, ad, bd, cd;
 
38
  int il, jl, kl, la, lb, lc, li, lk;
 
39
  int da, di, dj, dk;
 
40
  int Gad, Gdi, Gdj, Gdk, Glc, Gli, Glk,cnt;
 
41
  int nlinks;
 
42
  double value, F_val, t_val, E_val;
 
43
  double dijk, denom, *tvect, **Z;
 
44
  double value_ia, value_ka, denom_ia, denom_ka;
 
45
  dpdfile2 fIJ, fIJ2, fAB, fAB2, SIA_inc;
 
46
  dpdbuf4 SIjAb_inc, buf4_tmp;
 
47
  double ***W3, ***W3a;
 
48
 
 
49
  nirreps = CIjAb->params->nirreps;
 
50
 
 
51
  /* these are sent to T3 function */
 
52
  dpd_file2_init(&fIJ2, CC_OEI, 0, 0, 0, "fIJ");
 
53
  dpd_file2_init(&fAB2, CC_OEI, 0, 1, 1, "fAB");
 
54
  dpd_file2_mat_init(&fIJ2);
 
55
  dpd_file2_mat_init(&fAB2);
 
56
  dpd_file2_mat_rd(&fIJ2);
 
57
  dpd_file2_mat_rd(&fAB2);
 
58
 
 
59
  dpd_file2_init(&fIJ, CC_OEI, 0, 0, 0, "fIJ");
 
60
  dpd_file2_init(&fAB, CC_OEI, 0, 1, 1, "fAB");
 
61
  dpd_file2_mat_init(&fIJ);
 
62
  dpd_file2_mat_init(&fAB);
 
63
  dpd_file2_mat_rd(&fIJ);
 
64
  dpd_file2_mat_rd(&fAB);
 
65
  dpd_file2_mat_init(FME);
 
66
  dpd_file2_mat_rd(FME);
 
67
 
 
68
  GC = CIjAb->file.my_irrep;
 
69
  GWX3 = WAbEi->file.my_irrep;
 
70
  GX3 = GWX3^GC;
 
71
  GW = WmAEf->file.my_irrep;
 
72
  GS = SIjAb->file.my_irrep;
 
73
  if (GS != (GX3^GW)) {
 
74
    fprintf(outfile,"problem with irreps in cc3_sigma_RHF()\n"); 
 
75
    exit(1);
 
76
  }
 
77
  if (do_singles) {
 
78
    dpd_file2_init(&SIA_inc, CC_TMP0, GS, 0, 1, "CC3 SIA");  /* T3->S1 increment */
 
79
    dpd_file2_mat_init(&SIA_inc);
 
80
  }
 
81
  if (do_doubles) {
 
82
    dpd_buf4_init(&SIjAb_inc, CC_TMP0, GS, 0, 5, 0, 5, 0, "CC3 SIjAb");
 
83
    dpd_buf4_scm(&SIjAb_inc, 0);
 
84
  }
 
85
 
 
86
  for(h=0; h < nirreps; h++) {
 
87
    dpd_buf4_mat_irrep_init(CIjAb, h);
 
88
    dpd_buf4_mat_irrep_rd(CIjAb, h);
 
89
    dpd_buf4_mat_irrep_init(WMbIj, h);
 
90
    dpd_buf4_mat_irrep_rd(WMbIj, h);
 
91
 
 
92
    if (do_singles) {
 
93
      dpd_buf4_mat_irrep_init(Dints, h);
 
94
      dpd_buf4_mat_irrep_rd(Dints, h);
 
95
    }
 
96
    if (do_doubles) {
 
97
      dpd_buf4_mat_irrep_init(WMnIe, h);
 
98
      dpd_buf4_mat_irrep_rd(WMnIe, h);
 
99
    }
 
100
    dpd_buf4_mat_irrep_init(&SIjAb_inc, h);
 
101
  }
 
102
 
 
103
  W3 = (double ***) malloc(nirreps * sizeof(double **));
 
104
  W3a = (double ***) malloc(nirreps * sizeof(double **));
 
105
 
 
106
  for(Gi=0; Gi < nirreps; Gi++) {
 
107
    for(Gj=0; Gj < nirreps; Gj++) {
 
108
      Gij = Gji = Gi ^ Gj; 
 
109
      for(Gk=0; Gk < nirreps; Gk++) {
 
110
        Gkj = Gjk = Gk ^ Gj;
 
111
        Gik = Gki = Gi ^ Gk;
 
112
        Gijk = Gi ^ Gj ^ Gk;
 
113
        
 
114
        /* allocate memory for all irrep blocks of (ab,c) */
 
115
        for(Gab=0; Gab < nirreps; Gab++) {
 
116
          Gc = Gab ^ Gijk ^ GX3;
 
117
          W3[Gab] = dpd_block_matrix(WAbEi->params->coltot[Gab], virtpi[Gc]);
 
118
        }
 
119
        for(Ga=0; Ga < nirreps; Ga++) {
 
120
          Gbc = Ga ^ Gijk ^ GX3;
 
121
          W3a[Ga] = dpd_block_matrix(virtpi[Ga], WAbEi->params->coltot[Gbc]);
 
122
        }
 
123
        
 
124
        for(i=0; i < occpi[Gi]; i++) {
 
125
          I = occ_off[Gi] + i;
 
126
          for(j=0; j < occpi[Gj]; j++) {
 
127
            J = occ_off[Gj] + j;
 
128
            for(k=0; k < occpi[Gk]; k++) {
 
129
              K = occ_off[Gk] + k;
 
130
              
 
131
              ij = CIjAb->params->rowidx[I][J];
 
132
              ji = CIjAb->params->rowidx[J][I];
 
133
              ik = CIjAb->params->rowidx[I][K];
 
134
              ki = CIjAb->params->rowidx[K][I];
 
135
              jk = CIjAb->params->rowidx[J][K];
 
136
              kj = CIjAb->params->rowidx[K][J];
 
137
#ifdef T3_TIMER_ON              
 
138
  timer_on("T3_RHF");
 
139
#endif
 
140
              T3_RHF(W3, nirreps, I, Gi, J, Gj, K, Gk, CIjAb, WAbEi, WMbIj, 
 
141
                     &fIJ2, &fAB2, occpi, occ_off, virtpi, vir_off, omega);
 
142
#ifdef T3_TIMER_ON              
 
143
  timer_off("T3_RHF");
 
144
#endif
 
145
 
 
146
              /* do (Wmnie*X3(ab,c) --> SIjAb) contractions that use W3(ab,c) first */
 
147
              if (do_doubles) {
 
148
#ifdef T3_TIMER_ON              
 
149
  timer_on("X3*Wmnie");
 
150
#endif
 
151
                /* S_liba <-- -2.0 * t_ijkabc W_jklc */
 
152
                /* S_ilab <-- -2.0 * t_ijkabc W_jklc */
 
153
                for(Gl=0; Gl < nirreps; Gl++) {
 
154
                  Gli = Gl ^ Gi;
 
155
                  Gab = Gli ^ GS;
 
156
                  Gc = Gjk ^ Gl ^ GW;
 
157
 
 
158
                  lc = WMnIe->col_offset[Gjk][Gl];
 
159
                  nrows = SIjAb_inc.params->coltot[Gab];
 
160
                  ncols = occpi[Gl];
 
161
                  nlinks = virtpi[Gc];
 
162
 
 
163
                  if(nrows && ncols && nlinks) {
 
164
                    Z = dpd_block_matrix(nrows, ncols);
 
165
 
 
166
                    C_DGEMM('n', 't', nrows, ncols, nlinks, 1.0, W3[Gab][0], nlinks,
 
167
                        &(WMnIe->matrix[Gjk][jk][lc]), nlinks, 0.0, Z[0], ncols);
 
168
 
 
169
                    for(l=0; l < ncols; l++) {
 
170
                      L = occ_off[Gl] + l;
 
171
                      li = SIjAb_inc.params->rowidx[L][I];
 
172
                      il = SIjAb_inc.params->rowidx[I][L];
 
173
                      for(ab=0; ab < nrows; ab++) {
 
174
                        A = SIjAb_inc.params->colorb[Gab][ab][0];
 
175
                        B = SIjAb_inc.params->colorb[Gab][ab][1];
 
176
                        ba = SIjAb_inc.params->colidx[B][A];
 
177
                        SIjAb_inc.matrix[Gli][li][ba] -= 2.0 * Z[ab][l];
 
178
                        SIjAb_inc.matrix[Gli][il][ab] -= 2.0 * Z[ab][l];
 
179
                      }
 
180
                    }
 
181
                    dpd_free_block(Z, nrows, ncols);
 
182
                  }
 
183
                }
 
184
 
 
185
                /* S_lkba <-- + t_ijkabc W_jilc */
 
186
                /* S_klab <-- + t_ijkabc W_jilc */
 
187
                for(Gl=0; Gl < nirreps; Gl++) {
 
188
                  Glk = Gl ^ Gk;
 
189
                  Gab = Glk ^ GS;
 
190
                  Gc = Gji ^ Gl ^ GW;
 
191
 
 
192
                  lc = WMnIe->col_offset[Gji][Gl];
 
193
                  nrows = SIjAb_inc.params->coltot[Gab];
 
194
                  ncols = occpi[Gl];
 
195
                  nlinks = virtpi[Gc];
 
196
 
 
197
                  if(nrows && ncols && nlinks) {
 
198
                    Z = dpd_block_matrix(nrows, ncols);
 
199
 
 
200
                    C_DGEMM('n', 't', nrows, ncols, nlinks, 1.0, W3[Gab][0], nlinks,
 
201
                        &(WMnIe->matrix[Gji][ji][lc]), nlinks, 0.0, Z[0], ncols);
 
202
 
 
203
                    for(l=0; l < ncols; l++) {
 
204
                      L = occ_off[Gl] + l;
 
205
                      lk = SIjAb_inc.params->rowidx[L][K];
 
206
                      kl = SIjAb_inc.params->rowidx[K][L];
 
207
                      for(ab=0; ab < nrows; ab++) {
 
208
                        A = SIjAb_inc.params->colorb[Gab][ab][0];
 
209
                        B = SIjAb_inc.params->colorb[Gab][ab][1];
 
210
                        ba = SIjAb_inc.params->colidx[B][A];
 
211
                        SIjAb_inc.matrix[Glk][lk][ba] += Z[ab][l];
 
212
                        SIjAb_inc.matrix[Glk][kl][ab] += Z[ab][l];
 
213
                      }
 
214
                    }
 
215
                    dpd_free_block(Z, nrows, ncols);
 
216
                  }
 
217
                }
 
218
 
 
219
                /* S_liba <-- + t_ijkabc W_kjlc */
 
220
                /* S_ilab <-- + t_ijkabc W_kjlc */
 
221
                for(Gl=0; Gl < nirreps; Gl++) {
 
222
                  Gli = Gl ^ Gi;
 
223
                  Gab = Gli ^ GS;
 
224
                  Gc = Gjk ^ Gl ^ GW;
 
225
 
 
226
                  lc = WMnIe->col_offset[Gkj][Gl];
 
227
                  nrows = SIjAb_inc.params->coltot[Gab];
 
228
                  ncols = occpi[Gl];
 
229
                  nlinks = virtpi[Gc];
 
230
 
 
231
                  if(nrows && ncols && nlinks) {
 
232
                    Z = dpd_block_matrix(nrows, ncols);
 
233
 
 
234
                    C_DGEMM('n', 't', nrows, ncols, nlinks, 1.0, W3[Gab][0], nlinks,
 
235
                        &(WMnIe->matrix[Gkj][kj][lc]), nlinks, 0.0, Z[0], ncols);
 
236
 
 
237
                    for(l=0; l < ncols; l++) {
 
238
                      L = occ_off[Gl] + l;
 
239
                      li = SIjAb_inc.params->rowidx[L][I];
 
240
                      il = SIjAb_inc.params->rowidx[I][L];
 
241
                      for(ab=0; ab < nrows; ab++) {
 
242
                        A = SIjAb_inc.params->colorb[Gab][ab][0];
 
243
                        B = SIjAb_inc.params->colorb[Gab][ab][1];
 
244
                        ba = SIjAb_inc.params->colidx[B][A];
 
245
                        SIjAb_inc.matrix[Gli][li][ba] += Z[ab][l];
 
246
                        SIjAb_inc.matrix[Gli][il][ab] += Z[ab][l];
 
247
                      }
 
248
                    }
 
249
                    dpd_free_block(Z, nrows, ncols);
 
250
                  }
 
251
                }
 
252
#ifdef T3_TIMER_ON              
 
253
  timer_off("X3*Wmnie");
 
254
#endif
 
255
              } /* end Wmnie*X3 doubles contributions */
 
256
 
 
257
#ifdef T3_TIMER_ON              
 
258
  timer_on("X3_sort");
 
259
#endif
 
260
              /* sort W(ab,c) to W(a,bc) */
 
261
              for(Gab=0; Gab < nirreps; Gab++) {
 
262
                Gc = Gab ^ Gijk ^ GX3;
 
263
                for(ab=0; ab < WAbEi->params->coltot[Gab]; ab++ ){
 
264
                  A = WAbEi->params->colorb[Gab][ab][0];
 
265
                  B = WAbEi->params->colorb[Gab][ab][1];
 
266
                  Ga = WAbEi->params->rsym[A];
 
267
                  Gb = Gab^Ga;
 
268
                  a = A - vir_off[Ga];
 
269
                  for(c=0; c < virtpi[Gc]; c++) {
 
270
                    C = vir_off[Gc] + c;
 
271
                    bc = WAbEi->params->colidx[B][C];
 
272
                    W3a[Ga][a][bc] = W3[Gab][ab][c];
 
273
                  }
 
274
                }
 
275
              }
 
276
#ifdef T3_TIMER_ON              
 
277
  timer_off("X3_sort");
 
278
#endif
 
279
 
 
280
              /*** X3(a,bc)*Dints --> SIA Contributions ***/
 
281
              if (do_singles) {
 
282
                /* Note: the do_singles code has not been used where Dints!=A1
 
283
                   as this non-symmetric case does not arise for CC3 EOM energies */
 
284
 
 
285
                /* S_ia <-- t_ijkabc Djkbc */
 
286
                Ga = Gi ^ GS;
 
287
                Gbc = Gjk ^ GW;
 
288
                nrows = virtpi[Ga];
 
289
                ncols = Dints->params->coltot[Gbc];
 
290
 
 
291
                if(nrows && ncols)
 
292
                  C_DGEMV('n', nrows, ncols, 1.0, W3a[Ga][0], ncols,
 
293
                    &(Dints->matrix[Gjk][jk][0]), 1, 1.0, SIA_inc.matrix[Gi][i], 1);
 
294
 
 
295
                /* S_ka <-- tijkabc Djibc */
 
296
                Ga = Gk ^ GS;
 
297
                Gbc = Gji ^ GW;
 
298
                nrows = virtpi[Ga];
 
299
                ncols = Dints->params->coltot[Gbc];
 
300
 
 
301
                if(nrows && ncols)
 
302
                  C_DGEMV('n', nrows, ncols, -1.0, W3a[Ga][0], ncols,
 
303
                  &(Dints->matrix[Gji][ji][0]), 1, 1.0, SIA_inc.matrix[Gk][k], 1);
 
304
              }
 
305
 
 
306
              /*** X3(a,bc)*Wamef --> SIjAb Contributions ***/
 
307
              if (do_doubles) {
 
308
#ifdef T3_TIMER_ON              
 
309
  timer_on("X3*Wamef");
 
310
#endif
 
311
                /* S_IjAb <-- t_ijkabc F_kc */
 
312
                /* S_jIbA <-- t_ijkabc F_kc */
 
313
                Gc = Gk ^ GW;
 
314
                Gab = Gij ^ GS;
 
315
                nrows = SIjAb_inc.params->coltot[Gij^GS];
 
316
                ncols = virtpi[Gc];
 
317
                
 
318
                if(nrows && ncols) {
 
319
                  tvect = init_array(nrows);
 
320
                  C_DGEMV('n', nrows, ncols, 1.0, W3[Gab][0], ncols, FME->matrix[Gk][k], 1,
 
321
                    0.0, tvect, 1);
 
322
 
 
323
                  for(cnt=0; cnt<nrows; ++cnt) {
 
324
                    A = SIjAb_inc.params->colorb[Gab][cnt][0];
 
325
                    B = SIjAb_inc.params->colorb[Gab][cnt][1];
 
326
                    ba = SIjAb_inc.params->colidx[B][A];
 
327
                    SIjAb_inc.matrix[Gij][ij][cnt] += tvect[cnt];
 
328
                    SIjAb_inc.matrix[Gij][ji][ba] += tvect[cnt];
 
329
                  }
 
330
                  free(tvect);
 
331
                }
 
332
 
 
333
                /* S_KjAb <-- - t_ijkabc F_ic */
 
334
                /* S_jKbA <-- - t_ijkabc F_ic */
 
335
                Gc = Gi ^ GW;
 
336
                Gab = Gkj ^ GS;
 
337
                nrows = SIjAb_inc.params->coltot[Gkj^GS];
 
338
                ncols = virtpi[Gc];
 
339
                
 
340
                if(nrows && ncols) {
 
341
                  tvect = init_array(nrows);
 
342
                  C_DGEMV('n', nrows, ncols, 1.0, W3[Gab][0], ncols, FME->matrix[Gi][i], 1,
 
343
                      0.0, tvect, 1);
 
344
 
 
345
                  for(cnt=0; cnt<nrows; ++cnt) {
 
346
                    A = SIjAb_inc.params->colorb[Gab][cnt][0];
 
347
                    B = SIjAb_inc.params->colorb[Gab][cnt][1];
 
348
                    ba = SIjAb_inc.params->colidx[B][A];
 
349
                    SIjAb_inc.matrix[Gkj][kj][cnt] -= tvect[cnt];
 
350
                    SIjAb_inc.matrix[Gkj][jk][ba] -= tvect[cnt];
 
351
                  }
 
352
                  free(tvect);
 
353
                }
 
354
 
 
355
                /* S_ijad <-- 2.0 * t_ijkabc W_kdbc */
 
356
                /* S_jida <-- 2.0 * t_ijkabc W_kdbc */
 
357
                for(Gd=0; Gd < nirreps; Gd++) {
 
358
                  Gkd = Gk ^ Gd;
 
359
                  Ga = Gd ^ Gij ^ GS;
 
360
 
 
361
                  nrows = virtpi[Ga];
 
362
                  ncols = virtpi[Gd];
 
363
                  nlinks = WmAEf->params->coltot[Gkd^GW];
 
364
 
 
365
                  if(nrows && ncols && nlinks) {
 
366
                    kd = WmAEf->row_offset[Gkd][K];
 
367
                    WmAEf->matrix[Gkd] = dpd_block_matrix(virtpi[Gd], WmAEf->params->coltot[Gkd^GW]);
 
368
                    dpd_buf4_mat_irrep_rd_block(WmAEf, Gkd, kd, virtpi[Gd]);
 
369
 
 
370
                    Z = block_matrix(virtpi[Ga], virtpi[Gd]);
 
371
                    C_DGEMM('n', 't', nrows, ncols, nlinks, 1.0, W3a[Ga][0], nlinks,
 
372
                      WmAEf->matrix[Gkd][0], nlinks, 1.0, Z[0], ncols);
 
373
 
 
374
                    for(a=0; a < virtpi[Ga]; a++) {
 
375
                      A = vir_off[Ga] + a;
 
376
                      for(d=0; d < virtpi[Gd]; d++) {
 
377
                        D = vir_off[Gd] + d;
 
378
                        ad = SIjAb_inc.params->colidx[A][D];
 
379
                        da = SIjAb_inc.params->colidx[D][A];
 
380
                        SIjAb_inc.matrix[Gij][ij][ad] += 2.0 * Z[a][d];
 
381
                        SIjAb_inc.matrix[Gij][ji][da] += 2.0 * Z[a][d];
 
382
                      }
 
383
                    }
 
384
 
 
385
                    free_block(Z);
 
386
                    dpd_free_block(WmAEf->matrix[Gkd], virtpi[Gd], WmAEf->params->coltot[Gkd^GW]);
 
387
                  }
 
388
                }
 
389
 
 
390
                /* S_kjad <-- - t_ijkabc W_idbc */
 
391
                /* S_jkda <-- - t_ijkabc W_idbc */
 
392
                for(Gd=0; Gd < nirreps; Gd++) {
 
393
                  Gid = Gi ^ Gd;
 
394
                  Ga = Gd ^ Gjk ^ GS;
 
395
 
 
396
                  nrows = virtpi[Ga];
 
397
                  ncols = virtpi[Gd];
 
398
                  nlinks = WmAEf->params->coltot[Gid^GW];
 
399
 
 
400
                  if(nrows && ncols && nlinks) {
 
401
                    id = WmAEf->row_offset[Gid][I];
 
402
                    WmAEf->matrix[Gid] = dpd_block_matrix(virtpi[Gd], WmAEf->params->coltot[Gid^GW]);
 
403
                    dpd_buf4_mat_irrep_rd_block(WmAEf, Gid, id, virtpi[Gd]);
 
404
 
 
405
                    Z = block_matrix(virtpi[Ga], virtpi[Gd]);
 
406
                    C_DGEMM('n', 't', nrows, ncols, nlinks, 1.0, W3a[Ga][0], nlinks,
 
407
                      WmAEf->matrix[Gid][0], nlinks, 1.0, Z[0], ncols);
 
408
 
 
409
                    for(a=0; a < virtpi[Ga]; a++) {
 
410
                      A = vir_off[Ga] + a;
 
411
                      for(d=0; d < virtpi[Gd]; d++) {
 
412
                        D = vir_off[Gd] + d;
 
413
                        ad = SIjAb_inc.params->colidx[A][D];
 
414
                        da = SIjAb_inc.params->colidx[D][A];
 
415
                        SIjAb_inc.matrix[Gkj][kj][ad] -= Z[a][d];
 
416
                        SIjAb_inc.matrix[Gjk][jk][da] -= Z[a][d];
 
417
                      }
 
418
                    }
 
419
 
 
420
                    free_block(Z);
 
421
                    dpd_free_block(WmAEf->matrix[Gid], virtpi[Gd], WmAEf->params->coltot[Gid^GW]);
 
422
                  }
 
423
                }
 
424
 
 
425
                /* S_kida <-- - t_ijkabc W_jdbc */
 
426
                /* S_ikad <-- - t_ijkabc W_jdbc */
 
427
                for(Gd=0; Gd < nirreps; Gd++) {
 
428
                  Gjd = Gj ^ Gd;
 
429
                  Ga = Gd ^ Gik ^ GS;
 
430
 
 
431
                  nrows = virtpi[Ga];
 
432
                  ncols = virtpi[Gd];
 
433
                  nlinks = WmAEf->params->coltot[Gjd^GW];
 
434
 
 
435
                  if(nrows && ncols && nlinks) {
 
436
                    jd = WmAEf->row_offset[Gjd][J];
 
437
                    WmAEf->matrix[Gjd] = dpd_block_matrix(virtpi[Gd], WmAEf->params->coltot[Gjd^GW]);
 
438
                    dpd_buf4_mat_irrep_rd_block(WmAEf, Gjd, jd, virtpi[Gd]);
 
439
 
 
440
                    Z = block_matrix(virtpi[Ga], virtpi[Gd]);
 
441
                    C_DGEMM('n', 't', nrows, ncols, nlinks, 1.0, W3a[Ga][0], nlinks,
 
442
                      WmAEf->matrix[Gjd][0], nlinks, 1.0, Z[0], ncols);
 
443
 
 
444
                    for(a=0; a < virtpi[Ga]; a++) {
 
445
                      A = vir_off[Ga] + a;
 
446
                      for(d=0; d < virtpi[Gd]; d++) {
 
447
                        D = vir_off[Gd] + d;
 
448
                        ad = SIjAb_inc.params->colidx[A][D];
 
449
                        da = SIjAb_inc.params->colidx[D][A];
 
450
                        SIjAb_inc.matrix[Gki][ki][da] -= Z[a][d];
 
451
                        SIjAb_inc.matrix[Gik][ik][ad] -= Z[a][d];
 
452
                      }
 
453
                    }
 
454
                    free_block(Z);
 
455
                    dpd_free_block(WmAEf->matrix[Gjd], virtpi[Gd], WmAEf->params->coltot[Gjd^GW]);
 
456
                  }
 
457
                }
 
458
#ifdef T3_TIMER_ON              
 
459
timer_off("X3*Wamef");
 
460
#endif
 
461
              } /* end Wamef*X3->Sijab contributions */
 
462
            } /* k */
 
463
          } /* j */
 
464
         } /* i */
 
465
 
 
466
        for(Gab=0; Gab < nirreps; Gab++) {
 
467
          Gc = Gab ^ Gijk ^ GX3;
 
468
          dpd_free_block(W3[Gab], WAbEi->params->coltot[Gab], virtpi[Gc]);
 
469
        }
 
470
        for(Ga=0; Ga < nirreps; Ga++) {
 
471
          Gbc = Ga ^ Gijk ^ GX3;
 
472
          dpd_free_block(W3a[Ga], virtpi[Ga], WAbEi->params->coltot[Gbc]);
 
473
        }
 
474
      } /* Gk */
 
475
    } /* Gj */
 
476
  } /* Gi */
 
477
 
 
478
  /* close up files and update sigma vectors */
 
479
  dpd_file2_mat_close(&fIJ);
 
480
  dpd_file2_mat_close(&fAB);
 
481
  dpd_file2_mat_close(FME);
 
482
  dpd_file2_close(&fIJ);
 
483
  dpd_file2_close(&fAB);
 
484
 
 
485
  dpd_file2_mat_close(&fIJ2);
 
486
  dpd_file2_mat_close(&fAB2);
 
487
  dpd_file2_close(&fIJ2);
 
488
  dpd_file2_close(&fAB2);
 
489
 
 
490
  if (do_singles) {
 
491
    for(h=0; h < nirreps; h++)
 
492
      dpd_buf4_mat_irrep_close(Dints, h);
 
493
    dpd_file2_mat_wrt(&SIA_inc);
 
494
    dpd_file2_mat_close(&SIA_inc);
 
495
    dpd_file2_axpy(&SIA_inc, SIA, 1, 0);
 
496
    dpd_file2_close(&SIA_inc);
 
497
  }
 
498
 
 
499
  if (do_doubles) {
 
500
    for(h=0; h < nirreps; h++) {
 
501
      dpd_buf4_mat_irrep_close(WMnIe, h);
 
502
      dpd_buf4_mat_irrep_wrt(&SIjAb_inc, h);
 
503
      dpd_buf4_mat_irrep_close(&SIjAb_inc, h);
 
504
    }
 
505
    dpd_buf4_axpy(&SIjAb_inc, SIjAb, 1);
 
506
    dpd_buf4_close(&SIjAb_inc);
 
507
  }
 
508
}