1
#include <libdpd/dpd.h>
5
/* cc3_Wmbij(): Compute the Wmbij components of the
6
** T1-similarity-transformed Hamiltonian matrix, which is given in
9
** Wmbij = <mb||ij> + P(ij) t_i^e <mb||ej> - t_n^b Wmnij + t_i^e t_j^f <mb||ef>
11
** where the Wmnij intermediate is described in cc3_Wmnij.c.
16
void purge_Wmbij(void);
20
dpdbuf4 C, D, E, F, W, W1, Z, X, Z1;
21
dpdfile2 t1, tia, tIA;
23
if(params.ref == 0) { /** RHF **/
25
/* W(Mb,Ij) <-- <Mb|Ij> */
26
dpd_buf4_init(&E, CC_EINTS, 0, 0, 10, 0, 10, 0,"E <ij|ka>");
27
dpd_buf4_sort(&E, CC3_HET1, rspq, 10, 0, "CC3 WMbIj (Mb,Ij)");
30
dpd_file2_init(&t1, CC_OEI, 0, 0, 1, "tIA");
32
/* W(Mb,Ij) <-- - t(n,b) W(Mn,Ij) */
33
dpd_buf4_init(&W, CC3_HET1, 0, 10, 0, 10, 0, 0, "CC3 WMbIj (Mb,Ij)");
34
dpd_buf4_init(&W1, CC3_HET1, 0, 0, 0, 0, 0, 0, "CC3 WMnIj (Mn,Ij)");
35
dpd_contract424(&W1, &t1, &W, 1, 0, 1, -1, 1);
40
/* W(Mb,Ij) <-- + t(j,e) <Mb|Ie) */
41
dpd_buf4_init(&W, CC3_HET1, 0, 10, 0, 10, 0, 0, "CC3 WMbIj (Mb,Ij)");
42
dpd_buf4_init(&C, CC_CINTS, 0, 10, 10, 10, 10, 0, "C <ia|jb>");
43
dpd_contract424(&C, &t1, &W, 3, 1, 0, 1, 1);
47
/* Z(Mb,Ej) = <Mb|Ej> + t(j,f) * <Mb|Ef> */
48
dpd_buf4_init(&D, CC_DINTS, 0, 10, 11, 10, 11, 0, "D <ij|ab> (ib,aj)");
49
dpd_buf4_copy(&D, CC_TMP0, "CC3 ZMbEj (Mb,Ej)");
51
dpd_buf4_init(&Z, CC_TMP0, 0, 10, 11, 10, 11, 0, "CC3 ZMbEj (Mb,Ej)");
52
dpd_buf4_init(&F, CC_FINTS, 0, 10, 5, 10, 5, 0, "F <ia|bc>");
53
dpd_contract424(&F, &t1, &Z, 3, 1, 0, 1, 1);
57
/* W(Mb,Ij) <-- t(I,E) * Z(Mb,Ej) */
58
dpd_buf4_init(&W, CC3_HET1, 0, 10, 0, 10, 0, 0, "CC3 WMbIj (Mb,Ij)");
59
dpd_buf4_init(&Z, CC_TMP0, 0, 10, 11, 10, 11, 0, "CC3 ZMbEj (Mb,Ej)");
60
dpd_contract244(&t1, &Z, &W, 1, 2, 1, 1, 1);
62
dpd_buf4_sort(&W, CC3_HET1, rspq, 0, 10, "CC3 WMbIj (Ij,Mb)");
69
else if (params.ref == 1) {
70
/** W(MB,I>J) <--- <MB||IJ> **/
71
/** W(mb,i>j) <--- <mb||ij> **/
72
dpd_buf4_init(&E, CC_EINTS, 0, 2, 10, 2, 10, 0, "E <ij||ka> (i>j,ka)");
73
dpd_buf4_sort(&E, CC3_HET1, rspq, 10, 2, "CC3 WMBIJ (MB,I>J)");
74
dpd_buf4_sort(&E, CC3_HET1, rspq, 10, 2, "CC3 Wmbij (mb,i>j)");
77
/** W(Mb,Ij) <--- <Mb|Ij> **/
78
/** W(mB,iJ) <--- <mB|iJ> **/
79
dpd_buf4_init(&E, CC_EINTS, 0, 0, 10, 0, 10, 0, "E <ij|ka>");
80
dpd_buf4_sort(&E, CC3_HET1, rspq, 10, 0, "CC3 WMbIj (Mb,Ij)");
81
dpd_buf4_sort(&E, CC3_HET1, rspq, 10, 0, "CC3 WmBiJ (mB,iJ)");
86
dpd_file2_init(&tIA, CC_OEI, 0, 0, 1, "tIA");
87
dpd_file2_init(&tia, CC_OEI, 0, 0, 1, "tia");
89
/**** W(MB,I>J) <-- -ZMBJI <-- P(I/J)( -<JE||MB> * t1[I][E] ) ****/
90
dpd_buf4_init(&Z, CC_TMP0, 0, 10, 0, 10, 0, 0, "Z (MB,JI)");
91
dpd_buf4_init(&C, CC_CINTS, 0, 10, 10, 10, 10, 0, "C <ia||jb>");
92
dpd_contract424(&C, &tIA, &Z, 1, 1, 0, -1, 0);
95
dpd_buf4_sort(&Z, CC_TMP0, pqsr, 10, 0, "X (MB,IJ)");
96
dpd_buf4_init(&X, CC_TMP0, 0, 10, 0, 10, 0, 0, "X (MB,IJ)");
97
dpd_buf4_axpy(&Z, &X, -1);
99
dpd_buf4_init(&W, CC3_HET1, 0, 10, 0, 10, 2, 0, "CC3 WMBIJ (MB,I>J)");
100
dpd_buf4_axpy(&X, &W, 1);
104
/**** W(mb,i>j) <-- -Zmbji <-- P(i/j)( -<je||mb> * t1[i][e] ) ****/
105
dpd_buf4_init(&Z, CC_TMP0, 0, 10, 0, 10, 0, 0, "Z (mb,ji)");
106
dpd_buf4_init(&C, CC_CINTS, 0, 10, 10, 10, 10, 0, "C <ia||jb>");
107
dpd_contract424(&C, &tia, &Z, 1, 1, 0, -1, 0);
110
dpd_buf4_sort(&Z, CC_TMP0, pqsr, 10, 0, "X (mb,ij)");
111
dpd_buf4_init(&X, CC_TMP0, 0, 10, 0, 10, 0, 0, "X (mb,ij)");
112
dpd_buf4_axpy(&Z, &X, -1);
114
dpd_buf4_init(&W, CC3_HET1, 0, 10, 0, 10, 2, 0, "CC3 Wmbij (mb,i>j)");
115
dpd_buf4_axpy(&X, &W, 1);
119
/**** W(Mb,Ij) <-- ( <Mb|Ej> * t1[I][E] ) ****/
121
dpd_buf4_init(&W, CC3_HET1, 0, 10, 0, 10, 0, 0, "CC3 WMbIj (Mb,Ij)");
122
dpd_buf4_init(&D, CC_DINTS, 0, 10, 11, 10, 11, 0, "D <ij|ab> (ib,aj)");
123
dpd_contract244(&tIA, &D, &W, 1, 2, 1, 1, 1);
127
/**** W(Mb,Ij) <-- ( <Mb|Ie> * t1[j][e] ) ****/
129
dpd_buf4_init(&W, CC3_HET1, 0, 10, 0, 10, 0, 0, "CC3 WMbIj (Mb,Ij)");
130
dpd_buf4_init(&C, CC_CINTS, 0, 10, 10, 10, 10, 0, "C <ia|jb>");
131
dpd_contract424(&C, &tia, &W, 3, 1, 0, 1, 1);
135
/**** W(mB,iJ) <-- ( <mB|eJ> * t1[i][e] ) ****/
137
dpd_buf4_init(&W, CC3_HET1, 0, 10, 0, 10, 0, 0, "CC3 WmBiJ (mB,iJ)");
138
dpd_buf4_init(&D, CC_DINTS, 0, 10, 11, 10, 11, 0, "D <ij|ab> (ib,aj)");
139
dpd_contract244(&tia, &D, &W, 1, 2, 1, 1.0, 1.0);
143
/**** W(mB,iJ) <-- ( <mB|iE> * t1[J][E] ) ****/
145
dpd_buf4_init(&W, CC3_HET1, 0, 10, 0, 10, 0, 0, "CC3 WmBiJ (mB,iJ)");
146
dpd_buf4_init(&C, CC_CINTS, 0, 10, 10, 10, 10, 0, "C <ia|jb>");
147
dpd_contract424(&C, &tIA, &W, 3, 1, 0, 1.0, 1);
153
/**** W(MB,I>J) <-- ( t1[N][B] * W(MN,I>J) ****/
155
dpd_buf4_init(&W, CC3_HET1, 0, 10, 2, 10, 2, 0, "CC3 WMBIJ (MB,I>J)");
156
dpd_buf4_init(&Z, CC3_HET1, 0, 0, 2, 2, 2, 0, "CC3 WMNIJ (M>N,I>J)");
157
dpd_contract424(&Z, &tIA, &W, 1, 0, 1, -1, 1);
161
/**** W(mb,i>j) <-- ( t1[n][b] * W(mn,i>j) ****/
163
dpd_buf4_init(&W, CC3_HET1, 0, 10, 2, 10, 2, 0, "CC3 Wmbij (mb,i>j)");
164
dpd_buf4_init(&Z, CC3_HET1, 0, 0, 2, 2, 2, 0, "CC3 Wmnij (m>n,i>j)");
165
dpd_contract424(&Z, &tia, &W, 1, 0, 1, -1, 1);
169
/**** W(Mb,Ij) <-- ( t1[n][b] * W(Mn,Ij) ****/
171
dpd_buf4_init(&W, CC3_HET1, 0, 10, 0, 10, 0, 0, "CC3 WMbIj (Mb,Ij)");
172
dpd_buf4_init(&Z, CC3_HET1, 0, 0, 0, 0, 0, 0, "CC3 WMnIj (Mn,Ij)");
173
dpd_contract424(&Z, &tia, &W, 1, 0, 1, -1, 1);
177
/**** W(mB,iJ) <-- ( t1[N][B] * W(mN,iJ) ****/
179
dpd_buf4_init(&Z1, CC_TMP0, 0, 11, 0, 11, 0, 0, "Z (Bm,Ji)");
180
dpd_buf4_init(&Z, CC3_HET1, 0, 0, 0, 0, 0, 0, "CC3 WMnIj (Mn,Ij)");
181
dpd_contract244(&tIA, &Z, &Z1, 0, 0, 0, -1, 0);
183
dpd_buf4_sort_axpy(&Z1, CC3_HET1, qpsr, 10, 0, "CC3 WmBiJ (mB,iJ)", 1.0);
188
/**** W(MB,I>J) <-- 0.5*P(I/J)XMBIJ <--- ( t1[I][E] * ZMBEJ ) <-- <MB||EF> * t1[J][F] ****/
190
dpd_buf4_init(&Z, CC_TMP0, 0, 10, 11, 10, 11, 0, "Z (MB,EJ)");
191
dpd_buf4_init(&F, CC_FINTS, 0, 10, 5, 10, 5, 1, "F <ia|bc>");
192
dpd_contract424(&F, &tIA, &Z, 3, 1, 0, 1, 0);
195
dpd_buf4_init(&W, CC3_HET1, 0, 10, 0, 10, 2, 0, "CC3 WMBIJ (MB,I>J)");
196
dpd_contract244(&tIA, &Z, &W, 1, 2, 1, 1.0, 1.0);
200
/**** W(mb,i>j) <-- P(i/j) (Zmbif * t1[j][f]) <-- 0.5*( t1[i][e] * <mb||ef> ) ****/
202
dpd_buf4_init(&Z, CC_TMP0, 0, 10, 11, 10, 11, 0, "Z (mb,ej)");
203
dpd_buf4_init(&F, CC_FINTS, 0, 10, 5, 10, 5, 1, "F <ia|bc>");
204
dpd_contract424(&F, &tia, &Z, 3, 1, 0, 1, 0);
207
dpd_buf4_init(&W, CC3_HET1, 0, 10, 0, 10, 2, 0, "CC3 Wmbij (mb,i>j)");
208
dpd_contract244(&tia, &Z, &W, 1, 2, 1, 1.0, 1.0);
212
/**** W(Mb,Ij) <-- (ZIfMb * t1[j][f]) <-- t1[I][E] * <Mb|Ef> ****/
214
dpd_buf4_init(&Z, CC_TMP0, 0, 10, 10, 10, 10, 0, "Z (If,Mb)");
215
dpd_buf4_init(&F, CC_FINTS, 0, 10, 5, 10, 5, 0, "F <ia|bc>");
216
dpd_contract244(&tIA, &F, &Z, 1, 2, 0, 1, 0);
219
dpd_buf4_init(&W, CC3_HET1, 0, 10, 0, 10, 0, 0, "CC3 WMbIj (Mb,Ij)");
220
dpd_contract424(&Z, &tia, &W, 1, 1, 0, 1, 1);
224
/**** W(mB,iJ) <-- ZiFmB * t1[J][F] <-- t1[i][e] * <mB|eF> ****/
226
dpd_buf4_init(&Z, CC_TMP0, 0, 10, 10, 10, 10, 0, "Z (iF,mB)");
227
dpd_buf4_init(&F, CC_FINTS, 0, 10, 5, 10, 5, 0, "F <ia|bc>");
228
dpd_contract244(&tia, &F, &Z, 1, 2, 0, 1, 0);
231
dpd_buf4_init(&W, CC3_HET1, 0, 10, 0, 10, 0, 0, "CC3 WmBiJ (mB,iJ)");
232
dpd_contract424(&Z, &tIA, &W, 1, 1, 0, 1, 1);
236
/* do purge before sort */
239
/* do final sort to get (Ij,Mb) */
240
dpd_buf4_init(&W, CC3_HET1, 0, 10, 2, 10, 2, 0, "CC3 WMBIJ (MB,I>J)");
241
dpd_buf4_sort(&W, CC3_HET1, rspq, 2, 10, "CC3 WMBIJ (I>J,MB)");
243
dpd_buf4_init(&W, CC3_HET1, 0, 10, 2, 10, 2, 0, "CC3 Wmbij (mb,i>j)");
244
dpd_buf4_sort(&W, CC3_HET1, rspq, 2, 10, "CC3 Wmbij (i>j,mb)");
246
dpd_buf4_init(&W, CC3_HET1, 0, 10, 0, 10, 0, 0, "CC3 WMbIj (Mb,Ij)");
247
dpd_buf4_sort(&W, CC3_HET1, rspq, 0, 10, "CC3 WMbIj (Ij,Mb)");
249
dpd_buf4_init(&W, CC3_HET1, 0, 10, 0, 10, 0, 0, "CC3 WmBiJ (mB,iJ)");
250
dpd_buf4_sort(&W, CC3_HET1, rspq, 0, 10, "CC3 WmBiJ (iJ,mB)");
253
dpd_file2_close(&tIA);
254
dpd_file2_close(&tia);
257
else if (params.ref == 2) {
259
/** W(MB,I>J) <--- <MB||IJ> **/
260
dpd_buf4_init(&E, CC_EINTS, 0, 2, 20, 2, 20, 0, "E <IJ||KA> (I>J,KA)");
261
dpd_buf4_sort(&E, CC3_HET1, rspq, 20, 2, "CC3 WMBIJ (MB,I>J)");
264
/** W(mb,i>j) <--- <mb||ij> **/
265
dpd_buf4_init(&E, CC_EINTS, 0, 12, 30, 12, 30, 0, "E <ij||ka> (i>j,ka)");
266
dpd_buf4_sort(&E, CC3_HET1, rspq, 30, 12, "CC3 Wmbij (mb,i>j)");
269
/** W(Mb,Ij) <--- <Mb|Ij> **/
270
dpd_buf4_init(&E, CC_EINTS, 0, 22, 24, 22, 24, 0, "E <Ij|Ka>");
271
dpd_buf4_sort(&E, CC3_HET1, rspq, 24, 22, "CC3 WMbIj (Mb,Ij)");
274
/** W(mB,iJ) <--- <mB|iJ> **/
275
dpd_buf4_init(&E, CC_EINTS, 0, 23, 27, 23, 27, 0, "E <iJ|kA>");
276
dpd_buf4_sort(&E, CC3_HET1, rspq, 27, 23, "CC3 WmBiJ (mB,iJ)");
281
dpd_file2_init(&tIA, CC_OEI, 0, 0, 1, "tIA");
282
dpd_file2_init(&tia, CC_OEI, 0, 2, 3, "tia");
284
/**** W(MB,I>J) <-- -ZMBJI <-- P(I/J)( -<JE||MB> * t1[I][E] ) ****/
285
dpd_buf4_init(&Z, CC_TMP0, 0, 20, 0, 20, 0, 0, "Z (MB,JI)");
286
dpd_buf4_init(&C, CC_CINTS, 0, 20, 20, 20, 20, 0, "C <IA||JB>");
287
dpd_contract424(&C, &tIA, &Z, 1, 1, 0, -1, 0);
290
dpd_buf4_sort(&Z, CC_TMP0, pqsr, 20, 0, "X (MB,IJ)");
291
dpd_buf4_init(&X, CC_TMP0, 0, 20, 0, 20, 0, 0, "X (MB,IJ)");
292
dpd_buf4_axpy(&Z, &X, -1);
294
dpd_buf4_init(&W, CC3_HET1, 0, 20, 0, 20, 2, 0, "CC3 WMBIJ (MB,I>J)");
295
dpd_buf4_axpy(&X, &W, 1);
299
/**** W(mb,i>j) <-- -Zmbji <-- P(i/j)( -<je||mb> * t1[i][e] ) ****/
300
dpd_buf4_init(&Z, CC_TMP0, 0, 30, 10, 30, 10, 0, "Z (mb,ji)");
301
dpd_buf4_init(&C, CC_CINTS, 0, 30, 30, 30, 30, 0, "C <ia||jb>");
302
dpd_contract424(&C, &tia, &Z, 1, 1, 0, -1, 0);
305
dpd_buf4_sort(&Z, CC_TMP0, pqsr, 30, 10, "X (mb,ij)");
306
dpd_buf4_init(&X, CC_TMP0, 0, 30, 10, 30, 10, 0, "X (mb,ij)");
307
dpd_buf4_axpy(&Z, &X, -1);
309
dpd_buf4_init(&W, CC3_HET1, 0, 30, 10, 30, 12, 0, "CC3 Wmbij (mb,i>j)");
310
dpd_buf4_axpy(&X, &W, 1);
314
/**** W(Mb,Ij) <-- ( <Mb|Ej> * t1[I][E] ) ****/
316
dpd_buf4_init(&W, CC3_HET1, 0, 24, 22, 24, 22, 0, "CC3 WMbIj (Mb,Ij)");
317
dpd_buf4_init(&D, CC_DINTS, 0, 24, 26, 24, 26, 0, "D <Ij|Ab> (Ib,Aj)");
318
dpd_contract244(&tIA, &D, &W, 1, 2, 1, 1, 1);
322
/**** W(Mb,Ij) <-- ( <Mb|Ie> * t1[j][e] ) ****/
324
dpd_buf4_init(&W, CC3_HET1, 0, 24, 22, 24, 22, 0, "CC3 WMbIj (Mb,Ij)");
325
dpd_buf4_init(&C, CC_CINTS, 0, 24, 24, 24, 24, 0, "C <Ia|Jb>");
326
dpd_contract424(&C, &tia, &W, 3, 1, 0, 1, 1);
330
/**** W(mB,iJ) <-- ( <mB|eJ> * t1[i][e] ) ****/
332
dpd_buf4_init(&W, CC3_HET1, 0, 27, 23, 27, 23, 0, "CC3 WmBiJ (mB,iJ)");
333
dpd_buf4_init(&D, CC_DINTS, 0, 27, 25, 27, 25, 0, "D <iJ|aB> (iB,aJ)");
334
dpd_contract244(&tia, &D, &W, 1, 2, 1, 1.0, 1.0);
338
/**** W(mB,iJ) <-- ( <mB|iE> * t1[J][E] ) ****/
340
dpd_buf4_init(&W, CC3_HET1, 0, 27, 23, 27, 23, 0, "CC3 WmBiJ (mB,iJ)");
341
dpd_buf4_init(&C, CC_CINTS, 0, 27, 27, 27, 27, 0, "C <iA|jB>");
342
dpd_contract424(&C, &tIA, &W, 3, 1, 0, 1.0, 1);
348
/**** W(MB,I>J) <-- ( t1[N][B] * W(MN,I>J) ****/
350
dpd_buf4_init(&W, CC3_HET1, 0, 20, 2, 20, 2, 0, "CC3 WMBIJ (MB,I>J)");
351
dpd_buf4_init(&Z, CC3_HET1, 0, 0, 2, 2, 2, 0, "CC3 WMNIJ (M>N,I>J)");
352
dpd_contract424(&Z, &tIA, &W, 1, 0, 1, -1, 1);
356
/**** W(mb,i>j) <-- ( t1[n][b] * W(mn,i>j) ****/
358
dpd_buf4_init(&W, CC3_HET1, 0, 30, 12, 30, 12, 0, "CC3 Wmbij (mb,i>j)");
359
dpd_buf4_init(&Z, CC3_HET1, 0, 10, 12, 12, 12, 0, "CC3 Wmnij (m>n,i>j)");
360
dpd_contract424(&Z, &tia, &W, 1, 0, 1, -1, 1);
364
/**** W(Mb,Ij) <-- ( t1[n][b] * W(Mn,Ij) ****/
366
dpd_buf4_init(&W, CC3_HET1, 0, 24, 22, 24, 22, 0, "CC3 WMbIj (Mb,Ij)");
367
dpd_buf4_init(&Z, CC3_HET1, 0, 22, 22, 22, 22, 0, "CC3 WMnIj (Mn,Ij)");
368
dpd_contract424(&Z, &tia, &W, 1, 0, 1, -1, 1);
372
/**** W(mB,iJ) <-- ( t1[N][B] * W(mN,iJ) ****/
374
dpd_buf4_init(&Z1, CC_TMP0, 0, 26, 22, 26, 22, 0, "Z (Bm,Ji)");
375
dpd_buf4_init(&Z, CC3_HET1, 0, 22, 22, 22, 22, 0, "CC3 WMnIj (Mn,Ij)");
376
dpd_contract244(&tIA, &Z, &Z1, 0, 0, 0, -1, 0);
378
dpd_buf4_sort_axpy(&Z1, CC3_HET1, qpsr, 27, 23, "CC3 WmBiJ (mB,iJ)", 1.0);
383
/**** W(MB,I>J) <-- 0.5*P(I/J)XMBIJ <--- ( t1[I][E] * ZMBEJ ) <-- <MB||EF> * t1[J][F] ****/
385
dpd_buf4_init(&Z, CC_TMP0, 0, 20, 21, 20, 21, 0, "Z (MB,EJ)");
386
dpd_buf4_init(&F, CC_FINTS, 0, 20, 5, 20, 5, 1, "F <IA|BC>");
387
dpd_contract424(&F, &tIA, &Z, 3, 1, 0, 1, 0);
390
dpd_buf4_init(&W, CC3_HET1, 0, 20, 0, 20, 2, 0, "CC3 WMBIJ (MB,I>J)");
391
dpd_contract244(&tIA, &Z, &W, 1, 2, 1, 1.0, 1.0);
395
/**** W(mb,i>j) <-- P(i/j) (Zmbif * t1[j][f]) <-- 0.5*( t1[i][e] * <mb||ef> ) ****/
397
dpd_buf4_init(&Z, CC_TMP0, 0, 30, 31, 30, 31, 0, "Z (mb,ej)");
398
dpd_buf4_init(&F, CC_FINTS, 0, 30, 15, 30, 15, 1, "F <ia|bc>");
399
dpd_contract424(&F, &tia, &Z, 3, 1, 0, 1, 0);
402
dpd_buf4_init(&W, CC3_HET1, 0, 30, 10, 30, 12, 0, "CC3 Wmbij (mb,i>j)");
403
dpd_contract244(&tia, &Z, &W, 1, 2, 1, 1.0, 1.0);
407
/**** W(Mb,Ij) <-- (ZIfMb * t1[j][f]) <-- t1[I][E] * <Mb|Ef> ****/
409
dpd_buf4_init(&Z, CC_TMP0, 0, 24, 24, 24, 24, 0, "Z (If,Mb)");
410
dpd_buf4_init(&F, CC_FINTS, 0, 24, 28, 24, 28, 0, "F <Ia|Bc>");
411
dpd_contract244(&tIA, &F, &Z, 1, 2, 0, 1, 0);
414
dpd_buf4_init(&W, CC3_HET1, 0, 24, 22, 24, 22, 0, "CC3 WMbIj (Mb,Ij)");
415
dpd_contract424(&Z, &tia, &W, 1, 1, 0, 1, 1);
419
/**** W(mB,iJ) <-- ZiFmB * t1[J][F] <-- t1[i][e] * <mB|eF> ****/
421
dpd_buf4_init(&Z, CC_TMP0, 0, 27, 27, 27, 27, 0, "Z (iF,mB)");
422
dpd_buf4_init(&F, CC_FINTS, 0, 27, 29, 27, 29, 0, "F <iA|bC>");
423
dpd_contract244(&tia, &F, &Z, 1, 2, 0, 1, 0);
426
dpd_buf4_init(&W, CC3_HET1, 0, 27, 23, 27, 23, 0, "CC3 WmBiJ (mB,iJ)");
427
dpd_contract424(&Z, &tIA, &W, 1, 1, 0, 1, 1);
431
dpd_file2_close(&tIA);
432
dpd_file2_close(&tia);
434
dpd_buf4_init(&W, CC3_HET1, 0, 20, 2, 20, 2, 0, "CC3 WMBIJ (MB,I>J)");
435
dpd_buf4_sort(&W, CC3_HET1, rspq, 2, 20, "CC3 WMBIJ (I>J,MB)");
437
dpd_buf4_init(&W, CC3_HET1, 0, 30, 12, 30, 12, 0, "CC3 Wmbij (mb,i>j)");
438
dpd_buf4_sort(&W, CC3_HET1, rspq, 12, 30, "CC3 Wmbij (i>j,mb)");
440
dpd_buf4_init(&W, CC3_HET1, 0, 24, 22, 24, 22, 0, "CC3 WMbIj (Mb,Ij)");
441
dpd_buf4_sort(&W, CC3_HET1, rspq, 22, 24, "CC3 WMbIj (Ij,Mb)");
443
dpd_buf4_init(&W, CC3_HET1, 0, 27, 23, 27, 23, 0, "CC3 WmBiJ (mB,iJ)");
444
dpd_buf4_sort(&W, CC3_HET1, rspq, 23, 27, "CC3 WmBiJ (iJ,mB)");
449
void purge_Wmbij(void) {
452
int h, a, b, e, f, i, j, m, n;
453
int A, B, E, F, I, J, M, N;
454
int mn, ei, ma, ef, me, jb, mb, ij, ab;
455
int asym, bsym, esym, fsym, isym, jsym, msym, nsym;
456
int *occ_off, *vir_off;
457
int *occ_sym, *vir_sym;
458
int *openpi, nirreps;
460
nirreps = moinfo.nirreps;
461
occpi = moinfo.occpi; virtpi = moinfo.virtpi;
462
occ_off = moinfo.occ_off; vir_off = moinfo.vir_off;
463
occ_sym = moinfo.occ_sym; vir_sym = moinfo.vir_sym;
464
openpi = moinfo.openpi;
466
dpd_file4_init(&W, CC3_HET1, 0, 10, 2,"CC3 WMBIJ (MB,I>J)");
467
for(h=0; h < nirreps; h++) {
468
dpd_file4_mat_irrep_init(&W, h);
469
dpd_file4_mat_irrep_rd(&W, h);
470
for(mb=0; mb<W.params->rowtot[h]; mb++) {
471
b = W.params->roworb[h][mb][1];
472
bsym = W.params->qsym[b];
473
B = b - vir_off[bsym];
474
for(ij=0; ij<W.params->coltot[h]; ij++) {
475
if (B >= (virtpi[bsym] - openpi[bsym]))
476
W.matrix[h][mb][ij] = 0.0;
479
dpd_file4_mat_irrep_wrt(&W, h);
480
dpd_file4_mat_irrep_close(&W, h);
484
dpd_file4_init(&W, CC3_HET1, 0, 10, 2,"CC3 Wmbij (mb,i>j)");
485
for(h=0; h < nirreps; h++) {
486
dpd_file4_mat_irrep_init(&W, h);
487
dpd_file4_mat_irrep_rd(&W, h);
488
for(mb=0; mb<W.params->rowtot[h]; mb++) {
489
m = W.params->roworb[h][mb][0];
490
msym = W.params->psym[m];
491
M = m - occ_off[msym];
492
for(ij=0; ij<W.params->coltot[h]; ij++) {
493
i = W.params->colorb[h][ij][0];
494
j = W.params->colorb[h][ij][1];
495
isym = W.params->rsym[i];
496
jsym = W.params->ssym[j];
497
I = i - occ_off[isym];
498
J = j - occ_off[jsym];
499
if ((M >= (occpi[msym] - openpi[msym])) ||
500
(I >= (occpi[isym] - openpi[isym])) ||
501
(J >= (occpi[jsym] - openpi[jsym])) )
502
W.matrix[h][mb][ij] = 0.0;
505
dpd_file4_mat_irrep_wrt(&W, h);
506
dpd_file4_mat_irrep_close(&W, h);
510
dpd_file4_init(&W, CC3_HET1, 0, 10, 0,"CC3 WMbIj (Mb,Ij)");
511
for(h=0; h < nirreps; h++) {
512
dpd_file4_mat_irrep_init(&W, h);
513
dpd_file4_mat_irrep_rd(&W, h);
514
for(mb=0; mb<W.params->rowtot[h]; mb++) {
515
for(ij=0; ij<W.params->coltot[h]; ij++) {
516
j = W.params->colorb[h][ij][1];
517
jsym = W.params->ssym[j];
518
J = j - occ_off[jsym];
519
if (J >= (occpi[jsym] - openpi[jsym]))
520
W.matrix[h][mb][ij] = 0.0;
523
dpd_file4_mat_irrep_wrt(&W, h);
524
dpd_file4_mat_irrep_close(&W, h);
528
dpd_file4_init(&W, CC3_HET1, 0, 10, 0,"CC3 WmBiJ (mB,iJ)");
529
for(h=0; h < nirreps; h++) {
530
dpd_file4_mat_irrep_init(&W, h);
531
dpd_file4_mat_irrep_rd(&W, h);
532
for(mb=0; mb<W.params->rowtot[h]; mb++) {
533
m = W.params->roworb[h][mb][0];
534
b = W.params->roworb[h][mb][1];
535
msym = W.params->psym[m];
536
bsym = W.params->qsym[b];
537
M = m - occ_off[msym];
538
B = b - vir_off[bsym];
539
for(ij=0; ij<W.params->coltot[h]; ij++) {
540
i = W.params->colorb[h][ij][0];
541
isym = W.params->rsym[i];
542
I = i - occ_off[isym];
543
if ((M >= (occpi[msym] - openpi[msym])) ||
544
(B >= (virtpi[bsym] - openpi[bsym])) ||
545
(I >= (occpi[isym] - openpi[isym])) )
546
W.matrix[h][mb][ij] = 0.0;
549
dpd_file4_mat_irrep_wrt(&W, h);
550
dpd_file4_mat_irrep_close(&W, h);