← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/ccdensity/Gijab_ROHF.c

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2006-09-10 14:01:33 UTC
  • Revision ID: james.westby@ubuntu.com-20060910140133-ib2j86trekykfsfv
Tags: upstream-3.2.3
Import upstream version 3.2.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/ccdensity/Gijab_ROHF.c
 
1
#include <stdio.h>
 
2
#include <libdpd/dpd.h>
 
3
#define EXTERN
 
4
#include "globals.h"
 
5
 
 
6
void Gijab_ROHF(void)
 
7
{
 
8
  int h, nirreps, i, a, m, e, I, A, M, E, Isym, Asym, Msym, Esym, row, col;
 
9
  int j, b, J, B, Jsym, Bsym;
 
10
  double value;
 
11
  dpdfile2 T1, L1, g, ZZ, ZZ2, T1A, T1B;
 
12
  dpdbuf4 G, L, T, V, Z, Z1, Z2;
 
13
 
 
14
  nirreps = moinfo.nirreps;
 
15
 
 
16
  /* ( g(I,M) + L(M,E) T(I,E) ) --> Z(I,M)(TMP0)  */
 
17
  dpd_file2_init(&g, CC_OEI, 0, 0, 0, "GMI");
 
18
  dpd_file2_copy(&g, CC_TMP0, "Z(I,M)");
 
19
  dpd_file2_close(&g);
 
20
  dpd_file2_init(&ZZ, CC_TMP0, 0, 0, 0, "Z(I,M)");
 
21
  dpd_file2_init(&L1, CC_OEI, 0, 0, 1, "LIA");
 
22
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
23
  dpd_contract222(&T1, &L1, &ZZ, 0, 0, 1.0, 1.0);
 
24
  dpd_file2_close(&T1);
 
25
  dpd_file2_close(&L1);
 
26
  dpd_file2_close(&ZZ);
 
27
 
 
28
  /* ( g(i,m) + L(m,e) T(i,e) ) --> Z(i,m)(TMP1)  */
 
29
  dpd_file2_init(&g, CC_OEI, 0, 0, 0, "Gmi");
 
30
  dpd_file2_copy(&g, CC_TMP1, "Z(i,m)");
 
31
  dpd_file2_close(&g);
 
32
  dpd_file2_init(&ZZ, CC_TMP1, 0, 0, 0, "Z(i,m)");
 
33
  dpd_file2_init(&L1, CC_OEI, 0, 0, 1, "Lia");
 
34
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
 
35
  dpd_contract222(&T1, &L1, &ZZ, 0, 0, 1.0, 1.0);
 
36
  dpd_file2_close(&T1);
 
37
  dpd_file2_close(&L1);
 
38
  dpd_file2_close(&ZZ);
 
39
 
 
40
  /* ( g(E,A) - L(M,E) T(M,A) ) --> Z(E,A)(TMP2) */
 
41
  dpd_file2_init(&g, CC_OEI, 0, 1, 1, "GAE");
 
42
  dpd_file2_copy(&g, CC_TMP2, "Z(E,A)");
 
43
  dpd_file2_close(&g);
 
44
  dpd_file2_init(&ZZ, CC_TMP2, 0, 1, 1, "Z(E,A)");
 
45
  dpd_file2_init(&L1, CC_OEI, 0, 0, 1, "LIA");
 
46
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
47
  dpd_contract222(&L1, &T1, &ZZ, 1, 1, -1.0, 1.0);
 
48
  dpd_file2_close(&T1);
 
49
  dpd_file2_close(&L1);
 
50
  dpd_file2_close(&ZZ);
 
51
 
 
52
  /* ( g(e,a) - L(m,e) T(m,a) ) --> Z(e,a)(TMP3) */
 
53
  dpd_file2_init(&g, CC_OEI, 0, 1, 1, "Gae");
 
54
  dpd_file2_copy(&g, CC_TMP3, "Z(e,a)");
 
55
  dpd_file2_close(&g);
 
56
  dpd_file2_init(&ZZ, CC_TMP3, 0, 1, 1, "Z(e,a)");
 
57
  dpd_file2_init(&L1, CC_OEI, 0, 0, 1, "Lia");
 
58
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
 
59
  dpd_contract222(&L1, &T1, &ZZ, 1, 1, -1.0, 1.0);
 
60
  dpd_file2_close(&T1);
 
61
  dpd_file2_close(&L1);
 
62
  dpd_file2_close(&ZZ);
 
63
 
 
64
  dpd_file2_init(&T1A, CC_OEI, 0, 0, 1, "tIA");
 
65
  dpd_file2_mat_init(&T1A);
 
66
  dpd_file2_mat_rd(&T1A);
 
67
  dpd_file2_init(&T1B, CC_OEI, 0, 0, 1, "tia");
 
68
  dpd_file2_mat_init(&T1B);
 
69
  dpd_file2_mat_rd(&T1B);
 
70
 
 
71
  /* ( - T(IA,ME) + 2 * T(I,E) T(M,A) ) --> Z(IA,ME) */
 
72
  dpd_buf4_init(&T, CC_TAMPS, 0, 10, 10, 10, 10, 0, "tIAJB");
 
73
  dpd_buf4_copy(&T, CC_TMP4, "Z(IA,ME)");
 
74
  dpd_buf4_close(&T);
 
75
  dpd_buf4_init(&Z, CC_TMP4, 0, 10, 10, 10, 10, 0, "Z(IA,ME)");
 
76
  dpd_buf4_scm(&Z, -1.0);
 
77
  for(h=0; h < nirreps; h++) {
 
78
    dpd_buf4_mat_irrep_init(&Z, h); 0,
 
79
                                      dpd_buf4_mat_irrep_rd(&Z, h);
 
80
    for(row=0; row < Z.params->rowtot[h]; row++) {
 
81
      i = Z.params->roworb[h][row][0];
 
82
      a = Z.params->roworb[h][row][1];
 
83
      I = T1A.params->rowidx[i];  Isym = T1A.params->psym[i];
 
84
      A = T1A.params->colidx[a];  Asym = T1A.params->qsym[a];
 
85
          
 
86
      for(col=0; col < Z.params->coltot[h]; col++) {
 
87
        m = Z.params->colorb[h][col][0];
 
88
        e = Z.params->colorb[h][col][1];
 
89
        M = T1A.params->rowidx[m];  Msym = T1A.params->psym[m];
 
90
        E = T1A.params->colidx[e];  Esym = T1A.params->qsym[e];
 
91
 
 
92
        if((Isym==Esym) && (Msym==Asym))
 
93
          Z.matrix[h][row][col] += (2* T1A.matrix[Isym][I][E] *
 
94
                                    T1A.matrix[Msym][M][A]);
 
95
      }
 
96
    }
 
97
    dpd_buf4_mat_irrep_wrt(&Z, h);
 
98
    dpd_buf4_mat_irrep_close(&Z, h);
 
99
  }
 
100
  dpd_buf4_close(&Z);
 
101
 
 
102
  /* ( - T(ia,me) + 2 * T(i,e) T(m,a) ) --> Z(ia,me) */
 
103
  dpd_buf4_init(&T, CC_TAMPS, 0, 10, 10, 10, 10, 0, "tiajb");
 
104
  dpd_buf4_copy(&T, CC_TMP5, "Z(ia,me)");
 
105
  dpd_buf4_close(&T);
 
106
  dpd_buf4_init(&Z, CC_TMP5, 0, 10, 10, 10, 10, 0, "Z(ia,me)");
 
107
  dpd_buf4_scm(&Z, -1.0);
 
108
  for(h=0; h < nirreps; h++) {
 
109
    dpd_buf4_mat_irrep_init(&Z, h); 0,
 
110
                                      dpd_buf4_mat_irrep_rd(&Z, h);
 
111
    for(row=0; row < Z.params->rowtot[h]; row++) {
 
112
      i = Z.params->roworb[h][row][0];
 
113
      a = Z.params->roworb[h][row][1];
 
114
      I = T1B.params->rowidx[i];  Isym = T1B.params->psym[i];
 
115
      A = T1B.params->colidx[a];  Asym = T1B.params->qsym[a];
 
116
          
 
117
      for(col=0; col < Z.params->coltot[h]; col++) {
 
118
        m = Z.params->colorb[h][col][0];
 
119
        e = Z.params->colorb[h][col][1];
 
120
        M = T1B.params->rowidx[m];  Msym = T1B.params->psym[m];
 
121
        E = T1B.params->colidx[e];  Esym = T1B.params->qsym[e];
 
122
 
 
123
        if((Isym==Esym) && (Msym==Asym))
 
124
          Z.matrix[h][row][col] += (2* T1B.matrix[Isym][I][E] *
 
125
                                    T1B.matrix[Msym][M][A]);
 
126
      }
 
127
    }
 
128
    dpd_buf4_mat_irrep_wrt(&Z, h);
 
129
    dpd_buf4_mat_irrep_close(&Z, h);
 
130
  }
 
131
  dpd_buf4_close(&Z);
 
132
 
 
133
  /* ( - T(iA,Me) + 2 * T(i,e) T(M,A) ) --> Z(iA,Me) */
 
134
  dpd_buf4_init(&T, CC_TAMPS, 0, 10, 10, 10, 10, 0, "tjAIb");
 
135
  dpd_buf4_copy(&T, CC_TMP6, "Z(iA,Me)");
 
136
  dpd_buf4_close(&T);
 
137
  dpd_buf4_init(&Z, CC_TMP6, 0, 10, 10, 10, 10, 0, "Z(iA,Me)");
 
138
  for(h=0; h < nirreps; h++) {
 
139
    dpd_buf4_mat_irrep_init(&Z, h); 0,
 
140
                                      dpd_buf4_mat_irrep_rd(&Z, h);
 
141
    for(row=0; row < Z.params->rowtot[h]; row++) {
 
142
      i = Z.params->roworb[h][row][0];
 
143
      a = Z.params->roworb[h][row][1];
 
144
      I = T1B.params->rowidx[i];  Isym = T1B.params->psym[i];
 
145
      A = T1A.params->colidx[a];  Asym = T1A.params->qsym[a];
 
146
          
 
147
      for(col=0; col < Z.params->coltot[h]; col++) {
 
148
        m = Z.params->colorb[h][col][0];
 
149
        e = Z.params->colorb[h][col][1];
 
150
        M = T1A.params->rowidx[m];  Msym = T1A.params->psym[m];
 
151
        E = T1B.params->colidx[e];  Esym = T1B.params->qsym[e];
 
152
 
 
153
        if((Isym==Esym) && (Msym==Asym))
 
154
          Z.matrix[h][row][col] += (2* T1B.matrix[Isym][I][E] *
 
155
                                    T1A.matrix[Msym][M][A]);
 
156
      }
 
157
    }
 
158
    dpd_buf4_mat_irrep_wrt(&Z, h);
 
159
    dpd_buf4_mat_irrep_close(&Z, h);
 
160
  }
 
161
  dpd_buf4_close(&Z);
 
162
 
 
163
  /* ( - T(Ia,mE) + 2 * T(I,E) T(m,a) ) --> Z(Ia,mE) */
 
164
  dpd_buf4_init(&T, CC_TAMPS, 0, 10, 10, 10, 10, 0, "tIbjA");
 
165
  dpd_buf4_copy(&T, CC_TMP7, "Z(Ia,mE)");
 
166
  dpd_buf4_close(&T);
 
167
  dpd_buf4_init(&Z, CC_TMP7, 0, 10, 10, 10, 10, 0, "Z(Ia,mE)");
 
168
  for(h=0; h < nirreps; h++) {
 
169
    dpd_buf4_mat_irrep_init(&Z, h); 0,
 
170
                                      dpd_buf4_mat_irrep_rd(&Z, h);
 
171
    for(row=0; row < Z.params->rowtot[h]; row++) {
 
172
      i = Z.params->roworb[h][row][0];
 
173
      a = Z.params->roworb[h][row][1];
 
174
      I = T1A.params->rowidx[i];  Isym = T1A.params->psym[i];
 
175
      A = T1B.params->colidx[a];  Asym = T1B.params->qsym[a];
 
176
          
 
177
      for(col=0; col < Z.params->coltot[h]; col++) {
 
178
        m = Z.params->colorb[h][col][0];
 
179
        e = Z.params->colorb[h][col][1];
 
180
        M = T1B.params->rowidx[m];  Msym = T1B.params->psym[m];
 
181
        E = T1A.params->colidx[e];  Esym = T1A.params->qsym[e];
 
182
 
 
183
        if((Isym==Esym) && (Msym==Asym))
 
184
          Z.matrix[h][row][col] += (2* T1A.matrix[Isym][I][E] *
 
185
                                    T1B.matrix[Msym][M][A]);
 
186
      }
 
187
    }
 
188
    dpd_buf4_mat_irrep_wrt(&Z, h);
 
189
    dpd_buf4_mat_irrep_close(&Z, h);
 
190
  }
 
191
  dpd_buf4_close(&Z);
 
192
 
 
193
  dpd_file2_mat_close(&T1A);
 
194
  dpd_file2_close(&T1A);
 
195
  dpd_file2_mat_close(&T1B);
 
196
  dpd_file2_close(&T1B);
 
197
 
 
198
  /* L(IJ,AB) */
 
199
  dpd_buf4_init(&L, CC_LAMPS, 0, 2, 7, 2, 7, 0, "LIJAB");
 
200
  dpd_buf4_copy(&L, CC_GAMMA, "GIJAB");
 
201
  dpd_buf4_close(&L);
 
202
  dpd_buf4_init(&G, CC_GAMMA, 0, 2, 7, 2, 7, 0, "GIJAB");
 
203
  /* Tau(IJ,AB) */
 
204
  dpd_buf4_init(&T, CC_TAMPS, 0, 2, 7, 2, 7, 0, "tauIJAB");
 
205
  dpd_buf4_axpy(&T, &G, 1.0);
 
206
  dpd_buf4_close(&T);
 
207
  /* V(IJ,MN) Tau(MN,AB) */
 
208
  dpd_buf4_init(&T, CC_TAMPS, 0, 2, 7, 2, 7, 0, "tauIJAB");
 
209
  dpd_buf4_init(&V, CC_MISC, 0, 2, 2, 2, 2, 0, "VMNIJ");
 
210
  dpd_contract444(&V, &T, &G, 0, 1, 1.0, 1.0);
 
211
  dpd_buf4_close(&V);
 
212
  dpd_buf4_close(&T);
 
213
  dpd_buf4_close(&G);
 
214
  /* - ( Z(I,M) Tau(MJ,AB) - Z(J,M) Tau(MI,AB) ) */
 
215
  dpd_buf4_init(&Z1, CC_TMP8, 0, 0, 7, 0, 7, 0, "Z1(IJ,AB)");
 
216
  dpd_buf4_init(&T, CC_TAMPS, 0, 0, 7, 2, 7, 0, "tauIJAB");
 
217
  dpd_file2_init(&ZZ, CC_TMP0, 0, 0, 0, "Z(I,M)");
 
218
  dpd_contract244(&ZZ, &T, &Z1, 1, 0, 0, -1.0, 0.0);
 
219
  dpd_file2_close(&ZZ);
 
220
  dpd_buf4_close(&T);
 
221
  dpd_buf4_sort(&Z1, CC_TMP9, qprs, 0, 7, "Z2(JI,AB)");
 
222
  dpd_buf4_init(&Z2, CC_TMP9, 0, 0, 7, 0, 7, 0, "Z2(JI,AB)");
 
223
  dpd_buf4_axpy(&Z2, &Z1, -1.0);
 
224
  dpd_buf4_close(&Z2);
 
225
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 7, 2, 7, 0, "GIJAB");
 
226
  dpd_buf4_axpy(&Z1, &G, 1.0);
 
227
  dpd_buf4_close(&Z1);
 
228
  dpd_buf4_close(&G);
 
229
  /* - ( Z(E,A) Tau(IJ,BE) - Z(E,B) Tau(IJ,AE) ) */
 
230
  dpd_buf4_init(&Z1, CC_TMP8, 0, 2, 5, 2, 5, 0, "ZZ1(IJ,AB)");
 
231
  dpd_buf4_init(&T, CC_TAMPS, 0, 2, 5, 2, 7, 0, "tauIJAB");
 
232
  dpd_file2_init(&ZZ, CC_TMP2, 0, 1, 1, "Z(E,A)");
 
233
  dpd_contract424(&T, &ZZ, &Z1, 3, 0, 0, 1.0, 0.0);
 
234
  dpd_file2_close(&ZZ);
 
235
  dpd_buf4_close(&T);
 
236
  dpd_buf4_sort(&Z1, CC_TMP9, pqsr, 2, 5, "Z2(IJ,BA)");
 
237
  dpd_buf4_init(&Z2, CC_TMP9, 0, 2, 5, 2, 5, 0, "Z2(IJ,BA)");
 
238
  dpd_buf4_axpy(&Z2, &Z1, -1.0);
 
239
  dpd_buf4_close(&Z2);
 
240
  dpd_buf4_init(&G, CC_GAMMA, 0, 2, 5, 2, 7, 0, "GIJAB");
 
241
  dpd_buf4_axpy(&Z1, &G, 1.0);
 
242
  dpd_buf4_close(&Z1);
 
243
  dpd_buf4_close(&G);
 
244
  /* - P(IJ) P(AB) ( T'(IA,ME) (TMP4) V(JB,ME) + T(IA,me) (T2) V(JB,me) ) */
 
245
  dpd_buf4_init(&Z, CC_TMP8, 0, 10, 10, 10, 10, 0, "Z(IA,JB)");
 
246
  dpd_buf4_init(&T, CC_TMP4, 0, 10, 10, 10, 10, 0, "Z(IA,ME)");
 
247
  dpd_buf4_init(&V, CC_MISC, 0, 10, 10, 10, 10, 0, "VIAJB");
 
248
  dpd_contract444(&T, &V, &Z, 0, 0, 1.0, 0.0);
 
249
  dpd_buf4_close(&V);
 
250
  dpd_buf4_close(&T);
 
251
  dpd_buf4_init(&T, CC_TAMPS, 0, 10, 10, 10, 10, 0, "tIAjb");
 
252
  dpd_buf4_init(&V, CC_MISC, 0, 10, 10, 10, 10, 0, "VIAjb");
 
253
  dpd_contract444(&T, &V, &Z, 0, 0, -1.0, 1.0);
 
254
  dpd_buf4_close(&V);
 
255
  dpd_buf4_close(&T);
 
256
  dpd_buf4_sort(&Z, CC_TMP9, rqps, 10, 10, "Z(JA,IB)");
 
257
  dpd_buf4_sort(&Z, CC_TMP10, psrq, 10, 10, "Z(IB,JA)");
 
258
  dpd_buf4_sort(&Z, CC_TMP11, rspq, 10, 10, "Z(JB,IA)");
 
259
  dpd_buf4_init(&Z1, CC_TMP9, 0, 10, 10, 10, 10, 0, "Z(JA,IB)");
 
260
  dpd_buf4_axpy(&Z1, &Z, -1.0);
 
261
  dpd_buf4_close(&Z1);
 
262
  dpd_buf4_init(&Z1, CC_TMP10, 0, 10, 10, 10, 10, 0, "Z(IB,JA)");
 
263
  dpd_buf4_axpy(&Z1, &Z, -1.0);
 
264
  dpd_buf4_close(&Z1);
 
265
  dpd_buf4_init(&Z1, CC_TMP11, 0, 10, 10, 10, 10, 0, "Z(JB,IA)");
 
266
  dpd_buf4_axpy(&Z1, &Z, 1.0);
 
267
  dpd_buf4_close(&Z1);
 
268
  dpd_buf4_sort(&Z, CC_TMP9, prqs, 0, 5, "Z(IJ,AB)");
 
269
  dpd_buf4_close(&Z);
 
270
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 5, 2, 7, 0, "GIJAB");
 
271
  dpd_buf4_init(&Z, CC_TMP9, 0, 0, 5, 0, 5, 0, "Z(IJ,AB)");
 
272
  /* I don't understand this factor of 1/2 that shows up here */
 
273
  dpd_buf4_axpy(&Z, &G, -0.5);
 
274
  dpd_buf4_close(&Z);
 
275
  dpd_buf4_close(&G);
 
276
  /* T'(IA,ME) (TMP4) L(M,E) + T'(IA,me) (T2) L(m,e) --> ZZ(I,A) */
 
277
  dpd_file2_init(&ZZ, CC_TMP8, 0, 0, 1, "ZZ(I,A)");
 
278
  dpd_file2_init(&L1, CC_OEI, 0, 0, 1, "LIA");
 
279
  dpd_buf4_init(&T, CC_TMP4, 0, 10, 10, 10, 10, 0, "Z(IA,ME)");
 
280
  dpd_contract422(&T, &L1, &ZZ, 0, 0, 1.0, 0.0);
 
281
  dpd_buf4_close(&T);
 
282
  dpd_file2_close(&L1);
 
283
  dpd_file2_init(&L1, CC_OEI, 0, 0, 1, "Lia");
 
284
  dpd_buf4_init(&T, CC_TAMPS, 0, 10, 10, 10, 10, 0, "tIAjb");
 
285
  dpd_contract422(&T, &L1, &ZZ, 0, 0, -1.0, 1.0);
 
286
  dpd_buf4_close(&T);
 
287
  dpd_file2_close(&L1);
 
288
  dpd_file2_close(&ZZ);
 
289
  /* - P(IJ) P(AB) ZZ(I,A) T(J,B) --> G(IJ,AB) */
 
290
  dpd_file2_init(&ZZ, CC_TMP8, 0, 0, 1, "ZZ(I,A)");
 
291
  dpd_file2_mat_init(&ZZ);
 
292
  dpd_file2_mat_rd(&ZZ);
 
293
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
294
  dpd_file2_mat_init(&T1);
 
295
  dpd_file2_mat_rd(&T1);
 
296
 
 
297
  dpd_buf4_init(&G, CC_GAMMA, 0, 2, 7, 2, 7, 0, "GIJAB");
 
298
  for(h=0; h < nirreps; h++) {
 
299
    dpd_buf4_mat_irrep_init(&G, h); 0,
 
300
                                      dpd_buf4_mat_irrep_rd(&G, h);
 
301
    for(row=0; row < G.params->rowtot[h]; row++) {
 
302
      i = G.params->roworb[h][row][0];
 
303
      j = G.params->roworb[h][row][1];
 
304
          
 
305
      for(col=0; col < G.params->coltot[h]; col++) {
 
306
        a = G.params->colorb[h][col][0];
 
307
        b = G.params->colorb[h][col][1];
 
308
 
 
309
        value = 0.0;
 
310
 
 
311
        I = ZZ.params->rowidx[i]; Isym = ZZ.params->psym[i];
 
312
        J = T1.params->rowidx[j]; Jsym = T1.params->psym[j];
 
313
        A = ZZ.params->colidx[a]; Asym = ZZ.params->qsym[a];
 
314
        B = T1.params->colidx[b]; Bsym = T1.params->qsym[b];
 
315
 
 
316
        if((Isym==Asym) && (Jsym==Bsym))
 
317
          value += ZZ.matrix[Isym][I][A] * T1.matrix[Jsym][J][B];
 
318
 
 
319
        I = T1.params->rowidx[i]; Isym = T1.params->psym[i];
 
320
        J = ZZ.params->rowidx[j]; Jsym = ZZ.params->psym[j];
 
321
 
 
322
        if((Jsym==Asym) && (Isym==Bsym))
 
323
          value -= ZZ.matrix[Jsym][J][A] * T1.matrix[Isym][I][B];
 
324
 
 
325
        I = ZZ.params->rowidx[i]; Isym = ZZ.params->psym[i];
 
326
        J = T1.params->rowidx[j]; Jsym = T1.params->psym[j];
 
327
        A = T1.params->colidx[a]; Asym = T1.params->qsym[a];
 
328
        B = ZZ.params->colidx[b]; Bsym = ZZ.params->qsym[b];
 
329
 
 
330
        if((Isym==Bsym) && (Jsym==Asym))
 
331
          value -= ZZ.matrix[Isym][I][B] * T1.matrix[Jsym][J][A];
 
332
 
 
333
        I = T1.params->rowidx[i]; Isym = T1.params->psym[i];
 
334
        J = ZZ.params->rowidx[j]; Jsym = ZZ.params->psym[j];
 
335
 
 
336
        if((Isym==Asym) && (Jsym==Bsym))
 
337
          value += T1.matrix[Isym][I][A] * ZZ.matrix[Jsym][J][B];
 
338
 
 
339
        G.matrix[h][row][col] -= value;
 
340
              
 
341
      }
 
342
    }
 
343
    dpd_buf4_mat_irrep_wrt(&G, h);
 
344
    dpd_buf4_mat_irrep_close(&G, h);
 
345
  }
 
346
  dpd_buf4_close(&G);
 
347
 
 
348
  dpd_file2_mat_close(&T1);
 
349
  dpd_file2_close(&T1);
 
350
  dpd_file2_mat_close(&ZZ);
 
351
  dpd_file2_close(&ZZ);
 
352
 
 
353
  /* T(J,E) L(M,E) --> ZZ(J,M) */
 
354
  dpd_file2_init(&ZZ, CC_TMP8, 0, 0, 0, "Z(J,M)");
 
355
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
356
  dpd_file2_init(&L1, CC_OEI, 0, 0, 1, "LIA");
 
357
  dpd_contract222(&T1, &L1, &ZZ, 0, 0, 1.0, 0.0);
 
358
  dpd_file2_close(&L1);
 
359
  dpd_file2_close(&T1);
 
360
  /* ZZ(J,M) T(M,B) --> ZZ2(J,B) */
 
361
  dpd_file2_init(&ZZ2, CC_TMP9, 0, 0, 1, "ZZ2(J,B)");
 
362
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
363
  dpd_contract222(&ZZ, &T1, &ZZ2, 0, 1, 1.0, 0.0);
 
364
  dpd_file2_close(&T1);
 
365
  dpd_file2_close(&ZZ);
 
366
  dpd_file2_close(&ZZ2);
 
367
 
 
368
  /* 3 P(IJ) P(AB) T(I,A) ZZ(J,B) --> G(IJ,AB) */
 
369
  dpd_file2_init(&ZZ, CC_TMP9, 0, 0, 1, "ZZ2(J,B)");
 
370
  dpd_file2_mat_init(&ZZ);
 
371
  dpd_file2_mat_rd(&ZZ);
 
372
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
373
  dpd_file2_mat_init(&T1);
 
374
  dpd_file2_mat_rd(&T1);
 
375
 
 
376
  dpd_buf4_init(&G, CC_GAMMA, 0, 2, 7, 2, 7, 0, "GIJAB");
 
377
  for(h=0; h < nirreps; h++) {
 
378
    dpd_buf4_mat_irrep_init(&G, h); 0,
 
379
                                      dpd_buf4_mat_irrep_rd(&G, h);
 
380
    for(row=0; row < G.params->rowtot[h]; row++) {
 
381
      i = G.params->roworb[h][row][0];
 
382
      j = G.params->roworb[h][row][1];
 
383
          
 
384
      for(col=0; col < G.params->coltot[h]; col++) {
 
385
        a = G.params->colorb[h][col][0];
 
386
        b = G.params->colorb[h][col][1];
 
387
 
 
388
        value = 0.0;
 
389
 
 
390
        I = T1.params->rowidx[i]; Isym = T1.params->psym[i];
 
391
        J = ZZ.params->rowidx[j]; Jsym = ZZ.params->psym[j];
 
392
        A = T1.params->colidx[a]; Asym = T1.params->qsym[a];
 
393
        B = ZZ.params->colidx[b]; Bsym = ZZ.params->qsym[b];
 
394
 
 
395
        if((Isym==Asym) && (Jsym==Bsym))
 
396
          value += T1.matrix[Isym][I][A] * ZZ.matrix[Jsym][J][B];
 
397
 
 
398
        I = ZZ.params->rowidx[i]; Isym = ZZ.params->psym[i];
 
399
        J = T1.params->rowidx[j]; Jsym = T1.params->psym[j];
 
400
 
 
401
        if((Jsym==Asym) && (Isym==Bsym))
 
402
          value -= T1.matrix[Jsym][J][A] * ZZ.matrix[Isym][I][B];
 
403
 
 
404
        I = T1.params->rowidx[i]; Isym = T1.params->psym[i];
 
405
        J = ZZ.params->rowidx[j]; Jsym = ZZ.params->psym[j];
 
406
        A = ZZ.params->colidx[a]; Asym = ZZ.params->qsym[a];
 
407
        B = T1.params->colidx[b]; Bsym = T1.params->qsym[b];
 
408
 
 
409
        if((Isym==Bsym) && (Jsym==Asym))
 
410
          value -= T1.matrix[Isym][I][B] * ZZ.matrix[Jsym][J][A];
 
411
 
 
412
        I = ZZ.params->rowidx[i]; Isym = ZZ.params->psym[i];
 
413
        J = T1.params->rowidx[j]; Jsym = T1.params->psym[j];
 
414
 
 
415
        if((Isym==Asym) && (Jsym==Bsym))
 
416
          value += ZZ.matrix[Isym][I][A] * T1.matrix[Jsym][J][B];
 
417
 
 
418
        G.matrix[h][row][col] += 3.0 * value;
 
419
              
 
420
      }
 
421
    }
 
422
    dpd_buf4_mat_irrep_wrt(&G, h);
 
423
    dpd_buf4_mat_irrep_close(&G, h);
 
424
  }
 
425
  dpd_buf4_scm(&G, 0.5);
 
426
  dpd_buf4_close(&G);
 
427
 
 
428
  dpd_file2_mat_close(&T1);
 
429
  dpd_file2_close(&T1);
 
430
  dpd_file2_mat_close(&ZZ);
 
431
  dpd_file2_close(&ZZ);
 
432
 
 
433
 
 
434
 
 
435
  /* L(ij,ab) */
 
436
  dpd_buf4_init(&L, CC_LAMPS, 0, 2, 7, 2, 7, 0, "Lijab");
 
437
  dpd_buf4_copy(&L, CC_GAMMA, "Gijab");
 
438
  dpd_buf4_close(&L);
 
439
  dpd_buf4_init(&G, CC_GAMMA, 0, 2, 7, 2, 7, 0, "Gijab");
 
440
  /* Tau(ij,ab) */
 
441
  dpd_buf4_init(&T, CC_TAMPS, 0, 2, 7, 2, 7, 0, "tauijab");
 
442
  dpd_buf4_axpy(&T, &G, 1.0);
 
443
  dpd_buf4_close(&T);
 
444
  /* V(ij,mn) Tau(mn,ab) */
 
445
  dpd_buf4_init(&T, CC_TAMPS, 0, 2, 7, 2, 7, 0, "tauijab");
 
446
  dpd_buf4_init(&V, CC_MISC, 0, 2, 2, 2, 2, 0, "Vmnij");
 
447
  dpd_contract444(&V, &T, &G, 0, 1, 1.0, 1.0);
 
448
  dpd_buf4_close(&V);
 
449
  dpd_buf4_close(&T);
 
450
  dpd_buf4_close(&G);
 
451
  /* - ( Z(i,m) Tau(mj,ab) - Z(j,m) Tau(mi,ab) ) */
 
452
  dpd_buf4_init(&Z1, CC_TMP8, 0, 0, 7, 0, 7, 0, "Z1(ij,ab)");
 
453
  dpd_buf4_init(&T, CC_TAMPS, 0, 0, 7, 2, 7, 0, "tauijab");
 
454
  dpd_file2_init(&ZZ, CC_TMP1, 0, 0, 0, "Z(i,m)");
 
455
  dpd_contract244(&ZZ, &T, &Z1, 1, 0, 0, -1.0, 0.0);
 
456
  dpd_file2_close(&ZZ);
 
457
  dpd_buf4_close(&T);
 
458
  dpd_buf4_sort(&Z1, CC_TMP9, qprs, 0, 7, "Z2(ji,ab)");
 
459
  dpd_buf4_init(&Z2, CC_TMP9, 0, 0, 7, 0, 7, 0, "Z2(ji,ab)");
 
460
  dpd_buf4_axpy(&Z2, &Z1, -1.0);
 
461
  dpd_buf4_close(&Z2);
 
462
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 7, 2, 7, 0, "Gijab");
 
463
  dpd_buf4_axpy(&Z1, &G, 1.0);
 
464
  dpd_buf4_close(&Z1);
 
465
  dpd_buf4_close(&G);
 
466
  /* - ( Z(e,a) Tau(ij,be) - Z(e,b) Tau(ij,ae) ) */
 
467
  dpd_buf4_init(&Z1, CC_TMP8, 0, 2, 5, 2, 5, 0, "ZZ1(ij,ab)");
 
468
  dpd_buf4_init(&T, CC_TAMPS, 0, 2, 5, 2, 7, 0, "tauijab");
 
469
  dpd_file2_init(&ZZ, CC_TMP3, 0, 1, 1, "Z(e,a)");
 
470
  dpd_contract424(&T, &ZZ, &Z1, 3, 0, 0, 1.0, 0.0);
 
471
  dpd_file2_close(&ZZ);
 
472
  dpd_buf4_close(&T);
 
473
  dpd_buf4_sort(&Z1, CC_TMP9, pqsr, 2, 5, "Z2(ij,ba)");
 
474
  dpd_buf4_init(&Z2, CC_TMP9, 0, 2, 5, 2, 5, 0, "Z2(ij,ba)");
 
475
  dpd_buf4_axpy(&Z2, &Z1, -1.0);
 
476
  dpd_buf4_close(&Z2);
 
477
  dpd_buf4_init(&G, CC_GAMMA, 0, 2, 5, 2, 7, 0, "Gijab");
 
478
  dpd_buf4_axpy(&Z1, &G, 1.0);
 
479
  dpd_buf4_close(&Z1);
 
480
  dpd_buf4_close(&G);
 
481
  /* - P(ij) P(ab) ( T'(ia,me) (TMP5) V(jb,me) + T(ia,ME) (T2) V(jb,ME) ) */
 
482
  dpd_buf4_init(&Z, CC_TMP8, 0, 10, 10, 10, 10, 0, "Z(ia,jb)");
 
483
  dpd_buf4_init(&T, CC_TMP5, 0, 10, 10, 10, 10, 0, "Z(ia,me)");
 
484
  dpd_buf4_init(&V, CC_MISC, 0, 10, 10, 10, 10, 0, "Viajb");
 
485
  dpd_contract444(&T, &V, &Z, 0, 0, 1.0, 0.0);
 
486
  dpd_buf4_close(&V);
 
487
  dpd_buf4_close(&T);
 
488
  dpd_buf4_init(&T, CC_TAMPS, 0, 10, 10, 10, 10, 0, "tiaJB");
 
489
  dpd_buf4_init(&V, CC_MISC, 0, 10, 10, 10, 10, 0, "ViaJB");
 
490
  dpd_contract444(&T, &V, &Z, 0, 0, -1.0, 1.0);
 
491
  dpd_buf4_close(&V);
 
492
  dpd_buf4_close(&T);
 
493
  dpd_buf4_sort(&Z, CC_TMP9, rqps, 10, 10, "Z(ja,ib)");
 
494
  dpd_buf4_sort(&Z, CC_TMP10, psrq, 10, 10, "Z(ib,ja)");
 
495
  dpd_buf4_sort(&Z, CC_TMP11, rspq, 10, 10, "Z(jb,ia)");
 
496
  dpd_buf4_init(&Z1, CC_TMP9, 0, 10, 10, 10, 10, 0, "Z(ja,ib)");
 
497
  dpd_buf4_axpy(&Z1, &Z, -1.0);
 
498
  dpd_buf4_close(&Z1);
 
499
  dpd_buf4_init(&Z1, CC_TMP10, 0, 10, 10, 10, 10, 0, "Z(ib,ja)");
 
500
  dpd_buf4_axpy(&Z1, &Z, -1.0);
 
501
  dpd_buf4_close(&Z1);
 
502
  dpd_buf4_init(&Z1, CC_TMP11, 0, 10, 10, 10, 10, 0, "Z(jb,ia)");
 
503
  dpd_buf4_axpy(&Z1, &Z, 1.0);
 
504
  dpd_buf4_close(&Z1);
 
505
  dpd_buf4_sort(&Z, CC_TMP9, prqs, 0, 5, "Z(ij,ab)");
 
506
  dpd_buf4_close(&Z);
 
507
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 5, 2, 7, 0, "Gijab");
 
508
  dpd_buf4_init(&Z, CC_TMP9, 0, 0, 5, 0, 5, 0, "Z(ij,ab)");
 
509
  /* I don't understand this factor of 1/2 that shows up here */
 
510
  dpd_buf4_axpy(&Z, &G, -0.5);
 
511
  dpd_buf4_close(&Z);
 
512
  dpd_buf4_close(&G);
 
513
  /* T'(ia,me) (TMP5) L(m,e) + T'(ia,ME) (T2) L(M,E) --> ZZ(i,a) */
 
514
  dpd_file2_init(&ZZ, CC_TMP8, 0, 0, 1, "ZZ(i,a)");
 
515
  dpd_file2_init(&L1, CC_OEI, 0, 0, 1, "Lia");
 
516
  dpd_buf4_init(&T, CC_TMP5, 0, 10, 10, 10, 10, 0, "Z(ia,me)");
 
517
  dpd_contract422(&T, &L1, &ZZ, 0, 0, 1.0, 0.0);
 
518
  dpd_buf4_close(&T);
 
519
  dpd_file2_close(&L1);
 
520
  dpd_file2_init(&L1, CC_OEI, 0, 0, 1, "LIA");
 
521
  dpd_buf4_init(&T, CC_TAMPS, 0, 10, 10, 10, 10, 0, "tiaJB");
 
522
  dpd_contract422(&T, &L1, &ZZ, 0, 0, -1.0, 1.0);
 
523
  dpd_buf4_close(&T);
 
524
  dpd_file2_close(&L1);
 
525
  dpd_file2_close(&ZZ);
 
526
  /* - P(ij) P(ab) ZZ(i,a) T(j,b) --> G(ij,ab) */
 
527
  dpd_file2_init(&ZZ, CC_TMP8, 0, 0, 1, "ZZ(i,a)");
 
528
  dpd_file2_mat_init(&ZZ);
 
529
  dpd_file2_mat_rd(&ZZ);
 
530
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
 
531
  dpd_file2_mat_init(&T1);
 
532
  dpd_file2_mat_rd(&T1);
 
533
 
 
534
  dpd_buf4_init(&G, CC_GAMMA, 0, 2, 7, 2, 7, 0, "Gijab");
 
535
  for(h=0; h < nirreps; h++) {
 
536
    dpd_buf4_mat_irrep_init(&G, h); 0,
 
537
                                      dpd_buf4_mat_irrep_rd(&G, h);
 
538
    for(row=0; row < G.params->rowtot[h]; row++) {
 
539
      i = G.params->roworb[h][row][0];
 
540
      j = G.params->roworb[h][row][1];
 
541
          
 
542
      for(col=0; col < G.params->coltot[h]; col++) {
 
543
        a = G.params->colorb[h][col][0];
 
544
        b = G.params->colorb[h][col][1];
 
545
 
 
546
        value = 0.0;
 
547
 
 
548
        I = ZZ.params->rowidx[i]; Isym = ZZ.params->psym[i];
 
549
        J = T1.params->rowidx[j]; Jsym = T1.params->psym[j];
 
550
        A = ZZ.params->colidx[a]; Asym = ZZ.params->qsym[a];
 
551
        B = T1.params->colidx[b]; Bsym = T1.params->qsym[b];
 
552
 
 
553
        if((Isym==Asym) && (Jsym==Bsym))
 
554
          value += ZZ.matrix[Isym][I][A] * T1.matrix[Jsym][J][B];
 
555
 
 
556
        I = T1.params->rowidx[i]; Isym = T1.params->psym[i];
 
557
        J = ZZ.params->rowidx[j]; Jsym = ZZ.params->psym[j];
 
558
 
 
559
        if((Jsym==Asym) && (Isym==Bsym))
 
560
          value -= ZZ.matrix[Jsym][J][A] * T1.matrix[Isym][I][B];
 
561
 
 
562
        I = ZZ.params->rowidx[i]; Isym = ZZ.params->psym[i];
 
563
        J = T1.params->rowidx[j]; Jsym = T1.params->psym[j];
 
564
        A = T1.params->colidx[a]; Asym = T1.params->qsym[a];
 
565
        B = ZZ.params->colidx[b]; Bsym = ZZ.params->qsym[b];
 
566
 
 
567
        if((Isym==Bsym) && (Jsym==Asym))
 
568
          value -= ZZ.matrix[Isym][I][B] * T1.matrix[Jsym][J][A];
 
569
 
 
570
        I = T1.params->rowidx[i]; Isym = T1.params->psym[i];
 
571
        J = ZZ.params->rowidx[j]; Jsym = ZZ.params->psym[j];
 
572
 
 
573
        if((Isym==Asym) && (Jsym==Bsym))
 
574
          value += T1.matrix[Isym][I][A] * ZZ.matrix[Jsym][J][B];
 
575
 
 
576
        G.matrix[h][row][col] -= value;
 
577
              
 
578
      }
 
579
    }
 
580
    dpd_buf4_mat_irrep_wrt(&G, h);
 
581
    dpd_buf4_mat_irrep_close(&G, h);
 
582
  }
 
583
  dpd_buf4_close(&G);
 
584
 
 
585
  dpd_file2_mat_close(&T1);
 
586
  dpd_file2_close(&T1);
 
587
  dpd_file2_mat_close(&ZZ);
 
588
  dpd_file2_close(&ZZ);
 
589
 
 
590
  /* T(j,e) L(m,e) --> ZZ(j,m) */
 
591
  dpd_file2_init(&ZZ, CC_TMP8, 0, 0, 0, "Z(j,m)");
 
592
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
 
593
  dpd_file2_init(&L1, CC_OEI, 0, 0, 1, "Lia");
 
594
  dpd_contract222(&T1, &L1, &ZZ, 0, 0, 1.0, 0.0);
 
595
  dpd_file2_close(&L1);
 
596
  dpd_file2_close(&T1);
 
597
  /* ZZ(j,m) T(m,b) --> ZZ2(j,b) */
 
598
  dpd_file2_init(&ZZ2, CC_TMP9, 0, 0, 1, "ZZ2(j,b)");
 
599
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
 
600
  dpd_contract222(&ZZ, &T1, &ZZ2, 0, 1, 1.0, 0.0);
 
601
  dpd_file2_close(&T1);
 
602
  dpd_file2_close(&ZZ);
 
603
  dpd_file2_close(&ZZ2);
 
604
 
 
605
  /* 3 P(ij) P(ab) T(i,a) ZZ(j,b) --> G(ij,ab) */
 
606
  dpd_file2_init(&ZZ, CC_TMP9, 0, 0, 1, "ZZ2(j,b)");
 
607
  dpd_file2_mat_init(&ZZ);
 
608
  dpd_file2_mat_rd(&ZZ);
 
609
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
 
610
  dpd_file2_mat_init(&T1);
 
611
  dpd_file2_mat_rd(&T1);
 
612
 
 
613
  dpd_buf4_init(&G, CC_GAMMA, 0, 2, 7, 2, 7, 0, "Gijab");
 
614
  for(h=0; h < nirreps; h++) {
 
615
    dpd_buf4_mat_irrep_init(&G, h); 0,
 
616
                                      dpd_buf4_mat_irrep_rd(&G, h);
 
617
    for(row=0; row < G.params->rowtot[h]; row++) {
 
618
      i = G.params->roworb[h][row][0];
 
619
      j = G.params->roworb[h][row][1];
 
620
          
 
621
      for(col=0; col < G.params->coltot[h]; col++) {
 
622
        a = G.params->colorb[h][col][0];
 
623
        b = G.params->colorb[h][col][1];
 
624
 
 
625
        value = 0.0;
 
626
 
 
627
        I = T1.params->rowidx[i]; Isym = T1.params->psym[i];
 
628
        J = ZZ.params->rowidx[j]; Jsym = ZZ.params->psym[j];
 
629
        A = T1.params->colidx[a]; Asym = T1.params->qsym[a];
 
630
        B = ZZ.params->colidx[b]; Bsym = ZZ.params->qsym[b];
 
631
 
 
632
        if((Isym==Asym) && (Jsym==Bsym))
 
633
          value += T1.matrix[Isym][I][A] * ZZ.matrix[Jsym][J][B];
 
634
 
 
635
        I = ZZ.params->rowidx[i]; Isym = ZZ.params->psym[i];
 
636
        J = T1.params->rowidx[j]; Jsym = T1.params->psym[j];
 
637
 
 
638
        if((Jsym==Asym) && (Isym==Bsym))
 
639
          value -= T1.matrix[Jsym][J][A] * ZZ.matrix[Isym][I][B];
 
640
 
 
641
        I = T1.params->rowidx[i]; Isym = T1.params->psym[i];
 
642
        J = ZZ.params->rowidx[j]; Jsym = ZZ.params->psym[j];
 
643
        A = ZZ.params->colidx[a]; Asym = ZZ.params->qsym[a];
 
644
        B = T1.params->colidx[b]; Bsym = T1.params->qsym[b];
 
645
 
 
646
        if((Isym==Bsym) && (Jsym==Asym))
 
647
          value -= T1.matrix[Isym][I][B] * ZZ.matrix[Jsym][J][A];
 
648
 
 
649
        I = ZZ.params->rowidx[i]; Isym = ZZ.params->psym[i];
 
650
        J = T1.params->rowidx[j]; Jsym = T1.params->psym[j];
 
651
 
 
652
        if((Isym==Asym) && (Jsym==Bsym))
 
653
          value += ZZ.matrix[Isym][I][A] * T1.matrix[Jsym][J][B];
 
654
 
 
655
        G.matrix[h][row][col] += 3.0 * value;
 
656
              
 
657
      }
 
658
    }
 
659
    dpd_buf4_mat_irrep_wrt(&G, h);
 
660
    dpd_buf4_mat_irrep_close(&G, h);
 
661
  }
 
662
  dpd_buf4_scm(&G, 0.5);
 
663
  dpd_buf4_close(&G);
 
664
 
 
665
  dpd_file2_mat_close(&T1);
 
666
  dpd_file2_close(&T1);
 
667
  dpd_file2_mat_close(&ZZ);
 
668
  dpd_file2_close(&ZZ);
 
669
 
 
670
 
 
671
 
 
672
  /* L(Ij,Ab) */
 
673
  dpd_buf4_init(&L, CC_LAMPS, 0, 0, 5, 0, 5, 0, "LIjAb");
 
674
  dpd_buf4_copy(&L, CC_GAMMA, "GIjAb");
 
675
  dpd_buf4_close(&L);
 
676
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 5, 0, 5, 0, "GIjAb");
 
677
  /* Tau(Ij,Ab) */
 
678
  dpd_buf4_init(&T, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tauIjAb");
 
679
  dpd_buf4_axpy(&T, &G, 1.0);
 
680
  dpd_buf4_close(&T);
 
681
  /* V(Ij,Mn) Tau(Mn,Ab) */
 
682
  dpd_buf4_init(&T, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tauIjAb");
 
683
  dpd_buf4_init(&V, CC_MISC, 0, 0, 0, 0, 0, 0, "VMnIj");
 
684
  dpd_contract444(&V, &T, &G, 0, 1, 1.0, 1.0);
 
685
  dpd_buf4_close(&V);
 
686
  dpd_buf4_close(&T);
 
687
  dpd_buf4_close(&G);
 
688
  /* - ( Z(I,M) Tau(Mj,Ab) - Z(j,m) Tau(mI,bA) ) */
 
689
  dpd_buf4_init(&Z1, CC_TMP8, 0, 0, 5, 0, 5, 0, "Z1(Ij,Ab)");
 
690
  dpd_buf4_init(&T, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tauIjAb");
 
691
  dpd_file2_init(&ZZ, CC_TMP0, 0, 0, 0, "Z(I,M)");
 
692
  dpd_contract244(&ZZ, &T, &Z1, 1, 0, 0, 1.0, 0.0);
 
693
  dpd_file2_close(&ZZ);
 
694
  dpd_buf4_close(&T);
 
695
  dpd_buf4_init(&Z2, CC_TMP9, 0, 0, 5, 0, 5, 0, "Z2(jI,bA)");
 
696
  dpd_buf4_init(&T, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tauiJaB");
 
697
  dpd_file2_init(&ZZ, CC_TMP1, 0, 0, 0, "Z(i,m)");
 
698
  dpd_contract244(&ZZ, &T, &Z2, 1, 0, 0, 1.0, 0.0);
 
699
  dpd_file2_close(&ZZ);
 
700
  dpd_buf4_close(&T);
 
701
  dpd_buf4_sort(&Z2, CC_TMP10, qprs, 0, 5, "Z2(Ij,bA)");
 
702
  dpd_buf4_close(&Z2);
 
703
  dpd_buf4_init(&Z2, CC_TMP10, 0, 0, 5, 0, 5, 0, "Z2(Ij,bA)");
 
704
  dpd_buf4_sort(&Z2, CC_TMP9, pqsr, 0, 5, "Z2(Ij,Ab)");
 
705
  dpd_buf4_close(&Z2);
 
706
  dpd_buf4_init(&Z2, CC_TMP9, 0, 0, 5, 0, 5, 0, "Z2(Ij,Ab)");
 
707
  dpd_buf4_axpy(&Z2, &Z1, 1.0);
 
708
  dpd_buf4_close(&Z2);
 
709
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 5, 0, 5, 0, "GIjAb");
 
710
  dpd_buf4_axpy(&Z1, &G, -1.0);
 
711
  dpd_buf4_close(&Z1);
 
712
  dpd_buf4_close(&G);
 
713
  /* - ( Z(E,A) Tau(Ij,bE) - Z(e,b) Tau(Ij,Ae) ) */
 
714
  dpd_buf4_init(&Z1, CC_TMP8, 0, 0, 5, 0, 5, 0, "ZZ1(Ij,Ab)");
 
715
  dpd_buf4_init(&T, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tauIjAb");
 
716
  dpd_file2_init(&ZZ, CC_TMP3, 0, 1, 1, "Z(e,a)");
 
717
  dpd_contract424(&T, &ZZ, &Z1, 3, 0, 0, 1.0, 0.0);
 
718
  dpd_file2_close(&ZZ);
 
719
  dpd_buf4_close(&T);
 
720
  dpd_buf4_init(&Z2, CC_TMP9, 0, 0, 5, 0, 5, 0, "Z2(jI,bA)");
 
721
  dpd_buf4_init(&T, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tauiJaB");
 
722
  dpd_file2_init(&ZZ, CC_TMP2, 0, 1, 1, "Z(E,A)");
 
723
  dpd_contract424(&T, &ZZ, &Z2, 3, 0, 0, 1.0, 0.0);
 
724
  dpd_file2_close(&ZZ);
 
725
  dpd_buf4_close(&T);
 
726
  dpd_buf4_sort(&Z2, CC_TMP10, qprs, 0, 5, "Z2(Ij,bA)");
 
727
  dpd_buf4_close(&Z2);
 
728
  dpd_buf4_init(&Z2, CC_TMP10, 0, 0, 5, 0, 5, 0, "Z2(Ij,bA)");
 
729
  dpd_buf4_sort(&Z2, CC_TMP9, pqsr, 0, 5, "Z2(Ij,Ab)");
 
730
  dpd_buf4_close(&Z2);
 
731
  dpd_buf4_init(&Z2, CC_TMP9, 0, 0, 5, 0, 5, 0, "Z2(Ij,Ab)");
 
732
  dpd_buf4_axpy(&Z2, &Z1, 1.0);
 
733
  dpd_buf4_close(&Z2);
 
734
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 5, 0, 5, 0, "GIjAb");
 
735
  dpd_buf4_axpy(&Z1, &G, 1.0);
 
736
  dpd_buf4_close(&Z1);
 
737
  dpd_buf4_close(&G);
 
738
  /* - P(Ij) P(Ab) ( T'(IA,me) (T2) V(jb,me) + T'(IA,ME) (TMP4) V(jb,ME) ) */
 
739
  dpd_buf4_init(&Z, CC_TMP8, 0, 10, 10, 10, 10, 0, "Z(IA,jb)");
 
740
  dpd_buf4_init(&T, CC_TMP4, 0, 10, 10, 10, 10, 0, "Z(IA,ME)");
 
741
  dpd_buf4_init(&V, CC_MISC, 0, 10, 10, 10, 10, 0, "ViaJB");
 
742
  dpd_contract444(&T, &V, &Z, 0, 0, 1.0, 0.0);
 
743
  dpd_buf4_close(&V);
 
744
  dpd_buf4_close(&T);
 
745
  dpd_buf4_init(&T, CC_TAMPS, 0, 10, 10, 10, 10, 0, "tIAjb");
 
746
  dpd_buf4_init(&V, CC_MISC, 0, 10, 10, 10, 10, 0, "Viajb");
 
747
  dpd_contract444(&T, &V, &Z, 0, 0, -1.0, 1.0);
 
748
  dpd_buf4_close(&V);
 
749
  dpd_buf4_close(&T);
 
750
  /* T'(jA,Me) V(Ib,Me) */
 
751
  dpd_buf4_init(&Z1, CC_TMP9, 0, 10, 10, 10, 10, 0, "Z(jA,Ib)");
 
752
  dpd_buf4_init(&T, CC_TMP6, 0, 10, 10, 10, 10, 0, "Z(iA,Me)");
 
753
  dpd_buf4_init(&V, CC_MISC, 0, 10, 10, 10, 10, 0, "VIaJb");
 
754
  dpd_contract444(&T, &V, &Z1, 0, 0, 1.0, 0.0);
 
755
  dpd_buf4_close(&V);
 
756
  dpd_buf4_close(&T);
 
757
  dpd_buf4_sort(&Z1, CC_TMP10, rqps, 10, 10, "Z(IA,jb)");
 
758
  dpd_buf4_close(&Z1);
 
759
  dpd_buf4_init(&Z1, CC_TMP10, 0, 10, 10, 10, 10, 0, "Z(IA,jb)");
 
760
  dpd_buf4_axpy(&Z1, &Z, -1.0);
 
761
  dpd_buf4_close(&Z1);
 
762
  /* T'(Ib,mE) V(jA,mE) */
 
763
  dpd_buf4_init(&Z1, CC_TMP9, 0, 10, 10, 10, 10, 0, "Z(Ib,jA)");
 
764
  dpd_buf4_init(&T, CC_TMP7, 0, 10, 10, 10, 10, 0, "Z(Ia,mE)");
 
765
  dpd_buf4_init(&V, CC_MISC, 0, 10, 10, 10, 10, 0, "ViAjB");
 
766
  dpd_contract444(&T, &V, &Z1, 0, 0, 1.0, 0.0);
 
767
  dpd_buf4_close(&V);
 
768
  dpd_buf4_close(&T);
 
769
  dpd_buf4_sort(&Z1, CC_TMP10, psrq, 10, 10, "Z(IA,jb)");
 
770
  dpd_buf4_close(&Z1);
 
771
  dpd_buf4_init(&Z1, CC_TMP10, 0, 10, 10, 10, 10, 0, "Z(IA,jb)");
 
772
  dpd_buf4_axpy(&Z1, &Z, -1.0);
 
773
  dpd_buf4_close(&Z1);
 
774
  /* T'(jb,ME) (T2) V(IA,ME) + T'(jb,me) (TMP5) V(IA,me) */
 
775
  dpd_buf4_init(&Z1, CC_TMP9, 0, 10, 10, 10, 10, 0, "Z(IA,jb)");
 
776
  dpd_buf4_init(&T, CC_TMP5, 0, 10, 10, 10, 10, 0, "Z(ia,me)");
 
777
  dpd_buf4_init(&V, CC_MISC, 0, 10, 10, 10, 10, 0, "VIAjb");
 
778
  dpd_contract444(&V, &T, &Z1, 0, 0, 1.0, 0.0);
 
779
  dpd_buf4_close(&V);
 
780
  dpd_buf4_close(&T);
 
781
  dpd_buf4_init(&T, CC_TAMPS, 0, 10, 10, 10, 10, 0, "tiaJB");
 
782
  dpd_buf4_init(&V, CC_MISC, 0, 10, 10, 10, 10, 0, "VIAJB");
 
783
  dpd_contract444(&V, &T, &Z1, 0, 0, -1.0, 1.0);
 
784
  dpd_buf4_close(&V);
 
785
  dpd_buf4_close(&T);
 
786
  dpd_buf4_axpy(&Z1, &Z, 1.0);
 
787
  dpd_buf4_close(&Z1);
 
788
  /* - Z(IA,jb) --> G(Ij,Ab) */
 
789
  dpd_buf4_sort(&Z, CC_TMP9, prqs, 0, 5, "Z(Ij,Ab)");
 
790
  dpd_buf4_close(&Z);
 
791
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 5, 0, 5, 0, "GIjAb");
 
792
  dpd_buf4_init(&Z, CC_TMP9, 0, 0, 5, 0, 5, 0, "Z(Ij,Ab)");
 
793
  dpd_buf4_axpy(&Z, &G, -0.5);
 
794
  dpd_buf4_close(&Z);
 
795
  dpd_buf4_close(&G);
 
796
  /* T'(IA,me) (T2) L(m,e) + T'(IA,ME) (TMP4) L(M,E) --> ZZ(I,A) */
 
797
  dpd_file2_init(&ZZ, CC_TMP8, 0, 0, 1, "ZZ(I,A)");
 
798
  dpd_file2_init(&L1, CC_OEI, 0, 0, 1, "LIA");
 
799
  dpd_buf4_init(&T, CC_TMP4, 0, 10, 10, 10, 10, 0, "Z(IA,ME)");
 
800
  dpd_contract422(&T, &L1, &ZZ, 0, 0, 1.0, 0.0);
 
801
  dpd_buf4_close(&T);
 
802
  dpd_file2_close(&L1);
 
803
  dpd_file2_init(&L1, CC_OEI, 0, 0, 1, "Lia");
 
804
  dpd_buf4_init(&T, CC_TAMPS, 0, 10, 10, 10, 10, 0, "tIAjb");
 
805
  dpd_contract422(&T, &L1, &ZZ, 0, 0, -1.0, 1.0);
 
806
  dpd_buf4_close(&T);
 
807
  dpd_file2_close(&L1);
 
808
  dpd_file2_close(&ZZ);
 
809
  /* - ZZ(I,A) T(j,b) --> G(Ij,Ab) */
 
810
  dpd_file2_init(&ZZ, CC_TMP8, 0, 0, 1, "ZZ(I,A)");
 
811
  dpd_file2_mat_init(&ZZ);
 
812
  dpd_file2_mat_rd(&ZZ);
 
813
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
 
814
  dpd_file2_mat_init(&T1);
 
815
  dpd_file2_mat_rd(&T1);
 
816
 
 
817
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 5, 0, 5, 0, "GIjAb");
 
818
  for(h=0; h < nirreps; h++) {
 
819
    dpd_buf4_mat_irrep_init(&G, h); 0,
 
820
                                      dpd_buf4_mat_irrep_rd(&G, h);
 
821
    for(row=0; row < G.params->rowtot[h]; row++) {
 
822
      i = G.params->roworb[h][row][0];
 
823
      j = G.params->roworb[h][row][1];
 
824
          
 
825
      for(col=0; col < G.params->coltot[h]; col++) {
 
826
        a = G.params->colorb[h][col][0];
 
827
        b = G.params->colorb[h][col][1];
 
828
 
 
829
        value = 0.0;
 
830
 
 
831
        I = ZZ.params->rowidx[i]; Isym = ZZ.params->psym[i];
 
832
        J = T1.params->rowidx[j]; Jsym = T1.params->psym[j];
 
833
        A = ZZ.params->colidx[a]; Asym = ZZ.params->qsym[a];
 
834
        B = T1.params->colidx[b]; Bsym = T1.params->qsym[b];
 
835
 
 
836
        if((Isym==Asym) && (Jsym==Bsym))
 
837
          value += ZZ.matrix[Isym][I][A] * T1.matrix[Jsym][J][B];
 
838
 
 
839
        G.matrix[h][row][col] -= value;
 
840
              
 
841
      }
 
842
    }
 
843
    dpd_buf4_mat_irrep_wrt(&G, h);
 
844
    dpd_buf4_mat_irrep_close(&G, h);
 
845
  }
 
846
  dpd_buf4_close(&G);
 
847
 
 
848
  dpd_file2_mat_close(&T1);
 
849
  dpd_file2_close(&T1);
 
850
  dpd_file2_mat_close(&ZZ);
 
851
  dpd_file2_close(&ZZ);
 
852
 
 
853
  /* T'(jb,ME) (T2) L(M,E) + T'(jb,me) (TMP5) L(m,e) --> ZZ(j,b) */
 
854
  dpd_file2_init(&ZZ, CC_TMP8, 0, 0, 1, "ZZ(j,b)");
 
855
  dpd_file2_init(&L1, CC_OEI, 0, 0, 1, "LIA");
 
856
  dpd_buf4_init(&T, CC_TAMPS, 0, 10, 10, 10, 10, 0, "tiaJB");
 
857
  dpd_contract422(&T, &L1, &ZZ, 0, 0, -1.0, 0.0);
 
858
  dpd_buf4_close(&T);
 
859
  dpd_file2_close(&L1);
 
860
  dpd_file2_init(&L1, CC_OEI, 0, 0, 1, "Lia");
 
861
  dpd_buf4_init(&T, CC_TMP5, 0, 10, 10, 10, 10, 0, "Z(ia,me)");
 
862
  dpd_contract422(&T, &L1, &ZZ, 0, 0, 1.0, 1.0);
 
863
  dpd_buf4_close(&T);
 
864
  dpd_file2_close(&L1);
 
865
  dpd_file2_close(&ZZ);
 
866
  /* - ZZ(j,b) T(I,A) --> G(Ij,Ab) */
 
867
  dpd_file2_init(&ZZ, CC_TMP8, 0, 0, 1, "ZZ(j,b)");
 
868
  dpd_file2_mat_init(&ZZ);
 
869
  dpd_file2_mat_rd(&ZZ);
 
870
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
871
  dpd_file2_mat_init(&T1);
 
872
  dpd_file2_mat_rd(&T1);
 
873
 
 
874
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 5, 0, 5, 0, "GIjAb");
 
875
  for(h=0; h < nirreps; h++) {
 
876
    dpd_buf4_mat_irrep_init(&G, h); 0,
 
877
                                      dpd_buf4_mat_irrep_rd(&G, h);
 
878
    for(row=0; row < G.params->rowtot[h]; row++) {
 
879
      i = G.params->roworb[h][row][0];
 
880
      j = G.params->roworb[h][row][1];
 
881
          
 
882
      for(col=0; col < G.params->coltot[h]; col++) {
 
883
        a = G.params->colorb[h][col][0];
 
884
        b = G.params->colorb[h][col][1];
 
885
 
 
886
        value = 0.0;
 
887
 
 
888
        I = T1.params->rowidx[i]; Isym = T1.params->psym[i];
 
889
        J = ZZ.params->rowidx[j]; Jsym = ZZ.params->psym[j];
 
890
        A = T1.params->colidx[a]; Asym = T1.params->qsym[a];
 
891
        B = ZZ.params->colidx[b]; Bsym = ZZ.params->qsym[b];
 
892
 
 
893
        if((Isym==Asym) && (Jsym==Bsym))
 
894
          value += T1.matrix[Isym][I][A] * ZZ.matrix[Jsym][J][B];
 
895
 
 
896
        G.matrix[h][row][col] -= value;
 
897
              
 
898
      }
 
899
    }
 
900
    dpd_buf4_mat_irrep_wrt(&G, h);
 
901
    dpd_buf4_mat_irrep_close(&G, h);
 
902
  }
 
903
  dpd_buf4_close(&G);
 
904
 
 
905
  dpd_file2_mat_close(&T1);
 
906
  dpd_file2_close(&T1);
 
907
  dpd_file2_mat_close(&ZZ);
 
908
  dpd_file2_close(&ZZ);
 
909
 
 
910
  /* T(j,e) L(m,e) --> ZZ(j,m) */
 
911
  dpd_file2_init(&ZZ, CC_TMP8, 0, 0, 0, "Z(j,m)");
 
912
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
 
913
  dpd_file2_init(&L1, CC_OEI, 0, 0, 1, "Lia");
 
914
  dpd_contract222(&T1, &L1, &ZZ, 0, 0, 1.0, 0.0);
 
915
  dpd_file2_close(&L1);
 
916
  dpd_file2_close(&T1);
 
917
  /* ZZ(j,m) T(m,b) --> ZZ2(j,b) */
 
918
  dpd_file2_init(&ZZ2, CC_TMP9, 0, 0, 1, "ZZ2(j,b)");
 
919
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
 
920
  dpd_contract222(&ZZ, &T1, &ZZ2, 0, 1, 1.0, 0.0);
 
921
  dpd_file2_close(&T1);
 
922
  dpd_file2_close(&ZZ);
 
923
  dpd_file2_close(&ZZ2);
 
924
 
 
925
  /* 3 T(I,A) ZZ(j,b) --> G(Ij,Ab) */
 
926
  dpd_file2_init(&ZZ, CC_TMP9, 0, 0, 1, "ZZ2(j,b)");
 
927
  dpd_file2_mat_init(&ZZ);
 
928
  dpd_file2_mat_rd(&ZZ);
 
929
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
930
  dpd_file2_mat_init(&T1);
 
931
  dpd_file2_mat_rd(&T1);
 
932
 
 
933
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 5, 0, 5, 0, "GIjAb");
 
934
  for(h=0; h < nirreps; h++) {
 
935
    dpd_buf4_mat_irrep_init(&G, h); 0,
 
936
                                      dpd_buf4_mat_irrep_rd(&G, h);
 
937
    for(row=0; row < G.params->rowtot[h]; row++) {
 
938
      i = G.params->roworb[h][row][0];
 
939
      j = G.params->roworb[h][row][1];
 
940
          
 
941
      for(col=0; col < G.params->coltot[h]; col++) {
 
942
        a = G.params->colorb[h][col][0];
 
943
        b = G.params->colorb[h][col][1];
 
944
 
 
945
        value = 0.0;
 
946
 
 
947
        I = T1.params->rowidx[i]; Isym = T1.params->psym[i];
 
948
        J = ZZ.params->rowidx[j]; Jsym = ZZ.params->psym[j];
 
949
        A = T1.params->colidx[a]; Asym = T1.params->qsym[a];
 
950
        B = ZZ.params->colidx[b]; Bsym = ZZ.params->qsym[b];
 
951
 
 
952
        if((Isym==Asym) && (Jsym==Bsym))
 
953
          value += T1.matrix[Isym][I][A] * ZZ.matrix[Jsym][J][B];
 
954
 
 
955
        G.matrix[h][row][col] += 3.0 * value;
 
956
              
 
957
      }
 
958
    }
 
959
    dpd_buf4_mat_irrep_wrt(&G, h);
 
960
    dpd_buf4_mat_irrep_close(&G, h);
 
961
  }
 
962
  dpd_buf4_close(&G);
 
963
 
 
964
  dpd_file2_mat_close(&T1);
 
965
  dpd_file2_close(&T1);
 
966
  dpd_file2_mat_close(&ZZ);
 
967
  dpd_file2_close(&ZZ);
 
968
 
 
969
  /* T(I,E) L(M,E) --> ZZ(I,M) */
 
970
  dpd_file2_init(&ZZ, CC_TMP8, 0, 0, 0, "Z(I,M)");
 
971
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
972
  dpd_file2_init(&L1, CC_OEI, 0, 0, 1, "LIA");
 
973
  dpd_contract222(&T1, &L1, &ZZ, 0, 0, 1.0, 0.0);
 
974
  dpd_file2_close(&L1);
 
975
  dpd_file2_close(&T1);
 
976
 
 
977
  /* ZZ(I,M) T(M,A) --> ZZ2(I,A) */
 
978
  dpd_file2_init(&ZZ2, CC_TMP9, 0, 0, 1, "ZZ2(I,A)");
 
979
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
980
  dpd_contract222(&ZZ, &T1, &ZZ2, 0, 1, 1.0, 0.0);
 
981
  dpd_file2_close(&T1);
 
982
  dpd_file2_close(&ZZ);
 
983
  dpd_file2_close(&ZZ2);
 
984
 
 
985
  /* 3 T(j,b) ZZ(I,A) --> G(Ij,Ab) */
 
986
  dpd_file2_init(&ZZ, CC_TMP9, 0, 0, 1, "ZZ2(I,A)");
 
987
  dpd_file2_mat_init(&ZZ);
 
988
  dpd_file2_mat_rd(&ZZ);
 
989
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
 
990
  dpd_file2_mat_init(&T1);
 
991
  dpd_file2_mat_rd(&T1);
 
992
 
 
993
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 5, 0, 5, 0, "GIjAb");
 
994
  for(h=0; h < nirreps; h++) {
 
995
    dpd_buf4_mat_irrep_init(&G, h); 0,
 
996
                                      dpd_buf4_mat_irrep_rd(&G, h);
 
997
    for(row=0; row < G.params->rowtot[h]; row++) {
 
998
      i = G.params->roworb[h][row][0];
 
999
      j = G.params->roworb[h][row][1];
 
1000
          
 
1001
      for(col=0; col < G.params->coltot[h]; col++) {
 
1002
        a = G.params->colorb[h][col][0];
 
1003
        b = G.params->colorb[h][col][1];
 
1004
 
 
1005
        value = 0.0;
 
1006
 
 
1007
        I = ZZ.params->rowidx[i]; Isym = ZZ.params->psym[i];
 
1008
        J = T1.params->rowidx[j]; Jsym = T1.params->psym[j];
 
1009
        A = ZZ.params->colidx[a]; Asym = ZZ.params->qsym[a];
 
1010
        B = T1.params->colidx[b]; Bsym = T1.params->qsym[b];
 
1011
 
 
1012
        if((Isym==Asym) && (Jsym==Bsym))
 
1013
          value += ZZ.matrix[Isym][I][A] * T1.matrix[Jsym][J][B];
 
1014
 
 
1015
        G.matrix[h][row][col] += 3.0 * value;
 
1016
              
 
1017
      }
 
1018
    }
 
1019
    dpd_buf4_mat_irrep_wrt(&G, h);
 
1020
    dpd_buf4_mat_irrep_close(&G, h);
 
1021
  }
 
1022
  dpd_buf4_scm(&G, 0.5);
 
1023
  dpd_buf4_close(&G);
 
1024
 
 
1025
  dpd_file2_mat_close(&T1);
 
1026
  dpd_file2_close(&T1);
 
1027
  dpd_file2_mat_close(&ZZ);
 
1028
  dpd_file2_close(&ZZ);
 
1029
}