4
#include <libciomr/libciomr.h>
5
#include <libdpd/dpd.h>
9
double diagnostic(void)
11
int h, nirreps, Gi, Ga;
12
int i, a, I, A, row, col;
13
int num_elec, num_elec_a, num_elec_b;
15
int *occ_sym, *vir_sym;
18
double t1diag, t1diag_a, t1diag_b;
21
nirreps = moinfo.nirreps;
22
clsdpi = moinfo.clsdpi;
23
uoccpi = moinfo.uoccpi;
24
openpi = moinfo.openpi;
27
occpi = moinfo.occpi; virtpi = moinfo.virtpi;
28
occ_sym = moinfo.occ_sym; vir_sym = moinfo.vir_sym;
31
/* Compute the number of electrons */
32
for(h=0,num_elec_a=0,num_elec_b=0; h < nirreps; h++) {
33
num_elec_a += clsdpi[h] + openpi[h];
34
num_elec_b += clsdpi[h];
36
num_elec = num_elec_a + num_elec_b;
38
if(params.ref == 0) { /** RHF **/
40
dpd_file2_init(&T1A, CC_OEI, 0, 0, 1, "tIA");
41
t1diag = dpd_file2_dot_self(&T1A);
42
dpd_file2_close(&T1A);
45
t1diag = sqrt(t1diag);
48
else if(params.ref == 1) { /** ROHF **/
50
dpd_file2_init(&T1A, CC_OEI, 0, 0, 1, "tIA");
51
dpd_file2_mat_init(&T1A);
52
dpd_file2_mat_rd(&T1A);
53
dpd_file2_init(&T1B, CC_OEI, 0, 0, 1, "tia");
54
dpd_file2_mat_init(&T1B);
55
dpd_file2_mat_rd(&T1B);
58
for(h=0; h < nirreps; h++) {
60
for(i=0; i < (occpi[h] - openpi[h]); i++) {
61
for(a=0; a < (virtpi[h] - openpi[h]); a++) {
63
t1diag += (T1A.matrix[h][i][a] + T1B.matrix[h][i][a]) *
64
(T1A.matrix[h][i][a] + T1B.matrix[h][i][a]);
68
for(i=0; i < (occpi[h] - openpi[h]); i++) {
69
for(a=0; a < openpi[h]; a++) {
72
t1diag += 2 * T1B.matrix[h][i][A] * T1B.matrix[h][i][A];
76
for(i=0; i < openpi[h]; i++) {
78
for(a=0; a < (virtpi[h] - openpi[h]); a++) {
80
t1diag += 2 * T1A.matrix[h][I][a] * T1A.matrix[h][I][a];
86
t1diag = sqrt(t1diag);
89
dpd_file2_mat_close(&T1A);
90
dpd_file2_close(&T1A);
91
dpd_file2_mat_close(&T1B);
92
dpd_file2_close(&T1B);
95
else if(params.ref == 2) { /** UHF **/
97
dpd_file2_init(&T1A, CC_OEI, 0, 0, 1, "tIA");
98
dpd_file2_mat_init(&T1A);
99
dpd_file2_mat_rd(&T1A);
100
dpd_file2_init(&T1B, CC_OEI, 0, 2, 3, "tia");
101
dpd_file2_mat_init(&T1B);
102
dpd_file2_mat_rd(&T1B);
106
for(h=0; h < nirreps; h++) {
108
for(row=0; row < T1A.params->rowtot[h]; row++)
109
for(col=0; col < T1A.params->coltot[h]; col++)
110
t1diag_a += (T1A.matrix[h][row][col] * T1A.matrix[h][row][col]);
112
for(row=0; row < T1B.params->rowtot[h]; row++)
113
for(col=0; col < T1B.params->coltot[h]; col++)
114
t1diag_b += (T1B.matrix[h][row][col] * T1B.matrix[h][row][col]);
118
t1diag = sqrt((t1diag_a + t1diag_b)/(num_elec_a + num_elec_b));
120
dpd_file2_mat_close(&T1A);
121
dpd_file2_mat_close(&T1B);
122
dpd_file2_close(&T1A);
123
dpd_file2_close(&T1B);