4
#include <libciomr/libciomr.h>
5
#include <libdpd/dpd.h>
10
/* Computes a modified D1 diagnostic developed by T.J. Lee, but not yet
14
double d1diag_t1_rhf(void);
16
static double new_d1diag_t1_rohf(void)
19
int nclsd, nuocc, nopen;
20
double **T1_hp, **T1_hx, **T1_xp, **T1_sq;
22
double max_hp=0.0, max_xp=0.0, max_hx=0.0, max;
25
nirreps = moinfo.nirreps;
27
dpd_file2_init(&T1_a, CC_OEI, 0, 0, 1, "tIA");
28
dpd_file2_mat_init(&T1_a);
29
dpd_file2_mat_rd(&T1_a);
31
dpd_file2_init(&T1_b, CC_OEI, 0, 0, 1, "tia");
32
dpd_file2_mat_init(&T1_b);
33
dpd_file2_mat_rd(&T1_b);
35
for(h=0; h < nirreps; h++) {
36
nclsd = moinfo.clsdpi[h];
37
nuocc = moinfo.uoccpi[h];
38
nopen = moinfo.openpi[h];
41
T1_hp = block_matrix(nclsd, nuocc);
42
for (i=0; i < nclsd; i++)
43
for (j=0; j < nuocc; j++)
44
T1_hp[i][j] = (T1_a.matrix[h][i][j] + T1_b.matrix[h][i][j])/2.;
46
T1_sq = block_matrix(nclsd, nclsd);
47
C_DGEMM('n','t',nclsd,nclsd,nuocc,1.0,&(T1_hp[0][0]),nuocc,
48
&(T1_hp[0][0]),nuocc,0.0,&(T1_sq[0][0]),nclsd);
50
E = init_array(nclsd);
51
C = block_matrix(nclsd, nclsd);
52
sq_rsp(nclsd, nclsd, T1_sq, E, 0, C, 1e-12);
53
for(i=0; i < nclsd; i++) if(E[i] > max_hp) max_hp = E[i];
61
T1_hx = block_matrix(nclsd, nopen);
62
for (i=0; i < nclsd; i++)
63
for (j=0; j < nopen; j++)
64
T1_hx[i][j] = T1_b.matrix[h][i][nuocc+j]/sqrt(2.);
66
T1_sq = block_matrix(nclsd, nclsd);
67
C_DGEMM('n','t',nclsd,nclsd,nopen,1.0,&(T1_hx[0][0]),nopen,
68
&(T1_hx[0][0]),nopen,0.0,&(T1_sq[0][0]),nclsd);
70
E = init_array(nclsd);
71
C = block_matrix(nclsd, nclsd);
72
sq_rsp(nclsd, nclsd, T1_sq, E, 0, C, 1e-12);
73
for(i=0; i < nclsd; i++) if(E[i] > max_hx) max_hx = E[i];
81
T1_xp = block_matrix(nopen, nuocc);
82
for (i=0; i < nopen; i++)
83
for (j=0; j < nuocc; j++)
84
T1_xp[i][j] = T1_a.matrix[h][nclsd+i][j]/sqrt(2.);
86
T1_sq = block_matrix(nopen, nopen);
87
C_DGEMM('n','t',nopen,nopen,nuocc,1.0,&(T1_xp[0][0]),nuocc,
88
&(T1_xp[0][0]),nuocc,0.0,&(T1_sq[0][0]),nopen);
90
E = init_array(nopen);
91
C = block_matrix(nopen, nopen);
92
sq_rsp(nopen, nopen, T1_sq, E, 0, C, 1e-12);
93
for(i=0; i < nopen; i++) if(E[i] > max_xp) max_xp = E[i];
101
dpd_file2_mat_close(&T1_a);
102
dpd_file2_close(&T1_a);
104
dpd_file2_mat_close(&T1_b);
105
dpd_file2_close(&T1_b);
107
max_hp = sqrt(max_hp);
108
max_hx = sqrt(max_hx);
109
max_xp = sqrt(max_xp);
112
fprintf(outfile, "ND1: hp=%8.6f hx=%8.6f xp=%8.6f\n", max_hp, max_hx, max_xp);
116
if (max_hx > max) max = max_hx;
117
if (max_xp > max) max = max_xp;
122
double new_d1diag(void)
126
if(params.ref == 0) { /** RHF **/
127
norm = d1diag_t1_rhf();
129
else if (params.ref == 1) { /** ROHF **/
130
norm = new_d1diag_t1_rohf();