3
\brief Enter brief description of file here
6
#include <libdpd/dpd.h>
8
#include <libciomr/libciomr.h>
9
#include <libiwl/iwl.h>
16
namespace psi { namespace ccdensity {
18
/* SORTI_ROHF(): Place all the components of the ROHF Lagrangian into
19
** a large matrix, I (moinfo.I), which we also symmetrize by computing
20
** Ipq = 1/2 (Ipq + Iqp). This matrix is later written to disk in
21
** dump() for subsequent backtransformation. Note that some of the
22
** components of the Lagrangian computed into the IIJ, Iij, IIA, and
23
** Iia matrices remain non-symmetric (e.g., IIJ neq IJI). I re-used
24
** my sortone.c code here, so don't let some of the variable names
29
int h, nirreps, nmo, nfzv, nfzc, nclsd, nopen;
30
int row, col, i, j, I, J, a, b, A, B, p, q;
31
int *occpi, *virtpi, *occ_off, *vir_off;
32
int *occ_sym, *vir_sym, *openpi;
34
double **O, chksum, value;
42
nirreps = moinfo.nirreps;
43
occpi = moinfo.occpi; virtpi = moinfo.virtpi;
44
occ_off = moinfo.occ_off; vir_off = moinfo.vir_off;
45
occ_sym = moinfo.occ_sym; vir_sym = moinfo.vir_sym;
46
openpi = moinfo.openpi;
47
qt_occ = moinfo.qt_occ; qt_vir = moinfo.qt_vir;
49
O = block_matrix(nmo,nmo);
51
/* Sort alpha components first */
52
dpd_file2_init(&D, CC_OEI, 0, 0, 0, "I(I,J)");
53
dpd_file2_mat_init(&D);
55
for(h=0; h < nirreps; h++) {
56
for(i=0; i < occpi[h]; i++) {
57
I = qt_occ[occ_off[h] + i];
58
for(j=0; j < occpi[h]; j++) {
59
J = qt_occ[occ_off[h] + j];
60
O[I][J] += D.matrix[h][i][j];
64
dpd_file2_mat_close(&D);
67
dpd_file2_init(&D, CC_OEI, 0, 1, 1, "I'AB");
68
dpd_file2_mat_init(&D);
70
for(h=0; h < nirreps; h++) {
71
for(a=0; a < (virtpi[h] - openpi[h]); a++) {
72
A = qt_vir[vir_off[h] + a];
73
for(b=0; b < (virtpi[h] - openpi[h]); b++) {
74
B = qt_vir[vir_off[h] + b];
76
O[A][B] += D.matrix[h][a][b];
80
dpd_file2_mat_close(&D);
83
dpd_file2_init(&D, CC_OEI, 0, 0, 1, "I(I,A)");
84
dpd_file2_mat_init(&D);
86
for(h=0; h < nirreps; h++) {
87
for(i=0; i < occpi[h]; i++) {
88
I = qt_occ[occ_off[h] + i];
89
for(a=0; a < (virtpi[h] - openpi[h]); a++) {
90
A = qt_vir[vir_off[h] + a];
92
O[A][I] += D.matrix[h][i][a];
93
O[I][A] += D.matrix[h][i][a];
97
dpd_file2_mat_close(&D);
100
/* Sort beta components */
101
dpd_file2_init(&D, CC_OEI, 0, 0, 0, "I(i,j)");
102
dpd_file2_mat_init(&D);
103
dpd_file2_mat_rd(&D);
104
for(h=0; h < nirreps; h++) {
105
for(i=0; i < (occpi[h] - openpi[h]); i++) {
106
I = qt_occ[occ_off[h] + i];
107
for(j=0; j < (occpi[h] - openpi[h]); j++) {
108
J = qt_occ[occ_off[h] + j];
109
O[I][J] += D.matrix[h][i][j];
113
dpd_file2_mat_close(&D);
116
dpd_file2_init(&D, CC_OEI, 0, 1, 1, "I'ab");
117
dpd_file2_mat_init(&D);
118
dpd_file2_mat_rd(&D);
119
for(h=0; h < nirreps; h++) {
120
for(a=0; a < virtpi[h]; a++) {
121
A = qt_vir[vir_off[h] + a];
122
for(b=0; b < virtpi[h]; b++) {
123
B = qt_vir[vir_off[h] + b];
125
O[A][B] += D.matrix[h][a][b];
129
dpd_file2_mat_close(&D);
132
dpd_file2_init(&D, CC_OEI, 0, 0, 1, "I(i,a)");
133
dpd_file2_mat_init(&D);
134
dpd_file2_mat_rd(&D);
135
for(h=0; h < nirreps; h++) {
136
for(i=0; i < (occpi[h] - openpi[h]); i++) {
137
I = qt_occ[occ_off[h] + i];
138
for(a=0; a < virtpi[h]; a++) {
139
A = qt_vir[vir_off[h] + a];
141
O[A][I] += D.matrix[h][i][a];
142
O[I][A] += D.matrix[h][i][a];
146
dpd_file2_mat_close(&D);
149
/* Symmetrize the Lagrangian */
150
for(p=0; p < (nmo-nfzv); p++) {
151
for(q=0; q < p; q++) {
152
value = 0.5*(O[p][q] + O[q][p]);
153
O[p][q] = O[q][p] = value;
157
/* Multiply the Lagrangian by -2.0 for the final energy derivative
159
for(p=0; p < (nmo-nfzv); p++) {
160
for(q=0; q < (nmo-nfzv); q++) {
168
}} // namespace psi::ccdensity