1
#include <libdpd/dpd.h>
5
/* build_A(): Builds the CIS response matrix A, which is expressed in
8
** A(ai,bj) = delta_ij f_ab - delta_ab f_ij - <ja||ib>
10
** This routine will build each spin component of the entire matrix
11
** for later in-core diagonalization.
14
** A(AI,BJ) = delta_IJ f_AB - delta_AB f_IJ + 2 <IJ|AB> - <IA|JB>
16
** UHF and ROHF references:
17
** A(AI,BJ) = delta_IJ f_AB - delta_AB f_IJ - <JA||IB>
18
** A(ai,bj) = delta_ij f_ab - delta_ab f_ij - <ja||ib>
27
int a, b, i, j, ai, bj, A, B, I, J, Asym, Bsym, Isym, Jsym;
28
dpdbuf4 C, D, Amat, A_AA, A_BB, A_AB;
29
dpdfile2 fIJ, fij, fAB, fab;
31
nirreps = moinfo.nirreps;
33
if(params.ref == 0) { /** RHF **/
34
dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
35
dpd_buf4_sort(&D, CC_MISC, rpsq, 11, 11, "A(AI,BJ)");
37
dpd_buf4_init(&Amat, CC_MISC, 0, 11, 11, 11, 11, 0, "A(AI,BJ)");
38
dpd_buf4_scm(&Amat, 2);
39
dpd_buf4_close(&Amat);
40
dpd_buf4_init(&C, CC_CINTS, 0, 10, 10, 10, 10, 0, "C <ia|jb>");
41
dpd_buf4_sort_axpy(&C, CC_MISC, qpsr, 11, 11, "A(AI,BJ)", -1);
42
dpd_buf4_sort(&C, CC_MISC, qpsr, 11, 11, "A(AI,BJ) triplet");
44
dpd_buf4_init(&Amat, CC_MISC, 0, 11, 11, 11, 11, 0, "A(AI,BJ) triplet");
45
dpd_buf4_scm(&Amat, -1);
46
dpd_buf4_close(&Amat);
48
dpd_file2_init(&fIJ, CC_OEI, 0, 0, 0, "fIJ");
49
dpd_file2_mat_init(&fIJ);
50
dpd_file2_mat_rd(&fIJ);
51
dpd_file2_init(&fij, CC_OEI, 0, 0, 0, "fij");
52
dpd_file2_mat_init(&fij);
53
dpd_file2_mat_rd(&fij);
54
dpd_file2_init(&fAB, CC_OEI, 0, 1, 1, "fAB");
55
dpd_file2_mat_init(&fAB);
56
dpd_file2_mat_rd(&fAB);
57
dpd_file2_init(&fab, CC_OEI, 0, 1, 1, "fab");
58
dpd_file2_mat_init(&fab);
59
dpd_file2_mat_rd(&fab);
61
dpd_buf4_init(&Amat, CC_MISC, 0, 11, 11, 11, 11, 0, "A(AI,BJ)");
62
for(h=0; h < nirreps; h++) {
63
dpd_buf4_mat_irrep_init(&Amat, h);
64
dpd_buf4_mat_irrep_rd(&Amat, h);
65
for(ai=0; ai < Amat.params->rowtot[h]; ai++) {
66
a = Amat.params->roworb[h][ai][0];
67
i = Amat.params->roworb[h][ai][1];
68
A = fAB.params->rowidx[a];
69
I = fIJ.params->rowidx[i];
70
Asym = fAB.params->psym[a];
71
Isym = fIJ.params->psym[i];
72
for(bj=0; bj < Amat.params->coltot[h]; bj++) {
73
b = Amat.params->colorb[h][bj][0];
74
j = Amat.params->colorb[h][bj][1];
75
B = fAB.params->colidx[b];
76
J = fIJ.params->colidx[j];
77
Bsym = fAB.params->qsym[b];
78
Jsym = fIJ.params->qsym[j];
79
if((A==B) && (Isym==Jsym)) Amat.matrix[h][ai][bj] -= fIJ.matrix[Isym][I][J];
80
if((I==J) && (Asym==Bsym)) Amat.matrix[h][ai][bj] += fAB.matrix[Asym][A][B];
83
dpd_buf4_mat_irrep_wrt(&Amat, h);
84
dpd_buf4_mat_irrep_close(&Amat, h);
86
dpd_buf4_close(&Amat);
88
dpd_buf4_init(&Amat, CC_MISC, 0, 11, 11, 11, 11, 0, "A(AI,BJ) triplet");
89
for(h=0; h < nirreps; h++) {
90
dpd_buf4_mat_irrep_init(&Amat, h);
91
dpd_buf4_mat_irrep_rd(&Amat, h);
92
for(ai=0; ai < Amat.params->rowtot[h]; ai++) {
93
a = Amat.params->roworb[h][ai][0];
94
i = Amat.params->roworb[h][ai][1];
95
A = fAB.params->rowidx[a];
96
I = fIJ.params->rowidx[i];
97
Asym = fAB.params->psym[a];
98
Isym = fIJ.params->psym[i];
99
for(bj=0; bj < Amat.params->coltot[h]; bj++) {
100
b = Amat.params->colorb[h][bj][0];
101
j = Amat.params->colorb[h][bj][1];
102
B = fAB.params->colidx[b];
103
J = fIJ.params->colidx[j];
104
Bsym = fAB.params->qsym[b];
105
Jsym = fIJ.params->qsym[j];
106
if((A==B) && (Isym==Jsym)) Amat.matrix[h][ai][bj] -= fIJ.matrix[Isym][I][J];
107
if((I==J) && (Asym==Bsym)) Amat.matrix[h][ai][bj] += fAB.matrix[Asym][A][B];
110
dpd_buf4_mat_irrep_wrt(&Amat, h);
111
dpd_buf4_mat_irrep_close(&Amat, h);
113
dpd_buf4_close(&Amat);
115
dpd_file2_mat_close(&fab);
116
dpd_file2_close(&fab);
117
dpd_file2_mat_close(&fAB);
118
dpd_file2_close(&fAB);
119
dpd_file2_mat_close(&fij);
120
dpd_file2_close(&fij);
121
dpd_file2_mat_close(&fIJ);
122
dpd_file2_close(&fIJ);
124
else if(params.ref == 2) { /** UHF **/
126
/** <JA||IB> --> A(AI,BJ) **/
127
dpd_buf4_init(&C, CC_CINTS, 0, 20, 20, 20, 20, 0, "C <IA||JB>");
128
dpd_buf4_sort(&C, CC_MISC, qrsp, 21, 21, "A(AI,BJ)");
130
dpd_buf4_init(&Amat, CC_MISC, 0, 21, 21, 21, 21, 0, "A(AI,BJ)");
131
dpd_buf4_scm(&Amat, -1);
132
dpd_buf4_close(&Amat);
134
/** <ja||ib> --> A(ai,bj) **/
135
dpd_buf4_init(&C, CC_CINTS, 0, 30, 30, 30, 30, 0, "C <ia||jb>");
136
dpd_buf4_sort(&C, CC_MISC, qrsp, 31, 31, "A(ai,bj)");
138
dpd_buf4_init(&Amat, CC_MISC, 0, 31, 31, 31, 31, 0, "A(ai,bj)");
139
dpd_buf4_scm(&Amat, -1);
140
dpd_buf4_close(&Amat);
142
/** <Ij|Ab> --> A(AI,bj) **/
143
dpd_buf4_init(&D, CC_DINTS, 0, 22, 28, 22, 28, 0, "D <Ij|Ab>");
144
dpd_buf4_sort(&D, CC_MISC, rpsq, 21, 31, "A(AI,bj)");
147
dpd_file2_init(&fIJ, CC_OEI, 0, 0, 0, "fIJ");
148
dpd_file2_mat_init(&fIJ);
149
dpd_file2_mat_rd(&fIJ);
150
dpd_file2_init(&fij, CC_OEI, 0, 2, 2, "fij");
151
dpd_file2_mat_init(&fij);
152
dpd_file2_mat_rd(&fij);
153
dpd_file2_init(&fAB, CC_OEI, 0, 1, 1, "fAB");
154
dpd_file2_mat_init(&fAB);
155
dpd_file2_mat_rd(&fAB);
156
dpd_file2_init(&fab, CC_OEI, 0, 3, 3, "fab");
157
dpd_file2_mat_init(&fab);
158
dpd_file2_mat_rd(&fab);
160
dpd_buf4_init(&Amat, CC_MISC, 0, 21, 21, 21, 21, 0, "A(AI,BJ)");
161
for(h=0; h < nirreps; h++) {
162
dpd_buf4_mat_irrep_init(&Amat, h);
163
dpd_buf4_mat_irrep_rd(&Amat, h);
164
for(ai=0; ai < Amat.params->rowtot[h]; ai++) {
165
a = Amat.params->roworb[h][ai][0];
166
i = Amat.params->roworb[h][ai][1];
167
A = fAB.params->rowidx[a];
168
I = fIJ.params->rowidx[i];
169
Asym = fAB.params->psym[a];
170
Isym = fIJ.params->psym[i];
171
for(bj=0; bj < Amat.params->coltot[h]; bj++) {
172
b = Amat.params->colorb[h][bj][0];
173
j = Amat.params->colorb[h][bj][1];
174
B = fAB.params->colidx[b];
175
J = fIJ.params->colidx[j];
176
Bsym = fAB.params->qsym[b];
177
Jsym = fIJ.params->qsym[j];
178
if((A==B) && (Isym==Jsym)) Amat.matrix[h][ai][bj] -= fIJ.matrix[Isym][I][J];
179
if((I==J) && (Asym==Bsym)) Amat.matrix[h][ai][bj] += fAB.matrix[Asym][A][B];
182
dpd_buf4_mat_irrep_wrt(&Amat, h);
183
dpd_buf4_mat_irrep_close(&Amat, h);
185
dpd_buf4_close(&Amat);
187
dpd_buf4_init(&Amat, CC_MISC, 0, 31, 31, 31, 31, 0, "A(ai,bj)");
188
for(h=0; h < nirreps; h++) {
189
dpd_buf4_mat_irrep_init(&Amat, h);
190
dpd_buf4_mat_irrep_rd(&Amat, h);
191
for(ai=0; ai < Amat.params->rowtot[h]; ai++) {
192
a = Amat.params->roworb[h][ai][0];
193
i = Amat.params->roworb[h][ai][1];
194
A = fab.params->rowidx[a];
195
I = fij.params->rowidx[i];
196
Asym = fab.params->psym[a];
197
Isym = fij.params->psym[i];
198
for(bj=0; bj < Amat.params->coltot[h]; bj++) {
199
b = Amat.params->colorb[h][bj][0];
200
j = Amat.params->colorb[h][bj][1];
201
B = fab.params->colidx[b];
202
J = fij.params->colidx[j];
203
Bsym = fab.params->qsym[b];
204
Jsym = fij.params->qsym[j];
205
if((A==B) && (Isym==Jsym)) Amat.matrix[h][ai][bj] -= fij.matrix[Isym][I][J];
206
if((I==J) && (Asym==Bsym)) Amat.matrix[h][ai][bj] += fab.matrix[Asym][A][B];
209
dpd_buf4_mat_irrep_wrt(&Amat, h);
210
dpd_buf4_mat_irrep_close(&Amat, h);
212
dpd_buf4_close(&Amat);
214
dpd_file2_mat_close(&fab);
215
dpd_file2_close(&fab);
216
dpd_file2_mat_close(&fAB);
217
dpd_file2_close(&fAB);
218
dpd_file2_mat_close(&fij);
219
dpd_file2_close(&fij);
220
dpd_file2_mat_close(&fIJ);
221
dpd_file2_close(&fIJ);