← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/ccenergy/cc2_Wmnij.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/ccenergy/cc2_Wmnij.c
1
 
#include <libdpd/dpd.h>
2
 
#define EXTERN
3
 
#include "globals.h"
4
 
 
5
 
/* cc2_Wmnij(): Compute the Wmnij matrix from CC2 theory, which is
6
 
** given in spin-orbitals as:
7
 
**
8
 
** Wmnij = <mn||ij> + P(ij) t_j^e <mn||ie> + t_i^e t_j^f <mn||ef>
9
 
**
10
 
** TDC, Feb 2004
11
 
*/
12
 
 
13
 
void purge_cc2_Wmnij(void);
14
 
 
15
 
void cc2_Wmnij_build(void)
16
 
{
17
 
  dpdbuf4 A, E, D, Z, W, Z1, X;
18
 
  dpdfile2 t1, tIA, tia;
19
 
 
20
 
  timer_on("A->Wmnij");
21
 
  if(params.ref == 0) { /** RHF **/
22
 
    /* Wmnij <- <mn||ij> */
23
 
    dpd_buf4_init(&A, CC_AINTS, 0, 0, 0, 0, 0, 0, "A <ij|kl>");
24
 
    dpd_buf4_copy(&A, CC2_HET1, "CC2 WMnIj");
25
 
    dpd_buf4_close(&A);
26
 
  }
27
 
  else if (params.ref == 1) { /** ROHF **/
28
 
    /** W(M>N,I>J) <--- <MN||IJ> **/
29
 
    /** W(m>n,i>j) <--- <mn||ij> **/
30
 
    dpd_buf4_init(&A, CC_AINTS, 0, 2, 2, 0, 0, 1, "A <ij|kl>");
31
 
    dpd_buf4_copy(&A, CC2_HET1, "CC2 WMNIJ (M>N,I>J)");
32
 
    dpd_buf4_copy(&A, CC2_HET1, "CC2 Wmnij (m>n,i>j)");
33
 
    dpd_buf4_close(&A);
34
 
 
35
 
    /** W(Mn,Ij) <--- <Mn|Ij> **/
36
 
    dpd_buf4_init(&A, CC_AINTS, 0, 0, 0, 0, 0, 0, "A <ij|kl>");
37
 
    dpd_buf4_copy(&A, CC2_HET1, "CC2 WMnIj");
38
 
    dpd_buf4_close(&A);
39
 
  }
40
 
  else if (params.ref == 2) { /** UHF **/
41
 
    /** W(M>N,I>J) <--- <MN||IJ> **/
42
 
    dpd_buf4_init(&A, CC_AINTS, 0, 2, 2, 0, 0, 1, "A <IJ|KL>");
43
 
    dpd_buf4_copy(&A, CC2_HET1, "CC2 WMNIJ (M>N,I>J)");
44
 
    dpd_buf4_close(&A);
45
 
 
46
 
    /** W(m>n,i>j) <--- <mn||ij> **/
47
 
    dpd_buf4_init(&A, CC_AINTS, 0, 12, 12, 10, 10, 1, "A <ij|kl>");
48
 
    dpd_buf4_copy(&A, CC2_HET1, "CC2 Wmnij (m>n,i>j)");
49
 
    dpd_buf4_close(&A);
50
 
 
51
 
    /** W(Mn,Ij) <--- <Mn|Ij> **/
52
 
    dpd_buf4_init(&A, CC_AINTS, 0, 22, 22, 22, 22, 0, "A <Ij|Kl>");
53
 
    dpd_buf4_copy(&A, CC2_HET1, "CC2 WMnIj");
54
 
    dpd_buf4_close(&A);
55
 
  }
56
 
  timer_off("A->Wmnij");
57
 
 
58
 
  timer_on("E->Wmnij");
59
 
  if(params.ref == 0) { /** RHF **/
60
 
 
61
 
    dpd_file2_init(&t1, CC_OEI, 0, 0, 1, "tIA");
62
 
 
63
 
    /* Wmnij <- + P(ij) t(j,e) * <mn||ie> */
64
 
    dpd_buf4_init(&Z, CC_TMP0, 0, 0, 0, 0, 0, 0, "CC2 ZMnIj");
65
 
    dpd_buf4_init(&E, CC_EINTS, 0, 0, 10, 0, 10, 0, "E <ij|ka>");
66
 
    dpd_contract424(&E, &t1, &Z, 3, 1, 0, 1, 0);
67
 
    dpd_buf4_close(&E);
68
 
 
69
 
    dpd_buf4_init(&W, CC2_HET1, 0, 0, 0, 0, 0, 0, "CC2 WMnIj");
70
 
    dpd_buf4_axpy(&Z, &W, 1);
71
 
    dpd_buf4_close(&W);
72
 
    dpd_buf4_sort_axpy(&Z, CC2_HET1, qpsr, 0, 0, "CC2 WMnIj", 1);
73
 
    dpd_buf4_close(&Z);
74
 
 
75
 
    dpd_file2_close(&t1);
76
 
  }
77
 
  else if (params.ref == 1) { /** ROHF **/
78
 
 
79
 
    dpd_file2_init(&tIA, CC_OEI, 0, 0, 1, "tIA");
80
 
    dpd_file2_init(&tia, CC_OEI, 0, 0, 1, "tia");
81
 
 
82
 
    /**** W(M>N,I>J) <-- ZMNIJ <-- P(I/J)( <MN||IE> * t1[J][E] ) ****/
83
 
    dpd_buf4_init(&Z, CC_TMP0, 0, 2, 0, 2, 0, 0, "Z (M>N,IJ)");
84
 
    dpd_buf4_init(&E, CC_EINTS, 0, 2, 10, 2, 10, 0, "E <ij||ka> (i>j,ka)");
85
 
    dpd_contract424(&E, &tIA, &Z, 3, 1, 0, 1, 0);
86
 
    dpd_buf4_close(&E);
87
 
 
88
 
    dpd_buf4_sort(&Z, CC_TMP0, pqsr, 2, 0, "Z (M>N,JI)");
89
 
    dpd_buf4_init(&Z1, CC_TMP0, 0, 2, 0, 2, 0, 0, "Z (M>N,JI)");
90
 
    dpd_buf4_axpy(&Z1, &Z, -1);
91
 
    dpd_buf4_close(&Z1);
92
 
 
93
 
    dpd_buf4_init(&W, CC2_HET1, 0, 2, 0, 2, 2, 0, "CC2 WMNIJ (M>N,I>J)");
94
 
    dpd_buf4_axpy(&Z, &W, 1);
95
 
    dpd_buf4_close(&Z);
96
 
    dpd_buf4_close(&W);
97
 
 
98
 
    /**** W(m>n,i>j) <-- Zmnij <-- P(i/j)( <mn||ie> * t1[j][e] ) ****/
99
 
    dpd_buf4_init(&Z, CC_TMP0, 0, 2, 0, 2, 0, 0, "Z (m>n,ij)");
100
 
    dpd_buf4_init(&E, CC_EINTS, 0, 2, 10, 2, 10, 0, "E <ij||ka> (i>j,ka)");
101
 
    dpd_contract424(&E, &tia, &Z, 3, 1, 0, 1, 0);
102
 
    dpd_buf4_close(&E);
103
 
 
104
 
    dpd_buf4_sort(&Z, CC_TMP0, pqsr, 2, 0, "Z (m>n,ji)");
105
 
    dpd_buf4_init(&Z1, CC_TMP0, 0, 2, 0, 2, 0, 0, "Z (m>n,ji)");
106
 
    dpd_buf4_axpy(&Z1, &Z, -1);
107
 
    dpd_buf4_close(&Z1);
108
 
 
109
 
    dpd_buf4_init(&W, CC2_HET1, 0, 2, 0, 2, 2, 0, "CC2 Wmnij (m>n,i>j)");
110
 
    dpd_buf4_axpy(&Z, &W, 1);
111
 
    dpd_buf4_close(&Z);
112
 
    dpd_buf4_close(&W);
113
 
 
114
 
    /**** W(Mn,Ij) <-- <Mn|Ie> * t1[j][e] ****/
115
 
    dpd_buf4_init(&W, CC2_HET1, 0, 0, 0, 0, 0, 0, "CC2 WMnIj");
116
 
    dpd_buf4_init(&E, CC_EINTS, 0, 0, 10, 0, 10, 0, "E <ij|ka>");
117
 
    dpd_contract424(&E, &tia, &W, 3, 1, 0, 1, 1);
118
 
    dpd_buf4_close(&E);
119
 
    dpd_buf4_close(&W);
120
 
 
121
 
    /**** W(Mn,Ij) <-- <Mn|Ej> * t1[I][E] ****/
122
 
    dpd_buf4_init(&W, CC2_HET1, 0, 0, 0, 0, 0, 0, "CC2 WMnIj");
123
 
    dpd_buf4_init(&E, CC_EINTS, 0, 0, 11, 0, 11, 0, "E <ij|ak>");
124
 
    dpd_contract244(&tIA, &E, &W, 1, 2, 1, 1, 1);
125
 
    dpd_buf4_close(&E);
126
 
    dpd_buf4_close(&W);
127
 
 
128
 
    dpd_file2_close(&tIA);
129
 
    dpd_file2_close(&tia);
130
 
  }
131
 
  else if (params.ref == 2) { /** UHF **/
132
 
 
133
 
    dpd_file2_init(&tIA, CC_OEI, 0, 0, 1, "tIA");
134
 
    dpd_file2_init(&tia, CC_OEI, 0, 2, 3, "tia");
135
 
 
136
 
    /**** W(M>N,I>J) <-- ZMNIJ <-- P(I/J)( <MN||IE> * t1[J][E] ) ****/
137
 
    dpd_buf4_init(&Z, CC_TMP0, 0, 2, 0, 2, 0, 0, "Z (M>N,IJ)");
138
 
    dpd_buf4_init(&E, CC_EINTS, 0, 2, 20, 2, 20, 0, "E <IJ||KA> (I>J,KA)");
139
 
    dpd_contract424(&E, &tIA, &Z, 3, 1, 0, 1, 0);
140
 
    dpd_buf4_close(&E);
141
 
 
142
 
    dpd_buf4_sort(&Z, CC_TMP0, pqsr, 2, 0, "Z (M>N,JI)");
143
 
    dpd_buf4_init(&Z1, CC_TMP0, 0, 2, 0, 2, 0, 0, "Z (M>N,JI)");
144
 
    dpd_buf4_axpy(&Z1, &Z, -1);
145
 
    dpd_buf4_close(&Z1);
146
 
 
147
 
    dpd_buf4_init(&W, CC2_HET1, 0, 2, 0, 2, 2, 0, "CC2 WMNIJ (M>N,I>J)");
148
 
    dpd_buf4_axpy(&Z, &W, 1);
149
 
    dpd_buf4_close(&Z);
150
 
    dpd_buf4_close(&W);
151
 
 
152
 
    /**** W(m>n,i>j) <-- Zmnij <-- P(i/j)( <mn||ie> * t1[j][e] ) ****/
153
 
    dpd_buf4_init(&Z, CC_TMP0, 0, 12, 10, 12, 10, 0, "Z (m>n,ij)");
154
 
    dpd_buf4_init(&E, CC_EINTS, 0, 12, 30, 12, 30, 0, "E <ij||ka> (i>j,ka)");
155
 
    dpd_contract424(&E, &tia, &Z, 3, 1, 0, 1, 0);
156
 
    dpd_buf4_close(&E);
157
 
 
158
 
    dpd_buf4_sort(&Z, CC_TMP0, pqsr, 12, 10, "Z (m>n,ji)");
159
 
    dpd_buf4_init(&Z1, CC_TMP0, 0, 12, 10, 12, 10, 0, "Z (m>n,ji)");
160
 
    dpd_buf4_axpy(&Z1, &Z, -1);
161
 
    dpd_buf4_close(&Z1);
162
 
 
163
 
    dpd_buf4_init(&W, CC2_HET1, 0, 12, 10, 12, 12, 0, "CC2 Wmnij (m>n,i>j)");
164
 
    dpd_buf4_axpy(&Z, &W, 1);
165
 
    dpd_buf4_close(&Z);
166
 
    dpd_buf4_close(&W);
167
 
 
168
 
    /**** W(Mn,Ij) <-- <Mn|Ie> * t1[j][e] ****/
169
 
    dpd_buf4_init(&W, CC2_HET1, 0, 22, 22, 22, 22, 0, "CC2 WMnIj");
170
 
    dpd_buf4_init(&E, CC_EINTS, 0, 22, 24, 22, 24, 0, "E <Ij|Ka>");
171
 
    dpd_contract424(&E, &tia, &W, 3, 1, 0, 1, 1);
172
 
    dpd_buf4_close(&E);
173
 
    dpd_buf4_close(&W);
174
 
 
175
 
    /**** W(Mn,Ij) <-- <Mn|Ej> * t1[I][E] ****/
176
 
    dpd_buf4_init(&W, CC2_HET1, 0, 22, 22, 22, 22, 0, "CC2 WMnIj");
177
 
    dpd_buf4_init(&E, CC_EINTS, 0, 22, 26, 22, 26, 0, "E <Ij|Ak>");
178
 
    dpd_contract244(&tIA, &E, &W, 1, 2, 1, 1, 1);
179
 
    dpd_buf4_close(&E);
180
 
    dpd_buf4_close(&W);
181
 
 
182
 
    dpd_file2_close(&tIA);
183
 
    dpd_file2_close(&tia);
184
 
  }
185
 
  timer_off("E->Wmnij");
186
 
 
187
 
  timer_on("D->Wmnij");
188
 
  if(params.ref == 0) { /** RHF **/
189
 
 
190
 
    dpd_file2_init(&t1, CC_OEI, 0, 0, 1, "tIA");
191
 
 
192
 
    /* Wmnij<- +1/2 P(ij) t(i,e) t(j,f) * <mn||ef> */
193
 
    dpd_buf4_init(&Z, CC_TMP0, 0, 0, 10, 0, 10, 0, "CC2 ZMnIf (Mn,If)");
194
 
    dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
195
 
    dpd_contract244(&t1, &D, &Z, 1, 2, 1, 1, 0);
196
 
    dpd_buf4_close(&D);
197
 
 
198
 
    dpd_buf4_init(&Z1, CC_TMP0, 0, 0, 0, 0, 0, 0, "CC2 ZMnIj");
199
 
    dpd_contract424(&Z, &t1, &Z1, 3, 1, 0, 0.5, 0);
200
 
    dpd_buf4_close(&Z);
201
 
    dpd_buf4_init(&W, CC2_HET1, 0, 0, 0, 0, 0, 0, "CC2 WMnIj");
202
 
    dpd_buf4_axpy(&Z1, &W, 1);
203
 
    dpd_buf4_close(&W);
204
 
    dpd_buf4_sort_axpy(&Z1, CC2_HET1, qpsr, 0, 0, "CC2 WMnIj", 1);
205
 
    dpd_buf4_close(&Z1);
206
 
 
207
 
    dpd_file2_close(&t1);
208
 
  }
209
 
  else if (params.ref == 1) { /** ROHF **/
210
 
 
211
 
    dpd_file2_init(&tIA, CC_OEI, 0, 0, 1, "tIA");
212
 
    dpd_file2_init(&tia, CC_OEI, 0, 0, 1, "tia");
213
 
 
214
 
    /**** W(M>N,I>J) <-- tIE tJF <MN||EF> ****/
215
 
    dpd_buf4_init(&Z, CC_TMP0, 0, 2, 11, 2, 11, 0, "Z (M>N,EJ)");
216
 
    dpd_buf4_init(&D, CC_DINTS, 0, 2, 5, 2, 7, 0, "D <ij||ab> (i>j,a>b)");
217
 
    dpd_contract424(&D, &tIA, &Z, 3, 1, 0, 1, 0);
218
 
    dpd_buf4_close(&D);
219
 
 
220
 
    dpd_buf4_init(&W, CC2_HET1, 0, 2, 0, 2, 2, 0, "CC2 WMNIJ (M>N,I>J)");
221
 
    dpd_contract244(&tIA, &Z, &W, 1, 2, 1, 1, 1);
222
 
    dpd_buf4_close(&W);
223
 
    dpd_buf4_close(&Z);
224
 
 
225
 
    /**** W(M>N,I>J) <-- tIE tJF <MN||EF> ****/
226
 
    dpd_buf4_init(&Z, CC_TMP0, 0, 2, 11, 2, 11, 0, "Z (m>n,ej)");
227
 
    dpd_buf4_init(&D, CC_DINTS, 0, 2, 5, 2, 7, 0, "D <ij||ab> (i>j,a>b)");
228
 
    dpd_contract424(&D, &tia, &Z, 3, 1, 0, 1, 0);
229
 
    dpd_buf4_close(&D);
230
 
 
231
 
    dpd_buf4_init(&W, CC2_HET1, 0, 2, 0, 2, 2, 0, "CC2 Wmnij (m>n,i>j)");
232
 
    dpd_contract244(&tia, &Z, &W, 1, 2, 1, 1, 1);
233
 
    dpd_buf4_close(&W);
234
 
    dpd_buf4_close(&Z);
235
 
 
236
 
    /**** W(Mn,Ij) <-- tIE tjf <Mn|Ef> ****/
237
 
    dpd_buf4_init(&Z, CC_TMP0, 0, 0, 11, 0, 11, 0, "Z (Mn,Ej)");
238
 
    dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
239
 
    dpd_contract424(&D, &tia, &Z, 3, 1, 0, 1, 0);
240
 
    dpd_buf4_close(&D);
241
 
 
242
 
    dpd_buf4_init(&W, CC2_HET1, 0, 0, 0, 0, 0, 0, "CC2 WMnIj");
243
 
    dpd_contract244(&tIA, &Z, &W, 1, 2, 1, 1, 1);
244
 
    dpd_buf4_close(&W);
245
 
    dpd_buf4_close(&Z);
246
 
 
247
 
    dpd_file2_close(&tIA);
248
 
    dpd_file2_close(&tia);
249
 
 
250
 
    purge_cc2_Wmnij();
251
 
  }
252
 
  else if (params.ref == 2) { /** UHF **/
253
 
 
254
 
    dpd_file2_init(&tIA, CC_OEI, 0, 0, 1, "tIA");
255
 
    dpd_file2_init(&tia, CC_OEI, 0, 2, 3, "tia");
256
 
 
257
 
    /**** W(M>N,I>J) <-- tIE tJF <MN||EF> ****/
258
 
    dpd_buf4_init(&Z, CC_TMP0, 0, 2, 21, 2, 21, 0, "Z (M>N,EJ)");
259
 
    dpd_buf4_init(&D, CC_DINTS, 0, 2, 5, 2, 7, 0, "D <IJ||AB> (I>J,A>B)");
260
 
    dpd_contract424(&D, &tIA, &Z, 3, 1, 0, 1, 0);
261
 
    dpd_buf4_close(&D);
262
 
 
263
 
    dpd_buf4_init(&W, CC2_HET1, 0, 2, 0, 2, 2, 0, "CC2 WMNIJ (M>N,I>J)");
264
 
    dpd_contract244(&tIA, &Z, &W, 1, 2, 1, 1, 1);
265
 
    dpd_buf4_close(&W);
266
 
    dpd_buf4_close(&Z);
267
 
 
268
 
    /**** W(M>N,I>J) <-- tIE tJF <MN||EF> ****/
269
 
    dpd_buf4_init(&Z, CC_TMP0, 0, 12, 31, 12, 31, 0, "Z (m>n,ej)");
270
 
    dpd_buf4_init(&D, CC_DINTS, 0, 12, 15, 12, 17, 0, "D <ij||ab> (i>j,a>b)");
271
 
    dpd_contract424(&D, &tia, &Z, 3, 1, 0, 1, 0);
272
 
    dpd_buf4_close(&D);
273
 
 
274
 
    dpd_buf4_init(&W, CC2_HET1, 0, 12, 10, 12, 12, 0, "CC2 Wmnij (m>n,i>j)");
275
 
    dpd_contract244(&tia, &Z, &W, 1, 2, 1, 1, 1);
276
 
    dpd_buf4_close(&W);
277
 
    dpd_buf4_close(&Z);
278
 
 
279
 
    /**** W(Mn,Ij) <-- tIE tjf <Mn|Ef> ****/
280
 
    dpd_buf4_init(&Z, CC_TMP0, 0, 22, 26, 22, 26, 0, "Z (Mn,Ej)");
281
 
    dpd_buf4_init(&D, CC_DINTS, 0, 22, 28, 22, 28, 0, "D <Ij|Ab>");
282
 
    dpd_contract424(&D, &tia, &Z, 3, 1, 0, 1, 0);
283
 
    dpd_buf4_close(&D);
284
 
 
285
 
    dpd_buf4_init(&W, CC2_HET1, 0, 22, 22, 22, 22, 0, "CC2 WMnIj");
286
 
    dpd_contract244(&tIA, &Z, &W, 1, 2, 1, 1, 1);
287
 
    dpd_buf4_close(&W);
288
 
    dpd_buf4_close(&Z);
289
 
 
290
 
    dpd_file2_close(&tIA);
291
 
    dpd_file2_close(&tia);
292
 
  }
293
 
  timer_off("D->Wmnij");
294
 
 
295
 
}
296
 
 
297
 
 
298
 
void purge_cc2_Wmnij(void) {
299
 
  dpdfile2 FAE, Fmi, FME, Fme;
300
 
  dpdfile4 W;
301
 
  int *occpi, *virtpi;
302
 
  int h, a, b, e, f, i, j, m, n, omit;
303
 
  int    A, B, E, F, I, J, M, N;
304
 
  int mn, ei, ma, ef, me, jb, mb, ij, ab;
305
 
  int asym, bsym, esym, fsym, isym, jsym, msym, nsym;
306
 
  int *occ_off, *vir_off;
307
 
  int *occ_sym, *vir_sym;
308
 
  int *openpi, nirreps;
309
 
 
310
 
  nirreps = moinfo.nirreps;
311
 
  occpi = moinfo.occpi; virtpi = moinfo.virtpi;
312
 
  occ_off = moinfo.occ_off; vir_off = moinfo.vir_off;
313
 
  occ_sym = moinfo.occ_sym; vir_sym = moinfo.vir_sym;
314
 
  openpi = moinfo.openpi;
315
 
 
316
 
  /* Purge Wmnij matrix elements */
317
 
  dpd_file4_init(&W, CC2_HET1, 0, 2, 2,"CC2 Wmnij (m>n,i>j)");
318
 
  for(h=0; h < nirreps; h++) {
319
 
    dpd_file4_mat_irrep_init(&W, h);
320
 
    dpd_file4_mat_irrep_rd(&W, h);
321
 
    for(mn=0; mn < W.params->rowtot[h]; mn++) {
322
 
      m = W.params->roworb[h][mn][0];
323
 
      n = W.params->roworb[h][mn][1];
324
 
      msym = W.params->psym[m];
325
 
      nsym = W.params->qsym[n];
326
 
      M = m - occ_off[msym];
327
 
      N = n - occ_off[nsym];
328
 
      for(ij=0; ij < W.params->coltot[h]; ij++) {
329
 
        i = W.params->colorb[h][ij][0];
330
 
        j = W.params->colorb[h][ij][1];
331
 
        isym = W.params->rsym[i];
332
 
        jsym = W.params->ssym[j];
333
 
        I = i - occ_off[isym];
334
 
        J = j - occ_off[jsym];
335
 
        if ((I >= (occpi[isym] - openpi[isym])) ||
336
 
            (J >= (occpi[jsym] - openpi[jsym])) ||
337
 
            (M >= (occpi[msym] - openpi[msym])) ||
338
 
            (N >= (occpi[nsym] - openpi[nsym])) )
339
 
          W.matrix[h][mn][ij] = 0.0;
340
 
      }
341
 
    }
342
 
    dpd_file4_mat_irrep_wrt(&W, h);
343
 
    dpd_file4_mat_irrep_close(&W, h);
344
 
  }
345
 
  dpd_file4_close(&W);
346
 
 
347
 
  dpd_file4_init(&W, CC2_HET1, 0, 0, 0,"CC2 WMnIj");
348
 
  for(h=0; h < nirreps; h++) {
349
 
    dpd_file4_mat_irrep_init(&W, h);
350
 
    dpd_file4_mat_irrep_rd(&W, h);
351
 
    for(mn=0; mn < W.params->rowtot[h]; mn++) {
352
 
      n = W.params->roworb[h][mn][1];
353
 
      nsym = W.params->qsym[n];
354
 
      N = n - occ_off[nsym];
355
 
      for(ij=0; ij < W.params->coltot[h]; ij++) {
356
 
        j = W.params->colorb[h][ij][1];
357
 
        jsym = W.params->ssym[j];
358
 
        J = j - occ_off[jsym];
359
 
        if ((J >= (occpi[jsym] - openpi[jsym])) ||
360
 
            (N >= (occpi[nsym] - openpi[nsym])) )
361
 
          W.matrix[h][mn][ij] = 0.0;
362
 
      }
363
 
    }
364
 
    dpd_file4_mat_irrep_wrt(&W, h);
365
 
    dpd_file4_mat_irrep_close(&W, h);
366
 
  }
367
 
  dpd_file4_close(&W);
368
 
}