← Back to branch summary

~ubuntu-branches/ubuntu/karmic/psicode/karmic

« back to all changes in this revision

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

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2009-02-23 00:12:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20090223001202-rutldoy3dimfpesc
Tags: 3.4.0-1
* New upstream release.

[ Michael Banck ]
* debian/patches/01_DESTDIR.dpatch: Refreshed.
* debian/patches/02_FHS.dpatch: Removed, applied upstream.
* debian/patches/03_debian_docdir: Likewise.
* debian/patches/04_man.dpatch: Likewise.
* debian/patches/06_466828_fix_gcc_43_ftbfs.dpatch: Likewise.
* debian/patches/07_464867_move_executables: Fixed and refreshed.
* debian/patches/00list: Adjusted.
* debian/control: Improved description.
* debian/patches-held: Removed.
* debian/rules (install/psi3): Do not ship the ruby bindings for now.

[ Daniel Leidert ]
* debian/rules: Fix txtdir via DEB_MAKE_INSTALL_TARGET.
* debian/patches/01_DESTDIR.dpatch: Refreshed.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/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
 
}