← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cchbar/cc2_Wabei.cc

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

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

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

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/cchbar/cc2_Wabei.cc
 
1
/*! \file
 
2
    \ingroup CCHBAR
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#include <cstring>
 
6
#include <libdpd/dpd.h>
 
7
#include <libqt/qt.h>
 
8
#include "MOInfo.h"
 
9
#include "Params.h"
 
10
#define EXTERN
 
11
#include "globals.h"
 
12
 
 
13
namespace psi { namespace cchbar {
 
14
 
 
15
/* cc2_Wabei(): Compute the Wabei matrix from CC2 theory, which is
 
16
** given in spin orbitals as:
 
17
**
 
18
** Wabei = <ab||ei> - P(ab) t_m^a <mb||ei> + t_i^f <ab||ef> 
 
19
**         - P(ab) t_i^f t_m^b <am||ef> + t_m^a t_n^b <mn||ei>
 
20
**         + t_m^a t_i^f t_n^b <mn||ef>
 
21
**
 
22
** The basic strategy for this code is to generate two intermediate
 
23
** quantities, Z1(Ab,EI) and Z2(Ei,Ab), which are summed in the final
 
24
** step to give the complete W(Ei,Ab) intermediate.  This is sorted
 
25
** to W(iE,bA) storage for use in the triples equations.
 
26
**
 
27
** TDC, Feb 2004
 
28
*/
 
29
 
 
30
void purge_cc2_Wabei(void);
 
31
 
 
32
void cc2_Wabei_build(void)
 
33
{
 
34
  int omit = 0;
 
35
  int e, E;
 
36
  int Gef, Gab, Gei, Ge, Gf, Gi;
 
37
  int nrows, ncols, nlinks;
 
38
  dpdfile2 t1, tIA, tia;
 
39
  dpdbuf4 Z, Z1, Z2, Z3;
 
40
  dpdbuf4 B, C, D, F, W;
 
41
 
 
42
  timer_on("F->Wabei");
 
43
  if(params.ref == 0) { /** RHF **/
 
44
 
 
45
    dpd_buf4_init(&F, CC_FINTS, 0, 11, 5, 11, 5, 0, "F <ai|bc>");
 
46
    dpd_buf4_copy(&F, CC_TMP0, "CC2 WAbEi (Ei,Ab)");
 
47
    dpd_buf4_close(&F);
 
48
 
 
49
  }
 
50
 
 
51
  else if (params.ref == 1) { /* ROHF */
 
52
    dpd_file2_init(&tIA, CC_OEI, 0, 0, 1, "tIA");
 
53
    dpd_file2_init(&tia, CC_OEI, 0, 0, 1, "tia");
 
54
 
 
55
    /* W(A>B,EI) <--- <AB||EI> */
 
56
    /* W(a>b,ei) <--- <ab||ei> */
 
57
    dpd_buf4_init(&F, CC_FINTS, 0, 11, 7, 11, 5, 1, "F <ai|bc>");
 
58
    dpd_buf4_sort(&F, CC_TMP2, rspq, 7, 11, "CC2 WABEI (A>B,EI)");
 
59
    dpd_buf4_sort(&F, CC_TMP2, rspq, 7, 11, "CC2 Wabei (a>b,ei)");
 
60
    dpd_buf4_close(&F);
 
61
 
 
62
    /* W(Ab,Ei) <--- <Ab|Ei> */
 
63
    /* W(aB,eI) <--- <aB|eI> */
 
64
    dpd_buf4_init(&F, CC_FINTS, 0, 11, 5, 11, 5, 0, "F <ai|bc>");
 
65
    dpd_buf4_sort(&F, CC_TMP2, rspq, 5, 11, "CC2 WAbEi (Ab,Ei)");
 
66
    dpd_buf4_sort(&F, CC_TMP2, rspq, 5, 11, "CC2 WaBeI (aB,eI)");
 
67
    dpd_buf4_close(&F);
 
68
  }
 
69
 
 
70
  else if (params.ref == 2) { /* UHF */
 
71
    /* W(A>B,EI) <--- <AB||EI> */
 
72
    dpd_buf4_init(&F, CC_FINTS, 0, 21, 7, 21, 5, 1, "F <AI|BC>");
 
73
    dpd_buf4_sort(&F, CC_TMP0, rspq, 7, 21, "CC2 WABEI (A>B,EI)");
 
74
    dpd_buf4_close(&F);
 
75
 
 
76
    /* W(a>b,ei) <--- <ab||ei> */
 
77
    dpd_buf4_init(&F, CC_FINTS, 0, 31, 17, 31,15, 1, "F <ai|bc>");
 
78
    dpd_buf4_sort(&F, CC_TMP0, rspq, 17, 31, "CC2 Wabei (a>b,ei)");
 
79
    dpd_buf4_close(&F);
 
80
 
 
81
    /* W(Ab,Ei) <--- <Ab|Ei> */
 
82
    dpd_buf4_init(&F, CC_FINTS, 0, 28, 26, 28, 26, 0, "F <Ab|Ci>");
 
83
    dpd_buf4_copy(&F, CC_TMP0, "CC2 WAbEi (Ab,Ei)");
 
84
    dpd_buf4_close(&F);
 
85
 
 
86
    /* W(aB,eI) <--- <aB|eI> */
 
87
    dpd_buf4_init(&F, CC_FINTS, 0, 25, 29, 25, 29, 0, "F <aI|bC>");
 
88
    dpd_buf4_sort(&F, CC_TMP0, psrq, 29, 25, "CC2 WaBeI (aB,eI)");
 
89
    dpd_buf4_close(&F);
 
90
  }
 
91
  timer_off("F->Wabei");
 
92
 
 
93
  timer_on("B->Wabei");
 
94
  if(params.ref == 0) { /** RHF **/
 
95
 
 
96
    dpd_file2_init(&t1, CC_OEI, 0, 0, 1, "tIA");
 
97
    dpd_file2_mat_init(&t1);
 
98
    dpd_file2_mat_rd(&t1);
 
99
 
 
100
    /* WEbEi <-- <Ab|Ef> * t(i,f) */
 
101
 
 
102
    /* Plus Combination */
 
103
    dpd_buf4_init(&B, CC_BINTS, 0, 5, 8, 8, 8, 0, "B(+) <ab|cd> + <ab|dc>");
 
104
    dpd_buf4_init(&Z1, CC_TMP0, 0, 11, 8, 11, 8, 0, "Z1(ei,a>=b)");
 
105
 
 
106
    for(Gef=0; Gef < moinfo.nirreps; Gef++) {
 
107
      Gei = Gab = Gef; /* W and B are totally symmetric */
 
108
      for(Ge=0; Ge < moinfo.nirreps; Ge++) {
 
109
        Gf = Ge ^ Gef; Gi = Gf;  /* t1 is totally symmetric */
 
110
        B.matrix[Gef] = dpd_block_matrix(moinfo.virtpi[Gf],B.params->coltot[Gef]);
 
111
        Z1.matrix[Gef] = dpd_block_matrix(moinfo.occpi[Gi],Z1.params->coltot[Gei]);
 
112
        nrows = moinfo.occpi[Gi];
 
113
        ncols = Z1.params->coltot[Gef];
 
114
        nlinks = moinfo.virtpi[Gf];
 
115
        if(nrows && ncols && nlinks) {
 
116
          for(E=0; E < moinfo.virtpi[Ge]; E++) {
 
117
            e = moinfo.vir_off[Ge] + E;
 
118
            dpd_buf4_mat_irrep_rd_block(&B, Gef, B.row_offset[Gef][e], moinfo.virtpi[Gf]);
 
119
            C_DGEMM('n','n',nrows,ncols,nlinks,0.5,t1.matrix[Gi][0],nlinks,B.matrix[Gef][0],ncols,
 
120
                    0.0,Z1.matrix[Gei][0],ncols);
 
121
            dpd_buf4_mat_irrep_wrt_block(&Z1, Gei, Z1.row_offset[Gei][e], moinfo.occpi[Gi]);
 
122
          }
 
123
        }
 
124
        dpd_free_block(B.matrix[Gef], moinfo.virtpi[Gf], B.params->coltot[Gef]);
 
125
        dpd_free_block(Z1.matrix[Gef], moinfo.occpi[Gi], Z1.params->coltot[Gei]);
 
126
      }
 
127
    }
 
128
 
 
129
    dpd_buf4_close(&Z1);
 
130
    dpd_buf4_close(&B);
 
131
 
 
132
    /* Minus Combination */
 
133
    dpd_buf4_init(&B, CC_BINTS, 0, 5, 9, 9, 9, 0, "B(-) <ab|cd> - <ab|dc>");
 
134
    dpd_buf4_init(&Z2, CC_TMP0, 0, 11, 9, 11, 9, 0, "Z2(ei,a>=b)");
 
135
 
 
136
    for(Gef=0; Gef < moinfo.nirreps; Gef++) {
 
137
      Gei = Gab = Gef; /* W and B are totally symmetric */
 
138
      for(Ge=0; Ge < moinfo.nirreps; Ge++) {
 
139
        Gf = Ge ^ Gef; Gi = Gf;  /* t1 is totally symmetric */
 
140
        B.matrix[Gef] = dpd_block_matrix(moinfo.virtpi[Gf],B.params->coltot[Gef]);
 
141
        Z2.matrix[Gef] = dpd_block_matrix(moinfo.occpi[Gi],Z2.params->coltot[Gei]);
 
142
        nrows = moinfo.occpi[Gi];
 
143
        ncols = Z2.params->coltot[Gef];
 
144
        nlinks = moinfo.virtpi[Gf];
 
145
        if(nrows && ncols && nlinks) {
 
146
          for(E=0; E < moinfo.virtpi[Ge]; E++) {
 
147
            e = moinfo.vir_off[Ge] + E;
 
148
            dpd_buf4_mat_irrep_rd_block(&B, Gef, B.row_offset[Gef][e], moinfo.virtpi[Gf]);
 
149
            C_DGEMM('n','n',nrows,ncols,nlinks,0.5,t1.matrix[Gi][0],nlinks,B.matrix[Gef][0],ncols,
 
150
                    0.0,Z2.matrix[Gei][0],ncols);
 
151
            dpd_buf4_mat_irrep_wrt_block(&Z2, Gei, Z2.row_offset[Gei][e], moinfo.occpi[Gi]);
 
152
          }
 
153
        }
 
154
        dpd_free_block(B.matrix[Gef], moinfo.virtpi[Gf], B.params->coltot[Gef]);
 
155
        dpd_free_block(Z2.matrix[Gef], moinfo.occpi[Gi], Z2.params->coltot[Gei]);
 
156
      }
 
157
    }
 
158
 
 
159
    dpd_buf4_close(&Z2);
 
160
    dpd_buf4_close(&B);
 
161
 
 
162
    dpd_file2_mat_close(&t1);
 
163
    dpd_file2_close(&t1);
 
164
 
 
165
    dpd_buf4_init(&W, CC_TMP0, 0, 11, 5, 11, 5, 0, "CC2 WAbEi (Ei,Ab)");
 
166
    dpd_buf4_init(&Z1, CC_TMP0, 0, 11, 5, 11, 8, 0, "Z1(ei,a>=b)");
 
167
    dpd_buf4_axpy(&Z1, &W, 1);
 
168
    dpd_buf4_close(&Z1);
 
169
    dpd_buf4_init(&Z2, CC_TMP0, 0, 11, 5, 11, 9, 0, "Z2(ei,a>=b)");
 
170
    dpd_buf4_axpy(&Z2, &W, 1);
 
171
    dpd_buf4_close(&Z2);
 
172
    dpd_buf4_close(&W);
 
173
  }
 
174
  else if (params.ref == 1) { /* ROHF */
 
175
 
 
176
    /** W(A>B,EI) <--- <AB||EF> * t1[I][F] **/
 
177
    /** W(a>b,ei) <--- <ab||ef> * t1[i][f] **/
 
178
    dpd_buf4_init(&W, CC_TMP2, 0, 7, 11, 7, 11, 0, "CC2 WABEI (A>B,EI)");
 
179
    dpd_buf4_init(&B, CC_BINTS, 0, 7, 5, 5, 5, 1, "B <ab|cd>");
 
180
    dpd_contract424(&B, &tIA, &W, 3, 1, 0, 0.5, 1);
 
181
    dpd_buf4_close(&W);
 
182
 
 
183
    dpd_buf4_init(&W, CC_TMP2, 0, 7, 11, 7, 11, 0, "CC2 Wabei (a>b,ei)");
 
184
    dpd_contract424(&B, &tia, &W, 3, 1, 0, 0.5, 1);
 
185
    dpd_buf4_close(&B);
 
186
    dpd_buf4_close(&W);
 
187
 
 
188
    /** W(Ab,Ei) <--- <Ab|Ef> * t1[i][f] **/
 
189
    dpd_buf4_init(&W, CC_TMP2, 0, 5, 11, 5, 11, 0, "CC2 WAbEi (Ab,Ei)");
 
190
    dpd_buf4_init(&B, CC_BINTS, 0, 5, 5, 5, 5, 0, "B <ab|cd>");
 
191
    dpd_contract424(&B, &tia, &W, 3, 1, 0, 0.5, 1);
 
192
    dpd_buf4_close(&B);
 
193
    dpd_buf4_close(&W);
 
194
 
 
195
    /** W(aB,eI) <--- t1[I][F] * <aB|eF>  **/
 
196
    dpd_buf4_init(&W, CC_TMP2, 0, 5, 11, 5, 11, 0, "CC2 WaBeI (aB,eI)");
 
197
    dpd_buf4_init(&B, CC_BINTS, 0, 5, 5, 5, 5, 0, "B <ab|cd>");
 
198
    dpd_contract424(&B, &tIA, &W, 3, 1, 0, 0.5, 1.0);
 
199
    dpd_buf4_close(&B);
 
200
    dpd_buf4_close(&W);
 
201
 
 
202
    dpd_file2_close(&tIA);
 
203
    dpd_file2_close(&tia);
 
204
  }
 
205
 
 
206
  else if (params.ref == 2) { /* UHF */
 
207
    dpd_file2_init(&tIA, CC_OEI, 0, 0, 1, "tIA");
 
208
    dpd_file2_init(&tia, CC_OEI, 0, 2, 3, "tia");
 
209
 
 
210
    /** W(A>B,EI) <--- <AB||EF> * t1[I][F] **/
 
211
    dpd_buf4_init(&W, CC_TMP0, 0, 7, 21, 7, 21, 0, "CC2 WABEI (A>B,EI)");
 
212
    dpd_buf4_init(&B, CC_BINTS, 0, 7, 5, 5, 5, 1, "B <AB|CD>");
 
213
    dpd_contract424(&B, &tIA, &W, 3, 1, 0, 1, 1);
 
214
    dpd_buf4_close(&B);
 
215
    dpd_buf4_close(&W);
 
216
 
 
217
    /** W(a>b,ei) <--- <ab||ef> * t1[i][f] **/
 
218
    dpd_buf4_init(&W, CC_TMP0, 0, 17, 31, 17, 31, 0, "CC2 Wabei (a>b,ei)");
 
219
    dpd_buf4_init(&B, CC_BINTS, 0, 17, 15, 15, 15, 1, "B <ab|cd>");
 
220
    dpd_contract424(&B, &tia, &W, 3, 1, 0, 1, 1);
 
221
    dpd_buf4_close(&B);
 
222
    dpd_buf4_close(&W);
 
223
 
 
224
    /** W(Ab,Ei) <--- <Ab|Ef> * t1[i][f] **/
 
225
    dpd_buf4_init(&W, CC_TMP0, 0, 28, 26, 28, 26, 0, "CC2 WAbEi (Ab,Ei)");
 
226
    dpd_buf4_init(&B, CC_BINTS, 0, 28, 28, 28, 28, 0, "B <Ab|Cd>");
 
227
    dpd_contract424(&B, &tia, &W, 3, 1, 0, 1, 1);
 
228
    dpd_buf4_close(&B);
 
229
    dpd_buf4_close(&W);
 
230
 
 
231
    /** W(aB,eI) <--- t1[I][F] * <aB|eF>  **/
 
232
    dpd_buf4_init(&Z, CC_TMP0, 0, 28, 24, 28, 24, 0, "CC2 ZBaIe (Ba,Ie)");
 
233
    dpd_buf4_init(&B, CC_BINTS, 0, 28, 28, 28, 28, 0, "B <Ab|Cd>");
 
234
    dpd_contract244(&tIA, &B, &Z, 1, 2, 1, 1, 0);
 
235
    dpd_buf4_close(&B);
 
236
    dpd_buf4_sort_axpy(&Z, CC_TMP0, qpsr, 29, 25, "CC2 WaBeI (aB,eI)", 1);
 
237
    dpd_buf4_close(&Z);
 
238
 
 
239
    dpd_file2_close(&tIA);
 
240
    dpd_file2_close(&tia);
 
241
  }
 
242
  timer_off("B->Wabei");
 
243
 
 
244
  if(params.ref == 0) { /** RHF **/
 
245
    dpd_file2_init(&t1, CC_OEI, 0, 0, 1, "tIA");
 
246
 
 
247
    dpd_buf4_init(&W, CC_TMP0, 0, 11, 5, 11, 5, 0, "CC2 WAbEi (Ei,Ab)");
 
248
    dpd_buf4_init(&Z, CC_TMP0, 0, 11, 10, 11, 10, 0, "CC2 ZMbEj (Ej,Mb)");
 
249
    dpd_contract244(&t1, &Z, &W, 0, 2, 1, -1, 1);
 
250
    dpd_buf4_close(&Z);
 
251
    dpd_buf4_close(&W);
 
252
 
 
253
    dpd_buf4_init(&W, CC_TMP0, 0, 11, 5, 11, 5, 0, "CC2 WAbEi (Ei,Ab)");
 
254
    dpd_buf4_init(&Z, CC_TMP0, 0, 11, 11, 11, 11, 0, "CC2 ZMbeJ (eJ,bM)");
 
255
    dpd_contract424(&Z, &t1, &W, 3, 0, 0, -1, 1);
 
256
    dpd_buf4_close(&Z);
 
257
    dpd_buf4_close(&W);
 
258
 
 
259
    dpd_file2_close(&t1);
 
260
  }
 
261
  else if (params.ref == 2) { /* UHF */
 
262
    dpd_file2_init(&tIA, CC_OEI, 0, 0, 1, "tIA");
 
263
    dpd_file2_init(&tia, CC_OEI, 0, 2, 3, "tia");
 
264
 
 
265
    /** AAA ***/
 
266
    dpd_buf4_init(&Z1, CC_TMP0, 0, 5, 21, 5, 21, 0, "CC2 ZABEI (AB,EI)");
 
267
    dpd_buf4_init(&Z, CC_TMP0, 0, 20, 21, 20, 21, 0, "CC2 ZMBEJ");
 
268
    dpd_contract244(&tIA, &Z, &Z1, 0, 0, 0, -1, 0);
 
269
    dpd_buf4_close(&Z);
 
270
   
 
271
    dpd_buf4_sort(&Z1, CC_TMP0, qprs, 5, 21, "CC2 ZABEI (BA,EI)");
 
272
    dpd_buf4_init(&Z, CC_TMP0, 0, 5, 21, 5, 21, 0, "CC2 ZABEI (BA,EI)");
 
273
    dpd_buf4_axpy(&Z, &Z1, -1);
 
274
    dpd_buf4_close(&Z);
 
275
    dpd_buf4_init(&W, CC_TMP0, 0, 5, 21, 7, 21, 0, "CC2 WABEI (A>B,EI)");
 
276
    dpd_buf4_axpy(&Z1, &W, 1);
 
277
    dpd_buf4_close(&Z1);
 
278
    dpd_buf4_close(&W);
 
279
 
 
280
    /** ABAB **/
 
281
    dpd_buf4_init(&Z1, CC_TMP0, 0, 28, 26, 28, 26, 0, "CC2 ZAbEi (Ab,Ei)");
 
282
    dpd_buf4_init(&Z, CC_TMP0, 0, 24, 26, 24, 26, 0, "CC2 ZMbEj");
 
283
    dpd_contract244(&tIA, &Z, &Z1, 0, 0, 0, -1, 0);
 
284
    dpd_buf4_close(&Z);
 
285
    dpd_buf4_close(&Z1);
 
286
 
 
287
    dpd_buf4_init(&Z1, CC_TMP0, 0, 29, 26, 29, 26, 0, "CC2 ZAbEi (bA,Ei)");
 
288
    dpd_buf4_init(&Z, CC_TMP0, 0, 27, 26, 27, 26, 0, "CC2 ZmBEj");
 
289
    dpd_contract244(&tia, &Z, &Z1, 0, 0, 0, -1, 0);
 
290
    dpd_buf4_close(&Z);
 
291
    dpd_buf4_sort_axpy(&Z1, CC_TMP0, qprs, 28, 26, "CC2 ZAbEi (Ab,Ei)", 1);
 
292
    dpd_buf4_close(&Z1);
 
293
 
 
294
    dpd_buf4_init(&Z, CC_TMP0, 0, 28, 26, 28, 26, 0, "CC2 ZAbEi (Ab,Ei)");
 
295
    dpd_buf4_init(&W, CC_TMP0, 0, 28, 26, 28, 26, 0, "CC2 WAbEi (Ab,Ei)");
 
296
    dpd_buf4_axpy(&Z, &W, -1);
 
297
    dpd_buf4_close(&W);
 
298
    dpd_buf4_close(&Z);
 
299
 
 
300
    /** BABA **/
 
301
    dpd_buf4_init(&Z1, CC_TMP0, 0, 29, 25, 29, 25, 0, "CC2 ZaBeI (aB,eI)");
 
302
    dpd_buf4_init(&Z, CC_TMP0, 0, 27, 25, 27, 25, 0, "CC2 ZmBeJ");
 
303
    dpd_contract244(&tia, &Z, &Z1, 0, 0, 0, -1, 0);
 
304
    dpd_buf4_close(&Z);
 
305
    dpd_buf4_close(&Z1);
 
306
 
 
307
    dpd_buf4_init(&Z1, CC_TMP0, 0, 28, 25, 28, 25, 0, "CC2 ZaBeI (Ba,eI)");
 
308
    dpd_buf4_init(&Z, CC_TMP0, 0, 24, 25, 24, 25, 0, "CC2 ZMbeJ");
 
309
    dpd_contract244(&tIA, &Z, &Z1, 0, 0, 0, -1, 0);
 
310
    dpd_buf4_close(&Z);
 
311
    dpd_buf4_sort_axpy(&Z1, CC_TMP0, qprs, 29, 25, "CC2 ZaBeI (aB,eI)", 1);
 
312
    dpd_buf4_close(&Z1);
 
313
 
 
314
    dpd_buf4_init(&Z, CC_TMP0, 0, 29, 25, 29, 25, 0, "CC2 ZaBeI (aB,eI)");
 
315
    dpd_buf4_init(&W, CC_TMP0, 0, 29, 25, 29, 25, 0, "CC2 WaBeI (aB,eI)");
 
316
    dpd_buf4_axpy(&Z, &W, -1);
 
317
    dpd_buf4_close(&W);
 
318
    dpd_buf4_close(&Z);
 
319
 
 
320
    /** BBBB **/
 
321
    dpd_buf4_init(&Z1, CC_TMP0, 0, 15, 31, 15, 31, 0, "CC2 Zabei (ab,ei)");
 
322
    dpd_buf4_init(&Z, CC_TMP0, 0, 30, 31, 30, 31, 0, "CC2 Zmbej");
 
323
    dpd_contract244(&tia, &Z, &Z1, 0, 0, 0, -1, 0);
 
324
    dpd_buf4_close(&Z);
 
325
    dpd_buf4_sort(&Z1, CC_TMP0, qprs, 15, 31, "CC2 Zabei (ba,ei)");
 
326
    dpd_buf4_init(&Z, CC_TMP0, 0, 15, 31, 15, 31, 0, "CC2 Zabei (ba,ei)");
 
327
    dpd_buf4_axpy(&Z, &Z1, -1);
 
328
    dpd_buf4_close(&Z);
 
329
    dpd_buf4_init(&W, CC_TMP0, 0, 15, 31, 17, 31, 0, "CC2 Wabei (a>b,ei)");
 
330
    dpd_buf4_axpy(&Z1, &W, 1);
 
331
    dpd_buf4_close(&Z1);
 
332
    dpd_buf4_close(&W);
 
333
 
 
334
    dpd_file2_close(&tIA);
 
335
    dpd_file2_close(&tia);
 
336
  }
 
337
 
 
338
  timer_on("Wabei_sort");
 
339
  if (params.ref == 0) { /* RHF */
 
340
 
 
341
    dpd_buf4_init(&W, CC2_HET1, 0, 5, 11, 5, 11, 0, "CC2 WAbEi");
 
342
    dpd_buf4_scm(&W, 0);
 
343
    dpd_buf4_close(&W);
 
344
 
 
345
    dpd_buf4_init(&W, CC_TMP0, 0, 11, 5, 11, 5, 0, "CC2 WAbEi (Ei,Ab)");
 
346
    dpd_buf4_sort_axpy(&W, CC2_HET1, rspq, 5, 11, "CC2 WAbEi", 1);
 
347
    dpd_buf4_close(&W);
 
348
 
 
349
    if(!strcmp(params.wfn, "EOM_CC2")) {
 
350
      dpd_buf4_init(&W, CC_TMP0, 0, 11, 5, 11, 5, 0, "CC2 WAbEi (Ei,Ab)");
 
351
      dpd_buf4_copy(&W, CC2_HET1, "CC2 WAbEi (Ei,Ab)");
 
352
      dpd_buf4_close(&W);
 
353
    }
 
354
 
 
355
  }
 
356
  if (params.ref == 1) { /* ROHF */
 
357
 
 
358
    /* sort to Wabei (ei,ab) */
 
359
    dpd_buf4_init(&W, CC_TMP2, 0, 7, 11, 7, 11, 0, "CC2 WABEI (A>B,EI)");
 
360
    dpd_buf4_sort(&W, CC2_HET1, rspq, 11, 7, "CC2 WABEI (EI,A>B)");
 
361
    dpd_buf4_close(&W);
 
362
    dpd_buf4_init(&W, CC_TMP2, 0, 7, 11, 7, 11, 0, "CC2 Wabei (a>b,ei)");
 
363
    dpd_buf4_sort(&W, CC2_HET1, rspq, 11, 7, "CC2 Wabei (ei,a>b)");
 
364
    dpd_buf4_close(&W);
 
365
    dpd_buf4_init(&W, CC_TMP2, 0, 5, 11, 5, 11, 0, "CC2 WAbEi (Ab,Ei)");
 
366
    dpd_buf4_sort(&W, CC2_HET1, rspq, 11, 5, "CC2 WAbEi (Ei,Ab)");
 
367
    dpd_buf4_close(&W);
 
368
    dpd_buf4_init(&W, CC_TMP2, 0, 5, 11, 5, 11, 0, "CC2 WaBeI (aB,eI)");
 
369
    dpd_buf4_sort(&W, CC2_HET1, rspq, 11, 5, "CC2 WaBeI (eI,aB)");
 
370
    dpd_buf4_close(&W);
 
371
 
 
372
    /* purge before final sort */
 
373
 
 
374
    /*     purge_cc2_Wabei(); */
 
375
  }
 
376
  else if (params.ref == 2) { /* UHF */
 
377
 
 
378
    /* sort to Wabei (ei,ab) */
 
379
    dpd_buf4_init(&W, CC_TMP0, 0, 7, 21, 7, 21, 0, "CC2 WABEI (A>B,EI)");
 
380
    dpd_buf4_sort(&W, CC2_HET1, rspq, 21, 7, "CC2 WABEI (EI,A>B)");
 
381
    dpd_buf4_close(&W);
 
382
    dpd_buf4_init(&W, CC_TMP0, 0, 17, 31, 17, 31, 0, "CC2 Wabei (a>b,ei)");
 
383
    dpd_buf4_sort(&W, CC2_HET1, rspq, 31, 17, "CC2 Wabei (ei,a>b)");
 
384
    dpd_buf4_close(&W);
 
385
    dpd_buf4_init(&W, CC_TMP0, 0, 28, 26, 28, 26, 0, "CC2 WAbEi (Ab,Ei)");
 
386
    dpd_buf4_sort(&W, CC2_HET1, rspq, 26, 28, "CC2 WAbEi (Ei,Ab)");
 
387
    dpd_buf4_close(&W);
 
388
    dpd_buf4_init(&W, CC_TMP0, 0, 29, 25, 29, 25, 0, "CC2 WaBeI (aB,eI)");
 
389
    dpd_buf4_sort(&W, CC2_HET1, rspq, 25, 29, "CC2 WaBeI (eI,aB)");
 
390
    dpd_buf4_close(&W);
 
391
  }
 
392
  timer_off("Wabei_sort");
 
393
 
 
394
}
 
395
 
 
396
void purge_cc2_Wabei(void) {
 
397
  dpdfile4 W;
 
398
  int *occpi, *virtpi;
 
399
  int h, a, b, e, f, i, j, m, n;
 
400
  int    A, B, E, F, I, J, M, N;
 
401
  int mn, ei, ma, ef, me, jb, mb, ij, ab;
 
402
  int asym, bsym, esym, fsym, isym, jsym, msym, nsym;
 
403
  int *occ_off, *vir_off;
 
404
  int *occ_sym, *vir_sym;
 
405
  int *openpi, nirreps;
 
406
 
 
407
  nirreps = moinfo.nirreps;
 
408
  occpi = moinfo.occpi; virtpi = moinfo.virtpi;
 
409
  occ_off = moinfo.occ_off; vir_off = moinfo.vir_off;
 
410
  occ_sym = moinfo.occ_sym; vir_sym = moinfo.vir_sym;
 
411
  openpi = moinfo.openpi;
 
412
 
 
413
  /* Purge Wabei matrix elements */
 
414
  dpd_file4_init(&W, CC_TMP2, 0, 11, 7,"CC2 WABEI (EI,A>B)");
 
415
  for(h=0; h < nirreps; h++) {
 
416
    dpd_file4_mat_irrep_init(&W, h);
 
417
    dpd_file4_mat_irrep_rd(&W, h);
 
418
    for(ei=0; ei<W.params->rowtot[h]; ei++) {
 
419
      e = W.params->roworb[h][ei][0];
 
420
      esym = W.params->psym[e];
 
421
      E = e - vir_off[esym]; 
 
422
      for(ab=0; ab<W.params->coltot[h]; ab++) {
 
423
        a = W.params->colorb[h][ab][0];
 
424
        b = W.params->colorb[h][ab][1];
 
425
        asym = W.params->rsym[a];
 
426
        bsym = W.params->ssym[b]; 
 
427
        A = a - vir_off[asym];
 
428
        B = b - vir_off[bsym];
 
429
        if ((E >= (virtpi[esym] - openpi[esym])) ||
 
430
            (A >= (virtpi[asym] - openpi[asym])) ||
 
431
            (B >= (virtpi[bsym] - openpi[bsym])) )
 
432
          W.matrix[h][ei][ab] = 0.0;
 
433
      }
 
434
    }
 
435
    dpd_file4_mat_irrep_wrt(&W, h);
 
436
    dpd_file4_mat_irrep_close(&W, h);
 
437
  }
 
438
  dpd_file4_close(&W);
 
439
 
 
440
  dpd_file4_init(&W, CC_TMP2, 0, 11, 7,"CC2 Wabei (ei,a>b)");
 
441
  for(h=0; h < nirreps; h++) {
 
442
    dpd_file4_mat_irrep_init(&W, h);
 
443
    dpd_file4_mat_irrep_rd(&W, h);
 
444
    for(ei=0; ei<W.params->rowtot[h]; ei++) {
 
445
      i = W.params->roworb[h][ei][1];
 
446
      isym = W.params->qsym[i]; 
 
447
      I = i - occ_off[isym];
 
448
      for(ab=0; ab<W.params->coltot[h]; ab++) {
 
449
        if (I >= (occpi[isym] - openpi[isym]))
 
450
          W.matrix[h][ei][ab] = 0.0;
 
451
      }
 
452
    }
 
453
    dpd_file4_mat_irrep_wrt(&W, h);
 
454
    dpd_file4_mat_irrep_close(&W, h);
 
455
  }
 
456
  dpd_file4_close(&W);
 
457
 
 
458
  dpd_file4_init(&W, CC_TMP2, 0, 11, 5,"CC2 WAbEi (Ei,Ab)");
 
459
  for(h=0; h < nirreps; h++) {
 
460
    dpd_file4_mat_irrep_init(&W, h);
 
461
    dpd_file4_mat_irrep_rd(&W, h);
 
462
    for(ei=0; ei<W.params->rowtot[h]; ei++) {
 
463
      e = W.params->roworb[h][ei][0];
 
464
      i = W.params->roworb[h][ei][1];
 
465
      esym = W.params->psym[e];
 
466
      isym = W.params->qsym[i];
 
467
      E = e - vir_off[esym];
 
468
      I = i - occ_off[isym];
 
469
      for(ab=0; ab<W.params->coltot[h]; ab++) {
 
470
        a = W.params->colorb[h][ab][0];
 
471
        asym = W.params->rsym[a];
 
472
        bsym = W.params->ssym[b];
 
473
        A = a - vir_off[asym];
 
474
        if ((E >= (virtpi[esym] - openpi[esym])) ||
 
475
            (I >= (occpi[isym] - openpi[isym])) ||
 
476
            (A >= (virtpi[asym] - openpi[asym])) )
 
477
          W.matrix[h][ei][ab] = 0.0;
 
478
      }
 
479
    }
 
480
    dpd_file4_mat_irrep_wrt(&W, h);
 
481
    dpd_file4_mat_irrep_close(&W, h);
 
482
  }
 
483
  dpd_file4_close(&W);
 
484
 
 
485
  dpd_file4_init(&W, CC_TMP2, 0, 11, 5,"CC2 WaBeI (eI,aB)");
 
486
  for(h=0; h < nirreps; h++) {
 
487
    dpd_file4_mat_irrep_init(&W, h);
 
488
    dpd_file4_mat_irrep_rd(&W, h);
 
489
    for(ei=0; ei<W.params->rowtot[h]; ei++) {
 
490
      for(ab=0; ab<W.params->coltot[h]; ab++) {
 
491
        b = W.params->colorb[h][ab][1];
 
492
        bsym = W.params->ssym[b];
 
493
        B = b - vir_off[bsym];
 
494
        if (B >= (virtpi[bsym] - openpi[bsym]))
 
495
          W.matrix[h][ei][ab] = 0.0;
 
496
      }
 
497
    }
 
498
    dpd_file4_mat_irrep_wrt(&W, h);
 
499
    dpd_file4_mat_irrep_close(&W, h);
 
500
  }
 
501
  dpd_file4_close(&W);
 
502
 
 
503
}
 
504
 
 
505
 
 
506
}} // namespace psi::cchbar