← Back to branch summary

~ubuntu-branches/ubuntu/precise/psicode/precise

« back to all changes in this revision

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