1
#include <libdpd/dpd.h>
5
void rhf_sf_sort_I(void);
6
void uhf_sf_sort_I(void);
10
if(params.ref == 0) rhf_sf_sort_I();
11
else if(params.ref == 2) uhf_sf_sort_I();
17
int row, col, i, j, I, J, a, b, A, B, p, q;
18
int *occpi, *virpi, *occ_off, *vir_off;
19
int *occ_sym, *vir_sym;
36
O = block_matrix(nmo,nmo);
38
dpd_file2_init(&W, CC_OEI, 0, 0, 0, "WIJ");
39
dpd_file2_mat_init(&W);
41
for(h=0; h < nirreps; h++) {
42
for(i=0; i < occpi[h]; i++) {
43
I = qt_occ[occ_off[h] + i];
44
for(j=0; j < occpi[h]; j++) {
45
J = qt_occ[occ_off[h] + j];
46
O[I][J] += W.matrix[h][i][j];
50
dpd_file2_mat_close(&W);
53
dpd_file2_init(&W, CC_OEI, 0, 1, 1, "WAB");
54
dpd_file2_mat_init(&W);
56
for(h=0; h < nirreps; h++) {
57
for(a=0; a < virpi[h]; a++) {
58
A = qt_vir[vir_off[h] + a];
59
for(b=0; b < virpi[h]; b++) {
60
B = qt_vir[vir_off[h] + b];
61
O[A][B] += W.matrix[h][a][b];
65
dpd_file2_mat_close(&W);
68
dpd_file2_init(&W, CC_OEI, 0, 1, 0, "WAI");
69
dpd_file2_mat_init(&W);
71
for(h=0; h < nirreps; h++) {
72
for(a=0; a < virpi[h]; a++) {
73
A = qt_vir[vir_off[h] + a];
74
for(i=0; i < occpi[h]; i++) {
75
I = qt_occ[occ_off[h] + i];
76
O[A][I] += W.matrix[h][a][i];
80
dpd_file2_mat_close(&W);
83
dpd_file2_init(&W, CC_OEI, 0, 0, 1, "WIA");
84
dpd_file2_mat_init(&W);
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 < virpi[h]; a++) {
90
A = qt_vir[vir_off[h] + a];
91
O[I][A] += W.matrix[h][i][a];
95
dpd_file2_mat_close(&W);
98
for(p=0; p < nmo; p++) {
99
for(q=0; q < p; q++) {
100
value = 0.5 * (O[p][q] + O[q][p]);
101
O[p][q] = O[q][p] = value;
105
for(p=0; p < nmo; p++) {
106
for(q=0; q < nmo; q++) {
112
fprintf(outfile,"\n\tEnergy Weighted OPDM:\n");
113
print_mat(O,nmo,nmo,outfile);
116
psio_open(PSIF_MO_LAG, PSIO_OPEN_OLD);
117
psio_write_entry(PSIF_MO_LAG, "MO-basis Lagrangian", (char *)O[0],
118
sizeof(double)*nmo*nmo);
119
psio_close(PSIF_MO_LAG, 1);
124
void uhf_sort_W(void)
129
void rhf_sf_sort_I(void)
131
int h, nirreps, nmo, nfzv, nfzc, nclsd, nopen;
132
int row, col, i, j, I, J, a, b, A, B, p, q;
133
int *occpi, *virpi, *occ_off, *vir_off;
134
int *occ_sym, *vir_sym, *openpi;
135
int *qt_occ, *qt_vir;
136
double **O, chksum, value;
141
nirreps = mo.nirreps;
144
occ_off = mo.occ_off;
145
vir_off = mo.vir_off;
146
occ_sym = mo.occ_sym;
147
vir_sym = mo.vir_sym;
151
O = block_matrix(nmo,nmo);
153
/* Sort alpha components first */
154
dpd_file2_init(&D, CC_OEI, 0, 0, 0, "I(I,J)");
155
dpd_file2_mat_init(&D);
156
dpd_file2_mat_rd(&D);
157
for(h=0; h < nirreps; h++) {
158
for(i=0; i < occpi[h]; i++) {
159
I = qt_occ[occ_off[h] + i];
160
for(j=0; j < occpi[h]; j++) {
161
J = qt_occ[occ_off[h] + j];
162
O[I][J] += D.matrix[h][i][j];
166
dpd_file2_mat_close(&D);
169
dpd_file2_init(&D, CC_OEI, 0, 1, 1, "I'AB");
170
dpd_file2_mat_init(&D);
171
dpd_file2_mat_rd(&D);
172
for(h=0; h < nirreps; h++) {
173
for(a=0; a < virpi[h]; a++) {
174
A = qt_vir[vir_off[h] + a];
175
for(b=0; b < virpi[h]; b++) {
176
B = qt_vir[vir_off[h] + b];
177
O[A][B] += D.matrix[h][a][b];
181
dpd_file2_mat_close(&D);
184
dpd_file2_init(&D, CC_OEI, 0, 0, 1, "I(I,A)");
185
dpd_file2_mat_init(&D);
186
dpd_file2_mat_rd(&D);
187
for(h=0; h < nirreps; h++) {
188
for(i=0; i < occpi[h]; i++) {
189
I = qt_occ[occ_off[h] + i];
190
for(a=0; a < virpi[h]; a++) {
191
A = qt_vir[vir_off[h] + a];
192
O[A][I] += D.matrix[h][i][a];
193
O[I][A] += D.matrix[h][i][a];
197
dpd_file2_mat_close(&D);
200
/* Sort beta components */
201
dpd_file2_init(&D, CC_OEI, 0, 0, 0, "I(i,j)");
202
dpd_file2_mat_init(&D);
203
dpd_file2_mat_rd(&D);
204
for(h=0; h < nirreps; h++) {
205
for(i=0; i < occpi[h]; i++) {
206
I = qt_occ[occ_off[h] + i];
207
for(j=0; j < occpi[h]; j++) {
208
J = qt_occ[occ_off[h] + j];
209
O[I][J] += D.matrix[h][i][j];
213
dpd_file2_mat_close(&D);
216
dpd_file2_init(&D, CC_OEI, 0, 1, 1, "I'ab");
217
dpd_file2_mat_init(&D);
218
dpd_file2_mat_rd(&D);
219
for(h=0; h < nirreps; h++) {
220
for(a=0; a < virpi[h]; a++) {
221
A = qt_vir[vir_off[h] + a];
222
for(b=0; b < virpi[h]; b++) {
223
B = qt_vir[vir_off[h] + b];
224
O[A][B] += D.matrix[h][a][b];
228
dpd_file2_mat_close(&D);
231
dpd_file2_init(&D, CC_OEI, 0, 0, 1, "I(i,a)");
232
dpd_file2_mat_init(&D);
233
dpd_file2_mat_rd(&D);
234
for(h=0; h < nirreps; h++) {
235
for(i=0; i < occpi[h]; i++) {
236
I = qt_occ[occ_off[h] + i];
237
for(a=0; a < virpi[h]; a++) {
238
A = qt_vir[vir_off[h] + a];
239
O[A][I] += D.matrix[h][i][a];
240
O[I][A] += D.matrix[h][i][a];
244
dpd_file2_mat_close(&D);
247
/* Symmetrize the Lagrangian */
248
for(p=0; p < nmo; p++) {
249
for(q=0; q < p; q++) {
250
value = 0.5*(O[p][q] + O[q][p]);
251
O[p][q] = O[q][p] = value;
255
/* Multiply the Lagrangian by -2.0 for the final energy derivative
257
for(p=0; p < nmo; p++)
258
for(q=0; q < nmo; q++)
264
void uhf_sf_sort_I(void)