← Back to branch summary

~ubuntu-branches/ubuntu/intrepid/psicode/intrepid

« back to all changes in this revision

Viewing changes to src/bin/cceom/cc3_HC1.c

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2008-06-07 16:49:57 UTC
  • mfrom: (2.1.2 hardy)
  • Revision ID: james.westby@ubuntu.com-20080607164957-8pifvb133yjlkagn
Tags: 3.3.0-3
http://bugs.debian.org/464867
* debian/rules (DEB_MAKE_CHECK_TARGET): Do not abort test suite on
  failures.
* debian/rules (DEB_CONFIGURE_EXTRA_FLAGS): Set ${bindir} to /usr/lib/psi.
* debian/rules (install/psi3): Move psi3 file to /usr/bin.
* debian/patches/07_464867_move_executables.dpatch: New patch, add
  /usr/lib/psi to the $PATH, so that the moved executables are found.
  (closes: #464867)
* debian/patches/00list: Adjusted.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/cceom/cc3_HC1.c
 
1
#include <stdio.h>
 
2
#include <stdlib.h>
 
3
#include <libdpd/dpd.h>
 
4
#define EXTERN
 
5
#include "globals.h"
 
6
 
 
7
void HC1_F(int i, int C_irr);
 
8
void HC1_Wamef(int i, int C_irr);
 
9
void HC1_Wmnie(int i, int C_irr);
 
10
void HC1_Wmnij(int i, int C_irr);
 
11
void HC1_Wmbij(int i, int C_irr);
 
12
void HC1_Wmbej(int i, int C_irr);
 
13
void HC1_Wabei(int i, int C_irr); 
 
14
void purge_HC1(int C_irr);
 
15
void purge_Wmnie(int C_irr);
 
16
void purge_Wmbij(int C_irr);
 
17
void purge_Wabei(int C_irr);
 
18
 
 
19
void cc3_HC1 (int i, int C_irr) {
 
20
  HC1_F(i, C_irr);
 
21
  HC1_Wamef(i, C_irr);
 
22
  HC1_Wmnie(i, C_irr);
 
23
/*  HC1_Wmnij(i, C_irr); */
 
24
/*  HC1_Wmbij(i, C_irr); */
 
25
/*  HC1_Wmbej(i, C_irr); */
 
26
  HC1_Wabei(i, C_irr);   
 
27
 
 
28
  if (params.ref == 1) purge_HC1(C_irr);
 
29
 
 
30
  return;
 
31
}
 
32
 
 
33
/* constructs matrix elements of [H, C1] for CC3 EOM code */
 
34
 
 
35
void HC1_F(int i, int C_irr) {
 
36
 
 
37
  dpdfile2 FME, Fme, Cme, CME;
 
38
  dpdbuf4 D_anti, D;
 
39
  char CME_lbl[32], Cme_lbl[32];
 
40
  sprintf(CME_lbl, "%s %d", "CME", i);
 
41
  sprintf(Cme_lbl, "%s %d", "Cme", i);
 
42
        double tval;
 
43
 
 
44
  if(params.ref == 0) {
 
45
    /* HC1_F()  Fme = +C_n^f <mn||ef> */
 
46
 
 
47
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, CME_lbl);
 
48
 
 
49
    dpd_file2_init(&FME, CC3_HC1, C_irr, 0, 1, "HC1 FME");
 
50
    dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 0, "D 2<ij|ab> - <ij|ba>");
 
51
    dpd_dot13(&CME, &D, &FME, 0, 0, 1.0, 0.0);
 
52
    dpd_buf4_close(&D);
 
53
    dpd_file2_close(&FME);
 
54
 
 
55
    dpd_file2_close(&CME);
 
56
  }
 
57
 
 
58
  else if(params.ref == 1) { /** ROHF **/
 
59
    /* HC1_F()  Fme = +C_n^f <mn||ef> */
 
60
 
 
61
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, CME_lbl);
 
62
    dpd_file2_init(&Cme, EOM_Cme, C_irr, 0, 1, Cme_lbl);
 
63
 
 
64
    dpd_file2_init(&FME, CC3_HC1, C_irr, 0, 1, "HC1 FME");
 
65
    dpd_file2_init(&Fme, CC3_HC1, C_irr, 0, 1, "HC1 Fme");
 
66
  
 
67
    dpd_buf4_init(&D_anti, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij||ab>");
 
68
    dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
 
69
    dpd_dot13(&CME, &D_anti, &FME, 0, 0, 1.0, 0.0);
 
70
    dpd_dot13(&Cme, &D,      &FME, 0, 0, 1.0, 1.0);
 
71
    dpd_dot13(&Cme, &D_anti, &Fme, 0, 0, 1.0, 0.0);
 
72
    dpd_dot13(&CME, &D,      &Fme, 0, 0, 1.0, 1.0);
 
73
    dpd_buf4_close(&D_anti);
 
74
    dpd_buf4_close(&D);
 
75
 
 
76
    dpd_file2_close(&FME);
 
77
    dpd_file2_close(&Fme);
 
78
 
 
79
    dpd_file2_close(&CME);
 
80
    dpd_file2_close(&Cme);
 
81
 
 
82
  } /** RHF or ROHF **/
 
83
  else if(params.ref == 2) { /** UHF **/
 
84
    /* HC1_F()  Fme = +C_n^f <mn||ef> */
 
85
 
 
86
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, CME_lbl);
 
87
    dpd_file2_init(&Cme, EOM_Cme, C_irr, 2, 3, Cme_lbl);
 
88
 
 
89
    dpd_file2_init(&FME, CC3_HC1, C_irr, 0, 1, "HC1 FME");
 
90
    dpd_file2_init(&Fme, CC3_HC1, C_irr, 2, 3, "HC1 Fme");
 
91
 
 
92
    dpd_buf4_init(&D, CC_DINTS, 0, 20, 20, 20, 20, 0, "D <IJ||AB> (IA,JB)");
 
93
    dpd_contract422(&D, &CME, &FME, 0, 0, 1.0, 0.0);
 
94
    dpd_buf4_close(&D);
 
95
    dpd_buf4_init(&D, CC_DINTS, 0, 20, 30, 20, 30, 0, "D <Ij|Ab> (IA,jb)");
 
96
    dpd_contract422(&D, &Cme, &FME, 0, 0, 1.0, 1.0);
 
97
    dpd_buf4_close(&D);
 
98
    dpd_buf4_init(&D, CC_DINTS, 0, 30, 30, 30, 30, 0, "D <ij||ab> (ia,jb)");
 
99
    dpd_contract422(&D, &Cme, &Fme, 0, 0, 1.0, 0.0);
 
100
    dpd_buf4_close(&D);
 
101
    dpd_buf4_init(&D, CC_DINTS, 0, 30, 20, 30, 20, 0, "D <Ij|Ab> (ia,JB)");
 
102
    dpd_contract422(&D, &CME, &Fme, 0, 0, 1.0, 1.0);
 
103
    dpd_buf4_close(&D);
 
104
 
 
105
    dpd_file2_close(&FME);
 
106
    dpd_file2_close(&Fme);
 
107
 
 
108
    dpd_file2_close(&CME);
 
109
    dpd_file2_close(&Cme);
 
110
  }
 
111
 
 
112
  return;
 
113
}
 
114
 
 
115
 
 
116
void HC1_Wamef(int i, int C_irr) {
 
117
  dpdbuf4 Wamef, WAMEF, WAmEf, WaMeF, W, D_a, D;
 
118
  dpdfile2 Cme, CME;
 
119
  char CME_lbl[32], Cme_lbl[32];
 
120
  sprintf(CME_lbl, "%s %d", "CME", i);
 
121
  sprintf(Cme_lbl, "%s %d", "Cme", i);
 
122
 
 
123
  if(params.ref == 0) { /** RHF **/
 
124
    /* HC1_Wamef():  Wamef = -Cna <nm||ef> */
 
125
 
 
126
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, CME_lbl);
 
127
 
 
128
    dpd_buf4_init(&W, CC3_HC1, C_irr, 11, 5, 11, 5, 0, "HC1 WAmEf (Am,Ef)");
 
129
    dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
 
130
    dpd_contract244(&CME, &D, &W, 0, 0, 0, -1, 0.0);
 
131
    dpd_buf4_close(&D);
 
132
    dpd_buf4_sort(&W, CC3_HC1, qprs, 10, 5, "HC1 WAmEf (mA,Ef)");
 
133
    dpd_buf4_close(&W);
 
134
 
 
135
    dpd_file2_close(&CME);
 
136
  }
 
137
 
 
138
  else if(params.ref == 1) { /** ROHF **/
 
139
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, CME_lbl);
 
140
    dpd_file2_init(&Cme, EOM_Cme, C_irr, 0, 1, Cme_lbl);
 
141
 
 
142
    /* C(N,A) <NM||EF> --> W(AM,E>F) */
 
143
    dpd_buf4_init(&WAMEF, CC3_HC1, C_irr, 11, 7, 11, 7, 0, "HC1 WAMEF (AM,E>F)");
 
144
    dpd_buf4_init(&D, CC_DINTS, 0, 0, 7, 0, 7, 0, "D <ij||ab> (ij,a>b)");
 
145
    dpd_contract244(&CME, &D, &WAMEF, 0, 0, 0, -1, 0.0);
 
146
    dpd_buf4_close(&D);
 
147
    dpd_buf4_close(&WAMEF);  
 
148
 
 
149
    /* C(n,a) <nm||ef> --> W(am,e>f) */
 
150
    dpd_buf4_init(&Wamef, CC3_HC1, C_irr, 11, 7, 11, 7, 0, "HC1 Wamef (am,e>f)");
 
151
    dpd_buf4_init(&D, CC_DINTS, 0, 0, 7, 0, 7, 0, "D <ij||ab> (ij,a>b)");
 
152
    dpd_contract244(&Cme, &D, &Wamef, 0, 0, 0, -1, 0.0);
 
153
    dpd_buf4_close(&D);
 
154
    dpd_buf4_close(&Wamef); 
 
155
 
 
156
    /* C(N,A) <Nm|Ef> --> W(Am,Ef) */
 
157
    dpd_buf4_init(&WAmEf, CC3_HC1, C_irr, 11, 5, 11, 5, 0, "HC1 WAmEf (Am,Ef)");
 
158
    dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
 
159
    dpd_contract244(&CME, &D, &WAmEf, 0, 0, 0, -1, 0.0);
 
160
    dpd_buf4_close(&D);
 
161
    dpd_buf4_close(&WAmEf);
 
162
 
 
163
    /* C(n,a) <nM|eF> --> W(aM,eF) */
 
164
    dpd_buf4_init(&WaMeF, CC3_HC1, C_irr, 11, 5, 11, 5, 0, "HC1 WaMeF (aM,eF)");
 
165
    dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
 
166
    dpd_contract244(&Cme, &D, &WaMeF, 0, 0, 0, -1, 0.0);
 
167
    dpd_buf4_close(&D);
 
168
    dpd_buf4_close(&WaMeF);
 
169
 
 
170
    dpd_file2_close(&CME);
 
171
    dpd_file2_close(&Cme);
 
172
 
 
173
  } /** ROHF **/
 
174
  else if(params.ref == 2) { /** UHF **/
 
175
 
 
176
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, CME_lbl);
 
177
    dpd_file2_init(&Cme, EOM_Cme, C_irr, 2, 3, Cme_lbl);
 
178
 
 
179
    /* T(N,A) <NM||EF> --> W(AM,E>F) */
 
180
    dpd_buf4_init(&WAMEF, CC3_HC1, C_irr, 21, 7, 21, 7, 0, "HC1 WAMEF (AM,E>F)");
 
181
    dpd_buf4_init(&D, CC_DINTS, 0, 0, 7, 0, 7, 0, "D <IJ||AB> (IJ,A>B)");
 
182
    dpd_contract244(&CME, &D, &WAMEF, 0, 0, 0, -1, 0.0);
 
183
    dpd_buf4_close(&D);
 
184
    dpd_buf4_close(&WAMEF);  
 
185
 
 
186
    /* T(n,a) <nm||ef> --> W(am,e>f) */
 
187
    dpd_buf4_init(&Wamef, CC3_HC1, C_irr, 31, 17, 31, 17, 0, "HC1 Wamef (am,e>f)");
 
188
    dpd_buf4_init(&D, CC_DINTS, 0, 10, 17, 10, 17, 0, "D <ij||ab> (ij,a>b)");
 
189
    dpd_contract244(&Cme, &D, &Wamef, 0, 0, 0, -1, 0.0);
 
190
    dpd_buf4_close(&D);
 
191
    dpd_buf4_close(&Wamef); 
 
192
 
 
193
    /* T(N,A) <Nm|Ef> --> W(Am,Ef) */
 
194
    dpd_buf4_init(&WAmEf, CC3_HC1, C_irr, 26, 28, 26, 28, 0, "HC1 WAmEf (Am,Ef)");
 
195
    dpd_buf4_init(&D, CC_DINTS, 0, 22, 28, 22, 28, 0, "D <Ij|Ab>");
 
196
    dpd_contract244(&CME, &D, &WAmEf, 0, 0, 0, -1, 0.0);
 
197
    dpd_buf4_close(&D);
 
198
    dpd_buf4_close(&WAmEf);
 
199
 
 
200
    /* T(n,a) <nM|eF> --> W(aM,eF) */
 
201
    dpd_buf4_init(&WaMeF, CC3_HC1, C_irr, 25, 29, 25, 29, 0, "HC1 WaMeF (aM,eF)");
 
202
    dpd_buf4_init(&D, CC_DINTS, 0, 23, 29, 23, 29, 0, "D <iJ|aB>");
 
203
    dpd_contract244(&Cme, &D, &WaMeF, 0, 0, 0, -1, 0.0);
 
204
    dpd_buf4_close(&D);
 
205
    dpd_buf4_close(&WaMeF);
 
206
 
 
207
    /* Do final sort */
 
208
    dpd_buf4_init(&WAMEF, CC3_HC1, C_irr, 21, 7, 21, 7, 0, "HC1 WAMEF (AM,E>F)");
 
209
    dpd_buf4_sort(&WAMEF, CC3_HC1, qprs, 20, 7, "HC1 WAMEF (MA,F>E)");
 
210
    dpd_buf4_close(&WAMEF);  
 
211
    dpd_buf4_init(&WAMEF, CC3_HC1, C_irr, 20, 7, 20, 7, 0, "HC1 WAMEF (MA,F>E)");
 
212
    dpd_buf4_scm(&WAMEF, -1.0);
 
213
    dpd_buf4_close(&WAMEF);  
 
214
 
 
215
    dpd_buf4_init(&Wamef, CC3_HC1, C_irr, 31, 17, 31, 17, 0, "HC1 Wamef (am,e>f)");
 
216
    dpd_buf4_sort(&Wamef, CC3_HC1, qprs, 30, 17, "HC1 Wamef (ma,f>e)");
 
217
    dpd_buf4_close(&Wamef); 
 
218
    dpd_buf4_init(&Wamef, CC3_HC1, C_irr, 30, 17, 30, 17, 0, "HC1 Wamef (ma,f>e)");
 
219
    dpd_buf4_scm(&Wamef, -1.0);
 
220
    dpd_buf4_close(&Wamef);  
 
221
 
 
222
    dpd_buf4_init(&WAmEf, CC3_HC1, C_irr, 26, 28, 26, 28, 0, "HC1 WAmEf (Am,Ef)");
 
223
    dpd_buf4_sort(&WAmEf, CC3_HC1, qpsr, 27, 29, "HC1 WAmEf (mA,fE)");
 
224
    dpd_buf4_close(&WAmEf);
 
225
 
 
226
    dpd_buf4_init(&WaMeF, CC3_HC1, C_irr, 25, 29, 25, 29, 0, "HC1 WaMeF (aM,eF)");
 
227
    dpd_buf4_sort(&WaMeF, CC3_HC1, qpsr, 24, 28, "HC1 WaMeF (Ma,Fe)");
 
228
    dpd_buf4_close(&WaMeF);
 
229
 
 
230
    dpd_file2_close(&CME);
 
231
    dpd_file2_close(&Cme);
 
232
  } /** UHF **/
 
233
 
 
234
  return;
 
235
}
 
236
 
 
237
 
 
238
void HC1_Wmnie(int i, int C_irr) {
 
239
  dpdbuf4 W, Wmnie, WMNIE, WMnIe, WmNiE, WMniE, WmNIe;
 
240
  dpdbuf4 D, D_a, Z;
 
241
  dpdfile2 Cme, CME;
 
242
  char CME_lbl[32], Cme_lbl[32];
 
243
  sprintf(CME_lbl, "%s %d", "CME", i);
 
244
  sprintf(Cme_lbl, "%s %d", "Cme", i);
 
245
 
 
246
  if(params.ref == 0) { /** RHF **/
 
247
    /* HC1_Wmnie():  Wmnie = + Cif <mn||fe> */
 
248
 
 
249
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, CME_lbl);
 
250
 
 
251
    /* C(I,F) * D(Mn,Fe) --> W(Mn,Ie) */
 
252
    dpd_buf4_init(&WMnIe, CC3_HC1, C_irr, 0, 10, 0, 10, 0, "HC1 WMnIe (Mn,Ie)");
 
253
    dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
 
254
    dpd_contract244(&CME, &D, &WMnIe, 1, 2, 1, 1, 0.0);
 
255
    dpd_file2_close(&CME);
 
256
    dpd_buf4_close(&D);
 
257
    /* W(Mn,Ie) --> W(Mn,eI) */
 
258
    /* dpd_buf4_sort(&WMnIe, CC3_HC1, pqsr, 0, 11, "HC1 WMnIe (Mn,eI)"); */
 
259
    dpd_buf4_close(&WMnIe);
 
260
  }
 
261
 
 
262
  else if(params.ref == 1) { /** ROHF **/
 
263
    /* HC1_Wmnie():  Wmnie = + Cif <mn||fe> */
 
264
 
 
265
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, CME_lbl);
 
266
    dpd_file2_init(&Cme, EOM_Cme, C_irr, 0, 1, Cme_lbl);
 
267
 
 
268
    /* D(M>N,EF) * C(I,F) --> W(M>N,EI) */
 
269
    dpd_buf4_init(&WMNIE, CC3_HC1, C_irr, 2, 11, 2, 11, 0, "HC1 WMNIE (M>N,EI)");
 
270
    dpd_buf4_init(&D_a, CC_DINTS, 0, 2, 5, 2, 5,0, "D <ij||ab> (i>j,ab)");
 
271
    dpd_contract424(&D_a,&CME,&WMNIE, 3, 1, 0, -1, 0);
 
272
    dpd_buf4_close(&D_a);
 
273
    dpd_buf4_close(&WMNIE);
 
274
 
 
275
    /* D(m>n,ef) * C(i,f) --> W(m>n,ei) */
 
276
    dpd_buf4_init(&Wmnie, CC3_HC1, C_irr, 2, 11, 2, 11, 0, "HC1 Wmnie (m>n,ei)");
 
277
    dpd_buf4_init(&D_a, CC_DINTS, 0, 2, 5, 2, 5, 0, "D <ij||ab> (i>j,ab)");
 
278
    dpd_contract424(&D_a, &Cme, &Wmnie, 3, 1, 0, -1, 0);
 
279
    dpd_buf4_close(&D_a);
 
280
    dpd_buf4_close(&Wmnie);
 
281
 
 
282
    /* D(Mn,Fe) * C(I,F) --> W(Mn,Ie) */
 
283
    dpd_buf4_init(&WMnIe, CC_TMP0, C_irr, 0, 10, 0, 10, 0, "HC1 WMnIe (Mn,Ie)");
 
284
    dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
 
285
    dpd_contract244(&CME, &D, &WMnIe, 1, 2, 1, 1, 0);
 
286
    dpd_buf4_close(&D);
 
287
    /* W(Mn,Ie) --> W(Mn,eI) */
 
288
    dpd_buf4_sort(&WMnIe, CC3_HC1, pqsr, 0, 11, "HC1 WMnIe (Mn,eI)");
 
289
    dpd_buf4_close(&WMnIe);
 
290
 
 
291
    /* D(mN,fE) * C(i,f) --> W(mN,iE) */
 
292
    dpd_buf4_init(&WmNiE, CC_TMP1, C_irr, 0, 10, 0, 10, 0, "HC1 WmNiE (mN,iE)");
 
293
    dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
 
294
    dpd_contract244(&Cme,&D,&WmNiE, 1, 2, 1, 1, 0);
 
295
    dpd_buf4_close(&D);
 
296
    /* W(mN,iE) --> W(mN,Ei) */
 
297
    dpd_buf4_sort(&WmNiE, CC3_HC1, pqsr, 0, 11, "HC1 WmNiE (mN,Ei)");
 
298
    dpd_buf4_close(&WmNiE);
 
299
 
 
300
    dpd_file2_close(&CME);
 
301
    dpd_file2_close(&Cme);
 
302
 
 
303
    purge_Wmnie(C_irr); /* before sorting here */
 
304
 
 
305
    /* also put "normal" sorted versions in CC_HBAR */
 
306
    dpd_buf4_init(&WMNIE, CC3_HC1, C_irr, 2, 11, 2, 11, 0, "HC1 WMNIE (M>N,EI)");
 
307
    dpd_buf4_sort(&WMNIE, CC3_HC1, pqsr, 2, 10, "HC1 WMNIE (M>N,IE)");
 
308
    dpd_buf4_close(&WMNIE);
 
309
    dpd_buf4_init(&Wmnie, CC3_HC1, C_irr, 2, 11, 2, 11, 0, "HC1 Wmnie (m>n,ei)");
 
310
    dpd_buf4_sort(&Wmnie, CC3_HC1, pqsr, 2, 10, "HC1 Wmnie (m>n,ie)");
 
311
    dpd_buf4_close(&Wmnie);
 
312
    dpd_buf4_init(&WMnIe, CC3_HC1, C_irr, 0, 11, 0, 11, 0, "HC1 WMnIe (Mn,eI)");
 
313
    dpd_buf4_sort(&WMnIe, CC3_HC1, pqsr, 0, 10, "HC1 WMnIe (Mn,Ie)");
 
314
    dpd_buf4_close(&WMnIe);
 
315
    dpd_buf4_init(&WmNiE, CC3_HC1, C_irr, 0, 11, 0, 11, 0, "HC1 WmNiE (mN,Ei)");
 
316
    dpd_buf4_sort(&WmNiE, CC3_HC1, pqsr, 0, 10, "HC1 WmNiE (mN,iE)");
 
317
    dpd_buf4_close(&WmNiE);
 
318
  }
 
319
  else if(params.ref == 2) { /** UHF **/
 
320
 
 
321
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, CME_lbl);
 
322
    dpd_file2_init(&Cme, EOM_Cme, C_irr, 2, 3, Cme_lbl);
 
323
 
 
324
    /* <M>N||EF> T(I,F) --> W(M>N,EI) */
 
325
    dpd_buf4_init(&W, CC3_HC1, C_irr, 2, 21, 2, 21, 0, "HC1 WMNIE (M>N,EI)");
 
326
    dpd_buf4_init(&D, CC_DINTS, 0, 2, 5, 2, 5, 0, "D <IJ||AB> (I>J,AB)");
 
327
    dpd_contract424(&D, &CME, &W, 3, 1, 0, -1, 0);
 
328
    dpd_buf4_close(&D);
 
329
    dpd_buf4_close(&W);
 
330
 
 
331
    /* <m>n||ef> T(i,f) --> W(m>n,ei) */
 
332
    dpd_buf4_init(&W, CC3_HC1, C_irr, 12, 31, 12, 31, 0, "HC1 Wmnie (m>n,ei)");
 
333
    dpd_buf4_init(&D, CC_DINTS, 0, 12, 15, 12, 15, 0, "D <ij||ab> (i>j,ab)");
 
334
    dpd_contract424(&D, &Cme, &W, 3, 1, 0, -1, 0);
 
335
    dpd_buf4_close(&D);
 
336
    dpd_buf4_close(&W);
 
337
 
 
338
    /* Z(nM,eI) = <nM|eF> T(I,F) */
 
339
    dpd_buf4_init(&Z, CC_TMP1, C_irr, 23, 25, 23, 25, 0, "Z(nM,eI)");
 
340
    dpd_buf4_init(&D, CC_DINTS, 0, 23, 29, 23, 29, 0, "D <iJ|aB>");
 
341
    dpd_contract424(&D, &CME, &Z, 3, 1, 0, 1, 0);
 
342
    dpd_buf4_close(&D);
 
343
    /* Z(nM,eI) --> W(Mn,eI) */
 
344
    dpd_buf4_sort(&Z, CC3_HC1, qprs, 22, 25, "HC1 WMnIe (Mn,eI)");
 
345
    dpd_buf4_close(&Z);
 
346
 
 
347
    /* Z(Nm,Ei) = <Nm|Ef> T(i,f) */
 
348
    dpd_buf4_init(&Z, CC_TMP1, C_irr, 22, 26, 22, 26, 0, "Z(Nm,Ei)");
 
349
    dpd_buf4_init(&D, CC_DINTS, 0, 22, 28, 22, 28, 0, "D <Ij|Ab>");
 
350
    dpd_contract424(&D, &Cme, &Z, 3, 1, 0, 1, 0);
 
351
    dpd_buf4_close(&D);
 
352
    /* Z(Nm,Ei) --> W(mN,Ei) */
 
353
    dpd_buf4_sort(&Z, CC3_HC1, qprs, 23, 26, "HC1 WmNiE (mN,Ei)");
 
354
    dpd_buf4_close(&Z);
 
355
 
 
356
    dpd_file2_close(&CME);
 
357
    dpd_file2_close(&Cme);
 
358
 
 
359
    /* also put "normal" sorted versions in CC3_HC1 */
 
360
    dpd_buf4_init(&WMNIE, CC3_HC1, C_irr, 2, 21, 2, 21, 0, "HC1 WMNIE (M>N,EI)");
 
361
    dpd_buf4_sort(&WMNIE, CC3_HC1, pqsr, 2, 20, "HC1 WMNIE (M>N,IE)");
 
362
    dpd_buf4_close(&WMNIE);
 
363
    dpd_buf4_init(&Wmnie, CC3_HC1, C_irr, 12, 31, 12, 31, 0, "HC1 Wmnie (m>n,ei)");
 
364
    dpd_buf4_sort(&Wmnie, CC3_HC1, pqsr, 12, 30, "HC1 Wmnie (m>n,ie)");
 
365
    dpd_buf4_close(&Wmnie);
 
366
    dpd_buf4_init(&WMnIe, CC3_HC1, C_irr, 22, 25, 22, 25, 0, "HC1 WMnIe (Mn,eI)");
 
367
    dpd_buf4_sort(&WMnIe, CC3_HC1, pqsr, 22, 24, "HC1 WMnIe (Mn,Ie)");
 
368
    dpd_buf4_close(&WMnIe);
 
369
    dpd_buf4_init(&WmNiE, CC3_HC1, C_irr, 23, 26, 23, 26, 0, "HC1 WmNiE (mN,Ei)");
 
370
    dpd_buf4_sort(&WmNiE, CC3_HC1, pqsr, 23, 27, "HC1 WmNiE (mN,iE)");
 
371
    dpd_buf4_close(&WmNiE);
 
372
  }
 
373
 
 
374
  return;
 
375
}
 
376
 
 
377
void HC1_Wmnij(int i, int C_irr)
 
378
{
 
379
  dpdbuf4 WMNIJ, Wmnij, WMnIj, W; 
 
380
  dpdbuf4 Eijka, Eijka_anti, Eaijk, Eaijk_anti;
 
381
  dpdfile2 Cme, CME;
 
382
  char CME_lbl[32], Cme_lbl[32];
 
383
  sprintf(CME_lbl, "%s %d", "CME", i);
 
384
  sprintf(Cme_lbl, "%s %d", "Cme", i);
 
385
 
 
386
  if(params.ref == 0) { /** RHF **/
 
387
    /** HC1_Wmnij():  Wmnij = + P(ij) Cje <mn||ie> */
 
388
 
 
389
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, CME_lbl);
 
390
 
 
391
    dpd_buf4_init(&WMnIj, CC3_HC1, C_irr, 0, 0, 0, 0, 0, "HC1 WMnIj (Mn,Ij)");
 
392
    dpd_buf4_init(&Eaijk, CC_EINTS, 0, 11, 0, 11, 0, 0, "E <ai|jk>");
 
393
    dpd_contract244(&CME, &Eaijk, &WMnIj, 1, 0, 1, 1, 0.0);
 
394
    dpd_buf4_close(&Eaijk);
 
395
 
 
396
    dpd_buf4_init(&Eijka, CC_EINTS, 0, 0, 10, 0, 10, 0, "E <ij|ka>");
 
397
    dpd_contract424(&Eijka, &CME, &WMnIj, 3, 1, 0, 1, 1.0);
 
398
    dpd_buf4_close(&Eijka);
 
399
    dpd_buf4_close(&WMnIj);
 
400
 
 
401
    dpd_file2_close(&CME);
 
402
  }
 
403
 
 
404
  else if(params.ref == 1) { /** ROHF **/  
 
405
    /** HC1_Wmnij():  Wmnij = + P(ij) Cje <mn||ie> */
 
406
 
 
407
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, CME_lbl);
 
408
    dpd_file2_init(&Cme, EOM_Cme, C_irr, 0, 1, Cme_lbl);
 
409
 
 
410
    dpd_buf4_init(&Eijka_anti, CC_EINTS, 0, 2, 10, 2, 10, 0, "E <ij||ka> (i>j,ka)");
 
411
    dpd_buf4_init(&Eijka, CC_EINTS, 0, 0, 10, 0, 10, 0, "E <ij|ka>");
 
412
    dpd_buf4_init(&Eaijk_anti, CC_EINTS, 0, 11, 2, 11, 0, 1, "E <ai|jk>");
 
413
    dpd_buf4_init(&Eaijk, CC_EINTS, 0, 11, 0, 11, 0, 0, "E <ai|jk>");
 
414
 
 
415
    dpd_buf4_init(&WMNIJ, CC3_HC1, C_irr, 2, 0, 2, 2, 0, "HC1 WMNIJ (M>N,I>J)");
 
416
    dpd_contract424(&Eijka_anti, &CME, &WMNIJ, 3, 1, 0, 1, 0.0);
 
417
    dpd_contract244(&CME, &Eaijk_anti, &WMNIJ, 1, 0, 1, 1, 1);
 
418
    dpd_buf4_close(&WMNIJ);
 
419
 
 
420
    dpd_buf4_init(&Wmnij, CC3_HC1, C_irr, 2, 0, 2, 2, 0, "HC1 Wmnij (m>n,i>j)");
 
421
    dpd_contract424(&Eijka_anti, &Cme, &Wmnij, 3, 1, 0, 1, 0.0);
 
422
    dpd_contract244(&Cme, &Eaijk_anti, &Wmnij, 1, 0, 1, 1, 1);
 
423
    dpd_buf4_close(&Wmnij);
 
424
 
 
425
    dpd_buf4_init(&WMnIj, CC3_HC1, C_irr, 0, 0, 0, 0, 0, "HC1 WMnIj (Mn,Ij)");
 
426
    dpd_contract424(&Eijka, &Cme, &WMnIj, 3, 1, 0, 1, 0.0);
 
427
    dpd_contract244(&CME, &Eaijk, &WMnIj, 1, 0, 1, 1, 1);
 
428
    dpd_buf4_close(&WMnIj);
 
429
 
 
430
    dpd_buf4_close(&Eijka_anti);
 
431
    dpd_buf4_close(&Eijka);
 
432
    dpd_buf4_close(&Eaijk_anti);
 
433
    dpd_buf4_close(&Eaijk);
 
434
 
 
435
    dpd_file2_close(&CME);
 
436
    dpd_file2_close(&Cme);
 
437
  }
 
438
 
 
439
  else if(params.ref == 2) { /*** UHF ***/
 
440
    /** HC1_Wmnij():  Wmnij = + P(ij) Cje <mn||ie> */
 
441
 
 
442
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, CME_lbl);
 
443
    dpd_file2_init(&Cme, EOM_Cme, C_irr, 2, 3, Cme_lbl);
 
444
 
 
445
    dpd_buf4_init(&WMNIJ, CC3_HC1, C_irr, 2, 0, 2, 2, 0, "HC1 WMNIJ (M>N,I>J)");
 
446
    dpd_buf4_init(&Eijka, CC_EINTS, 0, 2, 20, 2, 20, 0, "E <IJ||KA> (I>J,KA)");
 
447
    dpd_buf4_init(&Eaijk, CC_EINTS, 0, 21, 2, 21, 0, 1, "E <AI|JK>");
 
448
    dpd_contract424(&Eijka, &CME, &WMNIJ, 3, 1, 0, 1, 0.0);
 
449
    dpd_contract244(&CME, &Eaijk, &WMNIJ, 1, 0, 1, 1, 1.0);
 
450
    dpd_buf4_close(&Eijka);
 
451
    dpd_buf4_close(&Eaijk);
 
452
    dpd_buf4_close(&WMNIJ);
 
453
 
 
454
    dpd_buf4_init(&Wmnij, CC3_HC1, C_irr, 12, 10, 12, 12, 0, "HC1 Wmnij (m>n,i>j)");
 
455
    dpd_buf4_init(&Eijka, CC_EINTS, 0, 12, 30, 12, 30, 0, "E <ij||ka> (i>j,ka)");
 
456
    dpd_buf4_init(&Eaijk, CC_EINTS, 0, 31, 12, 31, 10, 1, "E <ai|jk>");
 
457
    dpd_contract424(&Eijka, &Cme, &Wmnij, 3, 1, 0, 1, 0.0);
 
458
    dpd_contract244(&Cme, &Eaijk, &Wmnij, 1, 0, 1, 1, 1.0);
 
459
    dpd_buf4_close(&Eijka);
 
460
    dpd_buf4_close(&Eaijk);
 
461
    dpd_buf4_close(&Wmnij);
 
462
 
 
463
    dpd_buf4_init(&WMnIj, CC3_HC1, C_irr, 22, 22, 22, 22, 0, "HC1 WMnIj (Mn,Ij)");
 
464
    dpd_buf4_init(&Eijka, CC_EINTS, 0, 22, 24, 22, 24, 0, "E <Ij|Ka>");
 
465
    dpd_buf4_init(&Eaijk, CC_EINTS, 0, 26, 22, 26, 22, 0, "E <Ai|Jk>");
 
466
    dpd_contract424(&Eijka, &Cme, &WMnIj, 3, 1, 0, 1, 0.0);
 
467
    dpd_contract244(&CME, &Eaijk, &WMnIj, 1, 0, 1, 1, 1.0);
 
468
    dpd_buf4_close(&Eijka);
 
469
    dpd_buf4_close(&Eaijk);
 
470
    dpd_buf4_close(&WMnIj);
 
471
 
 
472
    dpd_file2_close(&CME);
 
473
    dpd_file2_close(&Cme);
 
474
  }
 
475
}
 
476
 
 
477
 
 
478
void HC1_Wmbij(int i, int C_irr)
 
479
{
 
480
  dpdbuf4 W, Wmnij, I, Z, Z1, Z2, C, D;
 
481
  dpdfile2 Cme, CME;
 
482
  char CME_lbl[32], Cme_lbl[32];
 
483
  sprintf(CME_lbl, "%s %d", "CME", i);
 
484
  sprintf(Cme_lbl, "%s %d", "Cme", i);
 
485
 
 
486
  if(params.ref == 0) { /** RHF **/
 
487
    /* HC1_Wmbij = +P(ij) Cie <mb||ej> - Cnb <mn||ij> */
 
488
 
 
489
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, CME_lbl);
 
490
    dpd_buf4_init(&W, CC3_HC1, C_irr, 10, 0, 10, 0, 0, "HC1 WMbIj (Mb,Ij)");
 
491
 
 
492
    /** - C_n^b <Mn|Ij> -> W(Mb,Ij) **/
 
493
    dpd_buf4_init(&I, CC_AINTS, 0, 0, 0, 0, 0, 0, "A <ij|kl>");
 
494
    dpd_contract424(&I, &CME, &W, 1, 0, 1, -1.0, 0.0);
 
495
    dpd_buf4_close(&I);
 
496
 
 
497
    /* <Mb||Ie> Cje -> W(Mb,Ij) */
 
498
    dpd_buf4_init(&C, CC_CINTS, 0, 10, 10, 10, 10, 0, "C <ia|jb>");
 
499
    dpd_contract424(&C, &CME, &W, 3, 1, 0, 1.0, 1.0);
 
500
    dpd_buf4_close(&C);
 
501
 
 
502
    /* CIE <Mb|Ej> -> W(Mb,Ij) */
 
503
    dpd_buf4_init(&D, CC_DINTS, 0, 10, 11, 10, 11, 0, "D <ij|ab> (ib,aj)");
 
504
    dpd_contract244(&CME, &D, &W, 1, 2, 1, 1.0, 1.0);
 
505
    dpd_buf4_close(&D);
 
506
 
 
507
    dpd_buf4_sort(&W, CC3_HC1, rspq, 0, 10, "HC1 WMbIj (Ij,Mb)");
 
508
 
 
509
    dpd_buf4_close(&W);
 
510
    dpd_file2_close(&CME);
 
511
  }
 
512
 
 
513
  else if(params.ref == 1) { /** ROHF **/
 
514
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, CME_lbl);
 
515
    dpd_file2_init(&Cme, EOM_Cme, C_irr, 0, 1, Cme_lbl);
 
516
 
 
517
    /** - C_N^B <MN||IJ> --> W(MB,IJ) **/
 
518
    dpd_buf4_init(&Wmnij, CC_AINTS, 0, 0, 2, 0, 0, 1, "A <ij|kl>");
 
519
    dpd_buf4_init(&W, CC3_HC1, C_irr, 10, 2, 10, 2, 0, "HC1 WMBIJ (MB,I>J)");
 
520
    dpd_contract424(&Wmnij, &CME, &W, 1, 0, 1, -1.0, 0.0);
 
521
    dpd_buf4_close(&W);
 
522
    dpd_buf4_close(&Wmnij);
 
523
 
 
524
    /** - C_n^b <mn||ij> **/
 
525
    dpd_buf4_init(&Wmnij, CC_AINTS, 0, 0, 2, 0, 0, 1, "A <ij|kl>");
 
526
    dpd_buf4_init(&W, CC3_HC1, C_irr, 10, 2, 10, 2, 0, "HC1 Wmbij (mb,i>j)");
 
527
    dpd_contract424(&Wmnij, &Cme, &W, 1, 0, 1, -1.0, 0.0);
 
528
    dpd_buf4_close(&W);
 
529
    dpd_buf4_close(&Wmnij);
 
530
 
 
531
    /** - C_n^b <Mn|Ij> **/
 
532
    dpd_buf4_init(&Wmnij, CC_AINTS, 0, 0, 0, 0, 0, 0, "A <ij|kl>");
 
533
    dpd_buf4_init(&W, CC3_HC1, C_irr, 10, 0, 10, 0, 0, "HC1 WMbIj (Mb,Ij)");
 
534
    dpd_contract424(&Wmnij, &Cme, &W, 1, 0, 1, -1.0, 0.0);
 
535
    dpd_buf4_close(&W);
 
536
    dpd_buf4_close(&Wmnij);
 
537
 
 
538
    /** - C_N^B <mN|iJ> **/
 
539
    dpd_buf4_init(&Wmnij, CC_AINTS, 0, 0, 0, 0, 0, 0, "A <ij|kl>");
 
540
    dpd_buf4_init(&W, CC3_HC1, C_irr, 10, 0, 10, 0, 0, "HC1 WmBiJ (mB,iJ)");
 
541
    dpd_contract424(&Wmnij, &CME, &W, 1, 0, 1, -1.0, 0.0);
 
542
    dpd_buf4_close(&W);
 
543
    dpd_buf4_close(&Wmnij);
 
544
 
 
545
    /* term 2: - P(ij) <mb||ei> Cje -> Wmbij */
 
546
 
 
547
    /** + P(IJ) C_J^E <MB||IE> -> WMBIJ **/
 
548
    dpd_buf4_init(&Z2, CC_TMP0, C_irr, 10, 0, 10, 0, 0, "Z1(MB,IJ)");
 
549
    dpd_buf4_init(&C, CC_CINTS, 0, 10, 10, 10, 10, 0, "C <ia||jb>");
 
550
    dpd_contract424(&C, &CME, &Z2, 3, 1, 0, 1.0, 0.0);
 
551
    dpd_buf4_close(&C);
 
552
 
 
553
    dpd_buf4_sort(&Z2, CC_TMP1, pqsr, 10, 0, "Z2(MB,JI)");
 
554
 
 
555
    dpd_buf4_init(&Z1, CC_TMP1, C_irr, 10, 0, 10, 0, 0, "Z2(MB,JI)");
 
556
    dpd_buf4_axpy(&Z1, &Z2, -1.0);
 
557
    dpd_buf4_close(&Z1);
 
558
    dpd_buf4_init(&W, CC3_HC1, C_irr, 10, 0, 10, 2, 0, "HC1 WMBIJ (MB,I>J)");
 
559
    dpd_buf4_axpy(&Z2, &W, 1.0);
 
560
    dpd_buf4_close(&Z2);
 
561
    dpd_buf4_close(&W);
 
562
 
 
563
    /** - P(ij) C_j^e ( <mb||ie> ) -> WMBIJ **/
 
564
 
 
565
    dpd_buf4_init(&Z2, CC_TMP0, C_irr, 10, 0, 10, 0, 0, "Z1(mb,ij)");
 
566
    dpd_buf4_init(&C, CC_CINTS, 0, 10, 10, 10, 10, 0, "C <ia||jb>");
 
567
    dpd_contract424(&C, &Cme, &Z2, 3, 1, 0, 1.0, 0.0);
 
568
    dpd_buf4_close(&C);
 
569
 
 
570
    dpd_buf4_sort(&Z2, CC_TMP1, pqsr, 10, 0, "Z2(mb,ji)");
 
571
 
 
572
    dpd_buf4_init(&Z1, CC_TMP1, C_irr, 10, 0, 10, 0, 0, "Z2(mb,ji)");
 
573
    dpd_buf4_axpy(&Z1, &Z2, -1.0);
 
574
    dpd_buf4_close(&Z1);
 
575
    dpd_buf4_init(&W, CC3_HC1, C_irr, 10, 0, 10, 2, 0, "HC1 Wmbij (mb,i>j)");
 
576
    dpd_buf4_axpy(&Z2, &W, 1.0);
 
577
    dpd_buf4_close(&Z2);
 
578
    dpd_buf4_close(&W);
 
579
 
 
580
    /* <Mb||Ie> Cje -> W(Mb,Ij) */
 
581
    dpd_buf4_init(&W, CC3_HC1, C_irr, 10, 0, 10, 0, 0, "HC1 WMbIj (Mb,Ij)");
 
582
    dpd_buf4_init(&C, CC_CINTS, 0, 10, 10, 10, 10, 0, "C <ia|jb>");
 
583
    dpd_contract424(&C, &Cme, &W, 3, 1, 0, 1.0, 1.0);
 
584
    dpd_buf4_close(&C); 
 
585
                    
 
586
    /* CIE <Mb|Ej> -> W(Mb,Ij) */
 
587
    dpd_buf4_init(&D, CC_DINTS, 0, 10, 11, 10, 11, 0, "D <ij|ab> (ib,aj)");
 
588
    dpd_contract244(&CME, &D, &W, 1, 2, 1, 1.0, 1.0);
 
589
    dpd_buf4_close(&D);
 
590
    dpd_buf4_close(&W);
 
591
 
 
592
    /** C_J^E <mB||iE> = Cje <mB|iE> **/
 
593
    dpd_buf4_init(&W, CC3_HC1, C_irr, 10, 0, 10, 0, 0, "HC1 WmBiJ (mB,iJ)");
 
594
    dpd_buf4_init(&C, CC_CINTS, 0, 10, 10, 10, 10, 0, "C <ia|jb>");
 
595
    dpd_contract424(&C, &CME, &W, 3, 1, 0, 1.0, 1.0);
 
596
    dpd_buf4_close(&C);
 
597
 
 
598
    /** -C_i^e <mB||Je> = +Cie <mB|eJ> = +Cie <mJ|eB> **/
 
599
    dpd_buf4_init(&D, CC_DINTS, 0, 10, 11, 10, 11, 0, "D <ij|ab> (ib,aj)");
 
600
    dpd_contract244(&Cme, &D, &W, 1, 2, 1, 1.0, 1.0);
 
601
    dpd_buf4_close(&D);
 
602
    dpd_buf4_close(&W);
 
603
 
 
604
    /* do purge before sort */
 
605
    purge_Wmbij(C_irr);
 
606
 
 
607
    /* do final sort to get (Ij,Mb) */
 
608
 
 
609
    dpd_buf4_init(&W, CC3_HC1, C_irr, 10, 2, 10, 2, 0, "HC1 WMBIJ (MB,I>J)");
 
610
    dpd_buf4_sort(&W, CC3_HC1, rspq, 2, 10, "HC1 WMBIJ (I>J,MB)");
 
611
    dpd_buf4_close(&W);
 
612
    dpd_buf4_init(&W, CC3_HC1, C_irr, 10, 2, 10, 2, 0, "HC1 Wmbij (mb,i>j)");
 
613
    dpd_buf4_sort(&W, CC3_HC1, rspq, 2, 10, "HC1 Wmbij (i>j,mb)");
 
614
    dpd_buf4_close(&W);
 
615
    dpd_buf4_init(&W, CC3_HC1, C_irr, 10, 0, 10, 0, 0, "HC1 WMbIj (Mb,Ij)");
 
616
    dpd_buf4_sort(&W, CC3_HC1, rspq, 0, 10, "HC1 WMbIj (Ij,Mb)");
 
617
    dpd_buf4_close(&W);
 
618
    dpd_buf4_init(&W, CC3_HC1, C_irr, 10, 0, 10, 0, 0, "HC1 WmBiJ (mB,iJ)");
 
619
    dpd_buf4_sort(&W, CC3_HC1, rspq, 0, 10, "HC1 WmBiJ (iJ,mB)");
 
620
    dpd_buf4_close(&W);
 
621
 
 
622
    dpd_file2_close(&CME);
 
623
    dpd_file2_close(&Cme);
 
624
  }
 
625
  else if(params.ref == 2) { /** UHF **/
 
626
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, CME_lbl);
 
627
    dpd_file2_init(&Cme, EOM_Cme, C_irr, 2, 3, Cme_lbl);
 
628
 
 
629
    /** - C_N^B W_MNIJ --> W(MB,IJ) **/
 
630
    dpd_buf4_init(&Wmnij, CC_AINTS, 0, 0, 2, 0, 0, 1, "A <IJ|KL>");
 
631
    dpd_buf4_init(&W, CC3_HC1, C_irr, 20, 2, 20, 2, 0, "HC1 WMBIJ (MB,I>J)");
 
632
    dpd_contract424(&Wmnij, &CME, &W, 1, 0, 1, -1.0, 0.0);
 
633
    dpd_buf4_close(&W);
 
634
    dpd_buf4_close(&Wmnij);
 
635
 
 
636
    /** - C_n^b W_mnij **/
 
637
    dpd_buf4_init(&Wmnij, CC_AINTS, 0, 10, 12, 10, 10, 1, "A <ij|kl>");
 
638
    dpd_buf4_init(&W, CC3_HC1, C_irr, 30, 12, 30, 12, 0, "HC1 Wmbij (mb,i>j)");
 
639
    dpd_contract424(&Wmnij, &Cme, &W, 1, 0, 1, -1.0, 0.0);
 
640
    dpd_buf4_close(&W);
 
641
    dpd_buf4_close(&Wmnij);
 
642
 
 
643
    /** - C_n^b W_MnIj **/
 
644
    dpd_buf4_init(&Wmnij, CC_AINTS, 0, 22, 22, 22, 22, 0, "A <Ij|Kl>");
 
645
    dpd_buf4_init(&W, CC3_HC1, C_irr, 24, 22, 24, 22, 0, "HC1 WMbIj (Mb,Ij)");
 
646
    dpd_contract424(&Wmnij, &Cme, &W, 1, 0, 1, -1.0, 0.0);
 
647
    dpd_buf4_close(&W);
 
648
    dpd_buf4_close(&Wmnij);
 
649
 
 
650
    /** - C_N^B W_mNiJ **/
 
651
    dpd_buf4_init(&Wmnij, CC_AINTS, 0, 22, 22, 22, 22, 0, "A <Ij|Kl>");
 
652
    dpd_buf4_sort(&Wmnij, CC_TMP0, qpsr, 23, 23, "A <iJ|kL>");
 
653
    dpd_buf4_close(&Wmnij);
 
654
    dpd_buf4_init(&Wmnij, CC_TMP0, 0, 23, 23, 23, 23, 0, "A <iJ|kL>");
 
655
    dpd_buf4_init(&W, CC3_HC1, C_irr, 27, 23, 27, 23, 0, "HC1 WmBiJ (mB,iJ)");
 
656
    dpd_contract424(&Wmnij, &CME, &W, 1, 0, 1, -1.0, 0.0);
 
657
    dpd_buf4_close(&W);
 
658
    dpd_buf4_close(&Wmnij);
 
659
 
 
660
    /* term 2: - P(ij) <mb||ei> Cje -> Wmbij */
 
661
 
 
662
    /** + P(IJ) C_J^E <MB||IE> -> WMBIJ **/
 
663
    dpd_buf4_init(&Z2, CC_TMP0, C_irr, 20, 0, 20, 0, 0, "Z1(MB,IJ)");
 
664
    dpd_buf4_init(&C, CC_CINTS, 0, 20, 20, 20, 20, 0, "C <IA||JB>");
 
665
    dpd_contract424(&C, &CME, &Z2, 3, 1, 0, 1.0, 0.0);
 
666
    dpd_buf4_close(&C);
 
667
 
 
668
    dpd_buf4_sort(&Z2, CC_TMP1, pqsr, 20, 0, "Z2(MB,JI)");
 
669
 
 
670
    dpd_buf4_init(&Z1, CC_TMP1, C_irr, 20, 0, 20, 0, 0, "Z2(MB,JI)");
 
671
    dpd_buf4_axpy(&Z1, &Z2, -1.0);
 
672
    dpd_buf4_close(&Z1);
 
673
    dpd_buf4_init(&W, CC3_HC1, C_irr, 20, 0, 20, 2, 0, "HC1 WMBIJ (MB,I>J)");
 
674
    dpd_buf4_axpy(&Z2, &W, 1.0);
 
675
    dpd_buf4_close(&Z2);
 
676
    dpd_buf4_close(&W);
 
677
 
 
678
    /** - P(ij) C_j^e ( <mb||ie> ) -> WMBIJ **/
 
679
 
 
680
    dpd_buf4_init(&Z2, CC_TMP0, C_irr, 30, 10, 30, 10, 0, "Z1(mb,ij)");
 
681
    dpd_buf4_init(&C, CC_CINTS, 0, 30, 30, 30, 30, 0, "C <ia||jb>");
 
682
    dpd_contract424(&C, &Cme, &Z2, 3, 1, 0, 1.0, 0.0);
 
683
    dpd_buf4_close(&C);
 
684
 
 
685
    dpd_buf4_sort(&Z2, CC_TMP1, pqsr, 30, 10, "Z2(mb,ji)");
 
686
 
 
687
    dpd_buf4_init(&Z1, CC_TMP1, C_irr, 30, 10, 30, 10, 0, "Z2(mb,ji)");
 
688
    dpd_buf4_axpy(&Z1, &Z2, -1.0);
 
689
    dpd_buf4_close(&Z1);
 
690
    dpd_buf4_init(&W, CC3_HC1, C_irr, 30, 10, 30, 12, 0, "HC1 Wmbij (mb,i>j)");
 
691
    dpd_buf4_axpy(&Z2, &W, 1.0);
 
692
    dpd_buf4_close(&Z2);
 
693
    dpd_buf4_close(&W);
 
694
 
 
695
    /* <Mb||Ie> Cje -> W(Mb,Ij) */
 
696
    dpd_buf4_init(&W, CC3_HC1, C_irr, 24, 22, 24, 22, 0, "HC1 WMbIj (Mb,Ij)");
 
697
    dpd_buf4_init(&C, CC_CINTS, 0, 24, 24, 24, 24, 0, "C <Ia|Jb>");
 
698
    dpd_contract424(&C, &Cme, &W, 3, 1, 0, 1.0, 1.0);
 
699
    dpd_buf4_close(&C); 
 
700
                    
 
701
    /* CIE <Mb|Ej> -> W(Mb,Ij) */
 
702
    dpd_buf4_init(&D, CC_DINTS, 0, 24, 26, 24, 26, 0, "D <Ij|Ab> (Ib,Aj)");
 
703
    dpd_contract244(&CME, &D, &W, 1, 2, 1, 1.0, 1.0);
 
704
    dpd_buf4_close(&D);
 
705
    dpd_buf4_close(&W);
 
706
 
 
707
    /** C_J^E <mB||iE> = C^J_E <mB|iE> **/
 
708
    dpd_buf4_init(&W, CC3_HC1, C_irr, 27, 23, 27, 23, 0, "HC1 WmBiJ (mB,iJ)");
 
709
    dpd_buf4_init(&C, CC_CINTS, 0, 27, 27, 27, 27, 0, "C <iA|jB>");
 
710
    dpd_contract424(&C, &CME, &W, 3, 1, 0, 1.0, 1.0);
 
711
    dpd_buf4_close(&C);
 
712
 
 
713
    /** -C_i^e <mB||Je> = +Cie <mB|eJ> = +Cie <mJ|eB> **/
 
714
    dpd_buf4_init(&D, CC_DINTS, 0, 27, 25, 27, 25, 0, "D <iJ|aB> (iB,aJ)");
 
715
    dpd_contract244(&Cme, &D, &W, 1, 2, 1, 1.0, 1.0);
 
716
    dpd_buf4_close(&D);
 
717
    dpd_buf4_close(&W);
 
718
 
 
719
    /** do final sort to (Ij,Mb) **/
 
720
    dpd_buf4_init(&W, CC3_HC1, C_irr, 20, 2, 20, 2, 0, "HC1 WMBIJ (MB,I>J)");
 
721
    dpd_buf4_sort(&W, CC3_HC1, rspq, 2, 20, "HC1 WMBIJ (I>J,MB)");
 
722
    dpd_buf4_close(&W);
 
723
 
 
724
    dpd_buf4_init(&W, CC3_HC1, C_irr, 30, 12, 30, 12, 0, "HC1 Wmbij (mb,i>j)");
 
725
    dpd_buf4_sort(&W, CC3_HC1, rspq, 12, 30, "HC1 Wmbij (i>j,mb)");
 
726
    dpd_buf4_close(&W);
 
727
 
 
728
    dpd_buf4_init(&W, CC3_HC1, C_irr, 24, 22, 24, 22, 0, "HC1 WMbIj (Mb,Ij)");
 
729
    dpd_buf4_sort(&W, CC3_HC1, rspq, 22, 24, "HC1 WMbIj (Ij,Mb)");
 
730
    dpd_buf4_close(&W);
 
731
 
 
732
    dpd_buf4_init(&W, CC3_HC1, C_irr, 27, 23, 27, 23, 0, "HC1 WmBiJ (mB,iJ)");
 
733
    dpd_buf4_sort(&W, CC3_HC1, rspq, 23, 27, "HC1 WmBiJ (iJ,mB)");
 
734
    dpd_buf4_close(&W);
 
735
 
 
736
    dpd_file2_close(&CME);
 
737
    dpd_file2_close(&Cme);
 
738
  }
 
739
 
 
740
}
 
741
 
 
742
 
 
743
void HC1_Wmbej(int i, int C_irr) {
 
744
  dpdbuf4 WMBEJ, Wmbej, WMbEj, WmBeJ, WmBEj, WMbeJ;
 
745
  dpdbuf4 D, C, F, E, X, Y, W, Z;
 
746
  dpdfile2 Cme, CME;
 
747
  char CME_lbl[32], Cme_lbl[32];
 
748
  sprintf(CME_lbl, "%s %d", "CME", i);
 
749
  sprintf(Cme_lbl, "%s %d", "Cme", i);
 
750
 
 
751
  if(params.ref == 0) { /** RHF **/
 
752
    /* Cjf <mb||ef> -> Wmbej */
 
753
 
 
754
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, CME_lbl);
 
755
 
 
756
    dpd_buf4_init(&WMbEj, CC_TMP0, C_irr, 10, 11, 10, 11, 0, "WMbEj");
 
757
    dpd_buf4_init(&F, CC_FINTS, 0, 10, 5, 10, 5, 0, "F <ia|bc>");
 
758
    dpd_contract424(&F, &CME, &WMbEj, 3, 1, 0, 1, 0.0);
 
759
    dpd_buf4_close(&F);
 
760
    dpd_buf4_close(&WMbEj);
 
761
 
 
762
    
 
763
    dpd_buf4_init(&Z, CC_TMP0, C_irr, 11, 11, 11, 11, 0, "Z(bM,eJ)");
 
764
    dpd_buf4_init(&F, CC_FINTS, 0, 11, 5, 11, 5, 0, "F <ai|bc>");
 
765
    dpd_contract424(&F, &CME, &Z, 3, 1, 0, -1, 0);
 
766
    dpd_buf4_close(&F);
 
767
 
 
768
    dpd_buf4_sort(&Z, CC_TMP0, qpsr, 10, 10, "WMbeJ"); /* (Mb,Je) */
 
769
    dpd_buf4_close(&Z);
 
770
 
 
771
    /* - Cnb <mn||ej> -> Wmbej */
 
772
 
 
773
    dpd_buf4_init(&E, CC_EINTS, 0, 11, 0, 11, 0, 0, "E <ai|jk>");
 
774
    dpd_buf4_init(&WMbEj, CC_TMP0, C_irr, 10, 11, 10, 11, 0, "WMbEj");
 
775
    dpd_contract424(&E, &CME, &WMbEj, 3, 0, 1, -1, 1.0);
 
776
    dpd_buf4_close(&WMbEj);
 
777
    dpd_buf4_close(&E);
 
778
 
 
779
    dpd_buf4_init(&E, CC_EINTS, 0, 0, 10, 0, 10, 0, "E <ij|ka>");
 
780
    dpd_buf4_init(&WMbeJ, CC_TMP0, C_irr, 10, 10, 10, 10, 0, "WMbeJ"); 
 
781
    dpd_contract424(&E, &CME, &WMbeJ, 1, 0, 1, 1, 1.0);
 
782
    dpd_buf4_close(&WMbeJ);
 
783
    dpd_buf4_close(&E);
 
784
 
 
785
    dpd_file2_close(&CME);
 
786
 
 
787
    /* Sort to (ME,JB) */
 
788
 
 
789
    dpd_buf4_init(&WMbEj, CC_TMP0, C_irr, 10, 11, 10, 11, 0, "WMbEj");
 
790
    dpd_buf4_sort(&WMbEj, CC3_HC1, prsq, 10, 10, "HC1 WMbEj (ME,jb)");
 
791
    dpd_buf4_close(&WMbEj);
 
792
 
 
793
    dpd_buf4_init(&WMbeJ, CC_TMP0, C_irr, 10, 10, 10, 10, 0, "WMbeJ");
 
794
    dpd_buf4_sort(&WMbeJ, CC3_HC1, psrq, 10, 10, "HC1 WMbeJ (Me,Jb)");
 
795
    dpd_buf4_close(&WMbeJ);
 
796
 
 
797
  }
 
798
  else if(params.ref == 1) { /** ROHF **/
 
799
 
 
800
 
 
801
    /* + Cjf <mb||ef> -> Wmbej*/
 
802
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, CME_lbl);
 
803
    dpd_file2_init(&Cme, EOM_Cme, C_irr, 0, 1, Cme_lbl);
 
804
 
 
805
    dpd_buf4_init(&F, CC_FINTS, 0, 10, 5, 10, 5, 1, "F <ia|bc>");
 
806
    dpd_buf4_init(&WMBEJ, CC_TMP0, C_irr, 10, 11, 10, 11, 0, "WMBEJ");
 
807
    dpd_contract424(&F, &CME, &WMBEJ, 3, 1, 0, 1, 0);
 
808
    dpd_buf4_close(&WMBEJ);
 
809
    dpd_buf4_init(&Wmbej, CC_TMP0, C_irr, 10, 11, 10, 11, 0, "Wmbej");
 
810
    dpd_contract424(&F, &Cme, &Wmbej, 3, 1, 0, 1, 0);
 
811
    dpd_buf4_close(&Wmbej);
 
812
    dpd_buf4_close(&F);
 
813
 
 
814
    dpd_buf4_init(&F, CC_FINTS, 0, 10, 5, 10, 5, 0, "F <ia|bc>");
 
815
    dpd_buf4_init(&WMbEj, CC_TMP0, C_irr, 10, 11, 10, 11, 0, "WMbEj");
 
816
    dpd_contract424(&F, &Cme, &WMbEj, 3, 1, 0, 1, 0);
 
817
    dpd_buf4_close(&WMbEj);
 
818
    dpd_buf4_init(&WmBeJ, CC_TMP0, C_irr, 10, 11, 10, 11, 0, "WmBeJ");
 
819
    dpd_contract424(&F, &CME, &WmBeJ, 3, 1, 0, 1, 0);
 
820
    dpd_buf4_close(&WmBeJ);
 
821
 
 
822
    dpd_buf4_init(&WMbeJ, CC_TMP0, C_irr, 10, 10, 10, 10, 0, "WMbeJ");
 
823
    dpd_contract244(&CME, &F, &WMbeJ, 1, 2, 1, -1, 0);
 
824
    dpd_buf4_close(&WMbeJ);
 
825
    dpd_buf4_init(&WmBEj, CC_TMP0, C_irr, 10, 10, 10, 10, 0, "WmBEj");
 
826
    dpd_contract244(&Cme, &F, &WmBEj, 1, 2, 1, -1, 0);
 
827
    dpd_buf4_close(&WmBEj);
 
828
    dpd_buf4_close(&F);
 
829
 
 
830
    /* - Cnb <mn||ej> -> Wmbej */
 
831
 
 
832
    dpd_buf4_init(&E, CC_EINTS, 0, 0, 11, 2, 11, 0, "E <ij||ka> (i>j,ak)");
 
833
    dpd_buf4_init(&WMBEJ, CC_TMP0, C_irr, 10, 11, 10, 11, 0, "WMBEJ");
 
834
    dpd_contract424(&E, &CME, &WMBEJ, 1, 0, 1, 1, 1);
 
835
    dpd_buf4_close(&WMBEJ);
 
836
    dpd_buf4_init(&Wmbej, CC_TMP0, C_irr, 10, 11, 10, 11, 0, "Wmbej");
 
837
    dpd_contract424(&E, &Cme, &Wmbej, 1, 0, 1, 1, 1);
 
838
    dpd_buf4_close(&Wmbej);
 
839
    dpd_buf4_close(&E);
 
840
 
 
841
    dpd_buf4_init(&E, CC_EINTS, 0, 11, 0, 11, 0, 0, "E <ai|jk>");
 
842
    dpd_buf4_init(&WMbEj, CC_TMP0, C_irr, 10, 11, 10, 11, 0, "WMbEj");
 
843
    dpd_contract424(&E, &Cme, &WMbEj, 3, 0, 1, -1, 1);
 
844
    dpd_buf4_close(&WMbEj);
 
845
    dpd_buf4_init(&WmBeJ, CC_TMP0, C_irr, 10, 11, 10, 11, 0, "WmBeJ");
 
846
    dpd_contract424(&E, &CME, &WmBeJ, 3, 0, 1, -1, 1);
 
847
    dpd_buf4_close(&WmBeJ);
 
848
    dpd_buf4_close(&E);
 
849
 
 
850
    dpd_buf4_init(&E, CC_EINTS, 0, 0, 10, 0, 10, 0, "E <ij|ka>");
 
851
    dpd_buf4_init(&WMbeJ, CC_TMP0, C_irr, 10, 10, 10, 10, 0, "WMbeJ");
 
852
    dpd_contract424(&E, &Cme, &WMbeJ, 1, 0, 1, 1, 1);
 
853
    dpd_buf4_close(&WMbeJ);
 
854
    dpd_buf4_init(&WmBEj, CC_TMP0, C_irr, 10, 10, 10, 10, 0, "WmBEj");
 
855
    dpd_contract424(&E, &CME, &WmBEj, 1, 0, 1, 1, 1);
 
856
    dpd_buf4_close(&WmBEj);
 
857
    dpd_buf4_close(&E);
 
858
 
 
859
    /* Convert to (ME,JB) for remaining terms */
 
860
 
 
861
    dpd_buf4_init(&WMBEJ, CC_TMP0, C_irr, 10, 11, 10, 11, 0, "WMBEJ");
 
862
    dpd_buf4_sort(&WMBEJ, CC3_HC1, prsq, 10, 10, "HC1 WMBEJ (ME,JB)");
 
863
    dpd_buf4_close(&WMBEJ);
 
864
 
 
865
    dpd_buf4_init(&Wmbej, CC_TMP0, C_irr, 10, 11, 10, 11, 0, "Wmbej");
 
866
    dpd_buf4_sort(&Wmbej, CC3_HC1, prsq, 10, 10, "HC1 Wmbej (me,jb)");
 
867
    dpd_buf4_close(&Wmbej);
 
868
 
 
869
    dpd_buf4_init(&WMbEj, CC_TMP0, C_irr, 10, 11, 10, 11, 0, "WMbEj");
 
870
    dpd_buf4_sort(&WMbEj, CC3_HC1, prsq, 10, 10, "HC1 WMbEj (ME,jb)");
 
871
    dpd_buf4_close(&WMbEj);
 
872
 
 
873
    dpd_buf4_init(&WmBeJ, CC_TMP0, C_irr, 10, 11, 10, 11, 0, "WmBeJ");
 
874
    dpd_buf4_sort(&WmBeJ, CC3_HC1, prsq, 10, 10, "HC1 WmBeJ (me,JB)");
 
875
    dpd_buf4_close(&WmBeJ);
 
876
 
 
877
    dpd_buf4_init(&WMbeJ, CC_TMP0, C_irr, 10, 10, 10, 10, 0, "WMbeJ");
 
878
    dpd_buf4_sort(&WMbeJ, CC3_HC1, psrq, 10, 10, "HC1 WMbeJ (Me,Jb)");
 
879
    dpd_buf4_close(&WMbeJ);
 
880
 
 
881
    dpd_buf4_init(&WmBEj, CC_TMP0, C_irr, 10, 10, 10, 10, 0, "WmBEj");
 
882
    dpd_buf4_sort(&WmBEj, CC3_HC1, psrq, 10, 10, "HC1 WmBEj (mE,jB)");
 
883
    dpd_buf4_close(&WmBEj);
 
884
 
 
885
  } /** ROHF **/
 
886
  else if(params.ref == 2) { /** UHF **/
 
887
 
 
888
    /* F -> Wmbej */ /* + Cjf <mb||ef> -> Wmbej*/
 
889
 
 
890
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, CME_lbl);
 
891
    dpd_file2_init(&Cme, EOM_Cme, C_irr, 2, 3, Cme_lbl);
 
892
 
 
893
    dpd_buf4_init(&W, CC_TMP0, C_irr, 20, 21, 20, 21, 0, "WMBEJ");
 
894
    dpd_buf4_init(&F, CC_FINTS, 0, 20, 5, 20, 5, 1, "F <IA|BC>");
 
895
    dpd_contract424(&F, &CME, &W, 3, 1, 0, 1, 0.0);
 
896
    dpd_buf4_close(&F);
 
897
    dpd_buf4_close(&W);
 
898
 
 
899
    dpd_buf4_init(&W, CC_TMP0, C_irr, 30, 31, 30, 31, 0, "Wmbej");
 
900
    dpd_buf4_init(&F, CC_FINTS, 0, 30, 15, 30, 15, 1, "F <ia|bc>");
 
901
    dpd_contract424(&F, &Cme, &W, 3, 1, 0, 1, 0.0);
 
902
    dpd_buf4_close(&F);
 
903
    dpd_buf4_close(&W);
 
904
 
 
905
    dpd_buf4_init(&W, CC_TMP0, C_irr, 24, 26, 24, 26, 0, "WMbEj");
 
906
    dpd_buf4_init(&F, CC_FINTS, 0, 24, 28, 24, 28, 0, "F <Ia|Bc>");
 
907
    dpd_contract424(&F, &Cme, &W, 3, 1, 0, 1, 0.0);
 
908
    dpd_buf4_close(&F);
 
909
    dpd_buf4_close(&W);
 
910
 
 
911
    dpd_buf4_init(&W, CC_TMP0, C_irr, 27, 25, 27, 25, 0, "WmBeJ");
 
912
    dpd_buf4_init(&F, CC_FINTS, 0, 27, 29, 27, 29, 0, "F <iA|bC>");
 
913
    dpd_contract424(&F, &CME, &W, 3, 1, 0, 1, 0.0);
 
914
    dpd_buf4_close(&F);
 
915
    dpd_buf4_close(&W);
 
916
 
 
917
    dpd_buf4_init(&W, CC_TMP0, C_irr, 24, 24, 24, 24, 0, "WMbeJ");
 
918
    dpd_buf4_init(&F, CC_FINTS, 0, 24, 28, 24, 28, 0, "F <Ia|Bc>");
 
919
    dpd_contract244(&CME, &F, &W, 1, 2, 1, -1, 0.0);
 
920
    dpd_buf4_close(&F);
 
921
    dpd_buf4_close(&W);
 
922
 
 
923
    dpd_buf4_init(&W, CC_TMP0, C_irr, 27, 27, 27, 27, 0, "WmBEj");
 
924
    dpd_buf4_init(&F, CC_FINTS, 0, 27, 29, 27, 29, 0, "F <iA|bC>");
 
925
    dpd_contract244(&Cme, &F, &W, 1, 2, 1, -1, 0.0);
 
926
    dpd_buf4_close(&F);
 
927
    dpd_buf4_close(&W);
 
928
 
 
929
    /* - Cnb <mn||ej> -> Wmbej */
 
930
 
 
931
    dpd_buf4_init(&W, CC_TMP0, C_irr, 20, 21, 20, 21, 0, "WMBEJ");
 
932
    dpd_buf4_init(&E, CC_EINTS, 0, 0, 21, 2, 21, 0, "E <IJ||KA> (I>J,AK)");
 
933
    dpd_contract424(&E, &CME, &W, 1, 0, 1, 1, 1);
 
934
    dpd_buf4_close(&E);
 
935
    dpd_buf4_close(&W);
 
936
 
 
937
    dpd_buf4_init(&W, CC_TMP0, C_irr, 30, 31, 30, 31, 0, "Wmbej");
 
938
    dpd_buf4_init(&E, CC_EINTS, 0, 10, 31, 12, 31, 0, "E <ij||ka> (i>j,ak)");
 
939
    dpd_contract424(&E, &Cme, &W, 1, 0, 1, 1, 1);
 
940
    dpd_buf4_close(&E);
 
941
    dpd_buf4_close(&W);
 
942
 
 
943
    dpd_buf4_init(&W, CC_TMP0, C_irr, 24, 26, 24, 26, 0, "WMbEj");
 
944
    dpd_buf4_init(&E, CC_EINTS, 0, 22, 26, 22, 26, 0, "E <Ij|Ak>");
 
945
    dpd_contract424(&E, &Cme, &W, 1, 0, 1, -1, 1);
 
946
    dpd_buf4_close(&E);
 
947
    dpd_buf4_close(&W);
 
948
 
 
949
    dpd_buf4_init(&W, CC_TMP0, C_irr, 27, 25, 27, 25, 0, "WmBeJ");
 
950
    dpd_buf4_init(&E, CC_EINTS, 0, 23, 25, 23, 25, 0, "E <iJ|aK>");
 
951
    dpd_contract424(&E, &CME, &W, 1, 0, 1, -1, 1);
 
952
    dpd_buf4_close(&E);
 
953
    dpd_buf4_close(&W);
 
954
 
 
955
    dpd_buf4_init(&W, CC_TMP0, C_irr, 24, 24, 24, 24, 0, "WMbeJ");
 
956
    dpd_buf4_init(&E, CC_EINTS, 0, 22, 24, 22, 24, 0, "E <Ij|Ka>");
 
957
    dpd_contract424(&E, &Cme, &W, 1, 0, 1, 1, 1);
 
958
    dpd_buf4_close(&E);
 
959
    dpd_buf4_close(&W);
 
960
 
 
961
    dpd_buf4_init(&W, CC_TMP0, C_irr, 27, 27, 27, 27, 0, "WmBEj");
 
962
    dpd_buf4_init(&E, CC_EINTS, 0, 23, 27, 23, 27, 0, "E <iJ|kA>");
 
963
    dpd_contract424(&E, &CME, &W, 1, 0, 1, 1, 1);
 
964
    dpd_buf4_close(&E);
 
965
    dpd_buf4_close(&W);
 
966
 
 
967
    /* Convert to (ME,JB) for remaining terms */
 
968
 
 
969
    dpd_buf4_init(&W, CC_TMP0, C_irr, 20, 21, 20, 21, 0, "WMBEJ");
 
970
    dpd_buf4_sort(&W, CC3_HC1, prsq, 20, 20, "HC1 WMBEJ (ME,JB)");
 
971
    dpd_buf4_close(&W);
 
972
 
 
973
    dpd_buf4_init(&W, CC_TMP0, C_irr, 30, 31, 30, 31, 0, "Wmbej");
 
974
    dpd_buf4_sort(&W, CC3_HC1, prsq, 30, 30, "HC1 Wmbej (me,jb)");
 
975
    dpd_buf4_close(&W);
 
976
 
 
977
    dpd_buf4_init(&W, CC_TMP0, C_irr, 24, 26, 24, 26, 0, "WMbEj");
 
978
    dpd_buf4_sort(&W, CC3_HC1, prsq, 20, 30, "HC1 WMbEj (ME,jb)");
 
979
    dpd_buf4_close(&W);
 
980
 
 
981
    dpd_buf4_init(&W, CC_TMP0, C_irr, 27, 25, 27, 25, 0, "WmBeJ");
 
982
    dpd_buf4_sort(&W, CC3_HC1, prsq, 30, 20, "HC1 WmBeJ (me,JB)");
 
983
    dpd_buf4_close(&W);
 
984
 
 
985
    dpd_buf4_init(&W, CC_TMP0, C_irr, 24, 24, 24, 24, 0, "WMbeJ");
 
986
    dpd_buf4_sort(&W, CC3_HC1, psrq, 24, 24, "HC1 WMbeJ (Me,Jb)");
 
987
    dpd_buf4_close(&W);
 
988
 
 
989
    dpd_buf4_init(&W, CC_TMP0, C_irr, 27, 27, 27, 27, 0, "WmBEj");
 
990
    dpd_buf4_sort(&W, CC3_HC1, psrq, 27, 27, "HC1 WmBEj (mE,jB)");
 
991
    dpd_buf4_close(&W);
 
992
  }
 
993
 
 
994
  return;
 
995
}
 
996
 
 
997
 
 
998
 
 
999
void HC1_Wabei(int i, int C_irr) { 
 
1000
  dpdbuf4 Z, Z1, Z2, Z3;
 
1001
  dpdbuf4 B, C, D, E, F, W;
 
1002
  dpdfile2 Cme, CME;
 
1003
  char CME_lbl[32], Cme_lbl[32];
 
1004
 
 
1005
  sprintf(CME_lbl, "%s %d", "CME", i);
 
1006
  sprintf(Cme_lbl, "%s %d", "Cme", i);
 
1007
 
 
1008
  if(params.ref == 0) { /** RHF **/
 
1009
 
 
1010
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, CME_lbl);
 
1011
 
 
1012
    dpd_buf4_init(&Z1, CC_TMP0, C_irr, 5, 11, 5, 11, 0, "CC3 Z(Ab,Ei)");
 
1013
    dpd_buf4_init(&Z2, CC_TMP0, C_irr, 11, 5, 11, 5, 0, "CC3 Z(Ei,Ab)");
 
1014
 
 
1015
    /* Z1(Ab,Ei) <-- <Ab|Ef> * C(i,f) */
 
1016
    dpd_buf4_init(&B, CC_BINTS, 0, 5, 5, 5, 5, 0, "B <ab|cd>");
 
1017
    dpd_contract424(&B, &CME, &Z1, 3, 1, 0, 1, 0);
 
1018
    dpd_buf4_close(&B);
 
1019
 
 
1020
    /* Z1(Ab,Ei) <--  - C(M,A) * <Mb|Ei> */
 
1021
    dpd_buf4_init(&D, CC_DINTS, 0, 10, 11, 10, 11, 0, "D <ij|ab> (ib,aj)");
 
1022
    dpd_contract244(&CME, &D, &Z1, 0, 0, 0, -1, 1);
 
1023
    dpd_buf4_close(&D);
 
1024
 
 
1025
    /* Z2(Ei,Ab) <-- - <mA,iE> C(m,b) */
 
1026
    dpd_buf4_init(&C, CC_CINTS, 0, 10, 10, 10, 10, 0, "C <ia|jb>");
 
1027
    dpd_buf4_sort(&C, CC_TMP0, qpsr, 11, 11, "CC3 Z(Ei,Am)");
 
1028
    dpd_buf4_close(&C);
 
1029
 
 
1030
    dpd_buf4_init(&Z, CC_TMP0, 0, 11, 11, 11, 11, 0, "CC3 Z(Ei,Am)");
 
1031
    dpd_contract424(&Z, &CME, &Z2, 3, 0, 0, -1, 0);
 
1032
    dpd_buf4_close(&Z);
 
1033
 
 
1034
    dpd_buf4_close(&Z2);
 
1035
 
 
1036
    /* W(Ab,Ei) = Z1(Ab,Ei) + Z2(Ei,Ab) */
 
1037
    dpd_buf4_sort_axpy(&Z1, CC_TMP0, rspq, 11, 5, "CC3 Z(Ei,Ab)", 1);
 
1038
    dpd_buf4_close(&Z1);
 
1039
    dpd_buf4_init(&Z2, CC_TMP0, C_irr, 11, 5, 11, 5, 0, "CC3 Z(Ei,Ab)");
 
1040
    dpd_buf4_sort(&Z2, CC3_HC1, qpsr, 10, 5, "CC3 WAbEi (Ie,Ab)");
 
1041
 
 
1042
    dpd_file2_close(&CME);
 
1043
  }
 
1044
 
 
1045
  else if (params.ref == 1) { /* ROHF */
 
1046
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, CME_lbl);
 
1047
    dpd_file2_init(&Cme, EOM_Cme, C_irr, 0, 1, Cme_lbl);
 
1048
 
 
1049
    dpd_buf4_init(&Z1, CC_TMP0, C_irr, 7, 11, 7, 11, 0, "WABEI (A>B,EI)");
 
1050
    dpd_buf4_init(&B, CC_BINTS, 0, 7, 5, 5, 5, 1, "B <ab|cd>");
 
1051
    dpd_contract424(&B, &CME, &Z1, 3, 1, 0, 1.0, 0.0);
 
1052
    dpd_buf4_close(&B);
 
1053
    dpd_buf4_close(&Z1);
 
1054
 
 
1055
    dpd_buf4_init(&Z1, CC_TMP0, C_irr, 7, 11, 7, 11, 0, "Wabei (a>b,ei)");
 
1056
    dpd_buf4_init(&B, CC_BINTS, 0, 7, 5, 5, 5, 1, "B <ab|cd>");
 
1057
    dpd_contract424(&B, &Cme, &Z1, 3, 1, 0, 1.0, 0.0);
 
1058
    dpd_buf4_close(&B);
 
1059
    dpd_buf4_close(&Z1);
 
1060
 
 
1061
    dpd_buf4_init(&Z1, CC_TMP0, C_irr, 5, 11, 5, 11, 0, "WAbEi (Ab,Ei)");
 
1062
    dpd_buf4_init(&B, CC_BINTS, 0, 5, 5, 5, 5, 0, "B <ab|cd>");
 
1063
    dpd_contract424(&B, &Cme, &Z1, 3, 1, 0, 1.0, 0.0);
 
1064
    dpd_buf4_close(&B);
 
1065
    dpd_buf4_close(&Z1);
 
1066
 
 
1067
    dpd_buf4_init(&Z1, CC_TMP0, C_irr, 5, 11, 5, 11, 0, "WaBeI (aB,eI)");
 
1068
    dpd_buf4_init(&B, CC_BINTS, 0, 5, 5, 5, 5, 0, "B <ab|cd>");
 
1069
    dpd_contract424(&B, &CME, &Z1, 3, 1, 0, 1.0, 0.0);
 
1070
    dpd_buf4_close(&B);
 
1071
    dpd_buf4_close(&Z1);
 
1072
 
 
1073
    /** -CMA <MI||EB> + CMB <MI||EA> **/
 
1074
 
 
1075
    dpd_buf4_init(&Z1, CC_TMP1, C_irr, 5, 11, 5, 11, 0, "Z (AB,EI)");
 
1076
 
 
1077
    dpd_buf4_init(&C, CC_CINTS, 0, 10, 11, 10, 11, 0, "C <ia||jb> (ia,bj)");
 
1078
    dpd_contract244(&CME, &C, &Z1, 0, 0, 0, 1.0, 0.0);
 
1079
    dpd_buf4_close(&C);
 
1080
 
 
1081
    dpd_buf4_sort(&Z1, CC_TMP1, qprs, 5, 11, "Z (BA,EI)");
 
1082
 
 
1083
    dpd_buf4_init(&W, CC_TMP0, C_irr, 5, 11, 7, 11, 0, "WABEI (A>B,EI)");
 
1084
    dpd_buf4_axpy(&Z1, &W, 1.0);
 
1085
    dpd_buf4_close(&Z1);
 
1086
    dpd_buf4_init(&Z1, CC_TMP1, C_irr, 5, 11, 5, 11, 0, "Z (BA,EI)");
 
1087
    dpd_buf4_axpy(&Z1, &W, -1.0);
 
1088
    dpd_buf4_close(&Z1);
 
1089
    dpd_buf4_close(&W);
 
1090
 
 
1091
    /** -Cma <mi||eb> + Cmb <mi||ea> **/
 
1092
 
 
1093
    dpd_buf4_init(&Z1, CC_TMP1, C_irr, 5, 11, 5, 11, 0, "Z (ab,ei)");
 
1094
 
 
1095
    dpd_buf4_init(&C, CC_CINTS, 0, 10, 11, 10, 11, 0, "C <ia||jb> (ia,bj)");
 
1096
    dpd_contract244(&Cme, &C, &Z1, 0, 0, 0, 1.0, 0.0);
 
1097
    dpd_buf4_close(&C);
 
1098
 
 
1099
    dpd_buf4_sort(&Z1, CC_TMP1, qprs, 5, 11, "Z (ba,ei)");
 
1100
 
 
1101
    dpd_buf4_init(&W, CC_TMP0, C_irr, 5, 11, 7, 11, 0, "Wabei (a>b,ei)");
 
1102
    dpd_buf4_axpy(&Z1, &W, 1.0);
 
1103
    dpd_buf4_close(&Z1);
 
1104
    dpd_buf4_init(&Z1, CC_TMP1, C_irr, 5, 11, 5, 11, 0, "Z (ba,ei)");
 
1105
    dpd_buf4_axpy(&Z1, &W, -1.0);
 
1106
    dpd_buf4_close(&Z1);
 
1107
    dpd_buf4_close(&W);
 
1108
 
 
1109
    /** -CMA <Mi||Eb> - Cmb <mA||iE> **/
 
1110
 
 
1111
    dpd_buf4_init(&Z1, CC_TMP0, C_irr, 5, 11, 5, 11, 0, "WAbEi (Ab,Ei)");
 
1112
    dpd_buf4_init(&D, CC_DINTS, 0, 10, 11, 10, 11, 0, "D <ij|ab> (ib,aj)");
 
1113
    dpd_contract244(&CME, &D, &Z1, 0, 0, 0, -1.0, 1.0);
 
1114
    dpd_buf4_close(&D);
 
1115
    dpd_buf4_close(&Z1);
 
1116
 
 
1117
    dpd_buf4_init(&Z1, CC_TMP1, C_irr, 5, 11, 5, 11, 0, "WAbEi (bA,Ei)");
 
1118
    dpd_buf4_init(&C, CC_CINTS, 0, 10, 11, 10, 11, 0, "C <ia|jb> (ia,bj)");
 
1119
    dpd_contract244(&Cme, &C, &Z1, 0, 0, 0, -1.0, 0.0);
 
1120
    dpd_buf4_close(&C);
 
1121
 
 
1122
    dpd_buf4_sort_axpy(&Z1, CC_TMP0, qprs, 5, 11, "WAbEi (Ab,Ei)", 1.0);
 
1123
    dpd_buf4_close(&Z1);
 
1124
 
 
1125
    /** **/
 
1126
    dpd_buf4_init(&Z1, CC_TMP0, C_irr, 5, 11, 5, 11, 0, "WaBeI (aB,eI)");
 
1127
    dpd_buf4_init(&D, CC_DINTS, 0, 10, 11, 10, 11, 0, "D <ij|ab> (ib,aj)");
 
1128
    dpd_contract244(&Cme, &D, &Z1, 0, 0, 0, -1.0, 1.0);
 
1129
    dpd_buf4_close(&D);
 
1130
    dpd_buf4_close(&Z1);
 
1131
 
 
1132
    dpd_buf4_init(&Z1, CC_TMP1, C_irr, 5, 11, 5, 11, 0, "WaBeI (Ba,eI)");
 
1133
    dpd_buf4_init(&C, CC_CINTS, 0, 10, 11, 10, 11, 0, "C <ia|jb> (ia,bj)");
 
1134
    dpd_contract244(&CME, &C, &Z1, 0, 0, 0, -1.0, 0.0);
 
1135
    dpd_buf4_close(&C);
 
1136
 
 
1137
    dpd_buf4_sort_axpy(&Z1, CC_TMP0, qprs, 5, 11, "WaBeI (aB,eI)", 1.0);
 
1138
    dpd_buf4_close(&Z1);
 
1139
 
 
1140
    dpd_file2_close(&CME);
 
1141
    dpd_file2_close(&Cme);
 
1142
 
 
1143
    /* final sort to (EI,AB) */
 
1144
 
 
1145
    dpd_buf4_init(&W, CC_TMP0, C_irr, 5, 11, 7, 11, 0, "WABEI (A>B,EI)");
 
1146
    dpd_buf4_sort(&W, CC_TMP0, rspq, 11, 7, "WABEI (EI,A>B)");
 
1147
    dpd_buf4_close(&W);
 
1148
    dpd_buf4_init(&W, CC_TMP0, C_irr, 5, 11, 7, 11, 0, "Wabei (a>b,ei)");
 
1149
    dpd_buf4_sort(&W, CC_TMP0, rspq, 11, 7, "Wabei (ei,a>b)");
 
1150
    dpd_buf4_close(&W);
 
1151
    dpd_buf4_init(&W, CC_TMP0, C_irr, 5, 11, 5, 11, 0, "WAbEi (Ab,Ei)");
 
1152
    dpd_buf4_sort(&W, CC_TMP0, rspq, 11, 5, "WAbEi (Ei,Ab)");
 
1153
    dpd_buf4_close(&W);
 
1154
    dpd_buf4_init(&W, CC_TMP0, C_irr, 5, 11, 5, 11, 0, "WaBeI (aB,eI)");
 
1155
    dpd_buf4_sort(&W, CC_TMP0, rspq, 11, 5, "WaBeI (eI,aB)");
 
1156
    dpd_buf4_close(&W);
 
1157
 
 
1158
    purge_Wabei(C_irr);
 
1159
 
 
1160
    /* final sort to Wabei (IE,AB) */
 
1161
    dpd_buf4_init(&W, CC_TMP0, C_irr, 11, 7, 11, 7, 0, "WABEI (EI,A>B)");
 
1162
    dpd_buf4_sort(&W, CC3_HC1, qprs, 10, 7, "HC1 WABEI (IE,A>B)");
 
1163
    dpd_buf4_close(&W);
 
1164
    dpd_buf4_init(&W, CC_TMP0, C_irr, 11, 7, 11, 7, 0, "Wabei (ei,a>b)");
 
1165
    dpd_buf4_sort(&W, CC3_HC1, qprs, 10, 7, "HC1 Wabei (ie,a>b)");
 
1166
    dpd_buf4_close(&W);
 
1167
    dpd_buf4_init(&W, CC_TMP0, C_irr, 11, 5, 11, 5, 0, "WAbEi (Ei,Ab)");
 
1168
    dpd_buf4_sort(&W, CC3_HC1, qprs, 10, 5, "HC1 WAbEi (iE,Ab)");
 
1169
    dpd_buf4_close(&W);
 
1170
    dpd_buf4_init(&W, CC_TMP0, C_irr, 11, 5, 11, 5, 0, "WaBeI (eI,aB)");
 
1171
    dpd_buf4_sort(&W, CC3_HC1, qprs, 10, 5, "HC1 WaBeI (Ie,aB)");
 
1172
    dpd_buf4_close(&W);
 
1173
  }
 
1174
  else if (params.ref == 2) { /* UHF */
 
1175
    dpd_file2_init(&CME, EOM_CME, C_irr, 0, 1, CME_lbl);
 
1176
    dpd_file2_init(&Cme, EOM_Cme, C_irr, 2, 3, Cme_lbl);
 
1177
 
 
1178
    /* term 1, Cif <ab||ef> */
 
1179
 
 
1180
    dpd_buf4_init(&Z1, CC_TMP0, C_irr, 7, 21, 7, 21, 0, "WABEI (A>B,EI)");
 
1181
    dpd_buf4_init(&B, CC_BINTS, 0, 7, 5, 5, 5, 1, "B <AB|CD>");
 
1182
    dpd_contract424(&B, &CME, &Z1, 3, 1, 0, 1.0, 0.0);
 
1183
    dpd_buf4_close(&B);
 
1184
    dpd_buf4_close(&Z1);
 
1185
 
 
1186
    dpd_buf4_init(&Z1, CC_TMP0, C_irr, 17, 31, 17, 31, 0, "Wabei (a>b,ei)");
 
1187
    dpd_buf4_init(&B, CC_BINTS, 0, 17, 15, 15, 15, 1, "B <ab|cd>");
 
1188
    dpd_contract424(&B, &Cme, &Z1, 3, 1, 0, 1.0, 0.0);
 
1189
    dpd_buf4_close(&B);
 
1190
    dpd_buf4_close(&Z1);
 
1191
 
 
1192
    dpd_buf4_init(&Z1, CC_TMP0, C_irr, 28, 26, 28, 26, 0, "WAbEi (Ab,Ei)");
 
1193
    dpd_buf4_init(&B, CC_BINTS, 0, 28, 28, 28, 28, 0, "B <Ab|Cd>");
 
1194
    dpd_contract424(&B, &Cme, &Z1, 3, 1, 0, 1.0, 0.0);
 
1195
    dpd_buf4_close(&B);
 
1196
    dpd_buf4_close(&Z1);
 
1197
 
 
1198
    dpd_buf4_init(&Z, CC_TMP0, C_irr, 24, 28, 24, 28, 0, "Z(Ie,Ba)");
 
1199
    dpd_buf4_init(&B, CC_BINTS, 0, 28, 28, 28, 28, 0, "B <Ab|Cd>");
 
1200
    dpd_contract244(&CME, &B, &Z, 1, 0, 0, 1.0, 0.0);
 
1201
    dpd_buf4_close(&B);
 
1202
    /** Z(Ie,Ba) --> W'(aB,eI) **/
 
1203
    /* srqp seems to have a bug
 
1204
     dpd_buf4_sort(&Z, CC_TMP0, srqp, 29, 25, "WaBeI (aB,eI)");
 
1205
    */
 
1206
    dpd_buf4_sort(&Z, CC_TMP0, rspq, 28, 24, "WaBeI (Ba,Ie) 1");
 
1207
    dpd_buf4_close(&Z);
 
1208
    dpd_buf4_init(&Z, CC_TMP0, C_irr, 28, 24, 28, 24, 0, "WaBeI (Ba,Ie)");
 
1209
    dpd_buf4_sort(&Z, CC_TMP0, qpsr, 29, 25, "WaBeI (aB,eI)");
 
1210
    dpd_buf4_close(&Z);
 
1211
 
 
1212
    /** UHF term 2 **/
 
1213
    /* -CMA <MI||EB> + CMB <MI||EA> **/
 
1214
 
 
1215
    dpd_buf4_init(&Z1, CC_TMP1, C_irr, 5, 21, 5, 21, 0, "Z (AB,EI)");
 
1216
 
 
1217
    dpd_buf4_init(&C, CC_CINTS, 0, 20, 21, 20, 21, 0, "C <IA||JB> (IA,BJ)");
 
1218
    dpd_contract244(&CME, &C, &Z1, 0, 0, 0, 1.0, 0.0);
 
1219
    dpd_buf4_close(&C);
 
1220
 
 
1221
    dpd_buf4_sort(&Z1, CC_TMP1, qprs, 5, 21, "Z (BA,EI)");
 
1222
 
 
1223
    dpd_buf4_init(&W, CC_TMP0, C_irr, 5, 21, 7, 21, 0, "WABEI (A>B,EI)");
 
1224
    dpd_buf4_axpy(&Z1, &W, 1.0);
 
1225
    dpd_buf4_close(&Z1);
 
1226
    dpd_buf4_init(&Z1, CC_TMP1, C_irr, 5, 21, 5, 21, 0, "Z (BA,EI)");
 
1227
    dpd_buf4_axpy(&Z1, &W, -1.0);
 
1228
    dpd_buf4_close(&Z1);
 
1229
 
 
1230
    /** -Cma <mi||eb> + Cmb <mi||ea> **/
 
1231
 
 
1232
    dpd_buf4_init(&Z1, CC_TMP1, C_irr, 15, 31, 15, 31, 0, "Z (ab,ei)");
 
1233
 
 
1234
    dpd_buf4_init(&C, CC_CINTS, 0, 30, 31, 30, 31, 0, "C <ia||jb> (ia,bj)");
 
1235
    dpd_contract244(&Cme, &C, &Z1, 0, 0, 0, 1.0, 0.0);
 
1236
    dpd_buf4_close(&C);
 
1237
 
 
1238
    dpd_buf4_sort(&Z1, CC_TMP1, qprs, 15, 31, "Z (ba,ei)");
 
1239
 
 
1240
    dpd_buf4_init(&W, CC_TMP0, C_irr, 15, 31, 17, 31, 0, "Wabei (a>b,ei)");
 
1241
    dpd_buf4_axpy(&Z1, &W, 1.0);
 
1242
    dpd_buf4_close(&Z1);
 
1243
    dpd_buf4_init(&Z1, CC_TMP1, C_irr, 15, 31, 15, 31, 0, "Z (ba,ei)");
 
1244
    dpd_buf4_axpy(&Z1, &W, -1.0);
 
1245
    dpd_buf4_close(&Z1);
 
1246
 
 
1247
    /** -CMA <Mi||Eb> - Cmb <mA|iE> **/
 
1248
 
 
1249
    dpd_buf4_init(&Z1, CC_TMP0, C_irr, 28, 26, 28, 26, 0, "WAbEi (Ab,Ei)");
 
1250
    dpd_buf4_init(&D, CC_DINTS, 0, 24, 26, 24, 26, 0, "D <Ij|Ab> (Ib,Aj)");
 
1251
    dpd_contract244(&CME, &D, &Z1, 0, 0, 0, -1.0, 1.0);
 
1252
    dpd_buf4_close(&D);
 
1253
    dpd_buf4_close(&Z1);
 
1254
 
 
1255
    dpd_buf4_init(&Z1, CC_TMP1, C_irr, 29, 26, 29, 26, 0, "WAbEi (bA,Ei)");
 
1256
    dpd_buf4_init(&C, CC_CINTS, 0, 27, 26, 27, 26, 0, "C <Ai|Bj> (iA,Bj)");
 
1257
    dpd_contract244(&Cme, &C, &Z1, 0, 0, 0, -1.0, 0.0);
 
1258
    dpd_buf4_close(&C);
 
1259
 
 
1260
    dpd_buf4_sort_axpy(&Z1, CC_TMP0, qprs, 28, 26, "WAbEi (Ab,Ei)", 1.0);
 
1261
    dpd_buf4_close(&Z1);
 
1262
 
 
1263
    /**  **/
 
1264
    dpd_buf4_init(&Z1, CC_TMP0, C_irr, 29, 25, 29, 25, 0, "WaBeI (aB,eI)");
 
1265
    dpd_buf4_init(&D, CC_DINTS, 0, 27, 25, 27, 25, 0, "D <iJ|aB> (iB,aJ)");
 
1266
    dpd_contract244(&Cme, &D, &Z1, 0, 0, 0, -1.0, 1.0);
 
1267
    dpd_buf4_close(&D);
 
1268
    dpd_buf4_close(&Z1);
 
1269
 
 
1270
    dpd_buf4_init(&Z1, CC_TMP1, C_irr, 28, 25, 28, 25, 0, "WaBeI (Ba,eI)");
 
1271
    dpd_buf4_init(&C, CC_CINTS, 0, 24, 25, 24, 25, 0, "C <Ia|Jb> (Ia,bJ)");
 
1272
    dpd_contract244(&CME, &C, &Z1, 0, 0, 0, -1.0, 0.0);
 
1273
    dpd_buf4_close(&C);
 
1274
 
 
1275
    dpd_buf4_sort_axpy(&Z1, CC_TMP0, qprs, 29, 25, "WaBeI (aB,eI)", 1.0);
 
1276
    dpd_buf4_close(&Z1);
 
1277
 
 
1278
    /* final sort and storage */
 
1279
    dpd_buf4_init(&W, CC_TMP0, C_irr, 7, 21, 7, 21, 0, "WABEI (A>B,EI)");
 
1280
    dpd_buf4_sort(&W, CC_TMP0, rspq, 21, 7, "WABEI (EI,A>B)");
 
1281
    dpd_buf4_close(&W);
 
1282
 
 
1283
    dpd_buf4_init(&W, CC_TMP0, C_irr, 21, 7, 21, 7, 0, "WABEI (EI,A>B)");
 
1284
    dpd_buf4_sort(&W, CC3_HC1, qprs, 20, 7, "HC1 WABEI (IE,B>A)");
 
1285
    dpd_buf4_close(&W);
 
1286
 
 
1287
    dpd_buf4_init(&W, CC3_HC1, C_irr, 20, 7, 20, 7, 0, "HC1 WABEI (IE,B>A)");
 
1288
    dpd_buf4_scm(&W, -1.0);
 
1289
    dpd_buf4_close(&W);
 
1290
 
 
1291
    /* final sort and storage */
 
1292
    dpd_buf4_init(&W, CC_TMP0, C_irr, 17, 31, 17, 31, 0, "Wabei (a>b,ei)");
 
1293
    dpd_buf4_sort(&W, CC_TMP0, rspq, 31, 17, "Wabei (ei,a>b)");
 
1294
    dpd_buf4_close(&W);
 
1295
 
 
1296
    dpd_buf4_init(&W, CC_TMP0, C_irr, 31, 17, 31, 17, 0, "Wabei (ei,a>b)");
 
1297
    dpd_buf4_sort(&W, CC3_HC1, qprs, 30, 17, "HC1 Wabei (ie,b>a)");
 
1298
    dpd_buf4_close(&W);
 
1299
 
 
1300
    dpd_buf4_init(&W, CC3_HC1, C_irr, 30, 17, 30, 17, 0, "HC1 Wabei (ie,b>a)");
 
1301
    dpd_buf4_scm(&W, -1.0);
 
1302
    dpd_buf4_close(&W);
 
1303
 
 
1304
    /* final sort and storage */
 
1305
    dpd_buf4_init(&W, CC_TMP0, C_irr, 28, 26, 28, 26, 0, "WAbEi (Ab,Ei)");
 
1306
    dpd_buf4_sort(&W, CC_TMP0, rspq, 26, 28, "WAbEi (Ei,Ab)");
 
1307
    dpd_buf4_close(&W);
 
1308
    dpd_buf4_init(&W, CC_TMP0, C_irr, 26, 28, 26, 28, 0, "WAbEi (Ei,Ab)");
 
1309
    dpd_buf4_sort(&W, CC3_HC1, qpsr, 27, 29, "HC1 WAbEi (iE,bA)");
 
1310
    dpd_buf4_close(&W);
 
1311
 
 
1312
    /* final sort and storage */
 
1313
    dpd_buf4_init(&W, CC_TMP0, C_irr, 29, 25, 29, 25, 0, "WaBeI (aB,eI)");
 
1314
    dpd_buf4_sort(&W, CC_TMP0, rspq, 25, 29, "WaBeI (eI,aB)");
 
1315
    dpd_buf4_close(&W);
 
1316
    dpd_buf4_init(&W, CC_TMP0, C_irr, 25, 29, 25, 29, 0, "WaBeI (eI,aB)");
 
1317
    dpd_buf4_sort(&W, CC3_HC1, qpsr, 24, 28, "HC1 WaBeI (Ie,Ba)");
 
1318
    dpd_buf4_close(&W);
 
1319
 
 
1320
    dpd_file2_close(&CME);
 
1321
    dpd_file2_close(&Cme);
 
1322
  }
 
1323
}
 
1324
 
 
1325
void purge_HC1(int C_irr) {
 
1326
  dpdfile2 FAE, Fmi, FME, Fme;
 
1327
  dpdfile4 W;
 
1328
  int *occpi, *virtpi;
 
1329
  int h, a, b, e, f, i, j, m, n;
 
1330
  int    A, B, E, F, I, J, M, N;
 
1331
  int mn, ei, ma, ef, me, jb, mb, ij, ab;
 
1332
  int asym, bsym, esym, fsym, isym, jsym, msym, nsym;
 
1333
  int *occ_off, *vir_off;
 
1334
  int *occ_sym, *vir_sym;
 
1335
  int *openpi, nirreps;
 
1336
 
 
1337
  nirreps = moinfo.nirreps;
 
1338
  occpi = moinfo.occpi; virtpi = moinfo.virtpi;
 
1339
  occ_off = moinfo.occ_off; vir_off = moinfo.vir_off;
 
1340
  occ_sym = moinfo.occ_sym; vir_sym = moinfo.vir_sym;
 
1341
  openpi = moinfo.openpi;
 
1342
 
 
1343
  /* Purge FME matrix elements */
 
1344
  dpd_file2_init(&FME, CC3_HC1, C_irr, 0, 1, "HC1 FME");
 
1345
  dpd_file2_mat_init(&FME);
 
1346
  dpd_file2_mat_rd(&FME);
 
1347
  for(h=0; h < nirreps; h++) {
 
1348
    for(m=0; m<occpi[h]; m++)
 
1349
      for(e=(virtpi[h]-openpi[h]); e<virtpi[h]; e++)
 
1350
        FME.matrix[h][m][e] = 0.0;
 
1351
  }
 
1352
  dpd_file2_mat_wrt(&FME);
 
1353
  dpd_file2_mat_close(&FME);
 
1354
  dpd_file2_close(&FME);
 
1355
 
 
1356
  /* Purge Fme matrix elements */
 
1357
  dpd_file2_init(&Fme, CC3_HC1, C_irr, 0, 1, "HC1 Fme");
 
1358
  dpd_file2_mat_init(&Fme);
 
1359
  dpd_file2_mat_rd(&Fme);
 
1360
  for(h=0; h < nirreps; h++) {
 
1361
    for(e=0; e<virtpi[h]; e++)
 
1362
      for(m=(occpi[h]-openpi[h]); m<occpi[h]; m++)
 
1363
        Fme.matrix[h][m][e] = 0.0;
 
1364
  }
 
1365
  dpd_file2_mat_wrt(&Fme);
 
1366
  dpd_file2_mat_close(&Fme);
 
1367
  dpd_file2_close(&Fme);
 
1368
 
 
1369
  /* Purge Wmnij matrix elements */
 
1370
  dpd_file4_init(&W, CC3_HC1, C_irr, 2, 2,"HC1 Wmnij (m>n,i>j)");
 
1371
  for(h=0; h < nirreps; h++) {
 
1372
    dpd_file4_mat_irrep_init(&W, h);
 
1373
    dpd_file4_mat_irrep_rd(&W, h);
 
1374
    for(mn=0; mn < W.params->rowtot[h]; mn++) {
 
1375
      m = W.params->roworb[h][mn][0];
 
1376
      n = W.params->roworb[h][mn][1];
 
1377
      msym = W.params->psym[m];
 
1378
      nsym = W.params->qsym[n];
 
1379
      M = m - occ_off[msym];
 
1380
      N = n - occ_off[nsym];
 
1381
      for(ij=0; ij < W.params->coltot[h]; ij++) {
 
1382
        i = W.params->colorb[h][ij][0];
 
1383
        j = W.params->colorb[h][ij][1];
 
1384
        isym = W.params->rsym[i];
 
1385
        jsym = W.params->ssym[j];
 
1386
        I = i - occ_off[isym];
 
1387
        J = j - occ_off[jsym];
 
1388
        if ((I >= (occpi[isym] - openpi[isym])) ||
 
1389
            (J >= (occpi[jsym] - openpi[jsym])) ||
 
1390
            (M >= (occpi[msym] - openpi[msym])) ||
 
1391
            (N >= (occpi[nsym] - openpi[nsym])) )
 
1392
          W.matrix[h][mn][ij] = 0.0;
 
1393
      }
 
1394
    }
 
1395
    dpd_file4_mat_irrep_wrt(&W, h);
 
1396
    dpd_file4_mat_irrep_close(&W, h);
 
1397
  }
 
1398
  dpd_file4_close(&W);
 
1399
 
 
1400
  dpd_file4_init(&W, CC3_HC1, C_irr, 0, 0,"HC1 WMnIj (Mn,Ij)");
 
1401
  for(h=0; h < nirreps; h++) {
 
1402
    dpd_file4_mat_irrep_init(&W, h);
 
1403
    dpd_file4_mat_irrep_rd(&W, h);
 
1404
    for(mn=0; mn < W.params->rowtot[h]; mn++) {
 
1405
      n = W.params->roworb[h][mn][1];
 
1406
      nsym = W.params->qsym[n];
 
1407
      N = n - occ_off[nsym];
 
1408
      for(ij=0; ij < W.params->coltot[h]; ij++) {
 
1409
        j = W.params->colorb[h][ij][1];
 
1410
        jsym = W.params->ssym[j];
 
1411
        J = j - occ_off[jsym];
 
1412
        if ((J >= (occpi[jsym] - openpi[jsym])) ||
 
1413
            (N >= (occpi[nsym] - openpi[nsym])) )
 
1414
          W.matrix[h][mn][ij] = 0.0;
 
1415
      }
 
1416
    }
 
1417
    dpd_file4_mat_irrep_wrt(&W, h);
 
1418
    dpd_file4_mat_irrep_close(&W, h);
 
1419
  }
 
1420
  dpd_file4_close(&W);
 
1421
 
 
1422
 
 
1423
  /* Purge Wmbej matrix elements */
 
1424
  dpd_file4_init(&W, CC3_HC1, C_irr, 10, 10,"HC1 WMBEJ (ME,JB)");
 
1425
  for(h=0; h < nirreps; h++) {
 
1426
    dpd_file4_mat_irrep_init(&W, h);
 
1427
    dpd_file4_mat_irrep_rd(&W, h);
 
1428
    for(me=0; me < W.params->rowtot[h]; me++) {
 
1429
      e = W.params->roworb[h][me][1];
 
1430
      esym = W.params->qsym[e];
 
1431
      E = e - vir_off[esym];
 
1432
      for(jb=0; jb < W.params->coltot[h]; jb++) {
 
1433
        b = W.params->colorb[h][jb][1];
 
1434
        bsym = W.params->ssym[b];
 
1435
        B = b - vir_off[bsym];
 
1436
        if ((E >= (virtpi[esym] - openpi[esym])) ||
 
1437
            (B >= (virtpi[bsym] - openpi[bsym])) )
 
1438
          W.matrix[h][me][jb] = 0.0;
 
1439
      }
 
1440
    }
 
1441
    dpd_file4_mat_irrep_wrt(&W, h);
 
1442
    dpd_file4_mat_irrep_close(&W, h);
 
1443
  }
 
1444
  dpd_file4_close(&W);
 
1445
 
 
1446
 
 
1447
  dpd_file4_init(&W, CC3_HC1, C_irr, 10, 10,"HC1 Wmbej (me,jb)");
 
1448
  for(h=0; h < nirreps; h++) {
 
1449
    dpd_file4_mat_irrep_init(&W, h);
 
1450
    dpd_file4_mat_irrep_rd(&W, h);
 
1451
    for(me=0; me < W.params->rowtot[h]; me++) {
 
1452
      m = W.params->roworb[h][me][0];
 
1453
      msym = W.params->psym[m];
 
1454
      M = m - occ_off[msym];
 
1455
      for(jb=0; jb < W.params->coltot[h]; jb++) {
 
1456
        j = W.params->colorb[h][jb][0];
 
1457
        jsym = W.params->rsym[j];
 
1458
        J = j - occ_off[jsym];
 
1459
        if ((M >= (occpi[msym] - openpi[msym])) ||
 
1460
            (J >= (occpi[jsym] - openpi[jsym])) )
 
1461
          W.matrix[h][me][jb] = 0.0;
 
1462
      }
 
1463
    }
 
1464
    dpd_file4_mat_irrep_wrt(&W, h);
 
1465
    dpd_file4_mat_irrep_close(&W, h);
 
1466
  }
 
1467
  dpd_file4_close(&W);
 
1468
 
 
1469
 
 
1470
  dpd_file4_init(&W, CC3_HC1, C_irr, 10, 10,"HC1 WMbEj (ME,jb)");
 
1471
  for(h=0; h < nirreps; h++) {
 
1472
    dpd_file4_mat_irrep_init(&W, h);
 
1473
    dpd_file4_mat_irrep_rd(&W, h);
 
1474
    for(me=0; me < W.params->rowtot[h]; me++) {
 
1475
      e = W.params->roworb[h][me][1];
 
1476
      esym = W.params->qsym[e];
 
1477
      E = e - vir_off[esym];
 
1478
      for(jb=0; jb < W.params->coltot[h]; jb++) {
 
1479
        j = W.params->colorb[h][jb][0];
 
1480
        jsym = W.params->rsym[j];
 
1481
        J = j - occ_off[jsym];
 
1482
        if ((E >= (virtpi[esym] - openpi[esym])) ||
 
1483
            (J >= (occpi[jsym] - openpi[jsym])) )
 
1484
          W.matrix[h][me][jb] = 0.0;
 
1485
      }
 
1486
    }
 
1487
    dpd_file4_mat_irrep_wrt(&W, h);
 
1488
    dpd_file4_mat_irrep_close(&W, h);
 
1489
  }
 
1490
  dpd_file4_close(&W);
 
1491
 
 
1492
 
 
1493
  dpd_file4_init(&W, CC3_HC1, C_irr, 10, 10,"HC1 WmBeJ (me,JB)");
 
1494
  for(h=0; h < nirreps; h++) {
 
1495
    dpd_file4_mat_irrep_init(&W, h);
 
1496
    dpd_file4_mat_irrep_rd(&W, h);
 
1497
    for(me=0; me < W.params->rowtot[h]; me++) {
 
1498
      m = W.params->roworb[h][me][0];
 
1499
      msym = W.params->psym[m];
 
1500
      M = m - occ_off[msym];
 
1501
      for(jb=0; jb < W.params->coltot[h]; jb++) {
 
1502
        b = W.params->colorb[h][jb][1];
 
1503
        bsym = W.params->ssym[b];
 
1504
        B = b - vir_off[bsym];
 
1505
        if ((M >= (occpi[msym] - openpi[msym])) ||
 
1506
            (B >= (virtpi[bsym] - openpi[bsym])) )
 
1507
          W.matrix[h][me][jb] = 0.0;
 
1508
      }
 
1509
    }
 
1510
    dpd_file4_mat_irrep_wrt(&W, h);
 
1511
    dpd_file4_mat_irrep_close(&W, h);
 
1512
  }
 
1513
  dpd_file4_close(&W);
 
1514
 
 
1515
 
 
1516
  dpd_file4_init(&W, CC3_HC1, C_irr, 10, 10,"HC1 WmBEj (mE,jB)");
 
1517
  for(h=0; h < nirreps; h++) {
 
1518
    dpd_file4_mat_irrep_init(&W, h);
 
1519
    dpd_file4_mat_irrep_rd(&W, h);
 
1520
    for(me=0; me < W.params->rowtot[h]; me++) {
 
1521
      m = W.params->roworb[h][me][0];
 
1522
      e = W.params->roworb[h][me][1];
 
1523
      msym = W.params->psym[m];
 
1524
      esym = W.params->qsym[e];
 
1525
      M = m - occ_off[msym];
 
1526
      E = e - vir_off[esym];
 
1527
      for(jb=0; jb < W.params->coltot[h]; jb++) {
 
1528
        j = W.params->colorb[h][jb][0];
 
1529
        b = W.params->colorb[h][jb][1];
 
1530
        jsym = W.params->rsym[j];
 
1531
        bsym = W.params->ssym[b];
 
1532
        J = j - occ_off[jsym];
 
1533
        B = b - vir_off[bsym];
 
1534
        if ((M >= (occpi[msym] - openpi[msym])) ||
 
1535
            (E >= (virtpi[esym] - openpi[esym])) ||
 
1536
            (J >= (occpi[jsym] - openpi[jsym])) ||
 
1537
            (B >= (virtpi[bsym] - openpi[bsym])) )
 
1538
          W.matrix[h][me][jb] = 0.0;
 
1539
      }
 
1540
    }
 
1541
    dpd_file4_mat_irrep_wrt(&W, h);
 
1542
    dpd_file4_mat_irrep_close(&W, h);
 
1543
  }
 
1544
  dpd_file4_close(&W);
 
1545
 
 
1546
 
 
1547
  /* Purge Wamef matrix elements */
 
1548
  dpd_file4_init(&W, CC3_HC1, C_irr, 11, 7,"HC1 WAMEF (AM,E>F)");
 
1549
  for(h=0; h < nirreps; h++) {
 
1550
    dpd_file4_mat_irrep_init(&W, h);
 
1551
    dpd_file4_mat_irrep_rd(&W, h);
 
1552
    for(ma=0; ma < W.params->rowtot[h]; ma++) {
 
1553
      a = W.params->roworb[h][ma][0];
 
1554
      asym = W.params->psym[a];
 
1555
      A = a - vir_off[asym];
 
1556
      for(ef=0; ef< W.params->coltot[h]; ef++) {
 
1557
        e = W.params->colorb[h][ef][0];
 
1558
        f = W.params->colorb[h][ef][1];
 
1559
        esym = W.params->rsym[e];
 
1560
        fsym = W.params->ssym[f];
 
1561
        E = e - vir_off[esym];
 
1562
        F = f - vir_off[fsym];
 
1563
        if ((A >= (virtpi[asym] - openpi[asym])) ||
 
1564
            (E >= (virtpi[esym] - openpi[esym])) ||
 
1565
            (F >= (virtpi[fsym] - openpi[fsym])) )
 
1566
          W.matrix[h][ma][ef] = 0.0;
 
1567
      }
 
1568
    }
 
1569
    dpd_file4_mat_irrep_wrt(&W, h);
 
1570
    dpd_file4_mat_irrep_close(&W, h);
 
1571
  }
 
1572
  dpd_file4_close(&W);
 
1573
 
 
1574
  dpd_file4_init(&W, CC3_HC1, C_irr, 11, 7,"HC1 Wamef (am,e>f)");
 
1575
  for(h=0; h < nirreps; h++) {
 
1576
    dpd_file4_mat_irrep_init(&W, h);
 
1577
    dpd_file4_mat_irrep_rd(&W, h);
 
1578
    for(ma=0; ma < W.params->rowtot[h]; ma++) {
 
1579
      m = W.params->roworb[h][ma][1];
 
1580
      msym = W.params->qsym[m];
 
1581
      M = m - occ_off[msym];
 
1582
      for(ef=0; ef< W.params->coltot[h]; ef++) {
 
1583
        if (M >=  (occpi[msym] - openpi[msym]))
 
1584
          W.matrix[h][ma][ef] = 0.0;
 
1585
      }
 
1586
    }
 
1587
    dpd_file4_mat_irrep_wrt(&W, h);
 
1588
    dpd_file4_mat_irrep_close(&W, h);
 
1589
  }
 
1590
  dpd_file4_close(&W);
 
1591
 
 
1592
  dpd_file4_init(&W, CC3_HC1, C_irr, 11, 5,"HC1 WAmEf (Am,Ef)");
 
1593
  for(h=0; h < nirreps; h++) {
 
1594
    dpd_file4_mat_irrep_init(&W, h);
 
1595
    dpd_file4_mat_irrep_rd(&W, h);
 
1596
    for(ma=0; ma < W.params->rowtot[h]; ma++) {
 
1597
      a = W.params->roworb[h][ma][0];
 
1598
      m = W.params->roworb[h][ma][1];
 
1599
      asym = W.params->psym[a];
 
1600
      msym = W.params->qsym[m];
 
1601
      M = m - occ_off[msym];
 
1602
      A = a - vir_off[asym];
 
1603
      for(ef=0; ef< W.params->coltot[h]; ef++) {
 
1604
        e = W.params->colorb[h][ef][0];
 
1605
        esym = W.params->rsym[e];
 
1606
        E = e - vir_off[esym];
 
1607
        if ((A >= (virtpi[asym] - openpi[asym])) ||
 
1608
            (M >=  (occpi[msym] - openpi[msym])) ||
 
1609
            (E >= (virtpi[esym] - openpi[esym])) )
 
1610
          W.matrix[h][ma][ef] = 0.0;
 
1611
      }
 
1612
    }
 
1613
    dpd_file4_mat_irrep_wrt(&W, h);
 
1614
    dpd_file4_mat_irrep_close(&W, h);
 
1615
  }
 
1616
  dpd_file4_close(&W);
 
1617
 
 
1618
  dpd_file4_init(&W, CC3_HC1, C_irr, 11, 5,"HC1 WaMeF (aM,eF)");
 
1619
  for(h=0; h < nirreps; h++) {
 
1620
    dpd_file4_mat_irrep_init(&W, h);
 
1621
    dpd_file4_mat_irrep_rd(&W, h);
 
1622
    for(ma=0; ma < W.params->rowtot[h]; ma++) {
 
1623
      for(ef=0; ef< W.params->coltot[h]; ef++) {
 
1624
        f = W.params->colorb[h][ef][1];
 
1625
        fsym = W.params->ssym[f];
 
1626
        F = f - vir_off[fsym];
 
1627
        if (F >= (virtpi[fsym] - openpi[fsym]))
 
1628
          W.matrix[h][ma][ef] = 0.0;
 
1629
      }
 
1630
    }
 
1631
    dpd_file4_mat_irrep_wrt(&W, h);
 
1632
    dpd_file4_mat_irrep_close(&W, h);
 
1633
  }
 
1634
  dpd_file4_close(&W);
 
1635
 
 
1636
 
 
1637
 
 
1638
  return;
 
1639
}
 
1640
 
 
1641
/* Purge Wmnie matrix elements */
 
1642
void purge_Wmnie(int C_irr) {
 
1643
  dpdfile4 W;
 
1644
  int *occpi, *virtpi;
 
1645
  int h, a, b, e, f, i, j, m, n;
 
1646
  int    A, B, E, F, I, J, M, N;
 
1647
  int mn, ei, ma, ef, me, jb, mb, ij, ab;
 
1648
  int asym, bsym, esym, fsym, isym, jsym, msym, nsym;
 
1649
  int *occ_off, *vir_off;
 
1650
  int *occ_sym, *vir_sym;
 
1651
  int *openpi, nirreps;
 
1652
 
 
1653
  nirreps = moinfo.nirreps;
 
1654
  occpi = moinfo.occpi; virtpi = moinfo.virtpi;
 
1655
  occ_off = moinfo.occ_off; vir_off = moinfo.vir_off;
 
1656
  occ_sym = moinfo.occ_sym; vir_sym = moinfo.vir_sym;
 
1657
  openpi = moinfo.openpi;
 
1658
 
 
1659
  dpd_file4_init(&W, CC3_HC1, C_irr, 0, 11,"HC1 WMnIe (Mn,eI)");
 
1660
  for(h=0; h < nirreps; h++) {
 
1661
    dpd_file4_mat_irrep_init(&W, h);
 
1662
    dpd_file4_mat_irrep_rd(&W, h);
 
1663
    for(mn=0; mn<W.params->rowtot[h]; mn++) {
 
1664
      n = W.params->roworb[h][mn][1];
 
1665
      nsym = W.params->qsym[n];
 
1666
      N = n - occ_off[nsym];
 
1667
      for(ei=0; ei<W.params->coltot[h]; ei++) {
 
1668
        if (N >= (occpi[nsym] - openpi[nsym]))
 
1669
          W.matrix[h][mn][ei] = 0.0;
 
1670
      }
 
1671
    }
 
1672
    dpd_file4_mat_irrep_wrt(&W, h);
 
1673
    dpd_file4_mat_irrep_close(&W, h);
 
1674
  }
 
1675
 
 
1676
  dpd_file4_init(&W, CC3_HC1, C_irr, 2, 11, "HC1 WMNIE (M>N,EI)");
 
1677
  for(h=0; h < W.params->nirreps; h++) {
 
1678
    dpd_file4_mat_irrep_init(&W, h);
 
1679
    dpd_file4_mat_irrep_rd(&W, h);
 
1680
    for(mn=0; mn<W.params->rowtot[h]; mn++) {
 
1681
      for(ei=0; ei<W.params->coltot[h]; ei++) {
 
1682
        e = W.params->colorb[h][ei][0];
 
1683
        esym = W.params->rsym[e];
 
1684
        E = e - vir_off[esym];
 
1685
        if (E >= (virtpi[esym] - openpi[esym]))
 
1686
          W.matrix[h][mn][ei] = 0.0;
 
1687
      }
 
1688
    }
 
1689
    dpd_file4_mat_irrep_wrt(&W, h);
 
1690
    dpd_file4_mat_irrep_close(&W, h);
 
1691
  }
 
1692
  dpd_file4_close(&W);
 
1693
 
 
1694
  dpd_file4_init(&W, CC3_HC1, C_irr, 2, 11,"HC1 Wmnie (m>n,ei)");
 
1695
  for(h=0; h < nirreps; h++) {
 
1696
    dpd_file4_mat_irrep_init(&W, h);
 
1697
    dpd_file4_mat_irrep_rd(&W, h);
 
1698
    for(mn=0; mn<W.params->rowtot[h]; mn++) {
 
1699
      m = W.params->roworb[h][mn][0];
 
1700
      n = W.params->roworb[h][mn][1];
 
1701
      msym = W.params->psym[m];
 
1702
      nsym = W.params->qsym[n];
 
1703
      M = m - occ_off[msym];
 
1704
      N = n - occ_off[nsym];
 
1705
      for(ei=0; ei<W.params->coltot[h]; ei++) {
 
1706
        i = W.params->colorb[h][ei][1];
 
1707
        isym = W.params->ssym[i];
 
1708
        I = i - occ_off[isym];
 
1709
        if ((M >= (occpi[msym] - openpi[msym])) ||
 
1710
          (N >= (occpi[nsym] - openpi[nsym])) ||
 
1711
          (I >= (occpi[isym] - openpi[isym])) )
 
1712
          W.matrix[h][mn][ei] = 0.0;
 
1713
      }
 
1714
    }
 
1715
    dpd_file4_mat_irrep_wrt(&W, h);
 
1716
    dpd_file4_mat_irrep_close(&W, h);
 
1717
  }
 
1718
  dpd_file4_close(&W);
 
1719
 
 
1720
  dpd_file4_init(&W, CC3_HC1, C_irr, 0, 11,"HC1 WmNiE (mN,Ei)");
 
1721
  for(h=0; h < nirreps; h++) {
 
1722
    dpd_file4_mat_irrep_init(&W, h);
 
1723
    dpd_file4_mat_irrep_rd(&W, h);
 
1724
    for(mn=0; mn<W.params->rowtot[h]; mn++) {
 
1725
      m = W.params->roworb[h][mn][0];
 
1726
      msym = W.params->psym[m];
 
1727
      M = m - occ_off[msym];
 
1728
      for(ei=0; ei<W.params->coltot[h]; ei++) {
 
1729
        e = W.params->colorb[h][ei][0];
 
1730
        i = W.params->colorb[h][ei][1];
 
1731
        esym = W.params->rsym[e];
 
1732
        isym = W.params->ssym[i];
 
1733
        E = e - vir_off[esym];
 
1734
        I = i - occ_off[isym];
 
1735
        if ((M >= (occpi[msym] - openpi[msym])) ||
 
1736
            (E >= (virtpi[esym] - openpi[esym])) ||
 
1737
            (I >= (occpi[isym] - openpi[isym])) )
 
1738
          W.matrix[h][mn][ei] = 0.0;
 
1739
      }
 
1740
    }
 
1741
    dpd_file4_mat_irrep_wrt(&W, h);
 
1742
    dpd_file4_mat_irrep_close(&W, h);
 
1743
  }
 
1744
  dpd_file4_close(&W);
 
1745
  return;
 
1746
}
 
1747
 
 
1748
/* Purge WMBIJ matrix elements */
 
1749
void purge_Wmbij(int C_irr) {
 
1750
  dpdfile4 W;
 
1751
  int *occpi, *virtpi;
 
1752
  int h, a, b, e, f, i, j, m, n;
 
1753
  int    A, B, E, F, I, J, M, N;
 
1754
  int mn, ei, ma, ef, me, jb, mb, ij, ab;
 
1755
  int asym, bsym, esym, fsym, isym, jsym, msym, nsym;
 
1756
  int *occ_off, *vir_off;
 
1757
  int *occ_sym, *vir_sym;
 
1758
  int *openpi, nirreps;
 
1759
 
 
1760
  nirreps = moinfo.nirreps;
 
1761
  occpi = moinfo.occpi; virtpi = moinfo.virtpi;
 
1762
  occ_off = moinfo.occ_off; vir_off = moinfo.vir_off;
 
1763
  occ_sym = moinfo.occ_sym; vir_sym = moinfo.vir_sym;
 
1764
  openpi = moinfo.openpi;
 
1765
 
 
1766
  dpd_file4_init(&W, CC3_HC1, C_irr, 10, 2,"HC1 WMBIJ (MB,I>J)");
 
1767
  for(h=0; h < nirreps; h++) {
 
1768
    dpd_file4_mat_irrep_init(&W, h);
 
1769
    dpd_file4_mat_irrep_rd(&W, h);
 
1770
    for(mb=0; mb<W.params->rowtot[h]; mb++) {
 
1771
      b = W.params->roworb[h][mb][1];
 
1772
      bsym = W.params->qsym[b];
 
1773
      B = b - vir_off[bsym];
 
1774
      for(ij=0; ij<W.params->coltot[h]; ij++) {
 
1775
        if (B >= (virtpi[bsym] - openpi[bsym]))
 
1776
          W.matrix[h][mb][ij] = 0.0;
 
1777
      }
 
1778
    }
 
1779
    dpd_file4_mat_irrep_wrt(&W, h);
 
1780
    dpd_file4_mat_irrep_close(&W, h);
 
1781
  }
 
1782
  dpd_file4_close(&W);
 
1783
 
 
1784
  dpd_file4_init(&W, CC3_HC1, C_irr, 10, 2,"HC1 Wmbij (mb,i>j)");
 
1785
  for(h=0; h < nirreps; h++) {
 
1786
    dpd_file4_mat_irrep_init(&W, h);
 
1787
    dpd_file4_mat_irrep_rd(&W, h);
 
1788
    for(mb=0; mb<W.params->rowtot[h]; mb++) {
 
1789
      m = W.params->roworb[h][mb][0];
 
1790
      msym = W.params->psym[m];
 
1791
      M = m - occ_off[msym];
 
1792
      for(ij=0; ij<W.params->coltot[h]; ij++) {
 
1793
        i = W.params->colorb[h][ij][0];
 
1794
        j = W.params->colorb[h][ij][1];
 
1795
        isym = W.params->rsym[i];
 
1796
        jsym = W.params->ssym[j];
 
1797
        I = i - occ_off[isym];
 
1798
        J = j - occ_off[jsym];
 
1799
        if ((M >= (occpi[msym] - openpi[msym])) ||
 
1800
            (I >= (occpi[isym] - openpi[isym])) ||
 
1801
            (J >= (occpi[jsym] - openpi[jsym])) )
 
1802
          W.matrix[h][mb][ij] = 0.0;
 
1803
      }
 
1804
    }
 
1805
    dpd_file4_mat_irrep_wrt(&W, h);
 
1806
    dpd_file4_mat_irrep_close(&W, h);
 
1807
  }
 
1808
  dpd_file4_close(&W);
 
1809
 
 
1810
  dpd_file4_init(&W, CC3_HC1, C_irr, 10, 0,"HC1 WMbIj (Mb,Ij)");
 
1811
  for(h=0; h < nirreps; h++) {
 
1812
    dpd_file4_mat_irrep_init(&W, h);
 
1813
    dpd_file4_mat_irrep_rd(&W, h);
 
1814
    for(mb=0; mb<W.params->rowtot[h]; mb++) {
 
1815
      for(ij=0; ij<W.params->coltot[h]; ij++) {
 
1816
        j = W.params->colorb[h][ij][1];
 
1817
        jsym = W.params->ssym[j];
 
1818
        J = j - occ_off[jsym];
 
1819
        if (J >= (occpi[jsym] - openpi[jsym]))
 
1820
          W.matrix[h][mb][ij] = 0.0;
 
1821
      }
 
1822
    }
 
1823
    dpd_file4_mat_irrep_wrt(&W, h);
 
1824
    dpd_file4_mat_irrep_close(&W, h);
 
1825
  }
 
1826
  dpd_file4_close(&W);
 
1827
 
 
1828
  dpd_file4_init(&W, CC3_HC1, C_irr, 10, 0,"HC1 WmBiJ (mB,iJ)");
 
1829
  for(h=0; h < nirreps; h++) {
 
1830
    dpd_file4_mat_irrep_init(&W, h);
 
1831
    dpd_file4_mat_irrep_rd(&W, h);
 
1832
    for(mb=0; mb<W.params->rowtot[h]; mb++) {
 
1833
      m = W.params->roworb[h][mb][0];
 
1834
      b = W.params->roworb[h][mb][1];
 
1835
      msym = W.params->psym[m];
 
1836
      bsym = W.params->qsym[b];
 
1837
      M = m - occ_off[msym];
 
1838
      B = b - vir_off[bsym];
 
1839
      for(ij=0; ij<W.params->coltot[h]; ij++) {
 
1840
        i = W.params->colorb[h][ij][0];
 
1841
        isym = W.params->rsym[i];
 
1842
        I = i - occ_off[isym];
 
1843
        if ((M >= (occpi[msym] - openpi[msym])) ||
 
1844
            (B >= (virtpi[bsym] - openpi[bsym])) ||
 
1845
            (I >= (occpi[isym] - openpi[isym])) )
 
1846
          W.matrix[h][mb][ij] = 0.0;
 
1847
      }
 
1848
    }
 
1849
    dpd_file4_mat_irrep_wrt(&W, h);
 
1850
    dpd_file4_mat_irrep_close(&W, h);
 
1851
  }
 
1852
  dpd_file4_close(&W);
 
1853
}
 
1854
 
 
1855
 
 
1856
/* Purge Wabei matrix elements */
 
1857
void purge_Wabei(int C_irr) {
 
1858
 
 
1859
  dpdfile4 W;
 
1860
  int *occpi, *virtpi;
 
1861
  int h, a, b, e, f, i, j, m, n;
 
1862
  int    A, B, E, F, I, J, M, N;
 
1863
  int mn, ei, ma, ef, me, jb, mb, ij, ab;
 
1864
  int asym, bsym, esym, fsym, isym, jsym, msym, nsym;
 
1865
  int *occ_off, *vir_off;
 
1866
  int *occ_sym, *vir_sym;
 
1867
  int *openpi, nirreps;
 
1868
 
 
1869
  nirreps = moinfo.nirreps;
 
1870
  occpi = moinfo.occpi; virtpi = moinfo.virtpi;
 
1871
  occ_off = moinfo.occ_off; vir_off = moinfo.vir_off;
 
1872
  occ_sym = moinfo.occ_sym; vir_sym = moinfo.vir_sym;
 
1873
  openpi = moinfo.openpi;
 
1874
 
 
1875
  dpd_file4_init(&W, CC_TMP0, C_irr, 11, 7,"WABEI (EI,A>B)");
 
1876
  for(h=0; h < nirreps; h++) {
 
1877
    dpd_file4_mat_irrep_init(&W, h);
 
1878
    dpd_file4_mat_irrep_rd(&W, h);
 
1879
    for(ei=0; ei<W.params->rowtot[h]; ei++) {
 
1880
      e = W.params->roworb[h][ei][0];
 
1881
      esym = W.params->psym[e];
 
1882
      E = e - vir_off[esym];
 
1883
      for(ab=0; ab<W.params->coltot[h]; ab++) {
 
1884
        a = W.params->colorb[h][ab][0];
 
1885
        b = W.params->colorb[h][ab][1];
 
1886
        asym = W.params->rsym[a];
 
1887
        bsym = W.params->ssym[b];
 
1888
        A = a - vir_off[asym];
 
1889
        B = b - vir_off[bsym];
 
1890
        if ((E >= (virtpi[esym] - openpi[esym])) ||
 
1891
            (A >= (virtpi[asym] - openpi[asym])) ||
 
1892
            (B >= (virtpi[bsym] - openpi[bsym])) )
 
1893
          W.matrix[h][ei][ab] = 0.0;
 
1894
      }
 
1895
    }
 
1896
    dpd_file4_mat_irrep_wrt(&W, h);
 
1897
    dpd_file4_mat_irrep_close(&W, h);
 
1898
  }
 
1899
  dpd_file4_close(&W);
 
1900
 
 
1901
  dpd_file4_init(&W, CC_TMP0, C_irr, 11, 7,"Wabei (ei,a>b)");
 
1902
  for(h=0; h < nirreps; h++) {
 
1903
    dpd_file4_mat_irrep_init(&W, h);
 
1904
    dpd_file4_mat_irrep_rd(&W, h);
 
1905
    for(ei=0; ei<W.params->rowtot[h]; ei++) {
 
1906
      i = W.params->roworb[h][ei][1];
 
1907
      isym = W.params->qsym[i];
 
1908
      I = i - occ_off[isym];
 
1909
      for(ab=0; ab<W.params->coltot[h]; ab++) {
 
1910
        if (I >= (occpi[isym] - openpi[isym]))
 
1911
          W.matrix[h][ei][ab] = 0.0;
 
1912
      }
 
1913
    }
 
1914
    dpd_file4_mat_irrep_wrt(&W, h);
 
1915
    dpd_file4_mat_irrep_close(&W, h);
 
1916
  }
 
1917
  dpd_file4_close(&W);
 
1918
 
 
1919
  dpd_file4_init(&W, CC_TMP0, C_irr, 11, 5,"WAbEi (Ei,Ab)");
 
1920
  for(h=0; h < nirreps; h++) {
 
1921
    dpd_file4_mat_irrep_init(&W, h);
 
1922
    dpd_file4_mat_irrep_rd(&W, h);
 
1923
    for(ei=0; ei<W.params->rowtot[h]; ei++) {
 
1924
      e = W.params->roworb[h][ei][0];
 
1925
      i = W.params->roworb[h][ei][1];
 
1926
      esym = W.params->psym[e];
 
1927
      isym = W.params->qsym[i];
 
1928
      E = e - vir_off[esym];
 
1929
      I = i - occ_off[isym];
 
1930
      for(ab=0; ab<W.params->coltot[h]; ab++) {
 
1931
        a = W.params->colorb[h][ab][0];
 
1932
        asym = W.params->rsym[a];
 
1933
        bsym = W.params->ssym[b];
 
1934
        A = a - vir_off[asym];
 
1935
        if ((E >= (virtpi[esym] - openpi[esym])) ||
 
1936
            (I >= (occpi[isym] - openpi[isym])) ||
 
1937
            (A >= (virtpi[asym] - openpi[asym])) )
 
1938
          W.matrix[h][ei][ab] = 0.0;
 
1939
      }
 
1940
    }
 
1941
    dpd_file4_mat_irrep_wrt(&W, h);
 
1942
    dpd_file4_mat_irrep_close(&W, h);
 
1943
  }
 
1944
  dpd_file4_close(&W);
 
1945
 
 
1946
  dpd_file4_init(&W, CC_TMP0, C_irr, 11, 5,"WaBeI (eI,aB)");
 
1947
  for(h=0; h < nirreps; h++) {
 
1948
    dpd_file4_mat_irrep_init(&W, h);
 
1949
    dpd_file4_mat_irrep_rd(&W, h);
 
1950
    for(ei=0; ei<W.params->rowtot[h]; ei++) {
 
1951
      for(ab=0; ab<W.params->coltot[h]; ab++) {
 
1952
        b = W.params->colorb[h][ab][1];
 
1953
        bsym = W.params->ssym[b];
 
1954
        B = b - vir_off[bsym];
 
1955
        if (B >= (virtpi[bsym] - openpi[bsym]))
 
1956
          W.matrix[h][ei][ab] = 0.0;
 
1957
      }
 
1958
    }
 
1959
    dpd_file4_mat_irrep_wrt(&W, h);
 
1960
    dpd_file4_mat_irrep_close(&W, h);
 
1961
  }
 
1962
  dpd_file4_close(&W);
 
1963
}