← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cclambda/cc3_l3l2.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/bin/cclambda/cc3_l3l2.cc
 
1
/*! \file
 
2
    \ingroup CCLAMBDA
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#include <cstdio>
 
6
#include <cstdlib>
 
7
#include <cmath>
 
8
#include <libdpd/dpd.h>
 
9
#include <libqt/qt.h>
 
10
#include "MOInfo.h"
 
11
#include "Params.h"
 
12
#define EXTERN
 
13
#include "globals.h"
 
14
 
 
15
namespace psi { namespace cclambda {
 
16
 
 
17
void cc3_l3l2_RHF_AAA(void);
 
18
void cc3_l3l2_RHF_AAB(void);
 
19
 
 
20
void L3_AAA(double ***W1, int nirreps, int I, int Gi, int J, int Gj, int K, int Gk, 
 
21
            dpdbuf4 *T2, dpdbuf4 *F, dpdbuf4 *E, dpdfile2 *fIJ, dpdfile2 *fAB, 
 
22
            dpdbuf4 *D, dpdbuf4 *LIJAB, dpdfile2 *LIA, dpdfile2 *FME,
 
23
            int *occpi, int *occ_off, int *virtpi, int *vir_off);
 
24
 
 
25
void L3_AAB(double ***W1, int nirreps, int I, int Gi, int J, int Gj, int K, int Gk, 
 
26
            dpdbuf4 *T2AA, dpdbuf4 *T2AB, dpdbuf4 *T2BA, dpdbuf4 *FAA, dpdbuf4 *FAB, dpdbuf4 *FBA,
 
27
            dpdbuf4 *EAA, dpdbuf4 *EAB, dpdbuf4 *EBA, dpdfile2 *fIJ, dpdfile2 *fij, 
 
28
            dpdfile2 *fAB, dpdfile2 *fab, dpdbuf4 *DAA, dpdbuf4 *DAB, dpdbuf4 *LIJAB, dpdbuf4 *LIjAb,
 
29
            dpdfile2 *LIA, dpdfile2 *Lia, dpdfile2 *FME, dpdfile2 *Fme,
 
30
            int *aoccpi, int *aocc_off, int *boccpi, int *bocc_off,
 
31
            int *avirtpi, int *avir_off, int *bvirtpi, int *bvir_off);
 
32
 
 
33
void cc3_l3l2(void)
 
34
{
 
35
  if(params.ref == 0) {
 
36
    cc3_l3l2_RHF_AAA();
 
37
    cc3_l3l2_RHF_AAB();
 
38
  }
 
39
}
 
40
 
 
41
void cc3_l3l2_RHF_AAA(void)
 
42
{
 
43
  int h, nirreps;
 
44
  int *occ_off, *occpi;
 
45
  int *vir_off, *virtpi;
 
46
  int Gijk;
 
47
  int Gi, Gj, Gk;
 
48
  int Ga, Gb, Gc;
 
49
  int Gab, ab;
 
50
  int i, j, k, I, J, K;
 
51
  int a, b, c, A, B, C;
 
52
  double ***W1, ***W2;
 
53
  dpdbuf4 L, E, F;
 
54
  dpdfile2 fIJ, fAB;
 
55
  dpdfile2 FME, LIA;
 
56
  dpdbuf4 Dints, LIJAB;
 
57
  dpdbuf4 WMAFE, WMNIE;
 
58
  dpdbuf4 ZIGDE, T2;
 
59
  dpdbuf4 ZDMAE;
 
60
  dpdbuf4 ZLMAO;
 
61
  dpdbuf4 ZIMLE;
 
62
  dpdbuf4 L2new, L2, D2;
 
63
  int Gjk, jk, Gid, id, Gik, ik;
 
64
  int Gd, d, DD;
 
65
  int cd, dc;
 
66
  int Gm, m, M;
 
67
  int Gmi, mi, im, mc;
 
68
  int Gjd, jd;
 
69
  int Gij, ij, Gmk, mk, am, Gbc, bc;
 
70
  int ac, mb;
 
71
  int nrows, ncols, nlinks;
 
72
  double **Z;
 
73
 
 
74
  nirreps = moinfo.nirreps;
 
75
  occpi = moinfo.occpi;
 
76
  occ_off = moinfo.occ_off;
 
77
  virtpi = moinfo.virtpi;
 
78
  vir_off = moinfo.vir_off;
 
79
 
 
80
  dpd_buf4_init(&WMAFE, CC3_HET1, 0, 10, 5, 10, 7, 0, "CC3 WABEI (IE,B>A)");
 
81
  dpd_buf4_init(&WMNIE, CC3_HET1, 0, 0, 10, 2, 10, 0, "CC3 WMBIJ (I>J,MB)");
 
82
  for(h=0; h < nirreps; h++) {
 
83
    dpd_buf4_mat_irrep_init(&WMNIE, h);
 
84
    dpd_buf4_mat_irrep_rd(&WMNIE, h);
 
85
  }
 
86
 
 
87
  dpd_buf4_init(&L2new, CC3_MISC, 0, 0, 5, 0, 5, 0, "CC3 LIJAB");
 
88
  for(h=0; h < nirreps; h++) dpd_buf4_mat_irrep_init(&L2new, h);
 
89
 
 
90
  dpd_buf4_init(&ZIGDE, CC3_MISC, 0, 10, 5, 10, 5, 0, "CC3 ZIGDE");
 
91
  dpd_buf4_scm(&ZIGDE, 0.0); /* must be cleared in each iteration */
 
92
 
 
93
  dpd_buf4_init(&ZDMAE, CC3_MISC, 0, 10, 5, 10, 5, 0, "CC3 ZDMAE (MD,AE)");
 
94
  dpd_buf4_scm(&ZDMAE, 0.0);
 
95
 
 
96
  dpd_buf4_init(&ZLMAO, CC3_MISC, 0, 0, 11, 0, 11, 0, "CC3 ZLMAO");
 
97
  for(h=0; h < nirreps; h++) dpd_buf4_mat_irrep_init(&ZLMAO, h);
 
98
 
 
99
  dpd_buf4_init(&ZIMLE, CC3_MISC, 0, 0, 10, 0, 10, 0, "CC3 ZIMLE");
 
100
  for(h=0; h < nirreps; h++) dpd_buf4_mat_irrep_init(&ZIMLE, h);
 
101
 
 
102
  dpd_buf4_init(&T2, CC_TAMPS, 0, 0, 5, 2, 7, 0, "tIJAB");
 
103
  for(h=0; h < nirreps; h++) {
 
104
    dpd_buf4_mat_irrep_init(&T2, h);
 
105
    dpd_buf4_mat_irrep_rd(&T2, h);
 
106
  }
 
107
 
 
108
  dpd_file2_init(&fIJ, CC_OEI, 0, 0, 0, "fIJ");
 
109
  dpd_file2_init(&fAB, CC_OEI, 0, 1, 1, "fAB");
 
110
 
 
111
  dpd_buf4_init(&L, CC_LAMBDA, 0, 0, 5, 2, 7, 0, "LIJAB");
 
112
  dpd_buf4_init(&F, CC3_HET1, 0, 10, 5, 10, 7, 0, "CC3 WAMEF (MA,F>E)");
 
113
  dpd_buf4_init(&E, CC3_HET1, 0, 0, 10, 2, 10, 0, "CC3 WMNIE (M>N,IE)");
 
114
 
 
115
  dpd_buf4_init(&Dints, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij||ab>");
 
116
  dpd_buf4_init(&LIJAB, CC_LAMBDA, 0, 0, 5, 2, 7, 0, "LIJAB");
 
117
  dpd_file2_init(&LIA, CC_LAMBDA, 0, 0, 1, "LIA");
 
118
  dpd_file2_init(&FME, CC_OEI, 0, 0, 1, "FME");
 
119
 
 
120
  W1 = (double ***) malloc(nirreps * sizeof(double **));
 
121
  W2 = (double ***) malloc(nirreps * sizeof(double **));
 
122
 
 
123
  for(Gi=0; Gi < nirreps; Gi++) {
 
124
    for(Gj=0; Gj < nirreps; Gj++) {
 
125
      Gij = Gi ^ Gj;
 
126
      for(Gk=0; Gk < nirreps; Gk++) {
 
127
        Gijk = Gi ^ Gj ^ Gk;
 
128
 
 
129
        Gjk = Gj ^ Gk;
 
130
        Gik = Gi ^ Gk;
 
131
 
 
132
        for(Gab=0; Gab < nirreps; Gab++) {
 
133
          Gc = Gab ^ Gijk; /* totally symmetric */
 
134
          W1[Gab] = dpd_block_matrix(F.params->coltot[Gab], virtpi[Gc]);
 
135
        }
 
136
        for(Ga=0; Ga < nirreps; Ga++) {
 
137
          Gbc = Ga ^ Gijk; /* totally symmetric */
 
138
          W2[Ga] = dpd_block_matrix(virtpi[Ga], F.params->coltot[Gbc]);
 
139
        }
 
140
 
 
141
        for(i=0; i < occpi[Gi]; i++) {
 
142
          I = occ_off[Gi] + i;
 
143
          for(j=0; j < occpi[Gj]; j++) {
 
144
            J = occ_off[Gj] + j;
 
145
            for(k=0; k < occpi[Gk]; k++) {
 
146
              K = occ_off[Gk] + k;
 
147
 
 
148
              L3_AAA(W1, nirreps, I, Gi, J, Gj, K, Gk, &L, &F, &E, &fIJ, &fAB,
 
149
                     &Dints, &LIJAB, &LIA, &FME, occpi, occ_off, virtpi, vir_off);
 
150
 
 
151
              /* L_JKDC <-- +1/2 t_IJKABC W_ABID (IDAB) */
 
152
              /* L_JKCD <-- -1/2 t_IJKABC W_ABID (IDAB) */
 
153
              jk = L2new.params->rowidx[J][K];
 
154
              for(Gd=0; Gd < nirreps; Gd++) {
 
155
                Gab = Gid = Gi ^ Gd;  /* assumes Wieab is totally symmetric */
 
156
                Gc = Gab ^ Gijk;  /* assumes T3 is totally symmetric */
 
157
 
 
158
                id = WMAFE.row_offset[Gid][I];
 
159
 
 
160
                Z = block_matrix(virtpi[Gc],virtpi[Gd]);
 
161
                WMAFE.matrix[Gid] = dpd_block_matrix(virtpi[Gd], WMAFE.params->coltot[Gid]);
 
162
                dpd_buf4_mat_irrep_rd_block(&WMAFE, Gid, id, virtpi[Gd]);
 
163
 
 
164
                nrows = virtpi[Gc];
 
165
                ncols = virtpi[Gd];
 
166
                nlinks = WMAFE.params->coltot[Gid];
 
167
 
 
168
                if(nrows && ncols && nlinks)
 
169
                  C_DGEMM('t', 't', nrows, ncols, nlinks, 0.5, W1[Gab][0], nrows, 
 
170
                          WMAFE.matrix[Gid][0], nlinks, 0.0, Z[0], ncols);
 
171
 
 
172
                for(c=0; c < virtpi[Gc]; c++) {
 
173
                  C = vir_off[Gc] + c;
 
174
                  for(d=0; d < virtpi[Gd]; d++) {
 
175
                    DD = vir_off[Gd] + d;
 
176
                    cd = L2new.params->colidx[C][DD];
 
177
                    dc = L2new.params->colidx[DD][C];
 
178
                    L2new.matrix[Gjk][jk][dc] += Z[c][d];
 
179
                    L2new.matrix[Gjk][jk][cd] += -Z[c][d];
 
180
                  }
 
181
                }
 
182
                dpd_free_block(WMAFE.matrix[Gid], virtpi[Gd], WMAFE.params->coltot[Gid]);
 
183
                free_block(Z);
 
184
              }
 
185
 
 
186
              /* t_MIAB <-- +1/2 t_IJKABC W_JKMC */
 
187
              /* t_IMAB <-- -1/2 t_IJKABC W_JKMC */
 
188
              jk = WMNIE.params->rowidx[J][K];
 
189
              for(Gm=0; Gm < nirreps; Gm++) {
 
190
                Gab = Gmi = Gm ^ Gi;  /* assumes totally symmetric */
 
191
                Gc = Gab ^ Gijk;      /* assumes totally symmetric */
 
192
 
 
193
                mc = WMNIE.col_offset[Gjk][Gm];
 
194
 
 
195
                nrows = F.params->coltot[Gab];
 
196
                ncols = occpi[Gm];
 
197
                nlinks = virtpi[Gc];
 
198
 
 
199
                Z = dpd_block_matrix(nrows, ncols);
 
200
 
 
201
                if(nrows && ncols && nlinks) 
 
202
                  C_DGEMM('n', 't', nrows, ncols, nlinks, 0.5, W1[Gab][0], nlinks,
 
203
                          &(WMNIE.matrix[Gjk][jk][mc]), nlinks, 0.0, Z[0], ncols);
 
204
 
 
205
                for(m=0; m < ncols; m++) {
 
206
                  M = occ_off[Gm] + m;
 
207
                  mi = L2new.params->rowidx[M][I];
 
208
                  im = L2new.params->rowidx[I][M];
 
209
                  for(ab=0; ab < nrows; ab++) {
 
210
                    L2new.matrix[Gmi][mi][ab] += Z[ab][m];
 
211
                    L2new.matrix[Gmi][im][ab] -= Z[ab][m];
 
212
                  }
 
213
                }
 
214
 
 
215
                dpd_free_block(Z, nrows, ncols);
 
216
              }
 
217
 
 
218
              /* Z_IDAB <-- 1/2 L_IJKABC t_JKDC */
 
219
 
 
220
              jk = T2.params->rowidx[J][K];
 
221
              for(Gab=0; Gab < nirreps; Gab++) {
 
222
                Gid = Gab; /* totally symmetric */
 
223
                Gc = Gab ^ Gijk; /* totally symmetric */
 
224
                Gd = Gi ^ Gid;
 
225
 
 
226
                nrows = virtpi[Gd];
 
227
                ncols = ZIGDE.params->coltot[Gid];
 
228
                nlinks = virtpi[Gc];
 
229
 
 
230
                dc = T2.col_offset[Gjk][Gd];
 
231
                id = ZIGDE.row_offset[Gid][I];
 
232
                ZIGDE.matrix[Gid] = dpd_block_matrix(nrows, ncols);
 
233
                dpd_buf4_mat_irrep_rd_block(&ZIGDE, Gid, id, nrows);
 
234
 
 
235
                if(nrows && ncols && nlinks)
 
236
                  C_DGEMM('n', 't', nrows, ncols, nlinks, 0.5, &(T2.matrix[Gjk][jk][dc]), nlinks,
 
237
                          W1[Gab][0], nlinks, 1.0, ZIGDE.matrix[Gid][0], ncols);
 
238
 
 
239
                dpd_buf4_mat_irrep_wrt_block(&ZIGDE, Gid, id, nrows);
 
240
                dpd_free_block(ZIGDE.matrix[Gid], nrows, ncols);
 
241
              }
 
242
 
 
243
              /* Z_JDAB <-- 1/2 L_IJKABC t_IKDC */
 
244
              ik = T2.params->rowidx[I][K];
 
245
              for(Gab=0; Gab < nirreps; Gab++) {
 
246
                Gjd = Gab; /* totally symmetric */
 
247
                Gc = Gab ^ Gijk; /* totally symmetric */
 
248
                Gd = Gj ^ Gjd;
 
249
 
 
250
                nrows = virtpi[Gd];
 
251
                ncols = ZDMAE.params->coltot[Gjd];
 
252
                nlinks = virtpi[Gc];
 
253
 
 
254
                dc = T2.col_offset[Gik][Gd];
 
255
                jd = ZDMAE.row_offset[Gjd][J];
 
256
                ZDMAE.matrix[Gjd] = dpd_block_matrix(nrows, ncols);
 
257
                dpd_buf4_mat_irrep_rd_block(&ZDMAE, Gjd, jd, nrows);
 
258
 
 
259
                if(nrows && ncols && nlinks) 
 
260
                  C_DGEMM('n', 't', nrows, ncols, nlinks, 0.5, &(T2.matrix[Gik][ik][dc]), nlinks,
 
261
                          W1[Gab][0], nlinks, 1.0, ZDMAE.matrix[Gjd][0], ncols);
 
262
 
 
263
                dpd_buf4_mat_irrep_wrt_block(&ZDMAE, Gjd, jd, nrows);
 
264
                dpd_free_block(ZDMAE.matrix[Gjd], nrows, ncols);
 
265
              }
 
266
 
 
267
              /* Z_IJAM <-- -1/2 L_IJKABC t_MKBC */
 
268
              /* sort W(AB,C) to W(A,BC) */
 
269
              for(Gab=0; Gab < nirreps; Gab++) {
 
270
                Gc = Gab ^ Gijk;
 
271
                for(ab=0; ab < F.params->coltot[Gab]; ab++) {
 
272
                  A = F.params->colorb[Gab][ab][0];
 
273
                  B = F.params->colorb[Gab][ab][1];
 
274
                  Ga = F.params->rsym[A];
 
275
                  a = A - vir_off[Ga];
 
276
                  for(c=0; c < virtpi[Gc]; c++) {
 
277
                    C = vir_off[Gc] + c;
 
278
                    bc = F.params->colidx[B][C];
 
279
                    W2[Ga][a][bc] = W1[Gab][ab][c];
 
280
                  }
 
281
                }
 
282
              }
 
283
 
 
284
              ij = ZLMAO.params->rowidx[I][J];
 
285
 
 
286
              for(Gm=0; Gm < nirreps; Gm++) {
 
287
                Gbc = Gmk = Gm ^ Gk; /* totally symmetric */
 
288
                Ga = Gij ^ Gm; /* totally symmetric */
 
289
 
 
290
                nrows = virtpi[Ga];
 
291
                ncols = T2.params->coltot[Gmk];
 
292
 
 
293
                for(m=0; m < occpi[Gm]; m++) {
 
294
                  M = occ_off[Gm] + m;
 
295
                  mk = T2.params->rowidx[M][K];
 
296
                  am = ZLMAO.col_offset[Gij][Ga] + m;
 
297
 
 
298
                  if(nrows && ncols)
 
299
                    C_DGEMV('n', nrows, ncols, -0.5, W2[Ga][0], ncols, T2.matrix[Gmk][mk], 1, 
 
300
                            1.0, &(ZLMAO.matrix[Gij][ij][am]), occpi[Gm]);
 
301
                }
 
302
              }
 
303
 
 
304
              /* Z_IJMB <-- -1/2 L_IJKABC t_MKAC */
 
305
              /* sort W(AB,C) to W(B,AC) */
 
306
              for(Gab=0; Gab < nirreps; Gab++) {
 
307
                Gc = Gab ^ Gijk;
 
308
                for(ab=0; ab < F.params->coltot[Gab]; ab++) {
 
309
                  A = F.params->colorb[Gab][ab][0];
 
310
                  B = F.params->colorb[Gab][ab][1];
 
311
                  Gb = F.params->ssym[B];
 
312
                  b = B - vir_off[Gb];
 
313
                  for(c=0; c < virtpi[Gc]; c++) {
 
314
                    C = vir_off[Gc] + c;
 
315
                    ac = F.params->colidx[A][C];
 
316
                    W2[Gb][b][ac] = W1[Gab][ab][c];
 
317
                  }
 
318
                }
 
319
              }
 
320
 
 
321
              ij = ZIMLE.params->rowidx[I][J];
 
322
 
 
323
              for(Gm=0; Gm < nirreps; Gm++) {
 
324
                Gb = Gm ^ Gij; /* totally symmetric */
 
325
                Gmk = Gm ^ Gk;
 
326
 
 
327
                nrows = virtpi[Gb];
 
328
                ncols = T2.params->coltot[Gmk];
 
329
 
 
330
                for(m=0; m < occpi[Gm]; m++) {
 
331
                  M = occ_off[Gm] + m;
 
332
                  mk = T2.params->rowidx[M][K];
 
333
                  mb = ZIMLE.col_offset[Gij][Gm] + m * virtpi[Gb];
 
334
 
 
335
                  if(nrows && ncols)
 
336
                    C_DGEMV('n', nrows, ncols, -0.5, W2[Gb][0], ncols, T2.matrix[Gmk][mk], 1,
 
337
                            1.0, &(ZIMLE.matrix[Gij][ij][mb]), 1);
 
338
                }
 
339
 
 
340
              }
 
341
 
 
342
 
 
343
            } /* k */
 
344
          } /* j */
 
345
        } /* i */
 
346
 
 
347
        for(Gab=0; Gab < nirreps; Gab++) {
 
348
          Gc = Gab ^ Gijk; /* totally symmetric */
 
349
          dpd_free_block(W1[Gab], F.params->coltot[Gab], virtpi[Gc]);
 
350
        }
 
351
        for(Ga=0; Ga < nirreps; Ga++) {
 
352
          Gbc = Ga ^ Gijk; /* totally symmetric */
 
353
          dpd_free_block(W2[Ga], virtpi[Ga], F.params->coltot[Gbc]);
 
354
        }
 
355
 
 
356
      } /* Gk */
 
357
    } /* Gj */
 
358
  } /* Gi */
 
359
 
 
360
  free(W1);
 
361
  free(W2);
 
362
 
 
363
  dpd_buf4_close(&E);
 
364
  dpd_buf4_close(&F);
 
365
  dpd_buf4_close(&L);
 
366
 
 
367
  dpd_file2_close(&fIJ);
 
368
  dpd_file2_close(&fAB);
 
369
 
 
370
  dpd_file2_close(&FME);
 
371
  dpd_file2_close(&LIA);
 
372
  dpd_buf4_close(&Dints);
 
373
  dpd_buf4_close(&LIJAB);
 
374
 
 
375
  dpd_buf4_close(&WMAFE);
 
376
  for(h=0; h < nirreps; h++) dpd_buf4_mat_irrep_close(&WMNIE, h);
 
377
  dpd_buf4_close(&WMNIE);
 
378
 
 
379
  dpd_buf4_close(&ZIGDE);
 
380
  dpd_buf4_close(&ZDMAE);
 
381
 
 
382
  for(h=0; h < nirreps; h++) {
 
383
    dpd_buf4_mat_irrep_wrt(&ZLMAO, h);
 
384
    dpd_buf4_mat_irrep_close(&ZLMAO, h);
 
385
  }
 
386
  dpd_buf4_close(&ZLMAO);
 
387
 
 
388
  for(h=0; h < nirreps; h++) {
 
389
    dpd_buf4_mat_irrep_wrt(&ZIMLE, h);
 
390
    dpd_buf4_mat_irrep_close(&ZIMLE, h);
 
391
  }
 
392
  dpd_buf4_close(&ZIMLE);
 
393
 
 
394
  for(h=0; h < nirreps; h++) dpd_buf4_mat_irrep_close(&T2, h);
 
395
  dpd_buf4_close(&T2);
 
396
 
 
397
  for(h=0; h < nirreps; h++) {
 
398
    dpd_buf4_mat_irrep_wrt(&L2new, h);
 
399
    dpd_buf4_mat_irrep_close(&L2new, h);
 
400
  }
 
401
  dpd_buf4_init(&D2, CC_DENOM, 0, 0, 5, 0, 5, 0, "dIjAb");
 
402
  dpd_buf4_dirprd(&D2, &L2new);
 
403
  dpd_buf4_close(&D2);
 
404
  dpd_buf4_init(&L2, CC_LAMBDA, 0, 0, 5, 2, 7, 0, "New LIJAB");
 
405
  dpd_buf4_axpy(&L2new, &L2, 1);
 
406
  dpd_buf4_close(&L2);
 
407
  dpd_buf4_close(&L2new);
 
408
 
 
409
}
 
410
 
 
411
void cc3_l3l2_RHF_AAB(void)
 
412
{
 
413
  int h, nirreps;
 
414
  int *occ_off, *occpi;
 
415
  int *vir_off, *virtpi;
 
416
  int Gi, Gj, Gk, Gijk;
 
417
  int Ga, Gb, Gc, Gab;
 
418
  int i, j, k, I, J, K;
 
419
  int a, b, c, A, B, C;
 
420
  int ab;
 
421
  double ***W1, ***W2;
 
422
  dpdbuf4 L2AA, L2AB, L2BA, EAA, EAB, EBA, FAA, FAB, FBA;
 
423
  dpdfile2 fIJ, fAB, fij, fab;
 
424
  dpdbuf4 DAAints, DABints, LIJAB, LIjAb;
 
425
  dpdfile2 LIA, Lia, FME, Fme;
 
426
  dpdbuf4 L2AAnew, L2ABnew, L2, D2;
 
427
  dpdbuf4 WmAfE, WMnIe, WMAFE, WMaFe, WMNIE, WmNiE;
 
428
  dpdbuf4 ZIGDE, T2AB, T2AA, ZIgDe;
 
429
  dpdbuf4 ZDMAE, ZDmAe, ZdMAe;
 
430
  dpdbuf4 ZLMAO, ZLmAo;
 
431
  dpdbuf4 ZIMLE, ZImLe, ZImlE;
 
432
  int nrows, ncols, nlinks;
 
433
  int Gcb, cb;
 
434
  int Gij, ij, Gji, ji, Gjk, jk, kj, Gkj;
 
435
  int Gd, d, DD, ad, da, Gkd, kd;
 
436
  int Gm, m, M, Gmi, mi, im, mc;
 
437
  int Gid, id, dc, cd;
 
438
  int Gac, ac, Gca, ca, bd, db;
 
439
  int Gbc, bc, Gmk, mk, km, Gim, ma, am;
 
440
  int Gik, ik, Gki, ki, Gjd, jd;
 
441
  int Gmj, mj, cm, Gjm, jm;
 
442
  int mb;
 
443
  double **Z;
 
444
 
 
445
  nirreps = moinfo.nirreps;
 
446
  occpi = moinfo.occpi;
 
447
  occ_off = moinfo.occ_off;
 
448
  virtpi = moinfo.virtpi;
 
449
  vir_off = moinfo.vir_off;
 
450
 
 
451
  dpd_buf4_init(&L2AAnew, CC3_MISC, 0, 0, 5, 0, 5, 0, "CC3 LIJAB");
 
452
  for(h=0; h < nirreps; h++) dpd_buf4_mat_irrep_init(&L2AAnew, h);
 
453
 
 
454
  dpd_buf4_init(&L2ABnew, CC3_MISC, 0, 0, 5, 0, 5, 0, "CC3 LIjAb");
 
455
  for(h=0; h < nirreps; h++) dpd_buf4_mat_irrep_init(&L2ABnew, h);
 
456
 
 
457
  dpd_buf4_init(&T2AB, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tIjAb");
 
458
  dpd_buf4_init(&T2AA, CC_TAMPS, 0, 0, 5, 2, 7, 0, "tIJAB");
 
459
  for(h=0; h < nirreps; h++) {
 
460
    dpd_buf4_mat_irrep_init(&T2AB, h);
 
461
    dpd_buf4_mat_irrep_rd(&T2AB, h);
 
462
 
 
463
    dpd_buf4_mat_irrep_init(&T2AA, h);
 
464
    dpd_buf4_mat_irrep_rd(&T2AA, h);
 
465
  }
 
466
 
 
467
  dpd_buf4_init(&ZIGDE, CC3_MISC, 0, 10, 5, 10, 5, 0, "CC3 ZIGDE");
 
468
  dpd_buf4_scm(&ZIGDE, 0.0); /* this must be cleared in each iteration */
 
469
  dpd_buf4_init(&ZIgDe, CC3_MISC, 0, 10, 5, 10, 5, 0, "CC3 ZIgDe");
 
470
  dpd_buf4_scm(&ZIgDe, 0.0); /* this must be cleared in each iteration */
 
471
 
 
472
  dpd_buf4_init(&ZDMAE, CC3_MISC, 0, 10, 5, 10, 5, 0, "CC3 ZDMAE (MD,AE)");
 
473
  dpd_buf4_scm(&ZDMAE, 0.0); /* must be cleared in each iteration */
 
474
  dpd_buf4_init(&ZDmAe, CC3_MISC, 0, 10, 5, 10, 5, 0, "CC3 ZDmAe (mD,Ae)");
 
475
  dpd_buf4_scm(&ZDmAe, 0.0); /* must be cleared in each iteration */
 
476
  dpd_buf4_init(&ZdMAe, CC3_MISC, 0, 10, 5, 10, 5, 0, "CC3 ZdMAe (Md,Ae)");
 
477
  dpd_buf4_scm(&ZdMAe, 0.0); /* must be cleared in each iteration */
 
478
 
 
479
  dpd_buf4_init(&ZLMAO, CC3_MISC, 0, 0, 11, 0, 11, 0, "CC3 ZLMAO");
 
480
  dpd_buf4_init(&ZLmAo, CC3_MISC, 0, 0, 11, 0, 11, 0, "CC3 ZLmAo");
 
481
  for(h=0; h < nirreps; h++) {
 
482
    dpd_buf4_mat_irrep_init(&ZLMAO, h);
 
483
    dpd_buf4_mat_irrep_rd(&ZLMAO, h);
 
484
 
 
485
    dpd_buf4_mat_irrep_init(&ZLmAo, h);
 
486
  }
 
487
 
 
488
  dpd_buf4_init(&ZIMLE, CC3_MISC, 0, 0, 10, 0, 10, 0, "CC3 ZIMLE");
 
489
  dpd_buf4_init(&ZImLe, CC3_MISC, 0, 0, 10, 0, 10, 0, "CC3 ZImLe");
 
490
  dpd_buf4_init(&ZImlE, CC3_MISC, 0, 0, 10, 0, 10, 0, "CC3 ZImlE");
 
491
  for(h=0; h < nirreps; h++) {
 
492
    dpd_buf4_mat_irrep_init(&ZIMLE, h);
 
493
    dpd_buf4_mat_irrep_rd(&ZIMLE, h);
 
494
 
 
495
    dpd_buf4_mat_irrep_init(&ZImLe, h);
 
496
    dpd_buf4_mat_irrep_init(&ZImlE, h);
 
497
  }
 
498
 
 
499
  dpd_buf4_init(&WmAfE, CC3_HET1, 0, 10, 5, 10, 5, 0, "CC3 WAbEi (iE,bA)");
 
500
  dpd_buf4_init(&WMAFE, CC3_HET1, 0, 10, 5, 10, 7, 0, "CC3 WABEI (IE,B>A)");
 
501
  dpd_buf4_init(&WMaFe, CC3_HET1, 0, 10, 5, 10, 5, 0, "CC3 WaBeI (Ie,Ba)");
 
502
 
 
503
  dpd_buf4_init(&WMnIe, CC3_HET1, 0, 0, 10, 0, 10, 0, "CC3 WMbIj (Ij,Mb)");
 
504
  dpd_buf4_init(&WMNIE, CC3_HET1, 0, 0, 10, 2, 10, 0, "CC3 WMBIJ (I>J,MB)");
 
505
  dpd_buf4_init(&WmNiE, CC3_HET1, 0, 0, 10, 0, 10, 0, "CC3 WmBiJ (iJ,mB)");
 
506
  for(h=0; h < nirreps; h++) {
 
507
    dpd_buf4_mat_irrep_init(&WMnIe, h);
 
508
    dpd_buf4_mat_irrep_rd(&WMnIe, h);
 
509
    dpd_buf4_mat_irrep_init(&WMNIE, h);
 
510
    dpd_buf4_mat_irrep_rd(&WMNIE, h);
 
511
    dpd_buf4_mat_irrep_init(&WmNiE, h);
 
512
    dpd_buf4_mat_irrep_rd(&WmNiE, h);
 
513
  }
 
514
 
 
515
  dpd_file2_init(&fIJ, CC_OEI, 0, 0, 0, "fIJ");
 
516
  dpd_file2_init(&fAB, CC_OEI, 0, 1, 1, "fAB");
 
517
  dpd_file2_init(&fij, CC_OEI, 0, 0, 0, "fij");
 
518
  dpd_file2_init(&fab, CC_OEI, 0, 1, 1, "fab");
 
519
 
 
520
  dpd_buf4_init(&L2AA, CC_LAMBDA, 0, 0, 5, 2, 7, 0, "LIJAB");
 
521
  dpd_buf4_init(&L2AB, CC_LAMBDA, 0, 0, 5, 0, 5, 0, "LIjAb");
 
522
  dpd_buf4_init(&L2BA, CC_LAMBDA, 0, 0, 5, 0, 5, 0, "LiJaB");
 
523
  dpd_buf4_init(&FAA, CC3_HET1, 0, 10, 5, 10, 7, 0, "CC3 WAMEF (MA,F>E)");
 
524
  dpd_buf4_init(&FAB, CC3_HET1, 0, 10, 5, 10, 5, 0, "CC3 WaMeF (Ma,Fe)");
 
525
  dpd_buf4_init(&FBA, CC3_HET1, 0, 10, 5, 10, 5, 0, "CC3 WAmEf (mA,fE)");
 
526
  dpd_buf4_init(&EAA, CC3_HET1, 0, 0, 10, 2, 10, 0, "CC3 WMNIE (M>N,IE)");
 
527
  dpd_buf4_init(&EAB, CC3_HET1, 0, 0, 10, 0, 10, 0, "CC3 WMnIe (Mn,Ie)");
 
528
  dpd_buf4_init(&EBA, CC3_HET1, 0, 0, 10, 0, 10, 0, "CC3 WmNiE (mN,iE)");
 
529
 
 
530
  dpd_buf4_init(&DAAints, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij||ab>");
 
531
  dpd_buf4_init(&DABints, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
 
532
  dpd_buf4_init(&LIJAB, CC_LAMBDA, 0, 0, 5, 2, 7, 0, "LIJAB");
 
533
  dpd_buf4_init(&LIjAb, CC_LAMBDA, 0, 0, 5, 0, 5, 0, "LIjAb");
 
534
  dpd_file2_init(&LIA, CC_LAMBDA, 0, 0, 1, "LIA");
 
535
  dpd_file2_init(&Lia, CC_LAMBDA, 0, 0, 1, "Lia");
 
536
  dpd_file2_init(&FME, CC_OEI, 0, 0, 1, "FME");
 
537
  dpd_file2_init(&Fme, CC_OEI, 0, 0, 1, "Fme");
 
538
 
 
539
  W1 = (double ***) malloc(nirreps * sizeof(double **));
 
540
  W2 = (double ***) malloc(nirreps * sizeof(double **));
 
541
 
 
542
  for(Gi=0; Gi < nirreps; Gi++) {
 
543
    for(Gj=0; Gj < nirreps; Gj++) {
 
544
      Gij = Gji = Gi ^ Gj;
 
545
      for(Gk=0; Gk < nirreps; Gk++) {
 
546
        Gijk = Gi ^ Gj ^ Gk;
 
547
 
 
548
        Gjk = Gkj = Gj ^ Gk;
 
549
        Gik = Gki = Gi ^ Gk;
 
550
 
 
551
        for(Gab=0; Gab < nirreps; Gab++) {
 
552
          Gc = Gab ^ Gijk; /* totally symmetric */
 
553
          W1[Gab] = dpd_block_matrix(FAA.params->coltot[Gab], virtpi[Gc]);
 
554
        }
 
555
        for(Ga=0; Ga < nirreps; Ga++) {
 
556
          Gcb = Ga ^ Gijk; /* assumes totally symmetric */
 
557
          W2[Ga] = dpd_block_matrix(virtpi[Ga], WmAfE.params->coltot[Gcb]);  /* alpha-beta-alpha */
 
558
        }
 
559
 
 
560
        for(i=0; i < occpi[Gi]; i++) {
 
561
          I = occ_off[Gi] + i;
 
562
          for(j=0; j < occpi[Gj]; j++) {
 
563
            J = occ_off[Gj] + j;
 
564
            for(k=0; k < occpi[Gk]; k++) {
 
565
              K = occ_off[Gk] + k;
 
566
 
 
567
              L3_AAB(W1, nirreps, I, Gi, J, Gj, K, Gk, &L2AA, &L2AB, &L2BA, 
 
568
                     &FAA, &FAB, &FBA, &EAA, &EAB, &EBA, &fIJ, &fij, &fAB, &fab,
 
569
                     &DAAints, &DABints, &LIJAB, &LIjAb, &LIA, &Lia, &FME, &Fme,
 
570
                     occpi, occ_off, occpi, occ_off, virtpi, vir_off, virtpi, vir_off);
 
571
 
 
572
              /* t_JIDA <-- t_IJkABc W_kDcB */
 
573
              /* sort W1(AB,c) to W2(A,cB) */
 
574
              for(Gab=0; Gab < nirreps; Gab++) {
 
575
                Gc = Gab ^ Gijk;
 
576
                for(ab=0; ab < FAA.params->coltot[Gab]; ab++) {
 
577
                  A = FAA.params->colorb[Gab][ab][0];
 
578
                  B = FAA.params->colorb[Gab][ab][1];
 
579
                  Ga = FAA.params->rsym[A];
 
580
                  a = A - vir_off[Ga];
 
581
                  for(c=0; c < virtpi[Gc]; c++) {
 
582
                    C = vir_off[Gc] + c;
 
583
                    cb = WmAfE.params->colidx[C][B];
 
584
                    W2[Ga][a][cb] = W1[Gab][ab][c];
 
585
                  }
 
586
                }
 
587
              }
 
588
 
 
589
              ji = L2AAnew.params->rowidx[J][I];
 
590
 
 
591
              for(Gd=0; Gd < nirreps; Gd++) {
 
592
                Gcb = Gkd = Gk ^ Gd; /* assumes totally symmetric */
 
593
                Ga = Gd ^ Gij;       /* assumes totally symmetric */
 
594
 
 
595
                kd = WmAfE.row_offset[Gkd][K];
 
596
                WmAfE.matrix[Gkd] = dpd_block_matrix(virtpi[Gd], WmAfE.params->coltot[Gkd]);
 
597
                dpd_buf4_mat_irrep_rd_block(&WmAfE, Gkd, kd, virtpi[Gd]);
 
598
                Z = block_matrix(virtpi[Ga], virtpi[Gd]);
 
599
 
 
600
                nrows = virtpi[Ga];
 
601
                ncols = virtpi[Gd];
 
602
                nlinks = WmAfE.params->coltot[Gkd];
 
603
 
 
604
                if(nrows && ncols && nlinks)
 
605
                  C_DGEMM('n', 't', nrows, ncols, nlinks, 1.0, W2[Ga][0], nlinks, 
 
606
                          WmAfE.matrix[Gkd][0], nlinks, 0.0, Z[0], ncols);
 
607
 
 
608
                for(a=0; a < virtpi[Ga]; a++) {
 
609
                  A = vir_off[Ga] + a;
 
610
                  for(d=0; d < virtpi[Gd]; d++) {
 
611
                    DD = vir_off[Gd] + d;
 
612
                    ad = L2AAnew.params->colidx[A][DD];
 
613
                    da = L2AAnew.params->colidx[DD][A];
 
614
                    L2AAnew.matrix[Gij][ji][ad] += -Z[a][d];
 
615
                    L2AAnew.matrix[Gij][ji][da] += Z[a][d];
 
616
                  }
 
617
                }
 
618
 
 
619
                dpd_free_block(WmAfE.matrix[Gkd], virtpi[Gd], WmAfE.params->coltot[Gkd]);
 
620
                free_block(Z);
 
621
              }
 
622
 
 
623
              /* t_MIAB <--- +t_IJkABc W_JkMc */
 
624
              /* t_IMAB <--- -t_IJkABc W_JkMc */
 
625
 
 
626
              jk = WMnIe.params->rowidx[J][K];
 
627
 
 
628
              for(Gm=0; Gm < nirreps; Gm++) {
 
629
                Gab = Gmi = Gm ^ Gi;  /* assumes totally symmetric */
 
630
                Gc = Gab ^ Gijk;      /* assumes totally symmetric */
 
631
 
 
632
                mc = WMnIe.col_offset[Gjk][Gm];
 
633
 
 
634
                nrows = FAA.params->coltot[Gab];
 
635
                ncols = occpi[Gm];
 
636
                nlinks = virtpi[Gc];
 
637
 
 
638
                Z = dpd_block_matrix(nrows, ncols);
 
639
 
 
640
                if(nrows && ncols && nlinks)
 
641
                  C_DGEMM('n', 't', nrows, ncols, nlinks, 1.0, W1[Gab][0], nlinks,
 
642
                          &(WMnIe.matrix[Gjk][jk][mc]), nlinks, 0.0, Z[0], ncols);
 
643
 
 
644
                for(m=0; m < ncols; m++) {
 
645
                  M = occ_off[Gm] + m;
 
646
                  mi = L2AAnew.params->rowidx[M][I];
 
647
                  im = L2AAnew.params->rowidx[I][M];
 
648
                  for(ab=0; ab < nrows; ab++) {
 
649
                    L2AAnew.matrix[Gmi][mi][ab] += Z[ab][m];
 
650
                    L2AAnew.matrix[Gmi][im][ab] -= Z[ab][m];
 
651
                  }
 
652
                }
 
653
 
 
654
                dpd_free_block(Z, nrows, ncols);
 
655
              }
 
656
 
 
657
              /* t_JkDc <-- 1/2 t_IJkABc W_IDAB */
 
658
              /* t_KjCd <-- 1/2 t_IJkABc W_IDAB */
 
659
 
 
660
              jk = L2ABnew.params->rowidx[J][K];
 
661
              kj = L2ABnew.params->rowidx[K][J];
 
662
 
 
663
              for(Gd=0; Gd < nirreps; Gd++) {
 
664
                Gab = Gid = Gi ^ Gd; /* assumes totally symmetric */
 
665
                Gc = Gab ^ Gijk;     /* assumes totally symmetric */
 
666
 
 
667
                id = WMAFE.row_offset[Gid][I];
 
668
                WMAFE.matrix[Gid] = dpd_block_matrix(virtpi[Gd], WMAFE.params->coltot[Gid]);
 
669
                dpd_buf4_mat_irrep_rd_block(&WMAFE, Gid, id, virtpi[Gd]);
 
670
                Z = block_matrix(virtpi[Gc], virtpi[Gd]);
 
671
 
 
672
                nrows = virtpi[Gc];
 
673
                ncols = virtpi[Gd];
 
674
                nlinks = WMAFE.params->coltot[Gid];
 
675
 
 
676
                if(nrows && ncols && nlinks)
 
677
                  C_DGEMM('t', 't', nrows, ncols, nlinks, 0.5, W1[Gab][0], nrows,
 
678
                          WMAFE.matrix[Gid][0], nlinks, 0.0, Z[0], ncols);
 
679
 
 
680
                for(c=0; c < virtpi[Gc]; c++) {
 
681
                  C = vir_off[Gc] + c;
 
682
                  for(d=0; d < virtpi[Gd]; d++) {
 
683
                    DD = vir_off[Gd] + d;
 
684
                    dc = L2ABnew.params->colidx[DD][C];
 
685
                    cd = L2ABnew.params->colidx[C][DD];
 
686
                    L2ABnew.matrix[Gjk][jk][dc] += Z[c][d];
 
687
                    L2ABnew.matrix[Gjk][kj][cd] += Z[c][d];
 
688
                  }
 
689
                }
 
690
 
 
691
                free_block(Z);
 
692
                dpd_free_block(WMAFE.matrix[Gid], virtpi[Gd], WMAFE.params->coltot[Gid]);
 
693
              }
 
694
 
 
695
              /* t_JkBd <-- t_IJkABc W_IdAc */
 
696
              /* t_KjBd <-- t_IJkABc W_IdAc */
 
697
              /* sort W1(AB,c) to W2(B,Ac) */
 
698
              for(Gab=0; Gab < nirreps; Gab++) {
 
699
                Gc = Gab ^ Gijk;
 
700
                for(ab=0; ab < FAA.params->coltot[Gab]; ab++) {
 
701
                  A = FAA.params->colorb[Gab][ab][0];
 
702
                  B = FAA.params->colorb[Gab][ab][1];
 
703
                  Gb = FAA.params->ssym[B];
 
704
                  b = B - vir_off[Gb];
 
705
                  for(c=0; c < virtpi[Gc]; c++) {
 
706
                    C = vir_off[Gc] + c;
 
707
                    ac = WMaFe.params->colidx[A][C];
 
708
                    W2[Gb][b][ac] = W1[Gab][ab][c];
 
709
                  }
 
710
                }
 
711
              }
 
712
 
 
713
              jk = L2ABnew.params->rowidx[J][K];
 
714
              kj = L2ABnew.params->rowidx[K][J];
 
715
 
 
716
              for(Gd=0; Gd < nirreps; Gd++) {
 
717
                Gac = Gid = Gi ^ Gd; /* assumes totally symmetric */
 
718
                Gb = Gac ^ Gijk;     /* assumes totally symmetric */
 
719
 
 
720
                id = WMaFe.row_offset[Gid][I]; 
 
721
                WMaFe.matrix[Gid] = dpd_block_matrix(virtpi[Gd], WMaFe.params->coltot[Gid]);
 
722
                dpd_buf4_mat_irrep_rd_block(&WMaFe, Gid, id, virtpi[Gd]);
 
723
                Z = block_matrix(virtpi[Gb], virtpi[Gd]);
 
724
 
 
725
                nrows = virtpi[Gb];
 
726
                ncols = virtpi[Gd];
 
727
                nlinks = WMaFe.params->coltot[Gid];
 
728
 
 
729
                if(nrows && ncols && nlinks)
 
730
                  C_DGEMM('n', 't', nrows, ncols, nlinks, 1.0, W2[Gb][0], nlinks,
 
731
                          WMaFe.matrix[Gid][0], nlinks, 0.0, Z[0], ncols);
 
732
 
 
733
                for(b=0; b < virtpi[Gb]; b++) {
 
734
                  B = vir_off[Gb] + b;
 
735
                  for(d=0; d < virtpi[Gd]; d++) {
 
736
                    DD = vir_off[Gd] + d;
 
737
                    bd = L2ABnew.params->colidx[B][DD];
 
738
                    db = L2ABnew.params->colidx[DD][B];
 
739
                    L2ABnew.matrix[Gjk][jk][bd] += Z[b][d];
 
740
                    L2ABnew.matrix[Gjk][kj][db] += Z[b][d];
 
741
                  }
 
742
                }
 
743
 
 
744
                dpd_free_block(WMaFe.matrix[Gid], virtpi[Gd], WMaFe.params->coltot[Gid]);
 
745
                free_block(Z);
 
746
              }
 
747
 
 
748
              /* t_MkBc <-- 1/2 t_IJkABc W_IJMA */
 
749
              /* sort W(AB,c) to W(A,Bc) */
 
750
              for(Gab=0; Gab < nirreps; Gab++) {
 
751
                Gc = Gab ^ Gijk;  /* assumes totally symmetric */
 
752
                for(ab=0; ab < FAA.params->coltot[Gab]; ab++ ){
 
753
                  A = FAA.params->colorb[Gab][ab][0];
 
754
                  B = FAA.params->colorb[Gab][ab][1];
 
755
                  Ga = FAA.params->rsym[A];
 
756
                  a = A - vir_off[Ga];
 
757
                  for(c=0; c < virtpi[Gc]; c++) {
 
758
                    C = vir_off[Gc] + c;
 
759
                    bc = L2ABnew.params->colidx[B][C];
 
760
                    W2[Ga][a][bc] = W1[Gab][ab][c];
 
761
                  }
 
762
                }
 
763
              }
 
764
 
 
765
              ij = WMNIE.params->rowidx[I][J];
 
766
 
 
767
              for(Gm=0; Gm < nirreps; Gm++) {
 
768
                Gbc = Gmk = Gm ^ Gk;  /* assumes totally symmetric */
 
769
                Ga = Gbc ^ Gijk;      /* assumes totally symmetric */
 
770
 
 
771
                ma = WMNIE.col_offset[Gij][Gm];
 
772
 
 
773
                nrows = L2ABnew.params->coltot[Gmk];
 
774
                ncols = occpi[Gm];
 
775
                nlinks = virtpi[Ga];
 
776
 
 
777
                Z = dpd_block_matrix(nrows, ncols);
 
778
 
 
779
                if(nrows && ncols && nlinks)
 
780
                  C_DGEMM('t', 't', nrows, ncols, nlinks, 0.5, W2[Ga][0], nrows,
 
781
                          &(WMNIE.matrix[Gij][ij][ma]), nlinks, 0.0, Z[0], ncols);
 
782
 
 
783
                for(m=0; m < occpi[Gm]; m++) {
 
784
                  M = occ_off[Gm] + m;
 
785
                  mk = L2ABnew.params->rowidx[M][K];
 
786
                  km = L2ABnew.params->rowidx[K][M];
 
787
                  for(Gb=0; Gb < nirreps; Gb++) {
 
788
                    Gc = Gbc ^ Gb;
 
789
                    for(b=0; b < virtpi[Gb]; b++) {
 
790
                      B = vir_off[Gb] + b;
 
791
                      for(c=0; c < virtpi[Gc]; c++) {
 
792
                        C = vir_off[Gc] + c;
 
793
                        bc = L2ABnew.params->colidx[B][C];
 
794
                        cb = L2ABnew.params->colidx[C][B];
 
795
                        L2ABnew.matrix[Gmk][mk][bc] += Z[bc][m];
 
796
                        L2ABnew.matrix[Gmk][km][cb] += Z[bc][m];
 
797
                      }
 
798
                    }
 
799
                  }
 
800
                }
 
801
 
 
802
                dpd_free_block(Z, nrows, ncols);
 
803
              }
 
804
 
 
805
              /* t_ImBc <-- t_IJkABc W_kJmA */
 
806
              /* sort W(AB,c) to W(A,Bc) */
 
807
              for(Gab=0; Gab < nirreps; Gab++) {
 
808
                Gc = Gab ^ Gijk;  /* assumes totally symmetric */
 
809
                for(ab=0; ab < FAA.params->coltot[Gab]; ab++ ){
 
810
                  A = FAA.params->colorb[Gab][ab][0];
 
811
                  B = FAA.params->colorb[Gab][ab][1];
 
812
                  Ga = FAA.params->rsym[A];
 
813
                  a = A - vir_off[Ga];
 
814
                  for(c=0; c < virtpi[Gc]; c++) {
 
815
                    C = vir_off[Gc] + c;
 
816
                    bc = L2ABnew.params->colidx[B][C];
 
817
                    W2[Ga][a][bc] = W1[Gab][ab][c];
 
818
                  }
 
819
                }
 
820
              }
 
821
 
 
822
              kj = WmNiE.params->rowidx[K][J];
 
823
 
 
824
              for(Gm=0; Gm < nirreps; Gm++) {
 
825
                Gbc = Gim = Gi ^ Gm;  /* assumes totally symmetric */
 
826
                Ga = Gbc ^ Gijk;      /* assumes totally symmetric */
 
827
 
 
828
                ma = WmNiE.col_offset[Gjk][Gm];
 
829
 
 
830
                nrows = L2ABnew.params->coltot[Gim];
 
831
                ncols = occpi[Gm];
 
832
                nlinks = virtpi[Ga];
 
833
 
 
834
                Z = dpd_block_matrix(nrows, ncols);
 
835
 
 
836
                if(nrows && ncols && nlinks)
 
837
                  C_DGEMM('t', 't', nrows, ncols, nlinks, 1.0, W2[Ga][0], nrows,
 
838
                          &(WmNiE.matrix[Gjk][kj][ma]), nlinks, 0.0, Z[0], ncols);
 
839
 
 
840
                for(m=0; m < occpi[Gm]; m++) {
 
841
                  M = occ_off[Gm] + m;
 
842
                  im = L2ABnew.params->rowidx[I][M];
 
843
                  mi = L2ABnew.params->rowidx[M][I];
 
844
                  for(Gb=0; Gb < nirreps; Gb++) {
 
845
                    Gc = Gbc ^ Gb;
 
846
                    for(b=0; b < virtpi[Gb]; b++) {
 
847
                      B = vir_off[Gb] + b;
 
848
                      for(c=0; c < virtpi[Gc]; c++) {
 
849
                        C = vir_off[Gc] + c;
 
850
                        bc = L2ABnew.params->colidx[B][C];
 
851
                        cb = L2ABnew.params->colidx[C][B];
 
852
                        L2ABnew.matrix[Gim][im][bc] += Z[bc][m];
 
853
                        L2ABnew.matrix[Gim][mi][cb] += Z[bc][m];
 
854
                      }
 
855
                    }
 
856
                  }
 
857
                }
 
858
 
 
859
                dpd_free_block(Z, nrows, ncols);
 
860
              }
 
861
 
 
862
              /* Z_IDAB <-- L_IJkABc t_JkDc */
 
863
 
 
864
              jk = T2AB.params->rowidx[J][K];
 
865
              for(Gab=0; Gab < nirreps; Gab++) {
 
866
                Gid = Gab; /* totally symmetric */
 
867
                Gc = Gab ^ Gijk; /* totally symmetric */
 
868
                Gd = Gi ^ Gid;
 
869
 
 
870
                nrows = virtpi[Gd];
 
871
                ncols = ZIGDE.params->coltot[Gid];
 
872
                nlinks = virtpi[Gc];
 
873
 
 
874
                dc = T2AB.col_offset[Gjk][Gd];
 
875
                id = ZIGDE.row_offset[Gid][I];
 
876
                ZIGDE.matrix[Gid] = dpd_block_matrix(nrows, ncols);
 
877
 
 
878
                if(nrows && ncols && nlinks) {
 
879
                  dpd_buf4_mat_irrep_rd_block(&ZIGDE, Gid, id, nrows);
 
880
 
 
881
                  C_DGEMM('n', 't', nrows, ncols, nlinks, 1.0, &(T2AB.matrix[Gjk][jk][dc]), nlinks,
 
882
                          W1[Gab][0], nlinks, 1.0, ZIGDE.matrix[Gid][0], ncols);
 
883
 
 
884
                  dpd_buf4_mat_irrep_wrt_block(&ZIGDE, Gid, id, nrows);
 
885
                }
 
886
 
 
887
                dpd_free_block(ZIGDE.matrix[Gid], nrows, ncols);
 
888
              }
 
889
 
 
890
              /* ZkDCa <-- 1/2 L_ijKabC t_ijdb */
 
891
 
 
892
              ij = T2AA.params->rowidx[I][J];
 
893
 
 
894
              /* sort W(ab,C) to W(b,Ca) */
 
895
              for(Gab=0; Gab < nirreps; Gab++) {
 
896
                Gc = Gab ^ Gijk;  /* assumes totally symmetric */
 
897
                for(ab=0; ab < FAA.params->coltot[Gab]; ab++ ){
 
898
                  A = FAA.params->colorb[Gab][ab][0];
 
899
                  B = FAA.params->colorb[Gab][ab][1];
 
900
                  Gb = FAA.params->ssym[B];
 
901
                  b = B - vir_off[Gb];
 
902
                  for(c=0; c < virtpi[Gc]; c++) {
 
903
                    C = vir_off[Gc] + c;
 
904
                    ca = ZIgDe.params->colidx[C][A];
 
905
                    W2[Gb][b][ca] = W1[Gab][ab][c];
 
906
                  }
 
907
                }
 
908
              }
 
909
 
 
910
              for(Gb=0; Gb < nirreps; Gb++) {
 
911
                Gd = Gb ^ Gij; /* totally symmetric */
 
912
                Gca = Gkd = Gk ^ Gd; /* totally symmetric */
 
913
 
 
914
                nrows = virtpi[Gd];
 
915
                ncols = ZIgDe.params->coltot[Gkd];
 
916
                nlinks = virtpi[Gb];
 
917
 
 
918
                db = T2AA.col_offset[Gij][Gd];
 
919
                kd = ZIgDe.row_offset[Gkd][K];
 
920
                ZIgDe.matrix[Gkd] = dpd_block_matrix(nrows, ncols);
 
921
                dpd_buf4_mat_irrep_rd_block(&ZIgDe, Gkd, kd, nrows);
 
922
 
 
923
                if(nrows && ncols && nlinks)
 
924
                  C_DGEMM('n', 'n', nrows, ncols, nlinks, 0.5, &(T2AA.matrix[Gij][ij][db]), nlinks,
 
925
                          W2[Gb][0], ncols, 1.0, ZIgDe.matrix[Gkd][0], ncols);
 
926
 
 
927
                dpd_buf4_mat_irrep_wrt_block(&ZIgDe, Gkd, kd, nrows);
 
928
                dpd_free_block(ZIgDe.matrix[Gkd], nrows, ncols);
 
929
              }
 
930
 
 
931
              /* Z_IdAc <-- L_IJkABc t_JkBd */
 
932
 
 
933
              jk = T2AB.params->rowidx[J][K];
 
934
 
 
935
              /* sort W(AB,c) to W(B,Ac) */
 
936
              for(Gab=0; Gab < nirreps; Gab++) {
 
937
                Gc = Gab ^ Gijk;  /* assumes totally symmetric */
 
938
                for(ab=0; ab < FAA.params->coltot[Gab]; ab++ ){
 
939
                  A = FAA.params->colorb[Gab][ab][0];
 
940
                  B = FAA.params->colorb[Gab][ab][1];
 
941
                  Gb = FAA.params->ssym[B];
 
942
                  b = B - vir_off[Gb];
 
943
                  for(c=0; c < virtpi[Gc]; c++) {
 
944
                    C = vir_off[Gc] + c;
 
945
                    ac = ZIgDe.params->colidx[A][C];
 
946
                    W2[Gb][b][ac] = W1[Gab][ab][c];
 
947
                  }
 
948
                }
 
949
              }
 
950
 
 
951
              for(Gb=0; Gb < nirreps; Gb++) {
 
952
                Gd = Gb ^ Gjk; /* totally symmetric */
 
953
                Gac = Gid = Gi ^ Gd; /* totally symmetric */
 
954
 
 
955
                nrows = virtpi[Gd]; 
 
956
                ncols = ZIgDe.params->coltot[Gid];
 
957
                nlinks = virtpi[Gb];
 
958
 
 
959
                bd = T2AB.col_offset[Gjk][Gb];
 
960
                id = ZIgDe.row_offset[Gid][I];
 
961
                ZIgDe.matrix[Gid] = dpd_block_matrix(nrows, ncols);
 
962
                dpd_buf4_mat_irrep_rd_block(&ZIgDe, Gid, id, nrows);
 
963
 
 
964
                if(nrows && ncols && nlinks)
 
965
                  C_DGEMM('t', 'n', nrows, ncols, nlinks, 1.0, &(T2AB.matrix[Gjk][jk][bd]), nrows,
 
966
                          W2[Gb][0], ncols, 1.0, ZIgDe.matrix[Gid][0], ncols);
 
967
 
 
968
                dpd_buf4_mat_irrep_wrt_block(&ZIgDe, Gid, id, nrows);
 
969
                dpd_free_block(ZIgDe.matrix[Gid], nrows, ncols);
 
970
              }
 
971
 
 
972
              /* Z_JDAB <-- 1/2 L_IJkABc t_IkDc */
 
973
              ik = T2AB.params->rowidx[I][K];
 
974
              for(Gab=0; Gab < nirreps; Gab++) {
 
975
                Gjd = Gab; /* totally symmetric */
 
976
                Gc = Gab ^ Gijk; /* totally symmetric */
 
977
                Gd = Gj ^ Gjd;
 
978
 
 
979
                nrows = virtpi[Gd];
 
980
                ncols = ZDMAE.params->coltot[Gjd];
 
981
                nlinks = virtpi[Gc];
 
982
 
 
983
                dc = T2AB.col_offset[Gik][Gd];
 
984
                jd = ZDMAE.row_offset[Gjd][J];
 
985
                ZDMAE.matrix[Gjd] = dpd_block_matrix(nrows, ncols);
 
986
                dpd_buf4_mat_irrep_rd_block(&ZDMAE, Gjd, jd, nrows);
 
987
 
 
988
                if(nrows && ncols && nlinks) 
 
989
                  C_DGEMM('n', 't', nrows, ncols, nlinks, 1.0, &(T2AB.matrix[Gik][ik][dc]), nlinks,
 
990
                          W1[Gab][0], nlinks, 1.0, ZDMAE.matrix[Gjd][0], ncols);
 
991
 
 
992
                dpd_buf4_mat_irrep_wrt_block(&ZDMAE, Gjd, jd, nrows);
 
993
                dpd_free_block(ZDMAE.matrix[Gjd], nrows, ncols);
 
994
              }
 
995
 
 
996
              /* Z_kDAc <-- 1/2 L_IJkABc t_IJDB */
 
997
              ij = T2AA.params->rowidx[I][J];
 
998
              /* sort W(AB,c) to W(B,Ac) */
 
999
              for(Gab=0; Gab < nirreps; Gab++) {
 
1000
                Gc = Gab ^ Gijk;  /* assumes totally symmetric */
 
1001
                for(ab=0; ab < FAA.params->coltot[Gab]; ab++ ){
 
1002
                  A = FAA.params->colorb[Gab][ab][0];
 
1003
                  B = FAA.params->colorb[Gab][ab][1];
 
1004
                  Gb = FAA.params->ssym[B];
 
1005
                  b = B - vir_off[Gb];
 
1006
                  for(c=0; c < virtpi[Gc]; c++) {
 
1007
                    C = vir_off[Gc] + c;
 
1008
                    ac = ZDmAe.params->colidx[A][C];
 
1009
                    W2[Gb][b][ac] = W1[Gab][ab][c];
 
1010
                  }
 
1011
                }
 
1012
              }
 
1013
 
 
1014
              for(Gb=0; Gb < nirreps; Gb++) {
 
1015
                Gd = Gb ^ Gij; /* totally symmetric */
 
1016
                Gac = Gkd = Gk ^ Gd; /* totally symmetric */
 
1017
 
 
1018
                nrows = virtpi[Gd]; 
 
1019
                ncols = ZDmAe.params->coltot[Gkd];
 
1020
                nlinks = virtpi[Gb];
 
1021
 
 
1022
                db = T2AA.col_offset[Gij][Gd];
 
1023
                kd = ZDmAe.row_offset[Gkd][K];
 
1024
                ZDmAe.matrix[Gkd] = dpd_block_matrix(nrows, ncols);
 
1025
                dpd_buf4_mat_irrep_rd_block(&ZDmAe, Gkd, kd, nrows);
 
1026
 
 
1027
                if(nrows && ncols && nlinks)
 
1028
                  C_DGEMM('n', 'n', nrows, ncols, nlinks, 0.5, &(T2AA.matrix[Gij][ij][db]), nlinks,
 
1029
                          W2[Gb][0], ncols, 1.0, ZDmAe.matrix[Gkd][0], ncols);
 
1030
 
 
1031
                dpd_buf4_mat_irrep_wrt_block(&ZDmAe, Gkd, kd, nrows);
 
1032
                dpd_free_block(ZDmAe.matrix[Gkd], nrows, ncols);
 
1033
              }
 
1034
 
 
1035
              /* Z_iDCa <-- L_ijKabC t_KjDb */
 
1036
              kj = T2AB.params->rowidx[K][J];
 
1037
              /* sort W(AB,c) to W(B,Ca) */
 
1038
              for(Gab=0; Gab < nirreps; Gab++) {
 
1039
                Gc = Gab ^ Gijk;  /* assumes totally symmetric */
 
1040
                for(ab=0; ab < FAA.params->coltot[Gab]; ab++ ){
 
1041
                  A = FAA.params->colorb[Gab][ab][0];
 
1042
                  B = FAA.params->colorb[Gab][ab][1];
 
1043
                  Gb = FAA.params->ssym[B];
 
1044
                  b = B - vir_off[Gb];
 
1045
                  for(c=0; c < virtpi[Gc]; c++) {
 
1046
                    C = vir_off[Gc] + c;
 
1047
                    ca = ZDmAe.params->colidx[C][A];
 
1048
                    W2[Gb][b][ca] = W1[Gab][ab][c];
 
1049
                  }
 
1050
                }
 
1051
              }
 
1052
 
 
1053
              for(Gb=0; Gb < nirreps; Gb++) {
 
1054
                Gd = Gb ^ Gkj; /* totally symmetric */
 
1055
                Gca = Gid = Gi ^ Gd; /* totally symmetric */
 
1056
 
 
1057
                nrows = virtpi[Gd]; 
 
1058
                ncols = ZDmAe.params->coltot[Gid];
 
1059
                nlinks = virtpi[Gb];
 
1060
 
 
1061
                db = T2AB.col_offset[Gkj][Gd];
 
1062
                id = ZDmAe.row_offset[Gid][I];
 
1063
                ZDmAe.matrix[Gid] = dpd_block_matrix(nrows, ncols);
 
1064
                dpd_buf4_mat_irrep_rd_block(&ZDmAe, Gid, id, nrows);
 
1065
 
 
1066
                if(nrows && ncols && nlinks)
 
1067
                  C_DGEMM('n', 'n', nrows, ncols, nlinks, 1.0, &(T2AB.matrix[Gkj][kj][db]), nlinks,
 
1068
                          W2[Gb][0], ncols, 1.0, ZDmAe.matrix[Gid][0], ncols);
 
1069
 
 
1070
                dpd_buf4_mat_irrep_wrt_block(&ZDmAe, Gid, id, nrows);
 
1071
                dpd_free_block(ZDmAe.matrix[Gid], nrows, ncols);
 
1072
              }
 
1073
 
 
1074
              /* Z_KdCa <-- -1/2 L_ijKabC t_ijdb */
 
1075
              ij = T2AA.params->rowidx[I][J];
 
1076
              /* sort W(AB,c) to W(B,Ca) */
 
1077
              for(Gab=0; Gab < nirreps; Gab++) {
 
1078
                Gc = Gab ^ Gijk;  /* assumes totally symmetric */
 
1079
                for(ab=0; ab < FAA.params->coltot[Gab]; ab++ ){
 
1080
                  A = FAA.params->colorb[Gab][ab][0];
 
1081
                  B = FAA.params->colorb[Gab][ab][1];
 
1082
                  Gb = FAA.params->ssym[B];
 
1083
                  b = B - vir_off[Gb];
 
1084
                  for(c=0; c < virtpi[Gc]; c++) {
 
1085
                    C = vir_off[Gc] + c;
 
1086
                    ca = ZdMAe.params->colidx[C][A];
 
1087
                    W2[Gb][b][ca] = W1[Gab][ab][c];
 
1088
                  }
 
1089
                }
 
1090
              }
 
1091
 
 
1092
              for(Gb=0; Gb < nirreps; Gb++) {
 
1093
                Gd = Gb ^ Gij; /* totally symmetric */
 
1094
                Gca = Gkd = Gk ^ Gd; /* totally symmetric */
 
1095
 
 
1096
                nrows = virtpi[Gd]; 
 
1097
                ncols = ZdMAe.params->coltot[Gkd];
 
1098
                nlinks = virtpi[Gb];
 
1099
 
 
1100
                db = T2AA.col_offset[Gij][Gd];
 
1101
                kd = ZdMAe.row_offset[Gkd][K];
 
1102
                ZdMAe.matrix[Gkd] = dpd_block_matrix(nrows, ncols);
 
1103
                dpd_buf4_mat_irrep_rd_block(&ZdMAe, Gkd, kd, nrows);
 
1104
 
 
1105
                if(nrows && ncols && nlinks)
 
1106
                  C_DGEMM('n', 'n', nrows, ncols, nlinks, -0.5, &(T2AA.matrix[Gij][ij][db]), nlinks,
 
1107
                          W2[Gb][0], ncols, 1.0, ZdMAe.matrix[Gkd][0], ncols);
 
1108
 
 
1109
                dpd_buf4_mat_irrep_wrt_block(&ZdMAe, Gkd, kd, nrows);
 
1110
                dpd_free_block(ZdMAe.matrix[Gkd], nrows, ncols);
 
1111
              }
 
1112
 
 
1113
              /* Z_JdAc <-- L_IJkABc t_IkBd */
 
1114
              ik = T2AB.params->rowidx[I][K];
 
1115
 
 
1116
              /* sort W(AB,c) to W(B,Ca) */
 
1117
              for(Gab=0; Gab < nirreps; Gab++) {
 
1118
                Gc = Gab ^ Gijk;  /* assumes totally symmetric */
 
1119
                for(ab=0; ab < FAA.params->coltot[Gab]; ab++ ){
 
1120
                  A = FAA.params->colorb[Gab][ab][0];
 
1121
                  B = FAA.params->colorb[Gab][ab][1];
 
1122
                  Gb = FAA.params->ssym[B];
 
1123
                  b = B - vir_off[Gb];
 
1124
                  for(c=0; c < virtpi[Gc]; c++) {
 
1125
                    C = vir_off[Gc] + c;
 
1126
                    ac = ZdMAe.params->colidx[A][C];
 
1127
                    W2[Gb][b][ac] = W1[Gab][ab][c];
 
1128
                  }
 
1129
                }
 
1130
              }
 
1131
 
 
1132
              for(Gb=0; Gb < nirreps; Gb++) {
 
1133
                Gd = Gb ^ Gik; /* totally symmetric */
 
1134
                Gca = Gjd = Gj ^ Gd; /* totally symmetric */
 
1135
 
 
1136
                nrows = virtpi[Gd]; 
 
1137
                ncols = ZdMAe.params->coltot[Gjd];
 
1138
                nlinks = virtpi[Gb];
 
1139
 
 
1140
                bd = T2AB.col_offset[Gik][Gb];
 
1141
                jd = ZdMAe.row_offset[Gjd][J];
 
1142
                ZdMAe.matrix[Gjd] = dpd_block_matrix(nrows, ncols);
 
1143
                dpd_buf4_mat_irrep_rd_block(&ZdMAe, Gjd, jd, nrows);
 
1144
 
 
1145
                if(nrows && ncols && nlinks)
 
1146
                  C_DGEMM('t', 'n', nrows, ncols, nlinks, 1.0, &(T2AB.matrix[Gik][ik][bd]), nrows,
 
1147
                          W2[Gb][0], ncols, 1.0, ZdMAe.matrix[Gjd][0], ncols);
 
1148
 
 
1149
                dpd_buf4_mat_irrep_wrt_block(&ZdMAe, Gjd, jd, nrows);
 
1150
                dpd_free_block(ZdMAe.matrix[Gjd], nrows, ncols);
 
1151
              }
 
1152
 
 
1153
              /* Z_IJAM <-- -1/2 L_IJkABc t_MkBc */
 
1154
              /* sort W(AB,C) to W(A,BC) */
 
1155
              for(Gab=0; Gab < nirreps; Gab++) {
 
1156
                Gc = Gab ^ Gijk;
 
1157
                for(ab=0; ab < FAA.params->coltot[Gab]; ab++) {
 
1158
                  A = FAA.params->colorb[Gab][ab][0];
 
1159
                  B = FAA.params->colorb[Gab][ab][1];
 
1160
                  Ga = FAA.params->rsym[A];
 
1161
                  a = A - vir_off[Ga];
 
1162
                  for(c=0; c < virtpi[Gc]; c++) {
 
1163
                    C = vir_off[Gc] + c;
 
1164
                    bc = FAA.params->colidx[B][C];
 
1165
                    W2[Ga][a][bc] = W1[Gab][ab][c];
 
1166
                  }
 
1167
                }
 
1168
              }
 
1169
 
 
1170
              ij = ZLMAO.params->rowidx[I][J];
 
1171
 
 
1172
              for(Gm=0; Gm < nirreps; Gm++) {
 
1173
                Gbc = Gmk = Gm ^ Gk; /* totally symmetric */
 
1174
                Ga = Gij ^ Gm; /* totally symmetric */
 
1175
 
 
1176
                nrows = virtpi[Ga];
 
1177
                ncols = T2AB.params->coltot[Gmk];
 
1178
 
 
1179
                for(m=0; m < occpi[Gm]; m++) {
 
1180
                  M = occ_off[Gm] + m;
 
1181
                  mk = T2AB.params->rowidx[M][K];
 
1182
                  am = ZLMAO.col_offset[Gij][Ga] + m;
 
1183
 
 
1184
                  if(nrows && ncols)
 
1185
                    C_DGEMV('n', nrows, ncols, -1.0, W2[Ga][0], ncols, T2AB.matrix[Gmk][mk], 1, 
 
1186
                            1.0, &(ZLMAO.matrix[Gij][ij][am]), occpi[Gm]);
 
1187
                }
 
1188
              }
 
1189
 
 
1190
              /* Z_KiCm <-- -1/2 L_ijKabC t_mjab */
 
1191
              ki = ZLmAo.params->rowidx[K][I];
 
1192
 
 
1193
              for(Gm=0; Gm < nirreps; Gm++) {
 
1194
                Gab = Gmj = Gm ^ Gj; /* totally symmetric */
 
1195
                Gc = Gm ^ Gki; /* totally symmetric */
 
1196
 
 
1197
                nrows = T2AA.params->coltot[Gmj];
 
1198
                ncols = virtpi[Gc];
 
1199
 
 
1200
                for(m=0; m < occpi[Gm]; m++) {
 
1201
                  M = occ_off[Gm] + m;
 
1202
                  mj = T2AA.params->rowidx[M][J];
 
1203
                  cm = ZLmAo.col_offset[Gki][Gc] + m;
 
1204
 
 
1205
                  if(nrows && ncols)
 
1206
                    C_DGEMV('t', nrows, ncols, -0.5, W1[Gab][0], ncols, T2AA.matrix[Gmj][mj], 1,
 
1207
                            1.0, &(ZLmAo.matrix[Gki][ki][cm]), occpi[Gm]);
 
1208
                }
 
1209
              }
 
1210
 
 
1211
              /* Z_IkAm <-- - L_IJkABc t_mJcB */
 
1212
              /* sort W(AB,C) to W(A,CB) */
 
1213
              for(Gab=0; Gab < nirreps; Gab++) {
 
1214
                Gc = Gab ^ Gijk;
 
1215
                for(ab=0; ab < FAA.params->coltot[Gab]; ab++) {
 
1216
                  A = FAA.params->colorb[Gab][ab][0];
 
1217
                  B = FAA.params->colorb[Gab][ab][1];
 
1218
                  Ga = FAA.params->rsym[A];
 
1219
                  a = A - vir_off[Ga];
 
1220
                  for(c=0; c < virtpi[Gc]; c++) {
 
1221
                    C = vir_off[Gc] + c;
 
1222
                    cb = FAA.params->colidx[C][B];
 
1223
                    W2[Ga][a][cb] = W1[Gab][ab][c];
 
1224
                  }
 
1225
                }
 
1226
              }
 
1227
 
 
1228
              ik = ZLmAo.params->rowidx[I][K];
 
1229
 
 
1230
              for(Gm=0; Gm < nirreps; Gm++) {
 
1231
                Gbc = Gmj = Gm ^ Gj; /* totally symmetric */
 
1232
                Ga = Gm ^ Gik; /* totally symmetric */
 
1233
 
 
1234
                nrows = virtpi[Ga];
 
1235
                ncols = T2AB.params->coltot[Gmj];
 
1236
 
 
1237
                for(m=0; m < occpi[Gm]; m++) {
 
1238
                  M = occ_off[Gm] + m;
 
1239
                  mj = T2AB.params->rowidx[M][J];
 
1240
                  am = ZLmAo.col_offset[Gik][Ga] + m;
 
1241
 
 
1242
                  if(nrows && ncols)
 
1243
                    C_DGEMV('n', nrows, ncols, -1.0, W2[Ga][0], ncols, T2AB.matrix[Gmj][mj], 1,
 
1244
                            1.0, &(ZLmAo.matrix[Gik][ik][am]), occpi[Gm]);
 
1245
                }
 
1246
              }
 
1247
 
 
1248
              /* Z_IJMB <-- - L_IJkABc t_MkAc */
 
1249
              /* sort W(AB,C) to W(B,AC) */
 
1250
              for(Gab=0; Gab < nirreps; Gab++) {
 
1251
                Gc = Gab ^ Gijk;
 
1252
                for(ab=0; ab < FAA.params->coltot[Gab]; ab++) {
 
1253
                  A = FAA.params->colorb[Gab][ab][0];
 
1254
                  B = FAA.params->colorb[Gab][ab][1];
 
1255
                  Gb = FAA.params->ssym[B];
 
1256
                  b = B - vir_off[Gb];
 
1257
                  for(c=0; c < virtpi[Gc]; c++) {
 
1258
                    C = vir_off[Gc] + c;
 
1259
                    ac = FAA.params->colidx[A][C];
 
1260
                    W2[Gb][b][ac] = W1[Gab][ab][c];
 
1261
                  }
 
1262
                }
 
1263
              }
 
1264
 
 
1265
              ij = ZIMLE.params->rowidx[I][J];
 
1266
 
 
1267
              for(Gm=0; Gm < nirreps; Gm++) {
 
1268
                Gb = Gm ^ Gij; /* totally symmetric */
 
1269
                Gmk = Gm ^ Gk;
 
1270
 
 
1271
                nrows = virtpi[Gb];
 
1272
                ncols = T2AB.params->coltot[Gmk];
 
1273
 
 
1274
                for(m=0; m < occpi[Gm]; m++) {
 
1275
                  M = occ_off[Gm] + m;
 
1276
                  mk = T2AB.params->rowidx[M][K];
 
1277
                  mb = ZIMLE.col_offset[Gij][Gm] + m * virtpi[Gb];
 
1278
 
 
1279
                  if(nrows && ncols)
 
1280
                    C_DGEMV('n', nrows, ncols, -1.0, W2[Gb][0], ncols, T2AB.matrix[Gmk][mk], 1,
 
1281
                            1.0, &(ZIMLE.matrix[Gij][ij][mb]), 1);
 
1282
                }
 
1283
 
 
1284
              }
 
1285
 
 
1286
              /* Z_IkMc <-- -1/2 L_IJkABc t_MJAB */ 
 
1287
              ik = ZImLe.params->rowidx[I][K];
 
1288
              for(Gm=0; Gm < nirreps; Gm++) {
 
1289
                Gc = Gm ^ Gik ; /* totally symmetric */
 
1290
                Gab = Gmj = Gm ^ Gj; /* totally symmetric */
 
1291
 
 
1292
                nrows = T2AA.params->coltot[Gmj];
 
1293
                ncols = virtpi[Gc];
 
1294
 
 
1295
                for(m=0; m < occpi[Gm]; m++) {
 
1296
                  M = occ_off[Gm] + m;
 
1297
                  mj = T2AA.params->rowidx[M][J];
 
1298
                  mc = ZImLe.col_offset[Gik][Gm] + m * virtpi[Gc];
 
1299
 
 
1300
                  if(nrows && ncols)
 
1301
                    C_DGEMV('t', nrows, ncols, -0.5, W1[Gab][0], ncols, T2AA.matrix[Gmj][mj], 1,
 
1302
                            1.0, &(ZImLe.matrix[Gik][ik][mc]), 1);
 
1303
                }
 
1304
              }
 
1305
 
 
1306
              /* Z_KiMa <-- - L_ijKabC t_MjCb */
 
1307
              /* sort W(AB,C) to W(A,CB) */
 
1308
              for(Gab=0; Gab < nirreps; Gab++) {
 
1309
                Gc = Gab ^ Gijk;
 
1310
                for(ab=0; ab < FAA.params->coltot[Gab]; ab++) {
 
1311
                  A = FAA.params->colorb[Gab][ab][0];
 
1312
                  B = FAA.params->colorb[Gab][ab][1];
 
1313
                  Ga = FAA.params->rsym[A];
 
1314
                  a = A - vir_off[Ga];
 
1315
                  for(c=0; c < virtpi[Gc]; c++) {
 
1316
                    C = vir_off[Gc] + c;
 
1317
                    cb = FAA.params->colidx[C][B];
 
1318
                    W2[Ga][a][cb] = W1[Gab][ab][c];
 
1319
                  }
 
1320
                }
 
1321
              }
 
1322
 
 
1323
              ki = ZImLe.params->rowidx[K][I];
 
1324
              for(Gm=0; Gm < nirreps; Gm++) {
 
1325
                Ga = Gm ^ Gki; /* totally symmetric */
 
1326
                Gmj = Gm ^ Gj;
 
1327
 
 
1328
                nrows = virtpi[Ga];
 
1329
                ncols = T2AB.params->coltot[Gmj];
 
1330
 
 
1331
                for(m=0; m < occpi[Gm]; m++) {
 
1332
                  M = occ_off[Gm] + m;
 
1333
                  mj = T2AB.params->rowidx[M][J];
 
1334
                  ma = ZImLe.col_offset[Gki][Gm] + m * virtpi[Ga];
 
1335
 
 
1336
                  if(nrows && ncols)
 
1337
                    C_DGEMV('n', nrows, ncols, -1.0, W2[Ga][0], ncols, T2AB.matrix[Gmj][mj], 1,
 
1338
                            1.0, &(ZImLe.matrix[Gki][ki][ma]), 1);
 
1339
                }
 
1340
              }
 
1341
 
 
1342
              /* Z_KimC <-- 1/2 L_ijKabC t_mjab */
 
1343
              ki = ZImlE.params->rowidx[K][I];
 
1344
              for(Gm=0; Gm < nirreps; Gm++) {
 
1345
                Gc = Gm ^ Gki; /* totally symmetric */
 
1346
                Gab = Gmj = Gm ^ Gj;
 
1347
 
 
1348
                nrows = T2AA.params->coltot[Gmj];
 
1349
                ncols = virtpi[Gc];
 
1350
 
 
1351
                for(m=0; m < occpi[Gm]; m++) {
 
1352
                  M = occ_off[Gm] + m;
 
1353
                  mj = T2AA.params->rowidx[M][J];
 
1354
                  mc = ZImlE.col_offset[Gki][Gm] + m * virtpi[Gc];
 
1355
 
 
1356
                  if(nrows && ncols)
 
1357
                    C_DGEMV('t', nrows, ncols, 0.5, W1[Gab][0], ncols, T2AA.matrix[Gmj][mj], 1,
 
1358
                            1.0, &(ZImlE.matrix[Gki][ki][mc]), 1);
 
1359
                }
 
1360
              }
 
1361
 
 
1362
              /* Z_IkmB <-- - l_IJkABc t_JmAc */
 
1363
              ik = ZImlE.params->rowidx[I][K];
 
1364
              /* sort W(AB,C) to W(B,AC) */
 
1365
              for(Gab=0; Gab < nirreps; Gab++) {
 
1366
                Gc = Gab ^ Gijk;
 
1367
                for(ab=0; ab < FAA.params->coltot[Gab]; ab++) {
 
1368
                  A = FAA.params->colorb[Gab][ab][0];
 
1369
                  B = FAA.params->colorb[Gab][ab][1];
 
1370
                  Gb = FAA.params->ssym[B];
 
1371
                  b = B - vir_off[Gb];
 
1372
                  for(c=0; c < virtpi[Gc]; c++) {
 
1373
                    C = vir_off[Gc] + c;
 
1374
                    ac = FAA.params->colidx[A][C];
 
1375
                    W2[Gb][b][ac] = W1[Gab][ab][c];
 
1376
                  }
 
1377
                }
 
1378
              }
 
1379
 
 
1380
              for(Gm=0; Gm < nirreps; Gm++) {
 
1381
                Gb = Gm ^ Gik; /* totally symmetric */
 
1382
                Gjm = Gm ^ Gj;
 
1383
 
 
1384
                nrows = virtpi[Gb];
 
1385
                ncols = T2AB.params->coltot[Gjm];
 
1386
 
 
1387
                for(m=0; m < occpi[Gm]; m++) {
 
1388
                  M = occ_off[Gm] + m;
 
1389
                  jm = T2AB.params->rowidx[J][M];
 
1390
                  mb = ZImlE.col_offset[Gki][Gm] + m * virtpi[Gb];
 
1391
 
 
1392
                  if(nrows && ncols)
 
1393
                    C_DGEMV('n', nrows, ncols, -1.0, W2[Gb][0], ncols, T2AB.matrix[Gjm][jm], 1,
 
1394
                            1.0, &(ZImlE.matrix[Gik][ik][mb]), 1);
 
1395
                }
 
1396
              }
 
1397
 
 
1398
            } /* k */
 
1399
          } /* j */
 
1400
        } /* i */
 
1401
 
 
1402
        for(Gab=0; Gab < nirreps; Gab++) {
 
1403
          Gc = Gab ^ Gijk; /* totally symmetric */
 
1404
          dpd_free_block(W1[Gab], FAA.params->coltot[Gab], virtpi[Gc]);
 
1405
        }
 
1406
        for(Ga=0; Ga < nirreps; Ga++) {
 
1407
          Gcb = Ga ^ Gijk; /* assumes totally symmetric */
 
1408
          dpd_free_block(W2[Ga], virtpi[Ga], WmAfE.params->coltot[Gcb]);
 
1409
        }
 
1410
 
 
1411
      } /* Gk */
 
1412
    } /* Gj */
 
1413
  } /* Gi */
 
1414
 
 
1415
  free(W1);
 
1416
  free(W2);
 
1417
 
 
1418
  dpd_buf4_close(&EAA);
 
1419
  dpd_buf4_close(&EAB);
 
1420
  dpd_buf4_close(&EBA);
 
1421
  dpd_buf4_close(&FAA);
 
1422
  dpd_buf4_close(&FAB);
 
1423
  dpd_buf4_close(&FBA);
 
1424
  dpd_buf4_close(&L2AA);
 
1425
  dpd_buf4_close(&L2AB);
 
1426
  dpd_buf4_close(&L2BA);
 
1427
 
 
1428
  dpd_file2_close(&fIJ);
 
1429
  dpd_file2_close(&fAB);
 
1430
  dpd_file2_close(&fij);
 
1431
  dpd_file2_close(&fab);
 
1432
 
 
1433
  dpd_file2_close(&FME);
 
1434
  dpd_file2_close(&Fme);
 
1435
  dpd_file2_close(&LIA);
 
1436
  dpd_file2_close(&Lia);
 
1437
 
 
1438
  dpd_buf4_close(&DAAints);
 
1439
  dpd_buf4_close(&DABints);
 
1440
  dpd_buf4_close(&LIJAB);
 
1441
  dpd_buf4_close(&LIjAb);
 
1442
 
 
1443
  dpd_buf4_close(&WmAfE);
 
1444
  dpd_buf4_close(&WMAFE);
 
1445
  dpd_buf4_close(&WMaFe);
 
1446
 
 
1447
  for(h=0; h < nirreps; h++) dpd_buf4_mat_irrep_close(&WMnIe, h);
 
1448
  for(h=0; h < nirreps; h++) dpd_buf4_mat_irrep_close(&WMNIE, h);
 
1449
  for(h=0; h < nirreps; h++) dpd_buf4_mat_irrep_close(&WmNiE, h);
 
1450
  dpd_buf4_close(&WMnIe);
 
1451
  dpd_buf4_close(&WMNIE);
 
1452
  dpd_buf4_close(&WmNiE);
 
1453
 
 
1454
  dpd_buf4_close(&ZIgDe);
 
1455
  dpd_buf4_close(&ZIGDE);
 
1456
 
 
1457
  dpd_buf4_close(&ZDMAE);
 
1458
  dpd_buf4_close(&ZDmAe);
 
1459
  dpd_buf4_close(&ZdMAe);
 
1460
 
 
1461
  for(h=0; h < nirreps; h++) {
 
1462
    dpd_buf4_mat_irrep_wrt(&ZLMAO, h);
 
1463
    dpd_buf4_mat_irrep_close(&ZLMAO, h);
 
1464
  }
 
1465
  dpd_buf4_close(&ZLMAO);
 
1466
  for(h=0; h < nirreps; h++) {
 
1467
    dpd_buf4_mat_irrep_wrt(&ZLmAo, h);
 
1468
    dpd_buf4_mat_irrep_close(&ZLmAo, h);
 
1469
  }
 
1470
  dpd_buf4_close(&ZLmAo);
 
1471
 
 
1472
  for(h=0; h < nirreps; h++) {
 
1473
    dpd_buf4_mat_irrep_wrt(&ZIMLE, h);
 
1474
    dpd_buf4_mat_irrep_close(&ZIMLE, h);
 
1475
  }
 
1476
  dpd_buf4_close(&ZIMLE);
 
1477
  for(h=0; h < nirreps; h++) {
 
1478
    dpd_buf4_mat_irrep_wrt(&ZImLe, h);
 
1479
    dpd_buf4_mat_irrep_close(&ZImLe, h);
 
1480
  }
 
1481
  dpd_buf4_close(&ZImLe);
 
1482
  for(h=0; h < nirreps; h++) {
 
1483
    dpd_buf4_mat_irrep_wrt(&ZImlE, h);
 
1484
    dpd_buf4_mat_irrep_close(&ZImlE, h);
 
1485
  }
 
1486
  dpd_buf4_close(&ZImlE);
 
1487
 
 
1488
  for(h=0; h < nirreps; h++) dpd_buf4_mat_irrep_close(&T2AB, h);
 
1489
  for(h=0; h < nirreps; h++) dpd_buf4_mat_irrep_close(&T2AA, h);
 
1490
  dpd_buf4_close(&T2AB);
 
1491
  dpd_buf4_close(&T2AA);
 
1492
 
 
1493
  for(h=0; h < nirreps; h++) {
 
1494
    dpd_buf4_mat_irrep_wrt(&L2AAnew, h);
 
1495
    dpd_buf4_mat_irrep_close(&L2AAnew, h);
 
1496
  }
 
1497
  dpd_buf4_init(&D2, CC_DENOM, 0, 0, 5, 0, 5, 0, "dIjAb");
 
1498
  dpd_buf4_dirprd(&D2, &L2AAnew);
 
1499
  dpd_buf4_close(&D2);
 
1500
  dpd_buf4_init(&L2, CC_LAMBDA, 0, 0, 5, 2, 7, 0, "New LIJAB");
 
1501
  dpd_buf4_axpy(&L2AAnew, &L2, 1);
 
1502
  dpd_buf4_close(&L2);
 
1503
  dpd_buf4_close(&L2AAnew);
 
1504
 
 
1505
  for(h=0; h < nirreps; h++) {
 
1506
    dpd_buf4_mat_irrep_wrt(&L2ABnew, h);
 
1507
    dpd_buf4_mat_irrep_close(&L2ABnew, h);
 
1508
  }
 
1509
  dpd_buf4_init(&D2, CC_DENOM, 0, 0, 5, 0, 5, 0, "dIjAb");
 
1510
  dpd_buf4_dirprd(&D2, &L2ABnew);
 
1511
  dpd_buf4_close(&D2);
 
1512
  dpd_buf4_init(&L2, CC_LAMBDA, 0, 0, 5, 0, 5, 0, "New LIjAb");
 
1513
  dpd_buf4_axpy(&L2ABnew, &L2, 1);
 
1514
  dpd_buf4_close(&L2);
 
1515
  dpd_buf4_close(&L2ABnew);
 
1516
 
 
1517
  /* Spin adaptation will remove this.  And yes, this means that all the above
 
1518
     calculations for LIJAB were pointless... -TDC */
 
1519
  dpd_buf4_init(&L2, CC_LAMBDA, 0, 2, 7, 0, 5, 1, "New LIjAb");
 
1520
  dpd_buf4_copy(&L2, CC_LAMBDA, "New LIJAB");
 
1521
  dpd_buf4_copy(&L2, CC_LAMBDA, "New Lijab");
 
1522
  dpd_buf4_close(&L2);
 
1523
 
 
1524
}
 
1525
 
 
1526
}} // namespace psi::cclambda