3
\brief Enter brief description of file here
13
namespace psi { namespace cceom {
15
/* This function orthogonalizes the r residual vector with the set of
16
numCs C vectors and adds the new vector to the C list if its norm is greater
17
than params.residual_tol */
19
extern double norm_C(dpdfile2 *CME, dpdfile2 *Cme,
20
dpdbuf4 *CMNEF, dpdbuf4 *Cmnef, dpdbuf4 *CMnEf);
22
extern void scm_C(dpdfile2 *CME, dpdfile2 *Cme, dpdbuf4 *CMNEF,
23
dpdbuf4 *Cmnef, dpdbuf4 *CMnEf, double a);
25
/* use for ROHF and UHF */
26
void schmidt_add(dpdfile2 *RIA, dpdfile2 *Ria,
27
dpdbuf4 *RIJAB, dpdbuf4 *Rijab, dpdbuf4 *RIjAb, int *numCs, int irrep)
32
dpdfile2 Cme, CME, Cme2, CME2;
33
dpdbuf4 CMNEF, Cmnef, CMnEf, CMNEF2, Cmnef2, CMnEf2;
35
char CME_lbl[32], Cme_lbl[32], CMNEF_lbl[32], Cmnef_lbl[32], CMnEf_lbl[32];
37
for (i=0; i<*numCs; i++) {
38
sprintf(CME_lbl, "%s %d", "CME", i);
39
sprintf(Cme_lbl, "%s %d", "Cme", i);
40
sprintf(CMNEF_lbl, "%s %d", "CMNEF", i);
41
sprintf(Cmnef_lbl, "%s %d", "Cmnef", i);
42
sprintf(CMnEf_lbl, "%s %d", "CMnEf", i);
44
dpd_file2_init(&CME, EOM_CME, irrep, 0, 1, CME_lbl);
45
dpd_buf4_init(&CMNEF, EOM_CMNEF, irrep, 2, 7, 2, 7, 0, CMNEF_lbl);
46
if (params.eom_ref == 1) {
47
dpd_file2_init(&Cme, EOM_Cme, irrep, 0, 1, Cme_lbl);
48
dpd_buf4_init(&Cmnef, EOM_Cmnef, irrep, 2, 7, 2, 7, 0, Cmnef_lbl);
49
dpd_buf4_init(&CMnEf, EOM_CMnEf, irrep, 0, 5, 0, 5, 0, CMnEf_lbl);
51
else if (params.eom_ref == 2) {
52
dpd_file2_init(&Cme, EOM_Cme, irrep, 2, 3, Cme_lbl);
53
dpd_buf4_init(&Cmnef, EOM_Cmnef, irrep, 12, 17, 12, 17, 0, Cmnef_lbl);
54
dpd_buf4_init(&CMnEf, EOM_CMnEf, irrep, 22, 28, 22, 28, 0, CMnEf_lbl);
57
dotval = dpd_file2_dot(RIA, &CME);
58
dotval += dpd_file2_dot(Ria, &Cme);
59
//fprintf(outfile, "OE Dotval for vector %d = %20.14f\n", i, dotval);
60
dotval += dpd_buf4_dot(RIJAB, &CMNEF);
61
dotval += dpd_buf4_dot(Rijab, &Cmnef);
62
dotval += dpd_buf4_dot(RIjAb, &CMnEf);
64
//fprintf(outfile, "Dotval for vector %d = %20.14f\n", i, dotval);
66
dpd_file2_axpy(&CME, RIA, -1.0*dotval, 0);
67
dpd_file2_axpy(&Cme, Ria, -1.0*dotval, 0);
68
dpd_buf4_axpy(&CMNEF, RIJAB, -1.0*dotval);
69
dpd_buf4_axpy(&Cmnef, Rijab, -1.0*dotval);
70
dpd_buf4_axpy(&CMnEf, RIjAb, -1.0*dotval);
72
dpd_file2_close(&CME);
73
dpd_file2_close(&Cme);
74
dpd_buf4_close(&CMNEF);
75
dpd_buf4_close(&Cmnef);
76
dpd_buf4_close(&CMnEf);
79
norm = norm_C(RIA, Ria, RIJAB, Rijab, RIjAb);
80
//fprintf(outfile, "Norm of residual (TDC) = %20.14f\n", norm);
82
if (norm < eom_params.schmidt_add_residual_tol) {
86
scm_C(RIA, Ria, RIJAB, Rijab, RIjAb, 1.0/norm);
87
sprintf(CME_lbl, "%s %d", "CME", *numCs);
88
sprintf(Cme_lbl, "%s %d", "Cme", *numCs);
89
sprintf(CMNEF_lbl, "%s %d", "CMNEF", *numCs);
90
sprintf(Cmnef_lbl, "%s %d", "Cmnef", *numCs);
91
sprintf(CMnEf_lbl, "%s %d", "CMnEf", *numCs);
93
dpd_file2_copy(RIA, EOM_CME, CME_lbl);
94
dpd_file2_copy(Ria, EOM_Cme, Cme_lbl);
95
dpd_buf4_copy(RIJAB, EOM_CMNEF, CMNEF_lbl);
96
dpd_buf4_copy(Rijab, EOM_Cmnef, Cmnef_lbl);
97
dpd_buf4_copy(RIjAb, EOM_CMnEf, CMnEf_lbl);
104
void schmidt_add_RHF(dpdfile2 *RIA, dpdbuf4 *RIjAb, int *numCs, int irrep)
106
double dotval, norm, R0, C0;
109
dpdbuf4 CMnEf, CAB1, CAB2;
112
char CME_lbl[32], Cme_lbl[32], CMNEF_lbl[32], Cmnef_lbl[32], CMnEf_lbl[32], C0_lbl[32];
114
if (params.full_matrix) psio_read_entry(EOM_R, "R0", (char *) &R0, sizeof(double));
116
for (i=0; i<*numCs; i++) {
117
/* Spin-adapt the residual */
118
dpd_buf4_copy(RIjAb, EOM_TMP, "RIjAb");
119
dpd_buf4_sort(RIjAb, EOM_TMP, pqsr, 0, 5, "RIjbA");
121
dpd_buf4_init(&R2a, EOM_TMP, irrep, 0, 5, 0, 5, 0, "RIjAb");
122
dpd_buf4_init(&R2b, EOM_TMP, irrep, 0, 5, 0, 5, 0, "RIjbA");
123
dpd_buf4_scm(&R2a, 2.0);
124
dpd_buf4_axpy(&R2b, &R2a, -1.0);
125
dpd_buf4_close(&R2b);
127
sprintf(CME_lbl, "%s %d", "CME", i);
128
sprintf(CMnEf_lbl, "%s %d", "CMnEf", i);
129
dpd_file2_init(&CME, EOM_CME, irrep, 0, 1, CME_lbl);
130
dpd_buf4_init(&CMnEf, EOM_CMnEf, irrep, 0, 5, 0, 5, 0, CMnEf_lbl);
131
dotval = 2.0 * dpd_file2_dot(RIA, &CME);
132
//fprintf(outfile, "OE Dotval for vector %d = %20.14f\n", i, dotval);
133
dotval += dpd_buf4_dot(&R2a, &CMnEf);
134
dpd_buf4_close(&R2a);
135
if (params.full_matrix) {
136
sprintf(C0_lbl, "%s %d", "C0", i);
137
psio_read_entry(EOM_CME, C0_lbl, (char *) &C0, sizeof(double));
141
//fprintf(outfile, "Dotval for vector %d = %20.14f\n", i, dotval);
142
R0 = R0 - 1.0 * dotval * C0;
143
dpd_file2_axpy(&CME, RIA, -1.0*dotval, 0);
144
dpd_buf4_axpy(&CMnEf, RIjAb, -1.0*dotval);
145
dpd_file2_close(&CME);
146
dpd_buf4_close(&CMnEf);
149
dpd_buf4_sort(RIjAb, EOM_TMP, pqsr, 0, 5, "RIjbA");
150
dpd_buf4_init(&R2b, EOM_TMP, irrep, 0, 5, 0, 5, 0, "RIjbA");
152
/* norm = norm_C_rhf(RIA, RIjAb, &R2b); */
153
norm = 2.0 * dpd_file2_dot_self(RIA);
154
norm += 2.0 * dpd_buf4_dot_self(RIjAb);
155
norm -= dpd_buf4_dot(RIjAb, &R2b);
156
if (params.full_matrix)
160
dpd_buf4_close(&R2b);
162
//fprintf(outfile, "Norm of residual (TDC) = %20.14f\n", norm);
164
if (norm < eom_params.schmidt_add_residual_tol) {
169
if (params.full_matrix) R0 *= 1.0/norm;
170
dpd_file2_scm(RIA, 1.0/norm);
171
dpd_buf4_scm(RIjAb, 1.0/norm);
174
dpd_buf4_sort(RIjAb, EOM_TMP, pqsr, 0, 5, "RIjbA");
175
dpd_buf4_init(&R2b, EOM_TMP, irrep, 0, 5, 0, 5, 0, "RIjbA");
176
norm = 2.0 * dpd_file2_dot_self(RIA);
177
norm += 2.0 * dpd_buf4_dot_self(RIjAb);
178
norm -= dpd_buf4_dot(RIjAb, &R2b);
179
if (params.full_matrix) norm += R0 * R0;
181
fprintf(outfile,"Norm of final new C in schmidt_add(): %20.15lf\n", norm);
182
dpd_buf4_close(&R2b);
185
sprintf(CME_lbl, "%s %d", "CME", *numCs);
186
sprintf(CMnEf_lbl, "%s %d", "CMnEf", *numCs);
188
dpd_file2_copy(RIA, EOM_CME, CME_lbl);
189
dpd_buf4_copy(RIjAb, EOM_CMnEf, CMnEf_lbl);
191
/* Generate AA and BB C2 vectors from AB vector */
192
/* C(IJ,AB) = C(ij,ab) = C(Ij,Ab) - C(Ij,bA) */
193
dpd_buf4_copy(RIjAb, EOM_TMP, "CMnEf");
194
dpd_buf4_sort(RIjAb, EOM_TMP, pqsr, 0, 5, "CMnfE");
196
dpd_buf4_init(&CAB1, EOM_TMP, irrep, 0, 5, 0, 5, 0, "CMnEf");
197
dpd_buf4_init(&CAB2, EOM_TMP, irrep, 0, 5, 0, 5, 0, "CMnfE");
198
dpd_buf4_axpy(&CAB2, &CAB1, -1.0);
199
dpd_buf4_close(&CAB2);
200
dpd_buf4_close(&CAB1);
202
if (params.full_matrix) {
203
sprintf(C0_lbl, "%s %d", "C0", *numCs);
204
psio_write_entry(EOM_CME, C0_lbl, (char *) &R0, sizeof(double));
211
}} // namespace psi::cceom