3
#include <libdpd/dpd.h>
12
int *vir_off, *virtpi;
19
dpdbuf4 T2, E, F, T2AA, T2AB, T2BA, EAA, EAB, EBA, FAA, FAB, FBA;
20
dpdfile2 fIJ, fAB, fij, fab;
23
int Gij, ij, Gbc, bc, Gjk, jk;
25
int **W_offset, offset;
27
nirreps = moinfo.nirreps;
29
occ_off = moinfo.occ_off;
30
virtpi = moinfo.virtpi;
31
vir_off = moinfo.vir_off;
33
W_offset = init_int_matrix(nirreps, nirreps);
34
for(Gab=0; Gab < nirreps; Gab++) {
35
for(Ga=0,offset=0; Ga < nirreps; Ga++) {
37
W_offset[Gab][Ga] = offset;
38
offset += virtpi[Ga] * virtpi[Gb];
42
dpd_file2_init(&XLD, CC3_MISC, 0, 0, 1, "CC3 XLD");
43
dpd_file2_mat_init(&XLD);
45
dpd_buf4_init(&L2, CC_LAMBDA, 0, 0, 5, 2, 7, 0, "LIJAB");
46
dpd_buf4_init(&L2AB, CC_LAMBDA, 0, 0, 5, 0, 5, 0, "LIjAb");
47
for(h=0; h < nirreps; h++) {
48
dpd_buf4_mat_irrep_init(&L2, h);
49
dpd_buf4_mat_irrep_rd(&L2, h);
51
dpd_buf4_mat_irrep_init(&L2AB, h);
52
dpd_buf4_mat_irrep_rd(&L2AB, h);
55
dpd_file2_init(&fIJ, CC_OEI, 0, 0, 0, "fIJ");
56
dpd_file2_init(&fAB, CC_OEI, 0, 1, 1, "fAB");
57
dpd_file2_init(&fij, CC_OEI, 0, 0, 0, "fij");
58
dpd_file2_init(&fab, CC_OEI, 0, 1, 1, "fab");
60
dpd_buf4_init(&T2, CC_TAMPS, 0, 0, 5, 2, 7, 0, "tIJAB");
61
dpd_buf4_init(&F, CC3_HET1, 0, 10, 5, 10, 7, 0, "CC3 WABEI (IE,B>A)");
62
dpd_buf4_init(&E, CC3_HET1, 0, 0, 10, 2, 10, 0, "CC3 WMBIJ (I>J,MB)");
65
dpd_buf4_init(&T2AB, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tIjAb");
66
dpd_buf4_init(&T2BA, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tiJaB");
68
dpd_buf4_init(&FAB, CC3_HET1, 0, 10, 5, 10, 5, 0, "CC3 WaBeI (Ie,Ba)");
69
dpd_buf4_init(&FBA, CC3_HET1, 0, 10, 5, 10, 5, 0, "CC3 WAbEi (iE,bA)");
71
dpd_buf4_init(&EAB, CC3_HET1, 0, 0, 10, 0, 10, 0, "CC3 WMbIj (Ij,Mb)");
72
dpd_buf4_init(&EBA, CC3_HET1, 0, 0, 10, 0, 10, 0, "CC3 WmBiJ (iJ,mB)");
74
/* target T3 amplitudes go in here */
75
W1 = (double ***) malloc(nirreps * sizeof(double **));
77
for(Gi=0; Gi < nirreps; Gi++) {
78
for(Gj=0; Gj < nirreps; Gj++) {
80
for(Gk=0; Gk < nirreps; Gk++) {
84
for(Gab=0; Gab < nirreps; Gab++) {
85
Gc = Gab ^ Gijk; /* totally symmetric */
86
W1[Gab] = dpd_block_matrix(F.params->coltot[Gab], virtpi[Gc]);
89
for(i=0; i < occpi[Gi]; i++) {
91
for(j=0; j < occpi[Gj]; j++) {
93
for(k=0; k < occpi[Gk]; k++) {
96
T3_AAA(W1, nirreps, I, Gi, J, Gj, K, Gk, &T2, &F, &E, &fIJ, &fAB,
97
occpi, occ_off, virtpi, vir_off, 0.0);
99
/* X_KC <-- 1/4 t_IJKABC <IJ||AB> */
101
Gc = Gk; /* assumes T1 is totally symmetric */
102
Gab = Gij; /* assumes <ij||ab> is totally symmetric */
104
ij = L2.params->rowidx[I][J];
106
nrows = L2.params->coltot[Gij];
110
C_DGEMV('t', nrows, ncols, 0.25, W1[Gab][0], ncols, L2.matrix[Gij][ij], 1,
111
1.0, XLD.matrix[Gk][k], 1);
113
T3_AAB(W1, nirreps, I, Gi, J, Gj, K, Gk, &T2AA, &T2AB, &T2BA,
114
&FAA, &FAB, &FBA, &EAA, &EAB, &EBA, &fIJ, &fij, &fAB, &fab,
115
occpi, occ_off, occpi, occ_off, virtpi, vir_off, virtpi, vir_off, 0.0);
117
/* t_IA <-- t_IJkABc <Jk|Bc> */
119
Ga = Gi; /* assumes T1 is totally symmetric */
120
Gbc = Gjk; /* assumes <jk|bc> is totally symmetric */
122
jk = L2AB.params->rowidx[J][K];
124
for(Gab=0; Gab < nirreps; Gab++) {
128
ab = W_offset[Gab][Ga];
129
bc = L2AB.col_offset[Gjk][Gb];
132
ncols = virtpi[Gb] * virtpi[Gc];
135
C_DGEMV('n', nrows, ncols, 1.0, W1[Gab][ab], ncols, &(L2AB.matrix[Gjk][jk][bc]), 1,
136
1.0, XLD.matrix[Gi][i], 1);
139
/* t_KC <-- 1/4 t_ijKabC <ij||ab> */
141
Gc = Gk; /* assumes T1 is totally symmetric */
142
Gab = Gij; /* assumes <ij||ab> is totally symmetric */
144
ij = L2.params->rowidx[I][J];
146
nrows = L2.params->coltot[Gij];
150
C_DGEMV('t', nrows, ncols, 0.25, W1[Gab][0], ncols, L2.matrix[Gij][ij], 1,
151
1.0, XLD.matrix[Gk][k], 1);
158
for(Gab=0; Gab < nirreps; Gab++) {
159
Gc = Gab ^ Gijk; /* totally symmetric */
160
dpd_free_block(W1[Gab], F.params->coltot[Gab], virtpi[Gc]);
171
dpd_file2_close(&fIJ);
172
dpd_file2_close(&fAB);
174
dpd_buf4_close(&EAB);
175
dpd_buf4_close(&EBA);
176
dpd_buf4_close(&FAB);
177
dpd_buf4_close(&FBA);
178
dpd_buf4_close(&T2AB);
179
dpd_buf4_close(&T2BA);
180
dpd_file2_close(&fij);
181
dpd_file2_close(&fab);
183
free_int_matrix(W_offset, nirreps);
185
for(h=0; h < nirreps; h++) {
186
dpd_buf4_mat_irrep_close(&L2, h);
187
dpd_buf4_mat_irrep_close(&L2AB, h);
190
dpd_buf4_close(&L2AB);
192
dpd_file2_mat_wrt(&XLD);
193
dpd_file2_close(&XLD);