← Back to branch summary

~ubuntu-branches/ubuntu/vivid/psicode/vivid

« back to all changes in this revision

Viewing changes to src/bin/ccenergy/cc2_Wmbij.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/ccenergy/cc2_Wmbij.cc
 
1
/*! \file
 
2
    \ingroup CCENERGY
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#include <cstdio>
 
6
#include <cstdlib>
 
7
#include <libciomr/libciomr.h>
 
8
#include <libdpd/dpd.h>
 
9
#include <libqt/qt.h>
 
10
#include "Params.h"
 
11
#include "MOInfo.h"
 
12
#define EXTERN
 
13
#include "globals.h"
 
14
 
 
15
namespace psi { namespace ccenergy {
 
16
 
 
17
void purge_cc2_Wmbij(void);
 
18
 
 
19
void cc2_Wmbij_build(void)
 
20
{
 
21
 
 
22
  dpdfile2 t1, tIA, tia;
 
23
  dpdbuf4 E, D, C, F;
 
24
  dpdbuf4 W, W1;
 
25
  dpdbuf4 Z, Z1;
 
26
  dpdbuf4 X;
 
27
 
 
28
  if(params.ref == 0) { /** RHF **/
 
29
    /* W(Mb,Ij) <-- <Mb|Ij> */
 
30
    dpd_buf4_init(&E, CC_EINTS, 0, 10, 0, 10, 0, 0, "E <ia|jk>");
 
31
    dpd_buf4_copy(&E, CC2_HET1, "CC2 WMbIj");
 
32
    dpd_buf4_close(&E);
 
33
  }
 
34
  else if(params.ref == 1) { /** ROHF **/  
 
35
    /** W(MB,I>J) <--- <MB||IJ> **/
 
36
    /** W(mb,i>j) <--- <mb||ij> **/
 
37
    dpd_buf4_init(&E, CC_EINTS, 0, 2, 10, 2, 10, 0, "E <ij||ka> (i>j,ka)");
 
38
    dpd_buf4_sort(&E, CC2_HET1, rspq, 10, 2, "CC2 WMBIJ (MB,I>J)");
 
39
    dpd_buf4_sort(&E, CC2_HET1, rspq, 10, 2, "CC2 Wmbij (mb,i>j)");
 
40
    dpd_buf4_close(&E);
 
41
 
 
42
    /** W(Mb,Ij) <--- <Mb|Ij> **/
 
43
    /** W(mB,iJ) <--- <mB|iJ> **/
 
44
    dpd_buf4_init(&E, CC_EINTS, 0, 0, 10, 0, 10, 0, "E <ij|ka>");
 
45
    dpd_buf4_sort(&E, CC2_HET1, rspq, 10, 0, "CC2 WMbIj");
 
46
    dpd_buf4_sort(&E, CC2_HET1, rspq, 10, 0, "CC2 WmBiJ (mB,iJ)");
 
47
    dpd_buf4_close(&E);
 
48
  }
 
49
  else if(params.ref == 2) { /*** UHF ***/
 
50
    /** W(MB,I>J) <--- <MB||IJ> **/
 
51
    dpd_buf4_init(&E, CC_EINTS, 0, 2, 20, 2, 20, 0, "E <IJ||KA> (I>J,KA)");
 
52
    dpd_buf4_sort(&E, CC2_HET1, rspq, 20, 2, "CC2 WMBIJ (MB,I>J)");
 
53
    dpd_buf4_close(&E);
 
54
 
 
55
    /** W(mb,i>j) <--- <mb||ij> **/
 
56
    dpd_buf4_init(&E, CC_EINTS, 0, 12, 30, 12, 30, 0, "E <ij||ka> (i>j,ka)");
 
57
    dpd_buf4_sort(&E, CC2_HET1, rspq, 30, 12, "CC2 Wmbij (mb,i>j)");
 
58
    dpd_buf4_close(&E);
 
59
 
 
60
    /** W(Mb,Ij) <--- <Mb|Ij> **/
 
61
    dpd_buf4_init(&E, CC_EINTS, 0, 22, 24, 22, 24, 0, "E <Ij|Ka>");
 
62
    dpd_buf4_sort(&E, CC2_HET1, rspq, 24, 22, "CC2 WMbIj");
 
63
    dpd_buf4_close(&E);
 
64
 
 
65
    /** W(mB,iJ) <--- <mB|iJ> **/
 
66
    dpd_buf4_init(&E, CC_EINTS, 0, 23, 27, 23, 27, 0, "E <iJ|kA>");
 
67
    dpd_buf4_sort(&E, CC2_HET1, rspq, 27, 23, "CC2 WmBiJ (mB,iJ)");
 
68
    dpd_buf4_close(&E);
 
69
  }
 
70
 
 
71
  if(params.ref == 0) { /** RHF **/
 
72
 
 
73
    dpd_file2_init(&t1, CC_OEI, 0, 0, 1, "tIA");
 
74
 
 
75
    /* W(Mb,Ij) <-- -t(n,b) * W(Mn,Ij) */
 
76
    dpd_buf4_init(&W, CC2_HET1, 0, 10, 0, 10, 0, 0, "CC2 WMbIj");
 
77
    dpd_buf4_init(&W1, CC2_HET1, 0, 0, 0, 0, 0, 0, "CC2 WMnIj");
 
78
    dpd_contract424(&W1, &t1, &W, 1, 0, 1, -0.5, 1);
 
79
    dpd_buf4_close(&W1);
 
80
    dpd_buf4_close(&W);
 
81
 
 
82
    dpd_file2_close(&t1);
 
83
  }
 
84
  else if(params.ref == 1) { /** ROHF **/  
 
85
 
 
86
    dpd_file2_init(&tIA, CC_OEI, 0, 0, 1, "tIA");
 
87
    dpd_file2_init(&tia, CC_OEI, 0, 0, 1, "tia");
 
88
 
 
89
    /**** W(MB,I>J) <-- -ZMBJI <-- P(I/J)( -<JE||MB> * t1[I][E] ) ****/
 
90
    dpd_buf4_init(&Z, CC_TMP0, 0, 10, 0, 10, 0, 0, "Z (MB,JI)");
 
91
    dpd_buf4_init(&C, CC_CINTS, 0, 10, 10, 10, 10, 0, "C <ia||jb>");
 
92
    dpd_contract424(&C, &tIA, &Z, 1, 1, 0, -1, 0);
 
93
    dpd_buf4_close(&C);
 
94
 
 
95
    dpd_buf4_sort(&Z, CC_TMP0, pqsr, 10, 0, "X (MB,IJ)");
 
96
    dpd_buf4_init(&X, CC_TMP0, 0, 10, 0, 10, 0, 0, "X (MB,IJ)");
 
97
    dpd_buf4_axpy(&Z, &X, -1);
 
98
    dpd_buf4_close(&Z);
 
99
    dpd_buf4_init(&W, CC2_HET1, 0, 10, 0, 10, 2, 0, "CC2 WMBIJ (MB,I>J)");
 
100
    dpd_buf4_axpy(&X, &W, 1);
 
101
    dpd_buf4_close(&X);
 
102
    dpd_buf4_close(&W);
 
103
 
 
104
    /**** W(mb,i>j) <-- -Zmbji <-- P(i/j)( -<je||mb> * t1[i][e] ) ****/
 
105
    dpd_buf4_init(&Z, CC_TMP0, 0, 10, 0, 10, 0, 0, "Z (mb,ji)");
 
106
    dpd_buf4_init(&C, CC_CINTS, 0, 10, 10, 10, 10, 0, "C <ia||jb>");
 
107
    dpd_contract424(&C, &tia, &Z, 1, 1, 0, -1, 0);
 
108
    dpd_buf4_close(&C);
 
109
 
 
110
    dpd_buf4_sort(&Z, CC_TMP0, pqsr, 10, 0, "X (mb,ij)");
 
111
    dpd_buf4_init(&X, CC_TMP0, 0, 10, 0, 10, 0, 0, "X (mb,ij)");
 
112
    dpd_buf4_axpy(&Z, &X, -1);
 
113
    dpd_buf4_close(&Z);
 
114
    dpd_buf4_init(&W, CC2_HET1, 0, 10, 0, 10, 2, 0, "CC2 Wmbij (mb,i>j)");
 
115
    dpd_buf4_axpy(&X, &W, 1);
 
116
    dpd_buf4_close(&X);
 
117
    dpd_buf4_close(&W);
 
118
 
 
119
    /**** W(Mb,Ij) <-- ( <Mb|Ej> * t1[I][E] ) ****/
 
120
    dpd_buf4_init(&W, CC2_HET1, 0, 10, 0, 10, 0, 0, "CC2 WMbIj");
 
121
    dpd_buf4_init(&D, CC_DINTS, 0, 10, 11, 10, 11, 0, "D <ij|ab> (ib,aj)");
 
122
    dpd_contract244(&tIA, &D, &W, 1, 2, 1, 1, 1);
 
123
    dpd_buf4_close(&D);
 
124
    dpd_buf4_close(&W);
 
125
 
 
126
    /**** W(Mb,Ij) <-- ( <Mb|Ie> * t1[j][e] ) ****/
 
127
    dpd_buf4_init(&W, CC2_HET1, 0, 10, 0, 10, 0, 0, "CC2 WMbIj");
 
128
    dpd_buf4_init(&C, CC_CINTS, 0, 10, 10, 10, 10, 0, "C <ia|jb>");
 
129
    dpd_contract424(&C, &tia, &W, 3, 1, 0, 1, 1);
 
130
    dpd_buf4_close(&C);
 
131
    dpd_buf4_close(&W);
 
132
 
 
133
    /**** W(mB,iJ) <-- ( <mB|eJ> * t1[i][e] ) ****/
 
134
    dpd_buf4_init(&W, CC2_HET1, 0, 10, 0, 10, 0, 0, "CC2 WmBiJ (mB,iJ)");
 
135
    dpd_buf4_init(&D, CC_DINTS, 0, 10, 11, 10, 11, 0, "D <ij|ab> (ib,aj)");
 
136
    dpd_contract244(&tia, &D, &W, 1, 2, 1, 1, 1);
 
137
    dpd_buf4_close(&D);
 
138
    dpd_buf4_close(&W);
 
139
 
 
140
    /**** W(mB,iJ) <-- ( <mB|iE> * t1[J][E] ) ****/
 
141
    dpd_buf4_init(&W, CC2_HET1, 0, 10, 0, 10, 0, 0, "CC2 WmBiJ (mB,iJ)");
 
142
    dpd_buf4_init(&C, CC_CINTS, 0, 10, 10, 10, 10, 0, "C <ia|jb>");
 
143
    dpd_contract424(&C, &tIA, &W, 3, 1, 0, 1, 1);
 
144
    dpd_buf4_close(&C);
 
145
    dpd_buf4_close(&W);
 
146
 
 
147
    dpd_file2_close(&tIA);
 
148
    dpd_file2_close(&tia);
 
149
  }
 
150
  else if(params.ref == 2) { /*** UHF ***/
 
151
 
 
152
    dpd_file2_init(&tIA, CC_OEI, 0, 0, 1, "tIA");
 
153
    dpd_file2_init(&tia, CC_OEI, 0, 2, 3, "tia");
 
154
 
 
155
    /**** W(MB,I>J) <-- -ZMBJI <-- P(I/J)( -<JE||MB> * t1[I][E] ) ****/
 
156
    dpd_buf4_init(&Z, CC_TMP0, 0, 20, 0, 20, 0, 0, "Z (MB,JI)");
 
157
    dpd_buf4_init(&C, CC_CINTS, 0, 20, 20, 20, 20, 0, "C <IA||JB>");
 
158
    dpd_contract424(&C, &tIA, &Z, 1, 1, 0, -1, 0);
 
159
    dpd_buf4_close(&C);
 
160
 
 
161
    dpd_buf4_sort(&Z, CC_TMP0, pqsr, 20, 0, "X (MB,IJ)");
 
162
    dpd_buf4_init(&X, CC_TMP0, 0, 20, 0, 20, 0, 0, "X (MB,IJ)");
 
163
    dpd_buf4_axpy(&Z, &X, -1);
 
164
    dpd_buf4_close(&Z);
 
165
    dpd_buf4_init(&W, CC2_HET1, 0, 20, 0, 20, 2, 0, "CC2 WMBIJ (MB,I>J)");
 
166
    dpd_buf4_axpy(&X, &W, 1);
 
167
    dpd_buf4_close(&X);
 
168
    dpd_buf4_close(&W);
 
169
 
 
170
    /**** W(mb,i>j) <-- -Zmbji <-- P(i/j)( -<je||mb> * t1[i][e] ) ****/
 
171
    dpd_buf4_init(&Z, CC_TMP0, 0, 30, 10, 30, 10, 0, "Z (mb,ji)");
 
172
    dpd_buf4_init(&C, CC_CINTS, 0, 30, 30, 30, 30, 0, "C <ia||jb>");
 
173
    dpd_contract424(&C, &tia, &Z, 1, 1, 0, -1, 0);
 
174
    dpd_buf4_close(&C);
 
175
 
 
176
    dpd_buf4_sort(&Z, CC_TMP0, pqsr, 30, 10, "X (mb,ij)");
 
177
    dpd_buf4_init(&X, CC_TMP0, 0, 30, 10, 30, 10, 0, "X (mb,ij)");
 
178
    dpd_buf4_axpy(&Z, &X, -1);
 
179
    dpd_buf4_close(&Z);
 
180
    dpd_buf4_init(&W, CC2_HET1, 0, 30, 10, 30, 12, 0, "CC2 Wmbij (mb,i>j)");
 
181
    dpd_buf4_axpy(&X, &W, 1);
 
182
    dpd_buf4_close(&X);
 
183
    dpd_buf4_close(&W);
 
184
 
 
185
    /**** W(Mb,Ij) <-- ( <Mb|Ej> * t1[I][E] ) ****/
 
186
    dpd_buf4_init(&W, CC2_HET1, 0, 24, 22, 24, 22, 0, "CC2 WMbIj");
 
187
    dpd_buf4_init(&D, CC_DINTS, 0, 24, 26, 24, 26, 0, "D <Ij|Ab> (Ib,Aj)");
 
188
    dpd_contract244(&tIA, &D, &W, 1, 2, 1, 1, 1);
 
189
    dpd_buf4_close(&D);
 
190
    dpd_buf4_close(&W);
 
191
 
 
192
    /**** W(Mb,Ij) <-- ( <Mb|Ie> * t1[j][e] ) ****/
 
193
    dpd_buf4_init(&W, CC2_HET1, 0, 24, 22, 24, 22, 0, "CC2 WMbIj");
 
194
    dpd_buf4_init(&C, CC_CINTS, 0, 24, 24, 24, 24, 0, "C <Ia|Jb>");
 
195
    dpd_contract424(&C, &tia, &W, 3, 1, 0, 1, 1);
 
196
    dpd_buf4_close(&C);
 
197
    dpd_buf4_close(&W);
 
198
 
 
199
    /**** W(mB,iJ) <-- ( <mB|eJ> * t1[i][e] ) ****/
 
200
    dpd_buf4_init(&W, CC2_HET1, 0, 27, 23, 27, 23, 0, "CC2 WmBiJ (mB,iJ)");
 
201
    dpd_buf4_init(&D, CC_DINTS, 0, 27, 25, 27, 25, 0, "D <iJ|aB> (iB,aJ)");
 
202
    dpd_contract244(&tia, &D, &W, 1, 2, 1, 1, 1);
 
203
    dpd_buf4_close(&D);
 
204
    dpd_buf4_close(&W);
 
205
 
 
206
    /**** W(mB,iJ) <-- ( <mB|iE> * t1[J][E] ) ****/
 
207
    dpd_buf4_init(&W, CC2_HET1, 0, 27, 23, 27, 23, 0, "CC2 WmBiJ (mB,iJ)");
 
208
    dpd_buf4_init(&C, CC_CINTS, 0, 27, 27, 27, 27, 0, "C <iA|jB>");
 
209
    dpd_contract424(&C, &tIA, &W, 3, 1, 0, 1, 1);
 
210
    dpd_buf4_close(&C);
 
211
    dpd_buf4_close(&W);
 
212
 
 
213
    dpd_file2_close(&tIA);
 
214
    dpd_file2_close(&tia);
 
215
  }
 
216
 
 
217
  if(params.ref == 0) { /** RHF **/
 
218
 
 
219
    dpd_file2_init(&t1, CC_OEI, 0, 0, 1, "tIA");
 
220
 
 
221
    /* W(Mb,Ij) <-- +P(ij) t(I,E) * <Mb||Ej> */
 
222
    dpd_buf4_init(&W, CC2_HET1, 0, 10, 0, 10, 0, 0, "CC2 WMbIj");
 
223
    dpd_buf4_init(&C, CC_CINTS, 0, 10, 10, 10, 10, 0, "C <ia|jb>");
 
224
    dpd_contract424(&C, &t1, &W, 3, 1, 0, 1, 1);
 
225
    dpd_buf4_close(&C);
 
226
    dpd_buf4_close(&W);
 
227
 
 
228
    dpd_buf4_init(&Z, CC_TMP0, 0, 0, 11, 0, 11, 0, "CC2 ZMbIj (jM,bI)");
 
229
    dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
 
230
    dpd_contract424(&D, &t1, &Z, 3, 1, 0, 1, 0);
 
231
    dpd_buf4_close(&D);
 
232
    dpd_buf4_sort_axpy(&Z, CC2_HET1, qrsp, 10, 0, "CC2 WMbIj", 1);
 
233
    dpd_buf4_close(&Z);
 
234
 
 
235
    dpd_file2_close(&t1);
 
236
  }
 
237
  else if(params.ref == 1) { /** ROHF **/  
 
238
 
 
239
    dpd_file2_init(&tIA, CC_OEI, 0, 0, 1, "tIA");
 
240
    dpd_file2_init(&tia, CC_OEI, 0, 0, 1, "tia");
 
241
 
 
242
    /**** W(MB,I>J) <-- ( t1[N][B] * W(MN,I>J) ****/
 
243
    dpd_buf4_init(&W, CC2_HET1, 0, 10, 2, 10, 2, 0, "CC2 WMBIJ (MB,I>J)");
 
244
    dpd_buf4_init(&Z, CC2_HET1, 0, 0, 2, 2, 2, 0, "CC2 WMNIJ (M>N,I>J)");
 
245
    dpd_contract424(&Z, &tIA, &W, 1, 0, 1, -0.5, 1);
 
246
    dpd_buf4_close(&Z);
 
247
    dpd_buf4_close(&W);
 
248
 
 
249
    /**** W(mb,i>j) <-- ( t1[n][b] * W(mn,i>j) ****/
 
250
    dpd_buf4_init(&W, CC2_HET1, 0, 10, 2, 10, 2, 0, "CC2 Wmbij (mb,i>j)");
 
251
    dpd_buf4_init(&Z, CC2_HET1, 0, 0, 2, 2, 2, 0, "CC2 Wmnij (m>n,i>j)");
 
252
    dpd_contract424(&Z, &tia, &W, 1, 0, 1, -0.5, 1);
 
253
    dpd_buf4_close(&Z);
 
254
    dpd_buf4_close(&W);
 
255
 
 
256
    /**** W(Mb,Ij) <-- ( t1[n][b] * W(Mn,Ij) ****/
 
257
    dpd_buf4_init(&W, CC2_HET1, 0, 10, 0, 10, 0, 0, "CC2 WMbIj");
 
258
    dpd_buf4_init(&Z, CC2_HET1, 0, 0, 0, 0, 0, 0, "CC2 WMnIj");
 
259
    dpd_contract424(&Z, &tia, &W, 1, 0, 1, -0.5, 1);
 
260
    dpd_buf4_close(&Z);
 
261
    dpd_buf4_close(&W);
 
262
 
 
263
    /**** W(mB,iJ) <-- ( t1[N][B] * W(mN,iJ) ****/
 
264
    dpd_buf4_init(&Z1, CC_TMP0, 0, 11, 0, 11, 0, 0, "Z (Bm,Ji)");
 
265
    dpd_buf4_init(&Z, CC2_HET1, 0, 0, 0, 0, 0, 0, "CC2 WMnIj");
 
266
    dpd_contract244(&tIA, &Z, &Z1, 0, 0, 0, -0.5, 0);
 
267
    dpd_buf4_close(&Z);
 
268
    dpd_buf4_sort_axpy(&Z1, CC2_HET1, qpsr, 10, 0, "CC2 WmBiJ (mB,iJ)", 1);
 
269
    dpd_buf4_close(&Z1);
 
270
 
 
271
    dpd_file2_close(&tIA);
 
272
    dpd_file2_close(&tia);
 
273
  }
 
274
  else if(params.ref == 2) { /*** UHF ***/
 
275
 
 
276
    dpd_file2_init(&tIA, CC_OEI, 0, 0, 1, "tIA");
 
277
    dpd_file2_init(&tia, CC_OEI, 0, 2, 3, "tia");
 
278
 
 
279
    /**** W(MB,I>J) <-- ( t1[N][B] * W(MN,I>J) ****/
 
280
    dpd_buf4_init(&W, CC2_HET1, 0, 20, 2, 20, 2, 0, "CC2 WMBIJ (MB,I>J)");
 
281
    dpd_buf4_init(&Z, CC2_HET1, 0, 0, 2, 2, 2, 0, "CC2 WMNIJ (M>N,I>J)");
 
282
    dpd_contract424(&Z, &tIA, &W, 1, 0, 1, -0.5, 1);
 
283
    dpd_buf4_close(&Z);
 
284
    dpd_buf4_close(&W);
 
285
 
 
286
    /**** W(mb,i>j) <-- ( t1[n][b] * W(mn,i>j) ****/
 
287
    dpd_buf4_init(&W, CC2_HET1, 0, 30, 12, 30, 12, 0, "CC2 Wmbij (mb,i>j)");
 
288
    dpd_buf4_init(&Z, CC2_HET1, 0, 10, 12, 12, 12, 0, "CC2 Wmnij (m>n,i>j)");
 
289
    dpd_contract424(&Z, &tia, &W, 1, 0, 1, -0.5, 1);
 
290
    dpd_buf4_close(&Z);
 
291
    dpd_buf4_close(&W);
 
292
 
 
293
    /**** W(Mb,Ij) <-- ( t1[n][b] * W(Mn,Ij) ****/
 
294
    dpd_buf4_init(&W, CC2_HET1, 0, 24, 22, 24, 22, 0, "CC2 WMbIj");
 
295
    dpd_buf4_init(&Z, CC2_HET1, 0, 22, 22, 22, 22, 0, "CC2 WMnIj");
 
296
    dpd_contract424(&Z, &tia, &W, 1, 0, 1, -0.5, 1);
 
297
    dpd_buf4_close(&Z);
 
298
    dpd_buf4_close(&W);
 
299
 
 
300
    /**** W(mB,iJ) <-- ( t1[N][B] * W(mN,iJ) ****/
 
301
    dpd_buf4_init(&Z1, CC_TMP0, 0, 26, 22, 26, 22, 0, "Z (Bm,Ji)");
 
302
    dpd_buf4_init(&Z, CC2_HET1, 0, 22, 22, 22, 22, 0, "CC2 WMnIj");
 
303
    dpd_contract244(&tIA, &Z, &Z1, 0, 0, 0, -0.5, 0);
 
304
    dpd_buf4_close(&Z);
 
305
    dpd_buf4_sort_axpy(&Z1, CC2_HET1, qpsr, 27, 23, "CC2 WmBiJ (mB,iJ)", 1);
 
306
    dpd_buf4_close(&Z1);
 
307
 
 
308
    dpd_file2_close(&tIA);
 
309
    dpd_file2_close(&tia);
 
310
  }
 
311
 
 
312
  if(params.ref == 0) { /** RHF **/
 
313
 
 
314
    dpd_file2_init(&t1, CC_OEI, 0, 0, 1, "tIA");
 
315
 
 
316
    /* Z(Mb,Ej) = <Mb|Ej> + t(j,f) * <Mb|Ef> */
 
317
    dpd_buf4_init(&Z, CC_TMP0, 0, 10, 11, 10, 11, 0, "CC2 ZMbEj (Mb,Ej)");
 
318
    dpd_buf4_init(&F, CC_FINTS, 0, 10, 5, 10, 5, 0, "F <ia|bc>");
 
319
    dpd_contract424(&F, &t1, &Z, 3, 1, 0, 1, 0);
 
320
    dpd_buf4_close(&F);
 
321
    dpd_buf4_close(&Z);
 
322
 
 
323
    /* W(Mb,Ij) <-- t(I,E) * Z(Mb,Ej) */
 
324
    dpd_buf4_init(&Z, CC_TMP0, 0, 10, 0, 10, 0, 0, "CC2 ZMbIj");
 
325
    dpd_buf4_init(&Z1, CC_TMP0, 0, 10, 11, 10, 11, 0, "CC2 ZMbEj (Mb,Ej)");
 
326
    dpd_contract244(&t1, &Z1, &Z, 1, 2, 1, 1, 0);
 
327
    dpd_buf4_close(&Z1);
 
328
    dpd_buf4_init(&W, CC2_HET1, 0, 10, 0, 10, 0, 0, "CC2 WMbIj");
 
329
    dpd_buf4_axpy(&Z, &W, 1);
 
330
    dpd_buf4_close(&W);
 
331
    dpd_buf4_close(&Z);
 
332
 
 
333
    dpd_file2_close(&t1);
 
334
  }
 
335
  else if(params.ref == 1) { /** ROHF **/  
 
336
 
 
337
    dpd_file2_init(&tIA, CC_OEI, 0, 0, 1, "tIA");
 
338
    dpd_file2_init(&tia, CC_OEI, 0, 0, 1, "tia");
 
339
 
 
340
    /**** W(MB,I>J) <-- 0.5*P(I/J)XMBIJ <--- ( t1[I][E] * ZMBEJ ) <--  <MB||EF> * t1[J][F] ****/
 
341
    dpd_buf4_init(&Z, CC_TMP0, 0, 10, 11, 10, 11, 0, "Z (MB,EJ)");
 
342
    dpd_buf4_init(&F, CC_FINTS, 0, 10, 5, 10, 5, 1, "F <ia|bc>");
 
343
    dpd_contract424(&F, &tIA, &Z, 3, 1, 0, 1, 0);
 
344
    dpd_buf4_close(&F);
 
345
 
 
346
    dpd_buf4_init(&W, CC2_HET1, 0, 10, 0, 10, 2, 0, "CC2 WMBIJ (MB,I>J)");
 
347
    dpd_contract244(&tIA, &Z, &W, 1, 2, 1, 1, 1);
 
348
    dpd_buf4_close(&Z);
 
349
    dpd_buf4_close(&W);
 
350
 
 
351
    /**** W(mb,i>j) <-- P(i/j) (Zmbif * t1[j][f]) <-- 0.5*( t1[i][e] * <mb||ef> ) ****/
 
352
    dpd_buf4_init(&Z, CC_TMP0, 0, 10, 11, 10, 11, 0, "Z (mb,ej)");
 
353
    dpd_buf4_init(&F, CC_FINTS, 0, 10, 5, 10, 5, 1, "F <ia|bc>");
 
354
    dpd_contract424(&F, &tia, &Z, 3, 1, 0, 1, 0);
 
355
    dpd_buf4_close(&F);
 
356
 
 
357
    dpd_buf4_init(&W, CC2_HET1, 0, 10, 0, 10, 2, 0, "CC2 Wmbij (mb,i>j)");
 
358
    dpd_contract244(&tia, &Z, &W, 1, 2, 1, 1, 1);
 
359
    dpd_buf4_close(&Z);
 
360
    dpd_buf4_close(&W);
 
361
 
 
362
    /**** W(Mb,Ij) <-- (ZIfMb * t1[j][f]) <--  t1[I][E] * <Mb|Ef> ****/
 
363
    dpd_buf4_init(&Z, CC_TMP0, 0, 10, 10, 10, 10, 0, "Z (If,Mb)");
 
364
    dpd_buf4_init(&F, CC_FINTS, 0, 10, 5, 10, 5, 0, "F <ia|bc>");
 
365
    dpd_contract244(&tIA, &F, &Z, 1, 2, 0, 1, 0);
 
366
    dpd_buf4_close(&F);
 
367
 
 
368
    dpd_buf4_init(&W, CC2_HET1, 0, 10, 0, 10, 0, 0, "CC2 WMbIj");
 
369
    dpd_contract424(&Z, &tia, &W, 1, 1, 0, 1, 1);
 
370
    dpd_buf4_close(&Z);
 
371
    dpd_buf4_close(&W);
 
372
 
 
373
    /**** W(mB,iJ) <-- ZiFmB * t1[J][F] <-- t1[i][e] * <mB|eF> ****/
 
374
    dpd_buf4_init(&Z, CC_TMP0, 0, 10, 10, 10, 10, 0, "Z (iF,mB)");
 
375
    dpd_buf4_init(&F, CC_FINTS, 0, 10, 5, 10, 5, 0, "F <ia|bc>");
 
376
    dpd_contract244(&tia, &F, &Z, 1, 2, 0, 1, 0);
 
377
    dpd_buf4_close(&F);
 
378
 
 
379
    dpd_buf4_init(&W, CC2_HET1, 0, 10, 0, 10, 0, 0, "CC2 WmBiJ (mB,iJ)");
 
380
    dpd_contract424(&Z, &tIA, &W, 1, 1, 0, 1, 1);
 
381
    dpd_buf4_close(&Z);
 
382
    dpd_buf4_close(&W);
 
383
 
 
384
    dpd_file2_close(&tIA);
 
385
    dpd_file2_close(&tia);
 
386
  }
 
387
  else if(params.ref == 2) { /*** UHF ***/
 
388
 
 
389
    dpd_file2_init(&tIA, CC_OEI, 0, 0, 1, "tIA");
 
390
    dpd_file2_init(&tia, CC_OEI, 0, 2, 3, "tia");
 
391
 
 
392
    /**** W(MB,I>J) <-- 0.5*P(I/J)XMBIJ <--- ( t1[I][E] * ZMBEJ ) <--  <MB||EF> * t1[J][F] ****/
 
393
    dpd_buf4_init(&Z, CC_TMP0, 0, 20, 21, 20, 21, 0, "Z (MB,EJ)");
 
394
    dpd_buf4_init(&F, CC_FINTS, 0, 20, 5, 20, 5, 1, "F <IA|BC>");
 
395
    dpd_contract424(&F, &tIA, &Z, 3, 1, 0, 1, 0);
 
396
    dpd_buf4_close(&F);
 
397
 
 
398
    dpd_buf4_init(&W, CC2_HET1, 0, 20, 0, 20, 2, 0, "CC2 WMBIJ (MB,I>J)");
 
399
    dpd_contract244(&tIA, &Z, &W, 1, 2, 1, 1, 1);
 
400
    dpd_buf4_close(&Z);
 
401
    dpd_buf4_close(&W);
 
402
 
 
403
    /**** W(mb,i>j) <-- P(i/j) (Zmbif * t1[j][f]) <-- 0.5*( t1[i][e] * <mb||ef> ) ****/
 
404
    dpd_buf4_init(&Z, CC_TMP0, 0, 30, 31, 30, 31, 0, "Z (mb,ej)");
 
405
    dpd_buf4_init(&F, CC_FINTS, 0, 30, 15, 30, 15, 1, "F <ia|bc>");
 
406
    dpd_contract424(&F, &tia, &Z, 3, 1, 0, 1, 0);
 
407
    dpd_buf4_close(&F);
 
408
 
 
409
    dpd_buf4_init(&W, CC2_HET1, 0, 30, 10, 30, 12, 0, "CC2 Wmbij (mb,i>j)");
 
410
    dpd_contract244(&tia, &Z, &W, 1, 2, 1, 1, 1);
 
411
    dpd_buf4_close(&Z);
 
412
    dpd_buf4_close(&W);
 
413
 
 
414
    /**** W(Mb,Ij) <-- (ZIfMb * t1[j][f]) <--  t1[I][E] * <Mb|Ef> ****/
 
415
    dpd_buf4_init(&Z, CC_TMP0, 0, 24, 24, 24, 24, 0, "Z (If,Mb)");
 
416
    dpd_buf4_init(&F, CC_FINTS, 0, 24, 28, 24, 28, 0, "F <Ia|Bc>");
 
417
    dpd_contract244(&tIA, &F, &Z, 1, 2, 0, 1, 0);
 
418
    dpd_buf4_close(&F);
 
419
 
 
420
    dpd_buf4_init(&W, CC2_HET1, 0, 24, 22, 24, 22, 0, "CC2 WMbIj");
 
421
    dpd_contract424(&Z, &tia, &W, 1, 1, 0, 1, 1);
 
422
    dpd_buf4_close(&Z);
 
423
    dpd_buf4_close(&W);
 
424
 
 
425
    /**** W(mB,iJ) <-- ZiFmB * t1[J][F] <-- t1[i][e] * <mB|eF> ****/
 
426
    dpd_buf4_init(&Z, CC_TMP0, 0, 27, 27, 27, 27, 0, "Z (iF,mB)");
 
427
    dpd_buf4_init(&F, CC_FINTS, 0, 27, 29, 27, 29, 0, "F <iA|bC>");
 
428
    dpd_contract244(&tia, &F, &Z, 1, 2, 0, 1, 0);
 
429
    dpd_buf4_close(&F);
 
430
 
 
431
    dpd_buf4_init(&W, CC2_HET1, 0, 27, 23, 27, 23, 0, "CC2 WmBiJ (mB,iJ)");
 
432
    dpd_contract424(&Z, &tIA, &W, 1, 1, 0, 1, 1);
 
433
    dpd_buf4_close(&Z);
 
434
    dpd_buf4_close(&W);
 
435
 
 
436
    dpd_file2_close(&tIA);
 
437
    dpd_file2_close(&tia);
 
438
  }
 
439
 
 
440
  if(params.ref == 1) { /** ROHF **/  
 
441
 
 
442
    /* do purge before sort */
 
443
    purge_cc2_Wmbij();
 
444
 
 
445
    /* do final sort to get (Ij,Mb) */
 
446
    dpd_buf4_init(&W, CC2_HET1, 0, 10, 2, 10, 2, 0, "CC2 WMBIJ (MB,I>J)");
 
447
    dpd_buf4_sort(&W, CC2_HET1, rspq, 2, 10, "CC2 WMBIJ (I>J,MB)");
 
448
    dpd_buf4_close(&W);
 
449
    dpd_buf4_init(&W, CC2_HET1, 0, 10, 2, 10, 2, 0, "CC2 Wmbij (mb,i>j)");
 
450
    dpd_buf4_sort(&W, CC2_HET1, rspq, 2, 10, "CC2 Wmbij (i>j,mb)");
 
451
    dpd_buf4_close(&W);
 
452
    dpd_buf4_init(&W, CC2_HET1, 0, 10, 0, 10, 0, 0, "CC2 WMbIj");
 
453
    dpd_buf4_sort(&W, CC2_HET1, rspq, 0, 10, "CC2 WMbIj (Ij,Mb)");
 
454
    dpd_buf4_close(&W);
 
455
    dpd_buf4_init(&W, CC2_HET1, 0, 10, 0, 10, 0, 0, "CC2 WmBiJ (mB,iJ)");
 
456
    dpd_buf4_sort(&W, CC2_HET1, rspq, 0, 10, "CC2 WmBiJ (iJ,mB)");
 
457
    dpd_buf4_close(&W); 
 
458
  }
 
459
  else if(params.ref == 2) { /*** UHF ***/
 
460
    dpd_buf4_init(&W, CC2_HET1, 0, 20, 2, 20, 2, 0, "CC2 WMBIJ (MB,I>J)");
 
461
    dpd_buf4_sort(&W, CC2_HET1, rspq, 2, 20, "CC2 WMBIJ (I>J,MB)");
 
462
    dpd_buf4_close(&W);
 
463
    dpd_buf4_init(&W, CC2_HET1, 0, 30, 12, 30, 12, 0, "CC2 Wmbij (mb,i>j)");
 
464
    dpd_buf4_sort(&W, CC2_HET1, rspq, 12, 30, "CC2 Wmbij (i>j,mb)");
 
465
    dpd_buf4_close(&W);
 
466
    dpd_buf4_init(&W, CC2_HET1, 0, 24, 22, 24, 22, 0, "CC2 WMbIj");
 
467
    dpd_buf4_sort(&W, CC2_HET1, rspq, 22, 24, "CC2 WMbIj (Ij,Mb)");
 
468
    dpd_buf4_close(&W);
 
469
    dpd_buf4_init(&W, CC2_HET1, 0, 27, 23, 27, 23, 0, "CC2 WmBiJ (mB,iJ)");
 
470
    dpd_buf4_sort(&W, CC2_HET1, rspq, 23, 27, "CC2 WmBiJ (iJ,mB)");
 
471
    dpd_buf4_close(&W);
 
472
  }
 
473
}
 
474
 
 
475
void purge_cc2_Wmbij(void) {
 
476
  dpdfile4 W;
 
477
  int *occpi, *virtpi;
 
478
  int h, a, b, e, f, i, j, m, n;
 
479
  int    A, B, E, F, I, J, M, N;
 
480
  int mn, ei, ma, ef, me, jb, mb, ij, ab;
 
481
  int asym, bsym, esym, fsym, isym, jsym, msym, nsym;
 
482
  int *occ_off, *vir_off;
 
483
  int *occ_sym, *vir_sym;
 
484
  int *openpi, nirreps;
 
485
 
 
486
  nirreps = moinfo.nirreps;
 
487
  occpi = moinfo.occpi; virtpi = moinfo.virtpi;
 
488
  occ_off = moinfo.occ_off; vir_off = moinfo.vir_off;
 
489
  occ_sym = moinfo.occ_sym; vir_sym = moinfo.vir_sym;
 
490
  openpi = moinfo.openpi;
 
491
 
 
492
  dpd_file4_init(&W, CC2_HET1, 0, 10, 2,"CC2 WMBIJ (MB,I>J)");
 
493
  for(h=0; h < nirreps; h++) {
 
494
    dpd_file4_mat_irrep_init(&W, h);
 
495
    dpd_file4_mat_irrep_rd(&W, h);
 
496
    for(mb=0; mb<W.params->rowtot[h]; mb++) {
 
497
      b = W.params->roworb[h][mb][1];
 
498
      bsym = W.params->qsym[b];
 
499
      B = b - vir_off[bsym];
 
500
      for(ij=0; ij<W.params->coltot[h]; ij++) {
 
501
        if (B >= (virtpi[bsym] - openpi[bsym]))
 
502
          W.matrix[h][mb][ij] = 0.0;
 
503
      }
 
504
    }
 
505
    dpd_file4_mat_irrep_wrt(&W, h);
 
506
    dpd_file4_mat_irrep_close(&W, h);
 
507
  }
 
508
  dpd_file4_close(&W);
 
509
 
 
510
  dpd_file4_init(&W, CC2_HET1, 0, 10, 2,"CC2 Wmbij (mb,i>j)");
 
511
  for(h=0; h < nirreps; h++) {
 
512
    dpd_file4_mat_irrep_init(&W, h);
 
513
    dpd_file4_mat_irrep_rd(&W, h);
 
514
    for(mb=0; mb<W.params->rowtot[h]; mb++) {
 
515
      m = W.params->roworb[h][mb][0];
 
516
      msym = W.params->psym[m];
 
517
      M = m - occ_off[msym];
 
518
      for(ij=0; ij<W.params->coltot[h]; ij++) {
 
519
        i = W.params->colorb[h][ij][0];
 
520
        j = W.params->colorb[h][ij][1];
 
521
        isym = W.params->rsym[i];
 
522
        jsym = W.params->ssym[j];
 
523
        I = i - occ_off[isym];
 
524
        J = j - occ_off[jsym];
 
525
        if ((M >= (occpi[msym] - openpi[msym])) ||
 
526
            (I >= (occpi[isym] - openpi[isym])) ||
 
527
            (J >= (occpi[jsym] - openpi[jsym])) )
 
528
          W.matrix[h][mb][ij] = 0.0;
 
529
      }
 
530
    }
 
531
    dpd_file4_mat_irrep_wrt(&W, h);
 
532
    dpd_file4_mat_irrep_close(&W, h);
 
533
  }
 
534
  dpd_file4_close(&W);
 
535
 
 
536
  dpd_file4_init(&W, CC2_HET1, 0, 10, 0,"CC2 WMbIj");
 
537
  for(h=0; h < nirreps; h++) {
 
538
    dpd_file4_mat_irrep_init(&W, h);
 
539
    dpd_file4_mat_irrep_rd(&W, h);
 
540
    for(mb=0; mb<W.params->rowtot[h]; mb++) {
 
541
      for(ij=0; ij<W.params->coltot[h]; ij++) {
 
542
        j = W.params->colorb[h][ij][1];
 
543
        jsym = W.params->ssym[j];
 
544
        J = j - occ_off[jsym];
 
545
        if (J >= (occpi[jsym] - openpi[jsym]))
 
546
          W.matrix[h][mb][ij] = 0.0;
 
547
      }
 
548
    }
 
549
    dpd_file4_mat_irrep_wrt(&W, h);
 
550
    dpd_file4_mat_irrep_close(&W, h);
 
551
  }
 
552
  dpd_file4_close(&W);
 
553
 
 
554
  dpd_file4_init(&W, CC2_HET1, 0, 10, 0,"CC2 WmBiJ (mB,iJ)");
 
555
  for(h=0; h < nirreps; h++) {
 
556
    dpd_file4_mat_irrep_init(&W, h);
 
557
    dpd_file4_mat_irrep_rd(&W, h);
 
558
    for(mb=0; mb<W.params->rowtot[h]; mb++) {
 
559
      m = W.params->roworb[h][mb][0];
 
560
      b = W.params->roworb[h][mb][1];
 
561
      msym = W.params->psym[m];
 
562
      bsym = W.params->qsym[b];
 
563
      M = m - occ_off[msym];
 
564
      B = b - vir_off[bsym];
 
565
      for(ij=0; ij<W.params->coltot[h]; ij++) {
 
566
        i = W.params->colorb[h][ij][0];
 
567
        isym = W.params->rsym[i];
 
568
        I = i - occ_off[isym];
 
569
        if ((M >= (occpi[msym] - openpi[msym])) ||
 
570
            (B >= (virtpi[bsym] - openpi[bsym])) ||
 
571
            (I >= (occpi[isym] - openpi[isym])) )
 
572
          W.matrix[h][mb][ij] = 0.0;
 
573
      }
 
574
    }
 
575
    dpd_file4_mat_irrep_wrt(&W, h);
 
576
    dpd_file4_mat_irrep_close(&W, h);
 
577
  }
 
578
  dpd_file4_close(&W);
 
579
}
 
580
}} // namespace psi::ccenergy