← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cclambda/L1.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/cclambda/L1.c
1
 
#include <stdio.h>
2
 
#include <stdlib.h>
3
 
#include <string.h>
4
 
#include <libdpd/dpd.h>
5
 
#include <libqt/qt.h>
6
 
#define EXTERN
7
 
#include "globals.h"
8
 
 
9
 
void local_filter_T1(dpdfile2 *);
10
 
 
11
 
void L1_build(struct L_Params L_params) {
12
 
  dpdfile2 newLIA, newLia, LIA, Lia;
13
 
  dpdfile2 dIA, dia, Fme, FME;
14
 
  dpdfile2 LFaet2, LFAEt2, LFmit2, LFMIt2;
15
 
  dpdfile2 GMI, Gmi, Gae, XIA, Xia;
16
 
  dpdfile2 GAE;
17
 
  dpdbuf4 WMBEJ, Wmbej, WMbEj, WmBeJ;
18
 
  dpdbuf4 WMBIJ, Wmbij, WMbIj, WmBiJ;
19
 
  dpdbuf4 LIJAB, Lijab, LIjAb, LiJaB, L2;
20
 
  dpdbuf4 WMNIE, Wmnie, WMnIe, WmNiE;
21
 
  dpdbuf4 WAMEF, Wamef, WAmEf, WaMeF, W;
22
 
  dpdbuf4 Z, D;
23
 
  dpdfile2 XLD;
24
 
  int Gim, Gi, Gm, Ga, Gam, nrows, ncols, A, a, am;
25
 
  int Gei, ei, e, i, Gef, Ge, Gf, E, I, af, fa, f;
26
 
  double *X;
27
 
  int L_irr;
28
 
  L_irr = L_params.irrep;
29
 
 
30
 
  /* ground state inhomogeneous term is Fme */
31
 
  if (L_params.ground) {
32
 
    if(params.ref == 0) {
33
 
      dpd_file2_init(&FME,CC_OEI, 0, 0, 1, "FME");
34
 
      dpd_file2_copy(&FME, CC_LAMBDA, "New LIA");
35
 
      dpd_file2_close(&FME);
36
 
    }
37
 
    else if(params.ref == 1) {
38
 
      dpd_file2_init(&Fme,CC_OEI, 0, 0, 1, "Fme");
39
 
      dpd_file2_init(&FME,CC_OEI, 0, 0, 1, "FME");
40
 
      dpd_file2_copy(&Fme, CC_LAMBDA, "New Lia");
41
 
      dpd_file2_copy(&FME, CC_LAMBDA, "New LIA");
42
 
      dpd_file2_close(&Fme);
43
 
      dpd_file2_close(&FME);
44
 
    }
45
 
    else if(params.ref == 2) {
46
 
      dpd_file2_init(&Fme,CC_OEI, 0, 2, 3, "Fme");
47
 
      dpd_file2_init(&FME,CC_OEI, 0, 0, 1, "FME");
48
 
      dpd_file2_copy(&Fme, CC_LAMBDA, "New Lia");
49
 
      dpd_file2_copy(&FME, CC_LAMBDA, "New LIA");
50
 
      dpd_file2_close(&Fme);
51
 
      dpd_file2_close(&FME);
52
 
    }
53
 
  }
54
 
  /* excited state - no inhomogenous term, first term is -energy*L*/
55
 
  else if (!params.zeta) {
56
 
    if (params.ref == 0 || params.ref == 1) {
57
 
      dpd_file2_init(&LIA, CC_LAMBDA, L_irr, 0, 1, "LIA");
58
 
      dpd_file2_init(&newLIA, CC_LAMBDA, L_irr, 0, 1, "New LIA");
59
 
      dpd_file2_init(&Lia, CC_LAMBDA, L_irr, 0, 1, "Lia");
60
 
      dpd_file2_init(&newLia, CC_LAMBDA, L_irr, 0, 1, "New Lia");
61
 
      dpd_file2_axpy(&LIA, &newLIA, -1.0 * L_params.cceom_energy,0.0);
62
 
      dpd_file2_axpy(&Lia, &newLia, -1.0 * L_params.cceom_energy,0.0);
63
 
      dpd_file2_close(&LIA);
64
 
      dpd_file2_close(&newLIA);
65
 
      dpd_file2_close(&Lia);
66
 
      dpd_file2_close(&newLia);
67
 
    }
68
 
    else if (params.ref == 2) {
69
 
      /* do nothing - TDC did not change to increments for the UHF case */
70
 
    }
71
 
  }
72
 
  /* solving zeta equations; inhomogeneous term is Xi */
73
 
  else {
74
 
    if (params.ref == 0) {
75
 
      dpd_file2_init(&XIA, EOM_XI, 0, 0, 1, "XIA");
76
 
      dpd_file2_copy(&XIA, CC_LAMBDA, "New LIA");
77
 
      dpd_file2_close(&XIA);
78
 
    }
79
 
    else if (params.ref == 1) {
80
 
      dpd_file2_init(&XIA, EOM_XI, 0, 0, 1, "XIA");
81
 
      dpd_file2_init(&Xia, EOM_XI, 0, 0, 1, "Xia");
82
 
      dpd_file2_copy(&XIA, CC_LAMBDA, "New LIA");
83
 
      dpd_file2_copy(&Xia, CC_LAMBDA, "New Lia");
84
 
      dpd_file2_close(&XIA);
85
 
      dpd_file2_close(&Xia);
86
 
    }
87
 
    else if(params.ref == 2) {
88
 
      dpd_file2_init(&XIA, EOM_XI, 0, 0, 1, "XIA");
89
 
      dpd_file2_init(&Xia, EOM_XI, 0, 2, 3, "Xia");
90
 
      dpd_file2_copy(&XIA, CC_LAMBDA, "New LIA");
91
 
      dpd_file2_copy(&Xia, CC_LAMBDA, "New Lia");
92
 
      dpd_file2_close(&XIA);
93
 
      dpd_file2_close(&Xia);
94
 
    }
95
 
  }
96
 
 
97
 
  if(params.ref == 0 || params.ref == 1) {
98
 
    dpd_file2_init(&newLIA, CC_LAMBDA, L_irr, 0, 1, "New LIA");
99
 
    dpd_file2_init(&newLia, CC_LAMBDA, L_irr, 0, 1, "New Lia");
100
 
  }
101
 
  else if(params.ref == 2) {
102
 
    dpd_file2_init(&newLIA, CC_LAMBDA, L_irr, 0, 1, "New LIA");
103
 
    dpd_file2_init(&newLia, CC_LAMBDA, L_irr, 2, 3, "New Lia");
104
 
  }
105
 
 
106
 
  if(params.ref == 0) { /** RHF **/
107
 
    dpd_file2_init(&LIA, CC_LAMBDA, L_irr, 0, 1, "LIA");
108
 
 
109
 
    /* L1 RHS += Lie*Fea */
110
 
    dpd_file2_init(&LFAEt2, CC_OEI, 0, 1, 1, "FAE");
111
 
    dpd_contract222(&LIA,&LFAEt2,&newLIA, 0, 1, 1.0, 1.0);
112
 
    dpd_file2_close(&LFAEt2);
113
 
 
114
 
    /* L1 RHS += -Lma*Fim */
115
 
    dpd_file2_init(&LFMIt2,CC_OEI, 0, 0, 0, "FMI");
116
 
    dpd_contract222(&LFMIt2,&LIA,&newLIA, 0, 1, -1.0, 1.0);
117
 
    dpd_file2_close(&LFMIt2);
118
 
 
119
 
    /* L1 RHS += Lme*Wieam */
120
 
    dpd_buf4_init(&W, CC_HBAR, 0, 10, 10, 10, 10, 0, "2 W(ME,jb) + W(Me,Jb)");
121
 
    dpd_contract422(&W, &LIA, &newLIA, 0, 0, 1.0, 1.0);
122
 
    dpd_buf4_close(&W);
123
 
 
124
 
    dpd_file2_close(&LIA);
125
 
  }
126
 
  else if(params.ref == 1) { /** ROHF **/
127
 
 
128
 
    /* L1 RHS += Lie*Fea */
129
 
    dpd_file2_init(&LIA, CC_LAMBDA, L_irr, 0, 1, "LIA");
130
 
    dpd_file2_init(&Lia, CC_LAMBDA, L_irr, 0, 1, "Lia");
131
 
 
132
 
    dpd_file2_init(&LFAEt2, CC_OEI, 0, 1, 1, "FAE");
133
 
    dpd_file2_init(&LFaet2, CC_OEI, 0, 1, 1, "Fae");
134
 
    dpd_contract222(&Lia,&LFaet2,&newLia, 0, 1, 1.0, 1.0);
135
 
    dpd_contract222(&LIA,&LFAEt2,&newLIA, 0, 1, 1.0, 1.0);
136
 
    dpd_file2_close(&LFaet2);
137
 
    dpd_file2_close(&LFAEt2);
138
 
 
139
 
    /* L1 RHS += -Lma*Fim */
140
 
    dpd_file2_init(&LFMIt2,CC_OEI, 0, 0, 0, "FMI");
141
 
    dpd_file2_init(&LFmit2,CC_OEI, 0, 0, 0, "Fmi");
142
 
    dpd_contract222(&LFmit2,&Lia,&newLia, 0, 1, -1.0, 1.0);
143
 
    dpd_contract222(&LFMIt2,&LIA,&newLIA, 0, 1, -1.0, 1.0);
144
 
    dpd_file2_close(&LFmit2);
145
 
    dpd_file2_close(&LFMIt2);
146
 
 
147
 
    /* L1 RHS += Lme*Wieam */
148
 
    dpd_buf4_init(&WMBEJ, CC_HBAR, 0, 10, 10, 10, 10, 0, "WMBEJ");
149
 
    dpd_contract422(&WMBEJ, &LIA, &newLIA, 0, 0, 1.0, 1.0);
150
 
    dpd_buf4_close(&WMBEJ);
151
 
 
152
 
    dpd_buf4_init(&WMbEj, CC_HBAR, 0, 10, 10, 10, 10, 0, "WMbEj");
153
 
    dpd_contract422(&WMbEj, &Lia, &newLIA, 0, 0, 1.0, 1.0);
154
 
    dpd_buf4_close(&WMbEj);
155
 
 
156
 
    dpd_buf4_init(&Wmbej, CC_HBAR, 0, 10, 10, 10, 10, 0, "Wmbej");
157
 
    dpd_contract422(&Wmbej, &Lia, &newLia, 0, 0, 1.0, 1.0);
158
 
    dpd_buf4_close(&Wmbej);
159
 
 
160
 
    dpd_buf4_init(&WmBeJ, CC_HBAR, 0, 10, 10, 10, 10, 0, "WmBeJ");
161
 
    dpd_contract422(&WmBeJ, &LIA, &newLia, 0, 0, 1.0, 1.0);
162
 
    dpd_buf4_close(&WmBeJ);
163
 
 
164
 
    dpd_file2_close(&LIA);
165
 
    dpd_file2_close(&Lia);
166
 
  }
167
 
  else if(params.ref == 2) { /** UHF **/
168
 
 
169
 
    /* L1 RHS += Lie*Fea */
170
 
    dpd_file2_init(&LIA, CC_LAMBDA, L_irr, 0, 1, "LIA");
171
 
    dpd_file2_init(&Lia, CC_LAMBDA, L_irr, 2, 3, "Lia");
172
 
 
173
 
    dpd_file2_init(&LFAEt2, CC_OEI, 0, 1, 1, "FAEt");
174
 
    dpd_file2_init(&LFaet2, CC_OEI, 0, 3, 3, "Faet");
175
 
    dpd_contract222(&Lia,&LFaet2,&newLia, 0, 1, 1, 1);
176
 
    dpd_contract222(&LIA,&LFAEt2,&newLIA, 0, 1, 1, 1);
177
 
    dpd_file2_close(&LFaet2);
178
 
    dpd_file2_close(&LFAEt2);
179
 
 
180
 
    /* L1 RHS += -Lma*Fim */
181
 
    dpd_file2_init(&LFMIt2,CC_OEI, 0, 0, 0, "FMIt");
182
 
    dpd_file2_init(&LFmit2,CC_OEI, 0, 2, 2, "Fmit");
183
 
    dpd_contract222(&LFmit2,&Lia,&newLia, 0, 1, -1, 1);
184
 
    dpd_contract222(&LFMIt2,&LIA,&newLIA, 0, 1, -1, 1);
185
 
    dpd_file2_close(&LFmit2);
186
 
    dpd_file2_close(&LFMIt2);
187
 
 
188
 
    /* L1 RHS += Lme*Wieam */
189
 
    dpd_buf4_init(&WMBEJ, CC_HBAR, 0, 20, 20, 20, 20, 0, "WMBEJ");
190
 
    dpd_contract422(&WMBEJ, &LIA, &newLIA, 0, 0, 1, 1);
191
 
    dpd_buf4_close(&WMBEJ);
192
 
 
193
 
    dpd_buf4_init(&WMbEj, CC_HBAR, 0, 20, 30, 20, 30, 0, "WMbEj");
194
 
    dpd_contract422(&WMbEj, &Lia, &newLIA, 0, 0, 1, 1);
195
 
    dpd_buf4_close(&WMbEj);
196
 
 
197
 
    dpd_buf4_init(&Wmbej, CC_HBAR, 0, 30, 30, 30, 30, 0, "Wmbej");
198
 
    dpd_contract422(&Wmbej, &Lia, &newLia, 0, 0, 1, 1);
199
 
    dpd_buf4_close(&Wmbej);
200
 
 
201
 
    dpd_buf4_init(&WmBeJ, CC_HBAR, 0, 30, 20, 30, 20, 0, "WmBeJ");
202
 
    dpd_contract422(&WmBeJ, &LIA, &newLia, 0, 0, 1, 1);
203
 
    dpd_buf4_close(&WmBeJ);
204
 
 
205
 
    dpd_file2_close(&LIA);
206
 
    dpd_file2_close(&Lia);
207
 
  }
208
 
 
209
 
  /* L1 RHS += 1/2 Limef*Wefam */
210
 
  /* L(i,a) += [ 2 L(im,ef) - L(im,fe) ] * W(am,ef) */
211
 
  /* Note: W(am,ef) is really Wabei (ei,ab) */
212
 
  if(params.ref == 0) { /** RHF **/
213
 
 
214
 
    dpd_buf4_init(&W, CC_HBAR, 0, 11, 5, 11, 5, 0, "WAbEi (Ei,Ab)");
215
 
    dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "2 LIjAb - LIjBa");
216
 
    /*    dpd_contract442(&L2, &W, &newLIA, 0, 0, 1.0, 1.0); */
217
 
    dpd_file2_mat_init(&newLIA);
218
 
    dpd_file2_mat_rd(&newLIA);
219
 
    for(Gam=0; Gam < moinfo.nirreps; Gam++) {
220
 
      Gef = Gam; /* W is totally symmetric */
221
 
      Gim = Gef ^ L_irr; 
222
 
      dpd_buf4_mat_irrep_init(&L2, Gim);
223
 
      dpd_buf4_mat_irrep_rd(&L2, Gim);
224
 
      dpd_buf4_mat_irrep_shift13(&L2, Gim);
225
 
 
226
 
      for(Gi=0; Gi < moinfo.nirreps; Gi++) {
227
 
        Ga = Gi ^ L_irr;
228
 
        Gm = Ga ^ Gam;
229
 
        W.matrix[Gam] = dpd_block_matrix(moinfo.occpi[Gm],W.params->coltot[Gam]);
230
 
 
231
 
        nrows = moinfo.occpi[Gi];
232
 
        ncols = moinfo.occpi[Gm] * W.params->coltot[Gam];
233
 
 
234
 
        for(A=0; A < moinfo.virtpi[Ga]; A++) {
235
 
          a = moinfo.vir_off[Ga] + A;
236
 
          am = W.row_offset[Gam][a];
237
 
 
238
 
          dpd_buf4_mat_irrep_rd_block(&W, Gam, am, moinfo.occpi[Gm]);
239
 
 
240
 
          if(nrows && ncols && moinfo.virtpi[Ga]) 
241
 
            C_DGEMV('n',nrows,ncols,1,L2.shift.matrix[Gim][Gi][0],ncols,W.matrix[Gam][0],1,
242
 
                    1, &(newLIA.matrix[Gi][0][A]), moinfo.virtpi[Ga]);
243
 
 
244
 
        }
245
 
 
246
 
        dpd_free_block(W.matrix[Gam], moinfo.occpi[Gm], W.params->coltot[Gam]);
247
 
      }
248
 
      dpd_buf4_mat_irrep_close(&L2, Gim);
249
 
    }
250
 
    dpd_file2_mat_wrt(&newLIA);
251
 
    dpd_file2_mat_close(&newLIA);
252
 
    dpd_buf4_close(&L2);
253
 
    dpd_buf4_close(&W);
254
 
 
255
 
  }
256
 
  else if(params.ref == 1) { /** ROHF **/
257
 
 
258
 
    dpd_buf4_init(&W, CC_HBAR, 0, 11, 7, 11, 7, 0, "WEIAB");
259
 
    dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 0, 7, 2, 7, 0, "LIJAB");
260
 
    dpd_contract442(&L2, &W, &newLIA, 0, 0, 1.0, 1.0);
261
 
    dpd_buf4_close(&W);
262
 
    dpd_buf4_close(&L2);
263
 
    dpd_buf4_init(&W, CC_HBAR, 0, 11, 5, 11, 5, 0, "WEiAb");
264
 
    dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "LIjAb");
265
 
    dpd_contract442(&L2, &W, &newLIA, 0, 0, 1.0, 1.0);
266
 
    dpd_buf4_close(&W);
267
 
    dpd_buf4_close(&L2);
268
 
 
269
 
    dpd_buf4_init(&W, CC_HBAR, 0, 11, 7, 11, 7, 0, "Weiab");
270
 
    dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 0, 7, 2, 7, 0, "Lijab");
271
 
    dpd_contract442(&L2, &W, &newLia, 0, 0, 1.0, 1.0);
272
 
    dpd_buf4_close(&W);
273
 
    dpd_buf4_close(&L2);
274
 
    dpd_buf4_init(&W, CC_HBAR, 0, 11, 5, 11, 5, 0, "WeIaB");
275
 
    dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "LiJaB");
276
 
    dpd_contract442(&L2, &W, &newLia, 0, 0, 1.0, 1.0);
277
 
    dpd_buf4_close(&W);
278
 
    dpd_buf4_close(&L2);
279
 
  }
280
 
  else if(params.ref == 2) {
281
 
 
282
 
    dpd_buf4_init(&W, CC_HBAR, 0, 21, 7, 21, 7, 0, "WEIAB");
283
 
    dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 0, 7, 2, 7, 0, "LIJAB");
284
 
    dpd_contract442(&L2, &W, &newLIA, 0, 0, 1, 1);
285
 
    dpd_buf4_close(&W);
286
 
    dpd_buf4_close(&L2);
287
 
    dpd_buf4_init(&W, CC_HBAR, 0, 26, 28, 26, 28, 0, "WEiAb");
288
 
    dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 22, 28, 22, 28, 0, "LIjAb");
289
 
    dpd_contract442(&L2, &W, &newLIA, 0, 0, 1, 1);
290
 
    dpd_buf4_close(&W);
291
 
    dpd_buf4_close(&L2);
292
 
 
293
 
    dpd_buf4_init(&W, CC_HBAR, 0, 31, 17, 31, 17, 0, "Weiab");
294
 
    dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 10, 17, 12, 17, 0, "Lijab");
295
 
    dpd_contract442(&L2, &W, &newLia, 0, 0, 1, 1);
296
 
    dpd_buf4_close(&W);
297
 
    dpd_buf4_close(&L2);
298
 
    dpd_buf4_init(&W, CC_HBAR, 0, 25, 29, 25, 29, 0, "WeIaB");
299
 
    dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 23, 29, 23, 29, 0, "LiJaB");
300
 
    dpd_contract442(&L2, &W, &newLia, 0, 0, 1, 1);
301
 
    dpd_buf4_close(&W);
302
 
    dpd_buf4_close(&L2);
303
 
 
304
 
  }
305
 
 
306
 
  /* L1 RHS += -1/2 Lmnae*Wiemn */
307
 
  if(params.ref == 0) {
308
 
    dpd_buf4_init(&WMbIj, CC_HBAR, 0, 10, 0, 10, 0, 0, "WMbIj");
309
 
    dpd_buf4_init(&LIjAb, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "2 LIjAb - LIjBa");
310
 
    dpd_contract442(&WMbIj, &LIjAb, &newLIA, 0, 2, -1.0, 1.0);
311
 
    dpd_buf4_close(&LIjAb);
312
 
    dpd_buf4_close(&WMbIj);
313
 
  }
314
 
  else if(params.ref == 1) {
315
 
 
316
 
    dpd_buf4_init(&WMBIJ, CC_HBAR, 0, 10, 2, 10, 2, 0, "WMBIJ");
317
 
    dpd_buf4_init(&LIJAB, CC_LAMBDA, L_irr, 2, 5, 2, 7, 0, "LIJAB");
318
 
    dpd_contract442(&WMBIJ, &LIJAB, &newLIA, 0, 2, -1.0, 1.0);
319
 
    dpd_buf4_close(&LIJAB);
320
 
    dpd_buf4_close(&WMBIJ);
321
 
 
322
 
    dpd_buf4_init(&WMbIj, CC_HBAR, 0, 10, 0, 10, 0, 0, "WMbIj");
323
 
    dpd_buf4_init(&LIjAb, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "LIjAb");
324
 
    dpd_contract442(&WMbIj, &LIjAb, &newLIA, 0, 2, -1.0, 1.0);
325
 
    dpd_buf4_close(&LIjAb);
326
 
    dpd_buf4_close(&WMbIj);
327
 
 
328
 
    dpd_buf4_init(&Wmbij, CC_HBAR, 0, 10, 2, 10, 2, 0, "Wmbij");
329
 
    dpd_buf4_init(&Lijab, CC_LAMBDA, L_irr, 2, 5, 2, 7, 0, "Lijab");
330
 
    dpd_contract442(&Wmbij, &Lijab, &newLia, 0, 2, -1.0, 1.0);
331
 
    dpd_buf4_close(&Lijab);
332
 
    dpd_buf4_close(&Wmbij);
333
 
 
334
 
    dpd_buf4_init(&WmBiJ, CC_HBAR, 0, 10, 0, 10, 0, 0, "WmBiJ");
335
 
    dpd_buf4_init(&LiJaB, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "LiJaB");
336
 
    dpd_contract442(&WmBiJ, &LiJaB, &newLia, 0, 2, -1.0, 1.0);
337
 
    dpd_buf4_close(&LiJaB);
338
 
    dpd_buf4_close(&WmBiJ);
339
 
  }
340
 
  else if(params.ref == 2) {
341
 
 
342
 
    dpd_buf4_init(&WMBIJ, CC_HBAR, 0, 20, 2, 20, 2, 0, "WMBIJ");
343
 
    dpd_buf4_init(&LIJAB, CC_LAMBDA, L_irr, 2, 5, 2, 7, 0, "LIJAB");
344
 
    dpd_contract442(&WMBIJ, &LIJAB, &newLIA, 0, 2, -1, 1);
345
 
    dpd_buf4_close(&LIJAB);
346
 
    dpd_buf4_close(&WMBIJ);
347
 
 
348
 
    dpd_buf4_init(&WMbIj, CC_HBAR, 0, 24, 22, 24, 22, 0, "WMbIj");
349
 
    dpd_buf4_init(&LIjAb, CC_LAMBDA, L_irr, 22, 28, 22, 28, 0, "LIjAb");
350
 
    dpd_contract442(&WMbIj, &LIjAb, &newLIA, 0, 2, -1, 1);
351
 
    dpd_buf4_close(&LIjAb);
352
 
    dpd_buf4_close(&WMbIj);
353
 
 
354
 
    dpd_buf4_init(&Wmbij, CC_HBAR, 0, 30, 12, 30, 12, 0, "Wmbij");
355
 
    dpd_buf4_init(&Lijab, CC_LAMBDA, L_irr, 12, 15, 12, 17, 0, "Lijab");
356
 
    dpd_contract442(&Wmbij, &Lijab, &newLia, 0, 2, -1, 1);
357
 
    dpd_buf4_close(&Lijab);
358
 
    dpd_buf4_close(&Wmbij);
359
 
 
360
 
    dpd_buf4_init(&WmBiJ, CC_HBAR, 0, 27, 23, 27, 23, 0, "WmBiJ");
361
 
    dpd_buf4_init(&LiJaB, CC_LAMBDA, L_irr, 23, 29, 23, 29, 0, "LiJaB");
362
 
    dpd_contract442(&WmBiJ, &LiJaB, &newLia, 0, 2, -1, 1);
363
 
    dpd_buf4_close(&LiJaB);
364
 
    dpd_buf4_close(&WmBiJ);
365
 
  }
366
 
 
367
 
 
368
 
  /* L1 RHS += -Gef*Weifa */
369
 
  if(params.ref == 0) {
370
 
 
371
 
/*     dpd_file2_init(&GAE, CC_LAMBDA, L_irr, 1, 1, "GAE"); */
372
 
/*     dpd_buf4_init(&WaMeF, CC_HBAR, 0, 11, 5, 11, 5, 0, "WAmEf 2(Am,Ef) - (Am,fE)"); */
373
 
/*     dpd_dot13(&GAE,&WaMeF,&newLIA, 0, 0, -1.0, 1.0); */
374
 
/*     dpd_buf4_close(&WaMeF); */
375
 
/*     dpd_file2_close(&GAE); */
376
 
 
377
 
    /* Above code replaced to remove disk-space and memory bottlenecks 7/26/05, -TDC */
378
 
    dpd_file2_init(&GAE, CC_LAMBDA, L_irr, 1, 1, "GAE");
379
 
    dpd_file2_mat_init(&GAE);
380
 
    dpd_file2_mat_rd(&GAE);
381
 
    dpd_file2_mat_init(&newLIA);
382
 
    dpd_file2_mat_rd(&newLIA);
383
 
    dpd_buf4_init(&W, CC_HBAR, 0, 11, 5, 11, 5, 0, "WAmEf");
384
 
    for(Gei=0; Gei < moinfo.nirreps; Gei++) {
385
 
      dpd_buf4_mat_irrep_row_init(&W, Gei);
386
 
      X = init_array(W.params->coltot[Gei]);
387
 
      for(ei=0; ei < W.params->rowtot[Gei]; ei++) {
388
 
        dpd_buf4_mat_irrep_row_rd(&W, Gei, ei);
389
 
        e = W.params->roworb[Gei][ei][0];
390
 
        i = W.params->roworb[Gei][ei][1];
391
 
        Ge = W.params->psym[e];
392
 
        Gf = Ge ^ L_irr;
393
 
        Gi = Ge ^ Gei;
394
 
        Ga = Gi ^ L_irr;
395
 
        E = e - moinfo.vir_off[Ge];
396
 
        I = i - moinfo.occ_off[Gi];
397
 
 
398
 
        zero_arr(X,W.params->coltot[Gei]);
399
 
 
400
 
        for(fa=0; fa < W.params->coltot[Gei]; fa++) {
401
 
          f = W.params->colorb[Gei][fa][0];
402
 
          a = W.params->colorb[Gei][fa][1];
403
 
          af = W.params->colidx[a][f];
404
 
          X[fa] = 2.0 * W.matrix[Gei][0][fa] - W.matrix[Gei][0][af];
405
 
        }
406
 
 
407
 
        nrows = moinfo.virtpi[Gf];
408
 
        ncols = moinfo.virtpi[Ga];
409
 
        if(nrows && ncols)
410
 
          C_DGEMV('t',nrows,ncols,-1,&X[W.col_offset[Gei][Gf]],ncols,
411
 
                  GAE.matrix[Ge][E],1,1,newLIA.matrix[Gi][I],1);
412
 
 
413
 
      }
414
 
      dpd_buf4_mat_irrep_row_close(&W, Gei);
415
 
      free(X);
416
 
    }
417
 
    dpd_buf4_close(&W);
418
 
    dpd_file2_mat_wrt(&newLIA);
419
 
    dpd_file2_mat_close(&newLIA);
420
 
    dpd_file2_mat_close(&GAE);
421
 
    dpd_file2_close(&GAE);
422
 
  }
423
 
  else if(params.ref == 1) {
424
 
 
425
 
    dpd_file2_init(&GAE, CC_LAMBDA, L_irr, 1, 1, "GAE");
426
 
    dpd_file2_init(&Gae, CC_LAMBDA, L_irr, 1, 1, "Gae");
427
 
 
428
 
    dpd_buf4_init(&WAMEF, CC_HBAR, 0, 11, 5, 11, 7, 0, "WAMEF");
429
 
    dpd_dot13(&GAE,&WAMEF,&newLIA, 0, 0, -1.0, 1.0);
430
 
    dpd_buf4_close(&WAMEF);
431
 
 
432
 
    dpd_buf4_init(&WaMeF, CC_HBAR, 0, 11, 5, 11, 5, 0, "WaMeF");
433
 
    dpd_dot13(&Gae,&WaMeF,&newLIA, 0, 0, -1.0, 1.0);
434
 
    dpd_buf4_close(&WaMeF);
435
 
 
436
 
    dpd_buf4_init(&Wamef, CC_HBAR, 0, 11, 5, 11, 7, 0, "Wamef");
437
 
    dpd_dot13(&Gae,&Wamef,&newLia, 0, 0, -1.0, 1.0);
438
 
    dpd_buf4_close(&Wamef);
439
 
 
440
 
    dpd_buf4_init(&WAmEf, CC_HBAR, 0, 11, 5, 11, 5, 0, "WAmEf");
441
 
    dpd_dot13(&GAE,&WAmEf,&newLia, 0, 0, -1.0, 1.0);
442
 
    dpd_buf4_close(&WAmEf);
443
 
 
444
 
    /*
445
 
      dpd_buf4_init(&WAMEF, CC_HBAR, 0, 10, 5, 10, 7, 0, "WAMEF");
446
 
      dpd_dot23(&GAE,&WAMEF,&newLIA, 0, 0, -1.0, 1.0);
447
 
      dpd_buf4_close(&WAMEF);
448
 
      dpd_buf4_init(&WaMeF, CC_HBAR, 0, 10, 5, 10, 5, 0, "WaMeF");
449
 
      dpd_dot23(&Gae,&WaMeF,&newLIA, 0, 0, -1.0, 1.0);
450
 
      dpd_buf4_close(&WaMeF);
451
 
      dpd_buf4_init(&Wamef, CC_HBAR, 0, 10, 5, 10, 7, 0, "Wamef");
452
 
      dpd_dot23(&Gae,&Wamef,&newLia, 0, 0, -1.0, 1.0);
453
 
      dpd_buf4_close(&Wamef);
454
 
      dpd_buf4_init(&WAmEf, CC_HBAR, 0, 10, 5, 10, 5, 0, "WAmEf");
455
 
      dpd_dot23(&GAE,&WAmEf,&newLia, 0, 0, -1.0, 1.0);
456
 
      dpd_buf4_close(&WAmEf);
457
 
    */
458
 
 
459
 
    dpd_file2_close(&Gae);
460
 
    dpd_file2_close(&GAE);
461
 
  }
462
 
  else if(params.ref == 2) {
463
 
 
464
 
    dpd_file2_init(&GAE, CC_LAMBDA, L_irr, 1, 1, "GAE");
465
 
    dpd_file2_init(&Gae, CC_LAMBDA, L_irr, 3, 3, "Gae");
466
 
 
467
 
    dpd_buf4_init(&W, CC_HBAR, 0, 21, 5, 21, 7, 0, "WAMEF");
468
 
    dpd_dot13(&GAE,&W,&newLIA, 0, 0, -1, 1);
469
 
    dpd_buf4_close(&W);
470
 
 
471
 
    dpd_buf4_init(&W, CC_HBAR, 0, 25, 29, 25, 29, 0, "WaMeF");
472
 
    dpd_dot13(&Gae,&W,&newLIA, 0, 0, -1, 1);
473
 
    dpd_buf4_close(&W);
474
 
 
475
 
    dpd_buf4_init(&W, CC_HBAR, 0, 31, 15, 31, 17, 0, "Wamef");
476
 
    dpd_dot13(&Gae,&W,&newLia, 0, 0, -1, 1);
477
 
    dpd_buf4_close(&W);
478
 
 
479
 
    dpd_buf4_init(&W, CC_HBAR, 0, 26, 28, 26, 28, 0, "WAmEf");
480
 
    dpd_dot13(&GAE,&W,&newLia, 0, 0, -1, 1);
481
 
    dpd_buf4_close(&W);
482
 
 
483
 
    dpd_file2_close(&Gae);
484
 
    dpd_file2_close(&GAE);
485
 
 
486
 
  }
487
 
 
488
 
  /* L1 RHS += -Gmn*Wmina */
489
 
  if(params.ref == 0) {
490
 
    dpd_file2_init(&GMI, CC_LAMBDA, L_irr, 0, 0, "GMI");
491
 
 
492
 
    dpd_buf4_init(&WmNiE, CC_HBAR, 0, 0, 11, 0, 11, 0, "2WMnIe - WnMIe (Mn,eI)");
493
 
    dpd_dot14(&GMI, &WmNiE, &newLIA, 0, 0, -1.0, 1.0);
494
 
    dpd_buf4_close(&WmNiE);
495
 
 
496
 
    dpd_file2_close(&GMI);
497
 
  }
498
 
  else if(params.ref == 1) {
499
 
 
500
 
    dpd_file2_init(&GMI, CC_LAMBDA, L_irr, 0, 0, "GMI");
501
 
    dpd_file2_init(&Gmi, CC_LAMBDA, L_irr, 0, 0, "Gmi");
502
 
 
503
 
    dpd_buf4_init(&WMNIE, CC_HBAR, 0, 0, 11, 2, 11, 0, "WMNIE (M>N,EI)");
504
 
    dpd_dot14(&GMI, &WMNIE, &newLIA, 0, 0, -1.0, 1.0); 
505
 
    dpd_buf4_close(&WMNIE);
506
 
 
507
 
    dpd_buf4_init(&WmNiE, CC_HBAR, 0, 0, 11, 0, 11, 0, "WmNiE (mN,Ei)");
508
 
    dpd_dot14(&Gmi, &WmNiE, &newLIA, 0, 0, -1.0, 1.0);
509
 
    dpd_buf4_close(&WmNiE);
510
 
 
511
 
    dpd_buf4_init(&Wmnie, CC_HBAR, 0, 0, 11, 2, 11, 0, "Wmnie (m>n,ei)");
512
 
    dpd_dot14(&Gmi, &Wmnie, &newLia, 0, 0, -1.0, 1.0);
513
 
    dpd_buf4_close(&Wmnie);
514
 
 
515
 
    dpd_buf4_init(&WMnIe, CC_HBAR, 0, 0, 11, 0, 11, 0, "WMnIe (Mn,eI)");
516
 
    dpd_dot14(&GMI, &WMnIe, &newLia, 0, 0, -1.0, 1.0);
517
 
    dpd_buf4_close(&WMnIe);
518
 
 
519
 
    dpd_file2_close(&Gmi);
520
 
    dpd_file2_close(&GMI);
521
 
 
522
 
  }
523
 
  else if(params.ref == 2) {
524
 
 
525
 
    dpd_file2_init(&GMI, CC_LAMBDA, L_irr, 0, 0, "GMI");
526
 
    dpd_file2_init(&Gmi, CC_LAMBDA, L_irr, 2, 2, "Gmi");
527
 
 
528
 
    dpd_buf4_init(&W, CC_HBAR, 0, 0, 21, 2, 21, 0, "WMNIE (M>N,EI)");
529
 
    dpd_dot14(&GMI, &W, &newLIA, 0, 0, -1, 1); 
530
 
    dpd_buf4_close(&W);
531
 
 
532
 
    dpd_buf4_init(&W, CC_HBAR, 0, 23, 26, 23, 26, 0, "WmNiE (mN,Ei)");
533
 
    dpd_dot14(&Gmi, &W, &newLIA, 0, 0, -1, 1);
534
 
    dpd_buf4_close(&W);
535
 
 
536
 
    dpd_buf4_init(&W, CC_HBAR, 0, 10, 31, 12, 31, 0, "Wmnie (m>n,ei)");
537
 
    dpd_dot14(&Gmi, &W, &newLia, 0, 0, -1, 1);
538
 
    dpd_buf4_close(&W);
539
 
 
540
 
    dpd_buf4_init(&W, CC_HBAR, 0, 22, 25, 22, 25, 0, "WMnIe (Mn,eI)");
541
 
    dpd_dot14(&GMI, &W, &newLia, 0, 0, -1, 1);
542
 
    dpd_buf4_close(&W);
543
 
 
544
 
    dpd_file2_close(&Gmi);
545
 
    dpd_file2_close(&GMI);
546
 
  }
547
 
 
548
 
  /* CC3 T3->L1 */
549
 
  if(!strcmp(params.wfn, "CC3")) {
550
 
    if(params.ref == 0) { 
551
 
 
552
 
      dpd_file2_init(&XLD, CC3_MISC, 0, 0, 1, "CC3 XLD");
553
 
      dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 0, "D 2<ij|ab> - <ij|ba>");
554
 
      dpd_dot24(&XLD, &D, &newLIA, 0, 0, 1, 1);
555
 
      dpd_buf4_close(&D);
556
 
      dpd_file2_close(&XLD);
557
 
 
558
 
      dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 0, 5, 2, 7, 0, "LIJAB");
559
 
      dpd_buf4_init(&Z, CC3_MISC, 0, 10, 0, 10, 0, 0, "CC3 ZIFLN");
560
 
      dpd_contract442(&Z, &L2, &newLIA, 0, 2, -0.5, 1);
561
 
      dpd_buf4_close(&Z);
562
 
      dpd_buf4_close(&L2);
563
 
 
564
 
      dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "LIjAb");
565
 
      dpd_buf4_init(&Z, CC3_MISC, 0, 10, 0, 10, 0, 0, "CC3 ZIfLn");
566
 
      dpd_contract442(&Z, &L2, &newLIA, 0, 2, -1, 1);
567
 
      dpd_buf4_close(&Z);
568
 
      dpd_buf4_close(&L2);
569
 
 
570
 
      dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 0, 5, 2, 7, 0, "LIJAB");
571
 
      dpd_buf4_init(&Z, CC3_MISC, 0, 11, 5, 11, 5, 0, "CC3 ZDFAN (AN,DF)");
572
 
      dpd_contract442(&L2, &Z, &newLIA, 0, 0, 0.5, 1);
573
 
      dpd_buf4_close(&Z);
574
 
      dpd_buf4_close(&L2);
575
 
 
576
 
      dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "LIjAb");
577
 
      dpd_buf4_init(&Z, CC3_MISC, 0, 11, 5, 11, 5, 0, "CC3 ZDfAn (An,Df)");
578
 
      dpd_contract442(&L2, &Z, &newLIA, 0, 0, 1.0, 1);
579
 
      dpd_buf4_close(&Z);
580
 
      dpd_buf4_close(&L2);
581
 
    }
582
 
  }
583
 
 
584
 
  dpd_file2_close(&newLIA);
585
 
  dpd_file2_close(&newLia);
586
 
 
587
 
  /* newLia * Dia */
588
 
  if(params.ref == 0) { /** RHF **/
589
 
 
590
 
    dpd_file2_init(&newLIA, CC_LAMBDA, L_irr, 0, 1, "New LIA");
591
 
    dpd_file2_copy(&newLIA, CC_LAMBDA, "New LIA Increment");
592
 
    dpd_file2_close(&newLIA);
593
 
 
594
 
    dpd_file2_init(&newLIA, CC_LAMBDA, L_irr, 0, 1, "New LIA Increment");
595
 
    if(params.local && local.filter_singles) local_filter_T1(&newLIA);
596
 
    else {
597
 
      dpd_file2_init(&dIA, CC_DENOM, L_irr, 0, 1, "dIA");
598
 
      dpd_file2_dirprd(&dIA, &newLIA);
599
 
      dpd_file2_close(&dIA);
600
 
    }
601
 
    dpd_file2_close(&newLIA);
602
 
 
603
 
    dpd_file2_init(&LIA, CC_LAMBDA, L_irr, 0, 1, "LIA");
604
 
    dpd_file2_copy(&LIA, CC_LAMBDA, "New LIA");
605
 
    dpd_file2_close(&LIA);
606
 
    dpd_file2_init(&newLIA, CC_LAMBDA, L_irr, 0, 1, "New LIA");
607
 
    dpd_file2_init(&LIA, CC_LAMBDA, L_irr, 0, 1, "New LIA Increment");
608
 
    dpd_file2_axpy(&LIA, &newLIA, 1, 0);
609
 
    dpd_file2_close(&LIA);
610
 
 
611
 
    dpd_file2_copy(&newLIA, CC_LAMBDA, "New Lia");  /* spin-adaptation for RHF */
612
 
    dpd_file2_close(&newLIA);
613
 
  }
614
 
  else if(params.ref == 1) { /** ROHF **/
615
 
 
616
 
    dpd_file2_init(&newLIA, CC_LAMBDA, L_irr, 0, 1, "New LIA");
617
 
    dpd_file2_copy(&newLIA, CC_LAMBDA, "New LIA Increment");
618
 
    dpd_file2_close(&newLIA);
619
 
 
620
 
    dpd_file2_init(&newLIA, CC_LAMBDA, L_irr, 0, 1, "New LIA Increment");
621
 
    dpd_file2_init(&dIA, CC_DENOM, L_irr, 0, 1, "dIA");
622
 
    dpd_file2_dirprd(&dIA, &newLIA);
623
 
    dpd_file2_close(&dIA);
624
 
    dpd_file2_close(&newLIA);
625
 
 
626
 
    dpd_file2_init(&LIA, CC_LAMBDA, L_irr, 0, 1, "LIA");
627
 
    dpd_file2_copy(&LIA, CC_LAMBDA, "New LIA");
628
 
    dpd_file2_close(&LIA);
629
 
    dpd_file2_init(&newLIA, CC_LAMBDA, L_irr, 0, 1, "New LIA");
630
 
    dpd_file2_init(&LIA, CC_LAMBDA, L_irr, 0, 1, "New LIA Increment");
631
 
    dpd_file2_axpy(&LIA, &newLIA, 1, 0);
632
 
    dpd_file2_close(&LIA);
633
 
    dpd_file2_close(&newLIA);
634
 
 
635
 
    dpd_file2_init(&newLia, CC_LAMBDA, L_irr, 0, 1, "New Lia");
636
 
    dpd_file2_copy(&newLia, CC_LAMBDA, "New Lia Increment");
637
 
    dpd_file2_close(&newLia);
638
 
 
639
 
    dpd_file2_init(&newLia, CC_LAMBDA, L_irr, 0, 1, "New Lia Increment");
640
 
    dpd_file2_init(&dia, CC_DENOM, L_irr, 0, 1, "dia");
641
 
    dpd_file2_dirprd(&dia, &newLia);
642
 
    dpd_file2_close(&dia);
643
 
    dpd_file2_close(&newLia);
644
 
 
645
 
    dpd_file2_init(&Lia, CC_LAMBDA, L_irr, 0, 1, "Lia");
646
 
    dpd_file2_copy(&Lia, CC_LAMBDA, "New Lia");
647
 
    dpd_file2_close(&Lia);
648
 
    dpd_file2_init(&newLia, CC_LAMBDA, L_irr, 0, 1, "New Lia");
649
 
    dpd_file2_init(&Lia, CC_LAMBDA, L_irr, 0, 1, "New Lia Increment");
650
 
    dpd_file2_axpy(&Lia, &newLia, 1, 0);
651
 
    dpd_file2_close(&Lia);
652
 
    dpd_file2_close(&newLia);
653
 
  }
654
 
  else if(params.ref == 2) {
655
 
 
656
 
    dpd_file2_init(&newLIA, CC_LAMBDA, L_irr, 0, 1, "New LIA");
657
 
    dpd_file2_init(&dIA, CC_DENOM, L_irr, 0, 1, "dIA");
658
 
    dpd_file2_dirprd(&dIA, &newLIA);
659
 
    dpd_file2_close(&dIA);
660
 
    dpd_file2_close(&newLIA);
661
 
 
662
 
    dpd_file2_init(&newLia, CC_LAMBDA, L_irr, 2, 3, "New Lia");
663
 
    dpd_file2_init(&dia, CC_DENOM, L_irr, 2, 3, "dia");
664
 
    dpd_file2_dirprd(&dia, &newLia);
665
 
    dpd_file2_close(&dia);
666
 
    dpd_file2_close(&newLia);
667
 
  }
668
 
 
669
 
#ifdef EOM_DEBUG
670
 
  check_sum("after L1 build",L_irr);
671
 
#endif
672
 
 
673
 
  return;
674
 
}
675
 
 
676