3
\brief Enter brief description of file here
5
#include <libdpd/dpd.h>
12
namespace psi { namespace cis {
14
void denom(int irrep, double root)
17
int ij, ab, i, j, a, b, I, J, A, B, isym, jsym, asym, bsym;
19
int *occpi, *virtpi, *occ_off, *vir_off;
20
int *aoccpi, *avirtpi, *aocc_off, *avir_off;
21
int *boccpi, *bvirtpi, *bocc_off, *bvir_off;
22
double fii, fjj, faa, fbb;
23
dpdfile2 fIJ, fij, fAB, fab;
27
nirreps = moinfo.nirreps;
29
if(params.ref == 0) { /** RHF **/
32
virtpi = moinfo.virtpi;
33
occ_off = moinfo.occ_off;
34
vir_off = moinfo.vir_off;
36
dpd_file2_init(&fIJ, CC_OEI, 0, 0, 0, "fIJ");
37
dpd_file2_mat_init(&fIJ);
38
dpd_file2_mat_rd(&fIJ);
40
dpd_file2_init(&fij, CC_OEI, 0, 0, 0, "fij");
41
dpd_file2_mat_init(&fij);
42
dpd_file2_mat_rd(&fij);
44
dpd_file2_init(&fAB, CC_OEI, 0, 1, 1, "fAB");
45
dpd_file2_mat_init(&fAB);
46
dpd_file2_mat_rd(&fAB);
48
dpd_file2_init(&fab, CC_OEI, 0, 1, 1, "fab");
49
dpd_file2_mat_init(&fab);
50
dpd_file2_mat_rd(&fab);
52
sprintf(lbl, "dIjAb[%d]", irrep);
53
dpd_buf4_init(&D, CC_MISC, irrep, 0, 5, 0, 5, 0, lbl);
54
for(Gij=0; Gij < nirreps; Gij++) {
56
dpd_buf4_mat_irrep_init(&D, Gij);
58
for(ij=0; ij < D.params->rowtot[Gij]; ij++) {
59
i = D.params->roworb[Gij][ij][0];
60
j = D.params->roworb[Gij][ij][1];
61
isym = D.params->psym[i];
62
jsym = D.params->qsym[j];
64
I = i - occ_off[isym];
65
J = j - occ_off[jsym];
67
fii = fIJ.matrix[isym][I][I];
68
fjj = fij.matrix[jsym][J][J];
70
for(ab=0; ab < D.params->coltot[Gij^irrep]; ab++) {
71
a = D.params->colorb[Gij^irrep][ab][0];
72
b = D.params->colorb[Gij^irrep][ab][1];
73
asym = D.params->rsym[a];
74
bsym = D.params->ssym[b];
76
A = a - vir_off[asym];
77
B = b - vir_off[bsym];
79
faa = fAB.matrix[asym][A][A];
80
fbb = fab.matrix[bsym][B][B];
82
D.matrix[Gij][ij][ab] = 1.0/(fii + fjj - faa - fbb + root);
86
dpd_buf4_mat_irrep_wrt(&D, Gij);
87
dpd_buf4_mat_irrep_close(&D, Gij);
94
else if(params.ref == 2) { /** UHF **/
96
aoccpi = moinfo.aoccpi;
97
boccpi = moinfo.boccpi;
98
avirtpi = moinfo.avirtpi;
99
bvirtpi = moinfo.bvirtpi;
100
aocc_off = moinfo.aocc_off;
101
bocc_off = moinfo.bocc_off;
102
avir_off = moinfo.avir_off;
103
bvir_off = moinfo.bvir_off;
105
dpd_file2_init(&fIJ, CC_OEI, 0, 0, 0, "fIJ");
106
dpd_file2_mat_init(&fIJ);
107
dpd_file2_mat_rd(&fIJ);
109
dpd_file2_init(&fij, CC_OEI, 0, 2, 2, "fij");
110
dpd_file2_mat_init(&fij);
111
dpd_file2_mat_rd(&fij);
113
dpd_file2_init(&fAB, CC_OEI, 0, 1, 1, "fAB");
114
dpd_file2_mat_init(&fAB);
115
dpd_file2_mat_rd(&fAB);
117
dpd_file2_init(&fab, CC_OEI, 0, 3, 3, "fab");
118
dpd_file2_mat_init(&fab);
119
dpd_file2_mat_rd(&fab);
121
sprintf(lbl, "dIJAB[%d]", irrep);
122
dpd_buf4_init(&D, CC_MISC, irrep, 1, 6, 1, 6, 0, lbl);
123
for(Gij=0; Gij < nirreps; Gij++) {
124
dpd_buf4_mat_irrep_init(&D, Gij);
125
for(ij=0; ij < D.params->rowtot[Gij]; ij++) {
126
i = D.params->roworb[Gij][ij][0];
127
j = D.params->roworb[Gij][ij][1];
128
isym = D.params->psym[i];
129
jsym = D.params->qsym[j];
130
I = i - aocc_off[isym];
131
J = j - aocc_off[jsym];
132
fii = fIJ.matrix[isym][I][I];
133
fjj = fIJ.matrix[jsym][J][J];
135
for(ab=0; ab < D.params->coltot[Gij^irrep]; ab++) {
136
a = D.params->colorb[Gij^irrep][ab][0];
137
b = D.params->colorb[Gij^irrep][ab][1];
138
asym = D.params->rsym[a];
139
bsym = D.params->ssym[b];
140
A = a - avir_off[asym];
141
B = b - avir_off[bsym];
142
faa = fAB.matrix[asym][A][A];
143
fbb = fAB.matrix[bsym][B][B];
145
D.matrix[Gij][ij][ab] = 1.0/(fii + fjj - faa - fbb + root);
148
dpd_buf4_mat_irrep_wrt(&D, Gij);
149
dpd_buf4_mat_irrep_close(&D, Gij);
153
sprintf(lbl, "dijab[%d]", irrep);
154
dpd_buf4_init(&D, CC_MISC, irrep, 11, 16, 11, 16, 0, lbl);
155
for(Gij=0; Gij < nirreps; Gij++) {
156
dpd_buf4_mat_irrep_init(&D, Gij);
157
for(ij=0; ij < D.params->rowtot[Gij]; ij++) {
158
i = D.params->roworb[Gij][ij][0];
159
j = D.params->roworb[Gij][ij][1];
160
isym = D.params->psym[i];
161
jsym = D.params->qsym[j];
162
I = i - bocc_off[isym];
163
J = j - bocc_off[jsym];
164
fii = fij.matrix[isym][I][I];
165
fjj = fij.matrix[jsym][J][J];
167
for(ab=0; ab < D.params->coltot[Gij^irrep]; ab++) {
168
a = D.params->colorb[Gij^irrep][ab][0];
169
b = D.params->colorb[Gij^irrep][ab][1];
170
asym = D.params->rsym[a];
171
bsym = D.params->ssym[b];
172
A = a - bvir_off[asym];
173
B = b - bvir_off[bsym];
174
faa = fab.matrix[asym][A][A];
175
fbb = fab.matrix[bsym][B][B];
177
D.matrix[Gij][ij][ab] = 1.0/(fii + fjj - faa - fbb + root);
180
dpd_buf4_mat_irrep_wrt(&D, Gij);
181
dpd_buf4_mat_irrep_close(&D, Gij);
185
sprintf(lbl, "dIjAb[%d]", irrep);
186
dpd_buf4_init(&D, CC_MISC, irrep, 22, 28, 22, 28, 0, lbl);
187
for(Gij=0; Gij < nirreps; Gij++) {
188
dpd_buf4_mat_irrep_init(&D, Gij);
189
for(ij=0; ij < D.params->rowtot[Gij]; ij++) {
190
i = D.params->roworb[Gij][ij][0];
191
j = D.params->roworb[Gij][ij][1];
192
isym = D.params->psym[i];
193
jsym = D.params->qsym[j];
194
I = i - aocc_off[isym];
195
J = j - bocc_off[jsym];
196
fii = fIJ.matrix[isym][I][I];
197
fjj = fij.matrix[jsym][J][J];
199
for(ab=0; ab < D.params->coltot[Gij^irrep]; ab++) {
200
a = D.params->colorb[Gij^irrep][ab][0];
201
b = D.params->colorb[Gij^irrep][ab][1];
202
asym = D.params->rsym[a];
203
bsym = D.params->ssym[b];
204
A = a - avir_off[asym];
205
B = b - bvir_off[bsym];
206
faa = fAB.matrix[asym][A][A];
207
fbb = fab.matrix[bsym][B][B];
209
D.matrix[Gij][ij][ab] = 1.0/(fii + fjj - faa - fbb + root);
212
dpd_buf4_mat_irrep_wrt(&D, Gij);
213
dpd_buf4_mat_irrep_close(&D, Gij);
217
dpd_file2_mat_close(&fIJ);
218
dpd_file2_mat_close(&fij);
219
dpd_file2_mat_close(&fAB);
220
dpd_file2_mat_close(&fab);
221
dpd_file2_close(&fIJ);
222
dpd_file2_close(&fij);
223
dpd_file2_close(&fAB);
224
dpd_file2_close(&fab);
228
}} // namespace psi::cis