← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/ccdensity/x_Gciab.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/ccdensity/x_Gciab.c
1
 
#include <stdio.h>
2
 
#include <libdpd/dpd.h>
3
 
#define EXTERN
4
 
#include "globals.h"
5
 
 
6
 
void x_Gciab_rohf(void);
7
 
void x_Gciab_6(void);
8
 
void x_Gciab_7(void);
9
 
void x_Gciab_8_rohf(void);
10
 
extern void x_Gciab_uhf(void);
11
 
 
12
 
/* This function computes the non-R0 parts of the 2pdm density matrix
13
 
   Gciab = 0.5 *(rho_abci + rho_ciab) */
14
 
void x_Gciab(void) {
15
 
  if (params.ref == 0 || params.ref == 1)
16
 
    x_Gciab_rohf();
17
 
  else
18
 
    x_Gciab_uhf();
19
 
}
20
 
 
21
 
void x_Gciab_rohf(void) { 
22
 
  int h, nirreps, i, j, k, a, I, J, K, A, Isym, Jsym, Ksym, Asym, row, col;
23
 
  int II,JJ,IIsym,JJsym;
24
 
  int L_irr, R_irr, G_irr;
25
 
  double value;
26
 
  dpdfile2 L1, T1, R1, I1;
27
 
  dpdbuf4 G, V, T, L, Z, Z2, R, Tau;
28
 
 
29
 
  L_irr = params.L_irr; R_irr = params.R_irr; G_irr = params.G_irr;
30
 
  nirreps = moinfo.nirreps;
31
 
 
32
 
  /* term 1, rho_abci += Lmiab * Rmc */
33
 
  dpd_buf4_init(&Z, EOM_TMP, G_irr, 7, 11, 7, 11, 0, "L2R1_VVOV(pqsr)");
34
 
  dpd_buf4_sort(&Z, EOM_TMP0, rspq, 11, 7, "GCIAB");
35
 
  dpd_buf4_close(&Z);
36
 
  dpd_buf4_init(&Z, EOM_TMP, G_irr, 7, 11, 7, 11, 0, "L2R1_vvov(pqsr)");
37
 
  dpd_buf4_sort(&Z, EOM_TMP0, rspq, 11, 7, "Gciab");
38
 
  dpd_buf4_close(&Z);
39
 
  dpd_buf4_init(&Z, EOM_TMP, G_irr, 5, 11, 5, 11, 0, "L2R1_VvoV(pqsr)");
40
 
  dpd_buf4_sort(&Z, EOM_TMP0, rspq, 11, 5, "GCiAb");
41
 
  dpd_buf4_close(&Z);
42
 
  dpd_buf4_init(&Z, EOM_TMP, G_irr, 5, 11, 5, 11, 0, "L2R1_VvOv(pqsr)");
43
 
  dpd_buf4_sort(&Z, EOM_TMP0, rsqp, 11, 5, "GcIaB");
44
 
  dpd_buf4_close(&Z);
45
 
  dpd_buf4_init(&Z, EOM_TMP0, G_irr, 11, 7, 11, 7, 0, "GCIAB");
46
 
  dpd_buf4_scm(&Z, -1.0);
47
 
  dpd_buf4_close(&Z);
48
 
  dpd_buf4_init(&Z, EOM_TMP0, G_irr, 11, 7, 11, 7, 0, "Gciab");
49
 
  dpd_buf4_scm(&Z, -1.0);
50
 
  dpd_buf4_close(&Z);
51
 
 
52
 
  /* term 2, rho_ciab += Rmiab * Lmc */
53
 
  dpd_buf4_init(&Z, EOM_TMP, G_irr, 7, 11, 7, 11, 0, "R2L1_VVOV(pqsr)");
54
 
  dpd_buf4_sort_axpy(&Z, EOM_TMP0, rspq, 11, 7, "GCIAB",-1.0);
55
 
  dpd_buf4_close(&Z);
56
 
  dpd_buf4_init(&Z, EOM_TMP, G_irr, 7, 11, 7, 11, 0, "R2L1_vvov(pqsr)");
57
 
  dpd_buf4_sort_axpy(&Z, EOM_TMP0, rspq, 11, 7, "Gciab",-1.0);
58
 
  dpd_buf4_close(&Z);
59
 
  dpd_buf4_init(&Z, EOM_TMP, G_irr, 5, 11, 5, 11, 0, "R2L1_VvoV(pqsr)");
60
 
  dpd_buf4_sort_axpy(&Z, EOM_TMP0, rspq, 11, 5, "GCiAb", 1.0);
61
 
  dpd_buf4_close(&Z);
62
 
  dpd_buf4_init(&Z, EOM_TMP, G_irr, 5, 11, 5, 11, 0, "R2L1_VvOv(pqsr)");
63
 
  dpd_buf4_sort_axpy(&Z, EOM_TMP0, rsqp, 11, 5, "GcIaB", 1.0);
64
 
  dpd_buf4_close(&Z);
65
 
  /* same in two-steps, perhaps suggestive of future out-of-core approach
66
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 7, 11, 7, 0, "GCIAB");
67
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 7, 11, 7, 0, "GCIAB");
68
 
  dpd_buf4_axpy(&Z, &G, -1.0);
69
 
  dpd_buf4_close(&Z);
70
 
  dpd_buf4_close(&G);
71
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 7, 11, 7, 0, "Gciab");
72
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 7, 11, 7, 0, "Gciab");
73
 
  dpd_buf4_axpy(&Z, &G, -1.0);
74
 
  dpd_buf4_close(&Z);
75
 
  dpd_buf4_close(&G);
76
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 5, 0, "GCiAb");
77
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 5, 11, 5, 0, "GCiAb");
78
 
  dpd_buf4_axpy(&Z, &G, 1.0);
79
 
  dpd_buf4_close(&Z);
80
 
  dpd_buf4_close(&G);
81
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 5, 0, "GcIaB");
82
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 5, 11, 5, 0, "GcIaB");
83
 
  dpd_buf4_axpy(&Z, &G, 1.0);
84
 
  dpd_buf4_close(&Z);
85
 
  dpd_buf4_close(&G);
86
 
  psio_close(EOM_TMP1,0);
87
 
  psio_open(EOM_TMP1, PSIO_OPEN_NEW);
88
 
  */
89
 
 
90
 
  /* term 3, rho_CIAB -= 0.5 LMNCE RMNAB tIE */
91
 
  dpd_buf4_init(&Z, EOM_TMP1, L_irr, 11, 2, 11, 2, 0, "Z(CI,MN)");
92
 
  dpd_buf4_init(&L, CC_GL, L_irr, 2, 5, 2, 7, 0, "LIJAB");
93
 
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
94
 
  dpd_contract424(&L, &T1, &Z, 3, 1, 1, 1.0, 0.0);
95
 
  dpd_file2_close(&T1);
96
 
  dpd_buf4_close(&L);
97
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 7, 11, 7, 0, "GCIAB");
98
 
  dpd_buf4_init(&R, CC_GR, R_irr, 2, 7, 2, 7, 0, "RIJAB");
99
 
  dpd_contract444(&Z, &R, &G, 0, 1, -1.0, 1.0);
100
 
  dpd_buf4_close(&R);
101
 
  dpd_buf4_close(&G);
102
 
  dpd_buf4_close(&Z);
103
 
  /* term 3, rho_ciab -= 0.5 Lmnce Rmnab tie */
104
 
  dpd_buf4_init(&Z, EOM_TMP1, L_irr, 11, 2, 11, 2, 0, "Z(ci,mn)");
105
 
  dpd_buf4_init(&L, CC_GL, L_irr, 2, 5, 2, 7, 0, "Lijab");
106
 
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
107
 
  dpd_contract424(&L, &T1, &Z, 3, 1, 1, 1.0, 0.0);
108
 
  dpd_file2_close(&T1);
109
 
  dpd_buf4_close(&L);
110
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 7, 11, 7, 0, "Gciab");
111
 
  dpd_buf4_init(&R, CC_GR, R_irr, 2, 7, 2, 7, 0, "Rijab");
112
 
  dpd_contract444(&Z, &R, &G, 0, 1, -1.0, 1.0);
113
 
  dpd_buf4_close(&R);
114
 
  dpd_buf4_close(&G);
115
 
  dpd_buf4_close(&Z);
116
 
  /* term 3, rho_CiAb -= 0.5 LMnCe RMnAb tie */
117
 
  dpd_buf4_init(&Z, EOM_TMP1, L_irr, 11, 0, 11, 0, 0, "Z(Ci,Mn)");
118
 
  dpd_buf4_init(&L, CC_GL, L_irr, 0, 5, 0, 5, 0, "LIjAb");
119
 
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
120
 
  dpd_contract424(&L, &T1, &Z, 3, 1, 1, 1.0, 0.0);
121
 
  dpd_file2_close(&T1);
122
 
  dpd_buf4_close(&L);
123
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 5, 0, "GCiAb");
124
 
  dpd_buf4_init(&R, CC_GR, R_irr, 0, 5, 0, 5, 0, "RIjAb");
125
 
  dpd_contract444(&Z, &R, &G, 0, 1, -1.0, 1.0);
126
 
  dpd_buf4_close(&R);
127
 
  dpd_buf4_close(&G);
128
 
  dpd_buf4_close(&Z);
129
 
  /* term 3, rho_cIaB -= 0.5 LmNcE RmNaB tIE */
130
 
  dpd_buf4_init(&Z, EOM_TMP1, L_irr, 11, 0, 11, 0, 0, "Z(cI,mN)");
131
 
  dpd_buf4_init(&L, CC_GL, L_irr, 0, 5, 0, 5, 0, "LiJaB");
132
 
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
133
 
  dpd_contract424(&L, &T1, &Z, 3, 1, 1, 1.0, 0.0);
134
 
  dpd_file2_close(&T1);
135
 
  dpd_buf4_close(&L);
136
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 5, 0, "GcIaB");
137
 
  dpd_buf4_init(&R, CC_GR, R_irr, 0, 5, 0, 5, 0, "RiJaB");
138
 
  dpd_contract444(&Z, &R, &G, 0, 1, -1.0, 1.0);
139
 
  dpd_buf4_close(&R);
140
 
  dpd_buf4_close(&G);
141
 
  dpd_buf4_close(&Z);
142
 
 
143
 
  psio_close(EOM_TMP1,0);
144
 
  psio_open(EOM_TMP1, PSIO_OPEN_NEW);
145
 
 
146
 
  /* term 4, rho_CIAB -= 0.5 LMNCE TauMNAB RIE */
147
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 2, 11, 2, 0, "Z(CI,MN)");
148
 
  dpd_buf4_init(&L, CC_GL, L_irr, 2, 5, 2, 7, 0, "LIJAB");
149
 
  dpd_file2_init(&R1, CC_GR, R_irr, 0, 1, "RIA");
150
 
  dpd_contract424(&L, &R1, &Z, 3, 1, 1, 1.0, 0.0);
151
 
  dpd_file2_close(&R1);
152
 
  dpd_buf4_close(&L);
153
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 7, 11, 7, 0, "GCIAB");
154
 
  dpd_buf4_init(&T, CC_TAMPS, 0, 2, 7, 2, 7, 0, "tauIJAB");
155
 
  dpd_contract444(&Z, &T, &G, 0, 1, -1.0, 1.0);
156
 
  dpd_buf4_close(&T);
157
 
  dpd_buf4_close(&G);
158
 
  dpd_buf4_close(&Z);
159
 
  /* term 4, rho_ciab -= 0.5 Lmnce Taumnab Rie */
160
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 2, 11, 2, 0, "Z(ci,mn)");
161
 
  dpd_buf4_init(&L, CC_GL, L_irr, 2, 5, 2, 7, 0, "Lijab");
162
 
  dpd_file2_init(&R1, CC_GR, R_irr, 0, 1, "Ria");
163
 
  dpd_contract424(&L, &R1, &Z, 3, 1, 1, 1.0, 0.0);
164
 
  dpd_file2_close(&R1);
165
 
  dpd_buf4_close(&L);
166
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 7, 11, 7, 0, "Gciab");
167
 
  dpd_buf4_init(&T, CC_TAMPS, 0, 2, 7, 2, 7, 0, "tauijab");
168
 
  dpd_contract444(&Z, &T, &G, 0, 1, -1.0, 1.0);
169
 
  dpd_buf4_close(&T);
170
 
  dpd_buf4_close(&G);
171
 
  dpd_buf4_close(&Z);
172
 
  /* term 4, rho_CiAb -= 0.5 LMnCe TauMnAb Rie */
173
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 0, 11, 0, 0, "Z(Ci,Mn)");
174
 
  dpd_buf4_init(&L, CC_GL, L_irr, 0, 5, 0, 5, 0, "LIjAb");
175
 
  dpd_file2_init(&R1, CC_GR, R_irr, 0, 1, "Ria");
176
 
  dpd_contract424(&L, &R1, &Z, 3, 1, 1, 1.0, 0.0);
177
 
  dpd_file2_close(&R1);
178
 
  dpd_buf4_close(&L);
179
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 5, 0, "GCiAb");
180
 
  dpd_buf4_init(&T, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tauIjAb");
181
 
  dpd_contract444(&Z, &T, &G, 0, 1, -1.0, 1.0);
182
 
  dpd_buf4_close(&T);
183
 
  dpd_buf4_close(&G);
184
 
  dpd_buf4_close(&Z);
185
 
  /* term 4, rho_cIaB -= 0.5 LmNcE TaumNaB RIE */
186
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 0, 11, 0, 0, "Z(cI,mN)");
187
 
  dpd_buf4_init(&L, CC_GL, L_irr, 0, 5, 0, 5, 0, "LiJaB");
188
 
  dpd_file2_init(&R1, CC_GR, R_irr, 0, 1, "RIA");
189
 
  dpd_contract424(&L, &R1, &Z, 3, 1, 1, 1.0, 0.0);
190
 
  dpd_file2_close(&R1);
191
 
  dpd_buf4_close(&L);
192
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 5, 0, "GcIaB");
193
 
  dpd_buf4_init(&T, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tauiJaB");
194
 
  dpd_contract444(&Z, &T, &G, 0, 1, -1.0, 1.0);
195
 
  dpd_buf4_close(&T);
196
 
  dpd_buf4_close(&G);
197
 
  dpd_buf4_close(&Z);
198
 
 
199
 
  /* term 5, rho_ciab += - (Lmnec Rme) (Tniab + P(ij) Tna Tib) */
200
 
  if (!params.connect_xi) {
201
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 7, 11, 7, 0, "GCIAB");
202
 
  dpd_buf4_init(&Tau, CC_TAMPS, 0, 0, 7, 2, 7, 0, "tauIJAB");
203
 
  dpd_file2_init(&I1, EOM_TMP, G_irr, 0, 1, "L2R1_OV");
204
 
  dpd_contract244(&I1, &Tau, &G, 0, 0, 0, 1.0, 1.0);
205
 
  dpd_file2_close(&I1);
206
 
  dpd_buf4_close(&Tau);
207
 
  dpd_buf4_close(&G);
208
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 7, 11, 7, 0, "Gciab");
209
 
  dpd_buf4_init(&Tau, CC_TAMPS, 0, 0, 7, 2, 7, 0, "tauijab");
210
 
  dpd_file2_init(&I1, EOM_TMP, G_irr, 0, 1, "L2R1_ov");
211
 
  dpd_contract244(&I1, &Tau, &G, 0, 0, 0, 1.0, 1.0);
212
 
  dpd_file2_close(&I1);
213
 
  dpd_buf4_close(&Tau);
214
 
  dpd_buf4_close(&G);
215
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 5, 0, "GCiAb");
216
 
  dpd_buf4_init(&Tau, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tauIjAb");
217
 
  dpd_file2_init(&I1, EOM_TMP, G_irr, 0, 1, "L2R1_OV");
218
 
  dpd_contract244(&I1, &Tau, &G, 0, 0, 0, 1.0, 1.0);
219
 
  dpd_file2_close(&I1);
220
 
  dpd_buf4_close(&Tau);
221
 
  dpd_buf4_close(&G);
222
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 5, 0, "GcIaB");
223
 
  dpd_buf4_init(&Tau, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tauiJaB");
224
 
  dpd_file2_init(&I1, EOM_TMP, G_irr, 0, 1, "L2R1_ov");
225
 
  dpd_contract244(&I1, &Tau, &G, 0, 0, 0, 1.0, 1.0);
226
 
  dpd_file2_close(&I1);
227
 
  dpd_buf4_close(&Tau);
228
 
  dpd_buf4_close(&G);
229
 
  }
230
 
 
231
 
  /* +P(ab) LR_VV(c,a) t(i,b) */
232
 
  x_Gciab_6();
233
 
 
234
 
  /* +P(ab) LR_TT(c,a) R(i,b) */
235
 
  x_Gciab_7();
236
 
 
237
 
  /* -P(ab) Lmnce Rinae Tmb, term 8 */
238
 
  /* -P(ab) Lmnce Tinae Rmb, term 9 */
239
 
  x_Gciab_8_rohf();
240
 
 
241
 
  psio_close(EOM_TMP1,0);
242
 
  psio_open(EOM_TMP1, PSIO_OPEN_NEW);
243
 
 
244
 
  /* term 10, -P(AB) LNMCE TIE TNA RMB */
245
 
  dpd_buf4_init(&Z, EOM_TMP1, L_irr, 2, 11, 2, 11, 0, "Z(NM,CI)");
246
 
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
247
 
  dpd_buf4_init(&L, CC_GL, L_irr, 2, 5, 2, 7, 0, "LIJAB");
248
 
  dpd_contract424(&L, &T1, &Z, 3, 1, 0, 1.0, 0.0);
249
 
  dpd_buf4_close(&L);
250
 
  dpd_buf4_close(&Z);
251
 
  dpd_buf4_init(&Z, EOM_TMP1, L_irr, 0, 11, 2, 11, 0, "Z(NM,CI)");
252
 
  dpd_buf4_init(&Z2, EOM_TMP1, L_irr, 11, 11, 11, 11, 0, "Z(CI,AM)");
253
 
  dpd_contract244(&T1, &Z, &Z2, 0, 0, 1, 1.0, 0.0);
254
 
  dpd_file2_close(&T1);
255
 
  dpd_buf4_close(&Z);
256
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 5, 11, 5, 0, "Z(CI,AB)");
257
 
  dpd_file2_init(&R1, CC_GR, R_irr, 0, 1, "RIA");
258
 
  dpd_contract424(&Z2, &R1, &Z, 3, 0, 0, 1.0, 0.0);
259
 
  dpd_file2_close(&R1);
260
 
  dpd_buf4_close(&Z2);
261
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 7, 0, "GCIAB");
262
 
  dpd_buf4_axpy(&Z, &G, -1.0);
263
 
  dpd_buf4_sort(&Z, EOM_TMP1, pqsr, 11, 5, "Z(CI,BA)");
264
 
  dpd_buf4_close(&Z);
265
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 5, 11, 5, 0, "Z(CI,BA)");
266
 
  dpd_buf4_axpy(&Z, &G, 1.0);
267
 
  dpd_buf4_close(&Z);
268
 
  dpd_buf4_close(&G);
269
 
 
270
 
  /* term 10, -P(ab) lmnce tie tna rmb */
271
 
  dpd_buf4_init(&Z, EOM_TMP1, L_irr, 2, 11, 2, 11, 0, "Z(nm,ci)");
272
 
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
273
 
  dpd_buf4_init(&L, CC_GL, L_irr, 2, 5, 2, 7, 0, "Lijab");
274
 
  dpd_contract424(&L, &T1, &Z, 3, 1, 0, 1.0, 0.0);
275
 
  dpd_buf4_close(&L);
276
 
  dpd_buf4_close(&Z);
277
 
  dpd_buf4_init(&Z, EOM_TMP1, L_irr, 0, 11, 2, 11, 0, "Z(nm,ci)");
278
 
  dpd_buf4_init(&Z2, EOM_TMP1, L_irr, 11, 11, 11, 11, 0, "Z(ci,am)");
279
 
  dpd_contract244(&T1, &Z, &Z2, 0, 0, 1, 1.0, 0.0);
280
 
  dpd_file2_close(&T1);
281
 
  dpd_buf4_close(&Z);
282
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 5, 11, 5, 0, "Z(ci,ab)");
283
 
  dpd_file2_init(&R1, CC_GR, R_irr, 0, 1, "Ria");
284
 
  dpd_contract424(&Z2, &R1, &Z, 3, 0, 0, 1.0, 0.0);
285
 
  dpd_file2_close(&R1);
286
 
  dpd_buf4_close(&Z2);
287
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 7, 0, "Gciab");
288
 
  dpd_buf4_axpy(&Z, &G, -1.0);
289
 
  dpd_buf4_sort(&Z, EOM_TMP1, pqsr, 11, 5, "Z(ci,ba)");
290
 
  dpd_buf4_close(&Z);
291
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 5, 11, 5, 0, "Z(ci,ba)");
292
 
  dpd_buf4_axpy(&Z, &G, 1.0);
293
 
  dpd_buf4_close(&Z);
294
 
  dpd_buf4_close(&G);
295
 
 
296
 
  /* term 10, GCiAb -= P(AB) LNmCe Tie TNA Rmb */
297
 
  dpd_buf4_init(&Z, EOM_TMP1, L_irr, 0, 11, 0, 11, 0, "Z(Nm,Ci)");
298
 
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
299
 
  dpd_buf4_init(&L, CC_GL, L_irr, 0, 5, 0, 5, 0, "LIjAb");
300
 
  dpd_contract424(&L, &T1, &Z, 3, 1, 0, 1.0, 0.0);
301
 
  dpd_buf4_close(&L);
302
 
  dpd_file2_close(&T1);
303
 
  dpd_buf4_init(&Z2, EOM_TMP1, L_irr, 11, 11, 11, 11, 0, "Z(Ci,Am)");
304
 
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
305
 
  dpd_contract244(&T1, &Z, &Z2, 0, 0, 1, 1.0, 0.0);
306
 
  dpd_file2_close(&T1);
307
 
  dpd_buf4_close(&Z);
308
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 5, 0, "GCiAb");
309
 
  dpd_file2_init(&R1, CC_GR, R_irr, 0, 1, "Ria");
310
 
  dpd_contract424(&Z2, &R1, &G, 3, 0, 0, -1.0, 1.0);
311
 
  dpd_file2_close(&R1);
312
 
  dpd_buf4_close(&Z2);
313
 
  dpd_buf4_close(&G);
314
 
 
315
 
  dpd_buf4_init(&Z, EOM_TMP1, L_irr, 0, 11, 0, 11, 0, "Z(nM,Ci)");
316
 
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
317
 
  dpd_buf4_init(&L, CC_GL, L_irr, 0, 5, 0, 5, 0, "LiJAb");
318
 
  dpd_contract424(&L, &T1, &Z, 3, 1, 0, 1.0, 0.0);
319
 
  dpd_buf4_close(&L);
320
 
  dpd_file2_close(&T1);
321
 
  dpd_buf4_init(&Z2, EOM_TMP1, L_irr, 11, 11, 11, 11, 0, "Z(Ci,bM)");
322
 
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
323
 
  dpd_contract244(&T1, &Z, &Z2, 0, 0, 1, 1.0, 0.0);
324
 
  dpd_file2_close(&T1);
325
 
  dpd_buf4_close(&Z);
326
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 5, 11, 5, 0, "Z(Ci,bA)");
327
 
  dpd_file2_init(&R1, CC_GR, R_irr, 0, 1, "RIA");
328
 
  dpd_contract424(&Z2, &R1, &Z, 3, 0, 0, 1.0, 0.0);
329
 
  dpd_file2_close(&R1);
330
 
  dpd_buf4_close(&Z2);
331
 
  dpd_buf4_sort(&Z, EOM_TMP1, pqsr, 11, 5, "Z(Ci,Ab)");
332
 
  dpd_buf4_close(&Z);
333
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 5, 11, 5, 0, "Z(Ci,Ab)");
334
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 5, 0, "GCiAb");
335
 
  dpd_buf4_axpy(&Z, &G, -1.0);
336
 
  dpd_buf4_close(&Z);
337
 
  dpd_buf4_close(&G);
338
 
 
339
 
  /* term 10, GcIaB - LnMcE TIE Tna RMB + LNmcE TIE TNB Rma */
340
 
  dpd_buf4_init(&Z, EOM_TMP1, L_irr, 0, 11, 0, 11, 0, "Z(nM,cI)");
341
 
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
342
 
  dpd_buf4_init(&L, CC_GL, L_irr, 0, 5, 0, 5, 0, "LiJaB");
343
 
  dpd_contract424(&L, &T1, &Z, 3, 1, 0, 1.0, 0.0);
344
 
  dpd_buf4_close(&L);
345
 
  dpd_file2_close(&T1);
346
 
  dpd_buf4_init(&Z2, EOM_TMP1, L_irr, 11, 11, 11, 11, 0, "Z(cI,aM)");
347
 
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tia");
348
 
  dpd_contract244(&T1, &Z, &Z2, 0, 0, 1, 1.0, 0.0);
349
 
  dpd_file2_close(&T1);
350
 
  dpd_buf4_close(&Z);
351
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 5, 0, "GcIaB");
352
 
  dpd_file2_init(&R1, CC_GR, R_irr, 0, 1, "RIA");
353
 
  dpd_contract424(&Z2, &R1, &G, 3, 0, 0, -1.0, 1.0);
354
 
  dpd_file2_close(&R1);
355
 
  dpd_buf4_close(&Z2);
356
 
  dpd_buf4_close(&G);
357
 
 
358
 
  dpd_buf4_init(&Z, EOM_TMP1, L_irr, 0, 11, 0, 11, 0, "Z(Nm,cI)");
359
 
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
360
 
  dpd_buf4_init(&L, CC_GL, L_irr, 0, 5, 0, 5, 0, "LIjaB");
361
 
  dpd_contract424(&L, &T1, &Z, 3, 1, 0, 1.0, 0.0);
362
 
  dpd_buf4_close(&L);
363
 
  dpd_file2_close(&T1);
364
 
  dpd_buf4_init(&Z2, EOM_TMP1, L_irr, 11, 11, 11, 11, 0, "Z(cI,Bm)");
365
 
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
366
 
  dpd_contract244(&T1, &Z, &Z2, 0, 0, 1, 1.0, 0.0);
367
 
  dpd_file2_close(&T1);
368
 
  dpd_buf4_close(&Z);
369
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 5, 11, 5, 0, "Z(cI,Ba)");
370
 
  dpd_file2_init(&R1, CC_GR, R_irr, 0, 1, "Ria");
371
 
  dpd_contract424(&Z2, &R1, &Z, 3, 0, 0, 1.0, 0.0);
372
 
  dpd_file2_close(&R1);
373
 
  dpd_buf4_close(&Z2);
374
 
  dpd_buf4_sort(&Z, EOM_TMP1, pqsr, 11, 5, "Z(cI,aB)");
375
 
  dpd_buf4_close(&Z);
376
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 5, 11, 5, 0, "Z(cI,aB)");
377
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 5, 0, "GcIaB");
378
 
  dpd_buf4_axpy(&Z, &G, -1.0);
379
 
  dpd_buf4_close(&Z);
380
 
  dpd_buf4_close(&G);
381
 
 
382
 
 
383
 
  /* add 1/2 to ground-state parts of density */
384
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 7, 11, 7, 0, "GCIAB");
385
 
  dpd_buf4_init(&V, CC_GAMMA, G_irr, 11, 7, 11, 7, 0, "GCIAB");
386
 
  dpd_buf4_axpy(&G, &V, 0.5);
387
 
  dpd_buf4_close(&V);
388
 
  dpd_buf4_close(&G);
389
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 7, 11, 7, 0, "Gciab");
390
 
  dpd_buf4_init(&V, CC_GAMMA, G_irr, 11, 7, 11, 7, 0, "Gciab");
391
 
  dpd_buf4_axpy(&G, &V, 0.5);
392
 
  dpd_buf4_close(&V);
393
 
  dpd_buf4_close(&G);
394
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 5, 0, "GCiAb");
395
 
  dpd_buf4_init(&V, CC_GAMMA, G_irr, 11, 5, 11, 5, 0, "GCiAb");
396
 
  dpd_buf4_axpy(&G, &V, 0.5);
397
 
  dpd_buf4_close(&V);
398
 
  dpd_buf4_close(&G);
399
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 5, 0, "GcIaB");
400
 
  dpd_buf4_init(&V, CC_GAMMA, G_irr, 11, 5, 11, 5, 0, "GcIaB");
401
 
  dpd_buf4_axpy(&G, &V, 0.5);
402
 
  dpd_buf4_close(&V);
403
 
  dpd_buf4_close(&G);
404
 
 
405
 
  /* clear out temporary files */
406
 
  psio_close(EOM_TMP0, 0);
407
 
  psio_open(EOM_TMP0, PSIO_OPEN_NEW);
408
 
 
409
 
  /*
410
 
  dpd_buf4_init(&V, CC_GAMMA, G_irr, 11, 7, 11, 7, 0, "GCIAB");
411
 
  value = dpd_buf4_dot_self(&V);
412
 
  dpd_buf4_close(&V);
413
 
  dpd_buf4_init(&V, CC_GAMMA, G_irr, 11, 7, 11, 7, 0, "Gciab");
414
 
  value += dpd_buf4_dot_self(&V);
415
 
  dpd_buf4_close(&V);
416
 
  dpd_buf4_init(&V, CC_GAMMA, G_irr, 11, 5, 11, 5, 0, "GCiAb");
417
 
  value += dpd_buf4_dot_self(&V);
418
 
  dpd_buf4_close(&V);
419
 
  dpd_buf4_init(&V, CC_GAMMA, G_irr, 11, 5, 11, 5, 0, "GcIaB");
420
 
  value += dpd_buf4_dot_self(&V);
421
 
  dpd_buf4_close(&V);
422
 
  fprintf(outfile,"<Gciab|Gciab> = %15.10lf\n",value);
423
 
  */
424
 
 
425
 
}
426
 
 
427
 
 
428
 
 
429
 
/* This function computes term 6,
430
 
   rho_ciab += P(ab) LR1_VV(c,a) T(i,b)
431
 
*/
432
 
 
433
 
void x_Gciab_6(void) { 
434
 
  int h, nirreps, c, i, a, b, C, I, A, B, Csym, Isym, Asym, Bsym, row, col;
435
 
  int AA, BB, AAsym, BBsym;
436
 
  int L_irr, R_irr, G_irr;
437
 
  dpdfile2 LR1A, LR1B, T1A, T1B;
438
 
  dpdbuf4 G;
439
 
 
440
 
  L_irr = params.L_irr; R_irr = params.R_irr; G_irr = params.G_irr;
441
 
  nirreps = moinfo.nirreps;
442
 
 
443
 
  /* open one-electron files for the nasty terms */
444
 
  dpd_file2_init(&LR1A, EOM_TMP, G_irr, 1, 1, "LR_VV");
445
 
  dpd_file2_init(&T1A, CC_OEI, 0, 0, 1, "tIA");
446
 
  if (params.ref == 0 || params.ref == 1) {
447
 
    dpd_file2_init(&LR1B, EOM_TMP, G_irr, 1, 1, "LR_vv");
448
 
    dpd_file2_init(&T1B, CC_OEI, 0, 0, 1, "tia");
449
 
  }
450
 
  else {
451
 
    dpd_file2_init(&LR1B, EOM_TMP, G_irr, 3, 3, "LR_vv");
452
 
    dpd_file2_init(&T1B, CC_OEI, 0, 2, 3, "tia");
453
 
  }
454
 
  dpd_file2_mat_init(&T1A);   dpd_file2_mat_init(&T1B);
455
 
  dpd_file2_mat_init(&LR1A);   dpd_file2_mat_init(&LR1B);
456
 
  dpd_file2_mat_rd(&T1A);     dpd_file2_mat_rd(&T1B);
457
 
  dpd_file2_mat_rd(&LR1A);     dpd_file2_mat_rd(&LR1B);
458
 
 
459
 
  /* rho_CIAB += LR1_VV(C,A) T(I,B) - LR1_VV(C,B) T(I,A) */
460
 
  if (params.ref == 0 || params.ref == 1)
461
 
    dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 7, 0, "GCIAB");
462
 
  else
463
 
    dpd_buf4_init(&G, EOM_TMP0, G_irr, 21, 5, 21, 7, 0, "GCIAB");
464
 
  for(h=0; h < nirreps; h++) {
465
 
    dpd_buf4_mat_irrep_init(&G, h);
466
 
    dpd_buf4_mat_irrep_rd(&G, h);
467
 
    for(row=0; row < G.params->rowtot[h]; row++) {
468
 
      c = G.params->roworb[h][row][0];
469
 
      i = G.params->roworb[h][row][1];
470
 
      C = LR1A.params->rowidx[c]; Csym = LR1A.params->psym[c];
471
 
      I = T1A.params->rowidx[i]; Isym = T1A.params->psym[i];
472
 
      for(col=0; col < G.params->coltot[h]; col++) {
473
 
        a = G.params->colorb[h^G_irr][col][0];
474
 
        b = G.params->colorb[h^G_irr][col][1];
475
 
        A = LR1A.params->colidx[a]; Asym = LR1A.params->qsym[a];
476
 
        B = T1A.params->colidx[b]; Bsym = T1A.params->qsym[b];
477
 
        AA = T1A.params->colidx[a]; AAsym = T1A.params->qsym[a];
478
 
        BB = LR1A.params->colidx[b]; BBsym = LR1A.params->qsym[b];
479
 
        if( ((Csym^Asym)==G_irr) && (Isym==Bsym))
480
 
          G.matrix[h][row][col] +=
481
 
            LR1A.matrix[Csym][C][A] * T1A.matrix[Isym][I][B];
482
 
        if( ((Csym^BBsym)==G_irr) && (Isym==AAsym))
483
 
          G.matrix[h][row][col] -=
484
 
            LR1A.matrix[Csym][C][BB] * T1A.matrix[Isym][I][AA];
485
 
      }
486
 
    }
487
 
    dpd_buf4_mat_irrep_wrt(&G, h);
488
 
    dpd_buf4_mat_irrep_close(&G, h);
489
 
  }
490
 
  dpd_buf4_close(&G);
491
 
  /* rho_ciab += LR1_vv(c,a) T(i,b) - LR1_vv(c,b) T(i,a) */
492
 
  if (params.ref == 0 || params.ref == 1)
493
 
    dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 7, 0, "Gciab");
494
 
  else
495
 
    dpd_buf4_init(&G, EOM_TMP0, G_irr, 31, 15, 31, 17, 0, "Gciab");
496
 
  for(h=0; h < nirreps; h++) {
497
 
    dpd_buf4_mat_irrep_init(&G, h);
498
 
    dpd_buf4_mat_irrep_rd(&G, h);
499
 
    for(row=0; row < G.params->rowtot[h]; row++) {
500
 
      c = G.params->roworb[h][row][0];
501
 
      i = G.params->roworb[h][row][1];
502
 
      C = LR1B.params->rowidx[c]; Csym = LR1B.params->psym[c];
503
 
      I = T1B.params->rowidx[i]; Isym = T1B.params->psym[i];
504
 
      for(col=0; col < G.params->coltot[h^G_irr]; col++) {
505
 
        a = G.params->colorb[h^G_irr][col][0];
506
 
        b = G.params->colorb[h^G_irr][col][1];
507
 
        A = LR1B.params->colidx[a]; Asym = LR1B.params->qsym[a];
508
 
        B = T1B.params->colidx[b]; Bsym = T1B.params->qsym[b];
509
 
        AA = T1B.params->colidx[a]; AAsym = T1B.params->qsym[a];
510
 
        BB = LR1B.params->colidx[b]; BBsym = LR1B.params->qsym[b];
511
 
        if( ((Csym^Asym)==G_irr) && (Isym==Bsym))
512
 
          G.matrix[h][row][col] +=
513
 
            LR1B.matrix[Csym][C][A] * T1B.matrix[Isym][I][B];
514
 
        if( ((Csym^BBsym)==G_irr) && (Isym==AAsym))
515
 
          G.matrix[h][row][col] -=
516
 
            LR1B.matrix[Csym][C][BB] * T1B.matrix[Isym][I][AA];
517
 
      }
518
 
    }
519
 
    dpd_buf4_mat_irrep_wrt(&G, h);
520
 
    dpd_buf4_mat_irrep_close(&G, h);
521
 
  }
522
 
  dpd_buf4_close(&G);
523
 
  /* rho_CiAb += LR1_VV(C,A) T(i,b) */
524
 
  if (params.ref == 0 || params.ref == 1)
525
 
    dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 5, 0, "GCiAb");
526
 
  else
527
 
    dpd_buf4_init(&G, EOM_TMP0, G_irr, 26, 28, 26, 28, 0, "GCiAb");
528
 
  for(h=0; h < nirreps; h++) {
529
 
    dpd_buf4_mat_irrep_init(&G, h);
530
 
    dpd_buf4_mat_irrep_rd(&G, h);
531
 
    for(row=0; row < G.params->rowtot[h]; row++) {
532
 
      c = G.params->roworb[h][row][0];
533
 
      i = G.params->roworb[h][row][1];
534
 
      C = LR1A.params->rowidx[c]; Csym = LR1A.params->psym[c];
535
 
      I = T1B.params->rowidx[i]; Isym = T1B.params->psym[i];
536
 
      for(col=0; col < G.params->coltot[h^G_irr]; col++) {
537
 
        a = G.params->colorb[h^G_irr][col][0];
538
 
        b = G.params->colorb[h^G_irr][col][1];
539
 
        A = LR1A.params->colidx[a]; Asym = LR1A.params->qsym[a];
540
 
        B = T1B.params->colidx[b]; Bsym = T1B.params->qsym[b];
541
 
        if( ((Csym^Asym)==G_irr) && (Isym==Bsym))
542
 
          G.matrix[h][row][col] +=
543
 
            LR1A.matrix[Csym][C][A] * T1B.matrix[Isym][I][B];
544
 
      }
545
 
    }
546
 
    dpd_buf4_mat_irrep_wrt(&G, h);
547
 
    dpd_buf4_mat_irrep_close(&G, h);
548
 
  }
549
 
  dpd_buf4_close(&G);
550
 
  /* rho_cIaB += LR1_vv(c,a) T(I,B) */
551
 
  if (params.ref == 0 || params.ref == 1)
552
 
    dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 5, 0, "GcIaB");
553
 
  else
554
 
    dpd_buf4_init(&G, EOM_TMP0, G_irr, 25, 29, 25, 29, 0, "GcIaB");
555
 
  for(h=0; h < nirreps; h++) {
556
 
    dpd_buf4_mat_irrep_init(&G, h);
557
 
    dpd_buf4_mat_irrep_rd(&G, h);
558
 
    for(row=0; row < G.params->rowtot[h]; row++) {
559
 
      c = G.params->roworb[h][row][0];
560
 
      i = G.params->roworb[h][row][1];
561
 
      C = LR1B.params->rowidx[c]; Csym = LR1B.params->psym[c];
562
 
      I = T1A.params->rowidx[i]; Isym = T1A.params->psym[i];
563
 
      for(col=0; col < G.params->coltot[h^G_irr]; col++) {
564
 
        a = G.params->colorb[h^G_irr][col][0];
565
 
        b = G.params->colorb[h^G_irr][col][1];
566
 
        A = LR1B.params->colidx[a]; Asym = LR1B.params->qsym[a];
567
 
        B = T1A.params->colidx[b]; Bsym = T1A.params->qsym[b];
568
 
        if( ((Csym^Asym)==G_irr) && (Isym==Bsym))
569
 
          G.matrix[h][row][col] +=
570
 
            LR1B.matrix[Csym][C][A] * T1A.matrix[Isym][I][B];
571
 
      }
572
 
    }
573
 
    dpd_buf4_mat_irrep_wrt(&G, h);
574
 
    dpd_buf4_mat_irrep_close(&G, h);
575
 
  }
576
 
  dpd_buf4_close(&G);
577
 
 
578
 
  dpd_file2_mat_close(&LR1A);
579
 
  dpd_file2_mat_close(&LR1B);
580
 
  dpd_file2_close(&LR1A);
581
 
  dpd_file2_close(&LR1B);
582
 
 
583
 
  dpd_file2_mat_close(&T1A);
584
 
  dpd_file2_mat_close(&T1B);
585
 
  dpd_file2_close(&T1A);
586
 
  dpd_file2_close(&T1B);
587
 
 
588
 
  return;
589
 
}
590
 
 
591
 
 
592
 
 
593
 
/* This function computes term 7,
594
 
   rho_ciab += P(ab) LT_VV(c,a) R(i,b)
595
 
*/
596
 
 
597
 
void x_Gciab_7(void) { 
598
 
  int h, nirreps, c, i, a, b, C, I, A, B, Csym, Isym, Asym, Bsym, row, col;
599
 
  int AA, BB, AAsym, BBsym;
600
 
  int L_irr, R_irr, G_irr;
601
 
  dpdfile2 LT1A, LT1B, R1A, R1B;
602
 
  dpdbuf4 G;
603
 
 
604
 
  L_irr = params.L_irr; R_irr = params.R_irr; G_irr = params.G_irr;
605
 
  nirreps = moinfo.nirreps;
606
 
 
607
 
  /* open one-electron files for the nasty terms */
608
 
  dpd_file2_init(&LT1A, EOM_TMP, L_irr, 1, 1, "LT_VV");
609
 
  dpd_file2_init(&R1A, CC_GR, R_irr, 0, 1, "RIA");
610
 
  if (params.ref == 0 || params.ref == 1) {
611
 
    dpd_file2_init(&LT1B, EOM_TMP, L_irr, 1, 1, "LT_vv");
612
 
    dpd_file2_init(&R1B, CC_GR, R_irr, 0, 1, "Ria");
613
 
  }
614
 
  else {
615
 
    dpd_file2_init(&LT1B, EOM_TMP, L_irr, 3, 3, "LT_vv");
616
 
    dpd_file2_init(&R1B, CC_GR, R_irr, 2, 3, "Ria");
617
 
  }
618
 
  dpd_file2_mat_init(&R1A);   dpd_file2_mat_init(&R1B);
619
 
  dpd_file2_mat_init(&LT1A);   dpd_file2_mat_init(&LT1B);
620
 
  dpd_file2_mat_rd(&R1A);     dpd_file2_mat_rd(&R1B);
621
 
  dpd_file2_mat_rd(&LT1A);     dpd_file2_mat_rd(&LT1B);
622
 
 
623
 
  /* rho_CIAB += LT_VV(C,A) R(I,B) - LT_VV(C,B) R(I,A) */
624
 
  if (params.ref == 0 || params.ref == 1)
625
 
    dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 7, 0, "GCIAB");
626
 
  else
627
 
    dpd_buf4_init(&G, EOM_TMP0, G_irr, 21, 5, 21, 7, 0, "GCIAB");
628
 
  for(h=0; h < nirreps; h++) {
629
 
    dpd_buf4_mat_irrep_init(&G, h);
630
 
    dpd_buf4_mat_irrep_rd(&G, h);
631
 
    for(row=0; row < G.params->rowtot[h]; row++) {
632
 
      c = G.params->roworb[h][row][0];
633
 
      i = G.params->roworb[h][row][1];
634
 
      C = LT1A.params->rowidx[c]; Csym = LT1A.params->psym[c];
635
 
      I = R1A.params->rowidx[i]; Isym = R1A.params->psym[i];
636
 
      for(col=0; col < G.params->coltot[h^G_irr]; col++) {
637
 
        a = G.params->colorb[h^G_irr][col][0];
638
 
        b = G.params->colorb[h^G_irr][col][1];
639
 
        A = LT1A.params->colidx[a]; Asym = LT1A.params->qsym[a];
640
 
        B = R1A.params->colidx[b]; Bsym = R1A.params->qsym[b];
641
 
        AA = R1A.params->colidx[a]; AAsym = R1A.params->qsym[a];
642
 
        BB = LT1A.params->colidx[b]; BBsym = LT1A.params->qsym[b];
643
 
        if( ((Csym^Asym)==L_irr) && ((Isym^Bsym)==R_irr) )
644
 
          G.matrix[h][row][col] +=
645
 
            LT1A.matrix[Csym][C][A] * R1A.matrix[Isym][I][B];
646
 
        if( ((Csym^BBsym)==L_irr) && ((Isym^AAsym)==R_irr) )
647
 
          G.matrix[h][row][col] -=
648
 
            LT1A.matrix[Csym][C][BB] * R1A.matrix[Isym][I][AA];
649
 
      }
650
 
    }
651
 
    dpd_buf4_mat_irrep_wrt(&G, h);
652
 
    dpd_buf4_mat_irrep_close(&G, h);
653
 
  }
654
 
  dpd_buf4_close(&G);
655
 
  /* rho_ciab += LT_vv(c,a) R(i,b) - LT_vv(c,b) R(i,a) */
656
 
  if (params.ref == 0 || params.ref == 1)
657
 
    dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 7, 0, "Gciab");
658
 
  else
659
 
    dpd_buf4_init(&G, EOM_TMP0, G_irr, 31, 15, 31, 17, 0, "Gciab");
660
 
  for(h=0; h < nirreps; h++) {
661
 
    dpd_buf4_mat_irrep_init(&G, h);
662
 
    dpd_buf4_mat_irrep_rd(&G, h);
663
 
    for(row=0; row < G.params->rowtot[h]; row++) {
664
 
      c = G.params->roworb[h][row][0];
665
 
      i = G.params->roworb[h][row][1];
666
 
      C = LT1B.params->rowidx[c]; Csym = LT1B.params->psym[c];
667
 
      I = R1B.params->rowidx[i]; Isym = R1B.params->psym[i];
668
 
      for(col=0; col < G.params->coltot[h^G_irr]; col++) {
669
 
        a = G.params->colorb[h^G_irr][col][0];
670
 
        b = G.params->colorb[h^G_irr][col][1];
671
 
        A = LT1B.params->colidx[a]; Asym = LT1B.params->qsym[a];
672
 
        B = R1B.params->colidx[b]; Bsym = R1B.params->qsym[b];
673
 
        AA = R1B.params->colidx[a]; AAsym = R1B.params->qsym[a];
674
 
        BB = LT1B.params->colidx[b]; BBsym = LT1B.params->qsym[b];
675
 
        if( ((Csym^Asym)==L_irr) && ((Isym^Bsym)==R_irr) )
676
 
          G.matrix[h][row][col] +=
677
 
            LT1B.matrix[Csym][C][A] * R1B.matrix[Isym][I][B];
678
 
        if( ((Csym^BBsym)==L_irr) && ((Isym^AAsym)==R_irr) )
679
 
          G.matrix[h][row][col] -=
680
 
            LT1B.matrix[Csym][C][BB] * R1B.matrix[Isym][I][AA];
681
 
      }
682
 
    }
683
 
    dpd_buf4_mat_irrep_wrt(&G, h);
684
 
    dpd_buf4_mat_irrep_close(&G, h);
685
 
  }
686
 
  dpd_buf4_close(&G);
687
 
  /* rho_CiAb += LT_VV(C,A) R(i,b) */
688
 
  if (params.ref == 0 || params.ref == 1)
689
 
    dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 5, 0, "GCiAb");
690
 
  else
691
 
    dpd_buf4_init(&G, EOM_TMP0, G_irr, 26, 28, 26, 28, 0, "GCiAb");
692
 
  for(h=0; h < nirreps; h++) {
693
 
    dpd_buf4_mat_irrep_init(&G, h);
694
 
    dpd_buf4_mat_irrep_rd(&G, h);
695
 
    for(row=0; row < G.params->rowtot[h]; row++) {
696
 
      c = G.params->roworb[h][row][0];
697
 
      i = G.params->roworb[h][row][1];
698
 
      C = LT1A.params->rowidx[c]; Csym = LT1A.params->psym[c];
699
 
      I = R1B.params->rowidx[i]; Isym = R1B.params->psym[i];
700
 
      for(col=0; col < G.params->coltot[h^G_irr]; col++) {
701
 
        a = G.params->colorb[h^G_irr][col][0];
702
 
        b = G.params->colorb[h^G_irr][col][1];
703
 
        A = LT1A.params->colidx[a]; Asym = LT1A.params->qsym[a];
704
 
        B = R1B.params->colidx[b]; Bsym = R1B.params->qsym[b];
705
 
        if( ((Csym^Asym)==L_irr) && ((Isym^Bsym)==R_irr) )
706
 
          G.matrix[h][row][col] +=
707
 
            LT1A.matrix[Csym][C][A] * R1B.matrix[Isym][I][B];
708
 
      }
709
 
    }
710
 
    dpd_buf4_mat_irrep_wrt(&G, h);
711
 
    dpd_buf4_mat_irrep_close(&G, h);
712
 
  }
713
 
  dpd_buf4_close(&G);
714
 
  /* rho_cIaB += LT_vv(c,a) R(I,B) */
715
 
  if (params.ref == 0 || params.ref == 1)
716
 
    dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 5, 0, "GcIaB");
717
 
  else
718
 
    dpd_buf4_init(&G, EOM_TMP0, G_irr, 25, 29, 25, 29, 0, "GcIaB");
719
 
  for(h=0; h < nirreps; h++) {
720
 
    dpd_buf4_mat_irrep_init(&G, h);
721
 
    dpd_buf4_mat_irrep_rd(&G, h);
722
 
    for(row=0; row < G.params->rowtot[h]; row++) {
723
 
      c = G.params->roworb[h][row][0];
724
 
      i = G.params->roworb[h][row][1];
725
 
      C = LT1B.params->rowidx[c]; Csym = LT1B.params->psym[c];
726
 
      I = R1A.params->rowidx[i]; Isym = R1A.params->psym[i];
727
 
      for(col=0; col < G.params->coltot[h^G_irr]; col++) {
728
 
        a = G.params->colorb[h^G_irr][col][0];
729
 
        b = G.params->colorb[h^G_irr][col][1];
730
 
        A = LT1B.params->colidx[a]; Asym = LT1B.params->qsym[a];
731
 
        B = R1A.params->colidx[b]; Bsym = R1A.params->qsym[b];
732
 
        if( ((Csym^Asym)==L_irr) && ((Isym^Bsym)==R_irr) )
733
 
          G.matrix[h][row][col] +=
734
 
            LT1B.matrix[Csym][C][A] * R1A.matrix[Isym][I][B];
735
 
      }
736
 
    }
737
 
    dpd_buf4_mat_irrep_wrt(&G, h);
738
 
    dpd_buf4_mat_irrep_close(&G, h);
739
 
  }
740
 
  dpd_buf4_close(&G);
741
 
 
742
 
  dpd_file2_mat_close(&LT1A);
743
 
  dpd_file2_mat_close(&LT1B);
744
 
  dpd_file2_close(&LT1A);
745
 
  dpd_file2_close(&LT1B);
746
 
 
747
 
  dpd_file2_mat_close(&R1A);
748
 
  dpd_file2_mat_close(&R1B);
749
 
  dpd_file2_close(&R1A);
750
 
  dpd_file2_close(&R1B);
751
 
 
752
 
  return;
753
 
}
754
 
 
755
 
 
756
 
 
757
 
/* This function computes terms of excited Gciab
758
 
   term 8  +P(ab) Lmnce Rinae Tmb
759
 
   term 9, +P(AB) LMNCE TINAE RMB
760
 
*/
761
 
 
762
 
void x_Gciab_8_rohf(void) { 
763
 
  int h, nirreps, i, j, k, a, I, J, K, A, Isym, Jsym, Ksym, Asym, row, col;
764
 
  int II,JJ,IIsym,JJsym;
765
 
  int L_irr, R_irr, G_irr;
766
 
  double value;
767
 
  dpdfile2 L1A, T1A, L1B, T1B, R1A, R1B, I1A, I1B;
768
 
  dpdbuf4 G, V, T, L, Z, Z1, Z2, Tau;
769
 
 
770
 
  L_irr = params.L_irr; R_irr = params.R_irr; G_irr = params.G_irr;
771
 
  nirreps = moinfo.nirreps;
772
 
 
773
 
  /* term 8, +P(AB) LMNCE RINAE TMB */
774
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 10, 5, 10, 5, 0, "Z(IA,BC)");
775
 
  dpd_buf4_init(&V, EOM_TMP, G_irr, 10, 10, 10, 10, 0, "R2L2_OVOV");
776
 
  dpd_file2_init(&T1A, CC_OEI, 0, 0, 1, "tIA");
777
 
  dpd_contract244(&T1A, &V, &Z, 0, 2, 1, 1.0, 0.0); 
778
 
  dpd_file2_close(&T1A);
779
 
  dpd_buf4_close(&V);
780
 
  dpd_buf4_sort(&Z, EOM_TMP1, sprq, 11, 5, "Z(CI,BA)");
781
 
  dpd_buf4_close(&Z);
782
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 5, 11, 5, 0, "Z(CI,BA)");
783
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 7, 0, "GCIAB");
784
 
  dpd_buf4_axpy(&Z, &G, 1.0);
785
 
  dpd_buf4_sort(&Z, EOM_TMP1, pqsr, 11, 5, "Z(CI,AB)");
786
 
  dpd_buf4_close(&Z);
787
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 5, 11, 5, 0, "Z(CI,AB)");
788
 
  dpd_buf4_axpy(&Z, &G, -1.0);
789
 
  dpd_buf4_close(&G);
790
 
  dpd_buf4_close(&Z);
791
 
  /* term 8, +P(ab) Lmnce Rinae Tmb */
792
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 10, 5, 10, 5, 0, "Z(ia,bc)");
793
 
  dpd_buf4_init(&V, EOM_TMP, G_irr, 10, 10, 10, 10, 0, "R2L2_ovov");
794
 
  dpd_file2_init(&T1A, CC_OEI, 0, 0, 1, "tia");
795
 
  dpd_contract244(&T1A, &V, &Z, 0, 2, 1, 1.0, 0.0); 
796
 
  dpd_file2_close(&T1A);
797
 
  dpd_buf4_close(&V);
798
 
  dpd_buf4_sort(&Z, EOM_TMP1, sprq, 11, 5, "Z(ci,ba)");
799
 
  dpd_buf4_close(&Z);
800
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 5, 11, 5, 0, "Z(ci,ba)");
801
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 7, 0, "Gciab");
802
 
  dpd_buf4_axpy(&Z, &G, 1.0);
803
 
  dpd_buf4_sort(&Z, EOM_TMP1, pqsr, 11, 5, "Z(ci,ab)");
804
 
  dpd_buf4_close(&Z);
805
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 5, 11, 5, 0, "Z(ci,ab)");
806
 
  dpd_buf4_axpy(&Z, &G, -1.0);
807
 
  dpd_buf4_close(&G);
808
 
  dpd_buf4_close(&Z);
809
 
  
810
 
  /* term 8, GCiAb -= LmNCe RiNAe tmb */
811
 
  dpd_buf4_init(&V, EOM_TMP, G_irr, 10, 10, 10, 10, 0, "R2L2_oVoV");
812
 
  dpd_buf4_sort(&V, EOM_TMP1, spqr, 11, 11, "Z(Ci,Am)");
813
 
  dpd_buf4_close(&V);
814
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 11, 11, 11, 0, "Z(Ci,Am)");
815
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 5, 0, "GCiAb");
816
 
  dpd_file2_init(&T1B, CC_OEI, 0, 0, 1, "tia");
817
 
  dpd_contract424(&Z, &T1B, &G, 3, 0, 0, -1.0, 1.0); 
818
 
  dpd_file2_close(&T1B);
819
 
  dpd_buf4_close(&G);
820
 
  dpd_buf4_close(&Z);
821
 
  /* term 8, GCiAb -= (LMnCe Rinbe + LMNCE RiNbE) TMA */
822
 
  dpd_buf4_init(&V, EOM_TMP, G_irr, 10, 10, 10, 10, 0, "R2L2_ovOV");
823
 
  dpd_buf4_sort(&V, EOM_TMP1, sprq, 11, 10, "Z(Ci,Mb)");
824
 
  dpd_buf4_close(&V);
825
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 10, 11, 10, 0, "Z(Ci,Mb)");
826
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 5, 0, "GCiAb");
827
 
  dpd_file2_init(&T1A, CC_OEI, 0, 0, 1, "tIA");
828
 
  dpd_contract244(&T1A, &Z, &G, 0, 2, 1, 1.0, 1.0); 
829
 
  dpd_file2_close(&T1A);
830
 
  dpd_buf4_close(&G);
831
 
  dpd_buf4_close(&Z);
832
 
 
833
 
  /* term 8, GcIaB -= LMncE RInaE tMB */
834
 
  dpd_buf4_init(&V, EOM_TMP, G_irr, 10, 10, 10, 10, 0, "R2L2_OvOv");
835
 
  dpd_buf4_sort(&V, EOM_TMP1, spqr, 11, 11, "Z(cI,aM)");
836
 
  dpd_buf4_close(&V);
837
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 11, 11, 11, 0, "Z(cI,aM)");
838
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 5, 0, "GcIaB");
839
 
  dpd_file2_init(&T1A, CC_OEI, 0, 0, 1, "tIA");
840
 
  dpd_contract424(&Z, &T1A, &G, 3, 0, 0, -1.0, 1.0); 
841
 
  dpd_file2_close(&T1A);
842
 
  dpd_buf4_close(&G);
843
 
  dpd_buf4_close(&Z);
844
 
  /* term 8, GcIaB -= (LmNcE RINBE + Lmnce RInBe) Tma */
845
 
  dpd_buf4_init(&V, EOM_TMP, G_irr, 10, 10, 10, 10, 0, "R2L2_OVov");
846
 
  dpd_buf4_sort(&V, EOM_TMP1, sprq, 11, 10, "Z(cI,mB)");
847
 
  dpd_buf4_close(&V);
848
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 10, 11, 10, 0, "Z(cI,mB)");
849
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 5, 0, "GcIaB");
850
 
  dpd_file2_init(&T1B, CC_OEI, 0, 0, 1, "tia");
851
 
  dpd_contract244(&T1B, &Z, &G, 0, 2, 1, 1.0, 1.0); 
852
 
  dpd_file2_close(&T1B);
853
 
  dpd_buf4_close(&G);
854
 
  dpd_buf4_close(&Z);
855
 
 
856
 
  psio_close(EOM_TMP1, 0);
857
 
  psio_open(EOM_TMP1, PSIO_OPEN_NEW);
858
 
 
859
 
  /* term 9, +P(AB) LMNCE TINAE RMB */
860
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 10, 5, 10, 5, 0, "Z(IA,BC)");
861
 
  dpd_buf4_init(&V, EOM_TMP, L_irr, 10, 10, 10, 10, 0, "VIAJB");
862
 
  dpd_file2_init(&R1A, CC_GR, R_irr, 0, 1, "RIA");
863
 
  dpd_contract244(&R1A, &V, &Z, 0, 2, 1, 1.0, 0.0); 
864
 
  dpd_file2_close(&R1A);
865
 
  dpd_buf4_close(&V);
866
 
  dpd_buf4_sort(&Z, EOM_TMP1, sprq, 11, 5, "Z(CI,BA)");
867
 
  dpd_buf4_close(&Z);
868
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 5, 11, 5, 0, "Z(CI,BA)");
869
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 7, 0, "GCIAB");
870
 
  dpd_buf4_axpy(&Z, &G, 1.0);
871
 
  dpd_buf4_sort(&Z, EOM_TMP1, pqsr, 11, 5, "Z(CI,AB)");
872
 
  dpd_buf4_close(&Z);
873
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 5, 11, 5, 0, "Z(CI,AB)");
874
 
  dpd_buf4_axpy(&Z, &G, -1.0);
875
 
  dpd_buf4_close(&G);
876
 
  dpd_buf4_close(&Z);
877
 
  /* term 9, +P(ab) Lmnce Tinae Rmb */
878
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 10, 5, 10, 5, 0, "Z(ia,bc)");
879
 
  dpd_buf4_init(&V, EOM_TMP, L_irr, 10, 10, 10, 10, 0, "Viajb");
880
 
  dpd_file2_init(&R1B, CC_GR, R_irr, 0, 1, "Ria");
881
 
  dpd_contract244(&R1B, &V, &Z, 0, 2, 1, 1.0, 0.0); 
882
 
  dpd_file2_close(&R1B);
883
 
  dpd_buf4_close(&V);
884
 
  dpd_buf4_sort(&Z, EOM_TMP1, sprq, 11, 5, "Z(ci,ba)");
885
 
  dpd_buf4_close(&Z);
886
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 5, 11, 5, 0, "Z(ci,ba)");
887
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 7, 0, "Gciab");
888
 
  dpd_buf4_axpy(&Z, &G, 1.0);
889
 
  dpd_buf4_sort(&Z, EOM_TMP1, pqsr, 11, 5, "Z(ci,ab)");
890
 
  dpd_buf4_close(&Z);
891
 
  dpd_buf4_init(&Z, EOM_TMP1, G_irr, 11, 5, 11, 5, 0, "Z(ci,ab)");
892
 
  dpd_buf4_axpy(&Z, &G, -1.0);
893
 
  dpd_buf4_close(&G);
894
 
  dpd_buf4_close(&Z);
895
 
  
896
 
  /* term 9, GCiAb -= LmNCe TiNAe Rmb */
897
 
  dpd_buf4_init(&V, EOM_TMP, L_irr, 10, 10, 10, 10, 0, "ViAjB");
898
 
  dpd_buf4_sort(&V, EOM_TMP1, spqr, 11, 11, "Z(Ci,Am)");
899
 
  dpd_buf4_close(&V);
900
 
  dpd_buf4_init(&Z, EOM_TMP1, L_irr, 11, 11, 11, 11, 0, "Z(Ci,Am)");
901
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 5, 0, "GCiAb");
902
 
  dpd_file2_init(&R1B, CC_GR, R_irr, 0, 1, "Ria");
903
 
  dpd_contract424(&Z, &R1B, &G, 3, 0, 0, -1.0, 1.0); 
904
 
  dpd_file2_close(&R1B);
905
 
  dpd_buf4_close(&G);
906
 
  dpd_buf4_close(&Z);
907
 
  /* term 9, GCiAb -= (LMnCe Tinbe + LMNCE TiNbE) RMA */
908
 
  dpd_buf4_init(&V, EOM_TMP, L_irr, 10, 10, 10, 10, 0, "ViaJB");
909
 
  dpd_buf4_sort(&V, EOM_TMP1, sprq, 11, 10, "Z(Ci,Mb)");
910
 
  dpd_buf4_close(&V);
911
 
  dpd_buf4_init(&Z, EOM_TMP1, L_irr, 11, 10, 11, 10, 0, "Z(Ci,Mb)");
912
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 5, 0, "GCiAb");
913
 
  dpd_file2_init(&R1A, CC_GR, R_irr, 0, 1, "RIA");
914
 
  dpd_contract244(&R1A, &Z, &G, 0, 2, 1, 1.0, 1.0); 
915
 
  dpd_file2_close(&R1A);
916
 
  dpd_buf4_close(&G);
917
 
  dpd_buf4_close(&Z);
918
 
 
919
 
  /* term 9, GcIaB -= LMncE TInaE RMB */
920
 
  dpd_buf4_init(&V, EOM_TMP, L_irr, 10, 10, 10, 10, 0, "VIaJb");
921
 
  dpd_buf4_sort(&V, EOM_TMP1, spqr, 11, 11, "Z(cI,aM)");
922
 
  dpd_buf4_close(&V);
923
 
  dpd_buf4_init(&Z, EOM_TMP1, L_irr, 11, 11, 11, 11, 0, "Z(cI,aM)");
924
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 5, 0, "GcIaB");
925
 
  dpd_file2_init(&R1A, CC_GR, R_irr, 0, 1, "RIA");
926
 
  dpd_contract424(&Z, &R1A, &G, 3, 0, 0, -1.0, 1.0); 
927
 
  dpd_file2_close(&R1A);
928
 
  dpd_buf4_close(&G);
929
 
  dpd_buf4_close(&Z);
930
 
  /* term 9, GcIaB -= (LmNcE TINBE + Lmnce TInBe) Rma */
931
 
  dpd_buf4_init(&V, EOM_TMP, L_irr, 10, 10, 10, 10, 0, "VIAjb");
932
 
  dpd_buf4_sort(&V, EOM_TMP1, sprq, 11, 10, "Z(cI,mB)");
933
 
  dpd_buf4_close(&V);
934
 
  dpd_buf4_init(&Z, EOM_TMP1, L_irr, 11, 10, 11, 10, 0, "Z(cI,mB)");
935
 
  dpd_buf4_init(&G, EOM_TMP0, G_irr, 11, 5, 11, 5, 0, "GcIaB");
936
 
  dpd_file2_init(&R1B, CC_GR, R_irr, 0, 1, "Ria");
937
 
  dpd_contract244(&R1B, &Z, &G, 0, 2, 1, 1.0, 1.0); 
938
 
  dpd_file2_close(&R1B);
939
 
  dpd_buf4_close(&G);
940
 
  dpd_buf4_close(&Z);
941
 
 
942
 
  return;
943
 
}