2
#include <libiwl/iwl.h>
3
#include <libdpd/dpd.h>
8
/* DEANTI_ROHF(): Convert the ROHF two-particle density from Dirac to
9
** Mulliken ordering. The original, Fock-adjusted density (see the
10
** comments in fock.c) corresponds to a two-electron energy (or energy
11
** derivative) expression of the form:
13
** E(TWO) = 1/4 sum_pqrs Gpqrs <pq||rs>
15
** However, the derivative two-electron integrals are produced in
16
** Mulliken-ordered, symmetric form rather than Dirac-ordered
17
** antisymmetric form. This code alters the two-particle density
18
** matrix ordering for the energy expression of the form:
20
** E(TWO) = 1/2 sum_pqrs Gpqrs <pq|rs>
22
** The final conversion to Mulliken ordering is taken care of in
25
** The second equation above may be derived via
27
** E(TWO) = 1/4 sum_pqrs Gpqrs (<pq|rs> - <pq|sr>)
28
** = 1/4 sum_pqrs Gpqrs <pq|rs> - 1/4 sum_pqrs Gpqrs <pq|sr>
29
** = 1/4 sum_pqrs Gpqrs <pq|rs> - 1/4 sum_pqrs Gpqsr <pq|rs>
30
** = 1/4 sum_pqrs (Gpqrs - Gpqsr) <pq|rs>
31
** = 1/2 sum_pqrs Gpqrs <pq|rs>
33
** Equations for the individual components are given explicitly in
37
void deanti_ROHF(struct RHO_Params rho_params)
41
double one_energy = 0.0, two_energy = 0.0, total_two_energy = 0.0;
42
dpdbuf4 A, B, C, DInts, E, FInts;
45
fprintf(outfile, "\n\tEnergies re-computed from Mulliken density:\n");
46
fprintf(outfile, "\t-------------------------------------------\n");
48
dpd_file2_init(&D, CC_OEI, 0, 0, 0, rho_params.DIJ_lbl);
49
dpd_file2_init(&F, CC_OEI, 0, 0, 0, "h(i,j)");
50
one_energy += dpd_file2_dot(&D, &F);
54
dpd_file2_init(&D, CC_OEI, 0, 0, 0, rho_params.Dij_lbl);
55
dpd_file2_init(&F, CC_OEI, 0, 0, 0, "h(i,j)");
56
one_energy += dpd_file2_dot(&D, &F);
60
dpd_file2_init(&D, CC_OEI, 0, 1, 1, rho_params.DAB_lbl);
61
dpd_file2_init(&F, CC_OEI, 0, 1, 1, "h(a,b)");
62
one_energy += dpd_file2_dot(&D, &F);
66
dpd_file2_init(&D, CC_OEI, 0, 1, 1, rho_params.Dab_lbl);
67
dpd_file2_init(&F, CC_OEI, 0, 1, 1, "h(a,b)");
68
one_energy += dpd_file2_dot(&D, &F);
72
dpd_file2_init(&D, CC_OEI, 0, 0, 1, rho_params.DIA_lbl);
73
dpd_file2_init(&F, CC_OEI, 0, 0, 1, "h(i,a)");
74
one_energy += dpd_file2_dot(&D, &F);
78
dpd_file2_init(&D, CC_OEI, 0, 0, 1, rho_params.Dia_lbl);
79
dpd_file2_init(&F, CC_OEI, 0, 0, 1, "h(i,a)");
80
one_energy += dpd_file2_dot(&D, &F);
84
dpd_file2_init(&D, CC_OEI, 0, 0, 1, rho_params.DAI_lbl);
85
dpd_file2_init(&F, CC_OEI, 0, 0, 1, "h(i,a)");
86
one_energy += dpd_file2_dot(&D, &F);
90
dpd_file2_init(&D, CC_OEI, 0, 0, 1, rho_params.Dai_lbl);
91
dpd_file2_init(&F, CC_OEI, 0, 0, 1, "h(i,a)");
92
one_energy += dpd_file2_dot(&D, &F);
96
fprintf(outfile, "\tOne-electron energy = %20.15f\n", one_energy);
100
/* G(Ij,Kl) <-- 1/2 G(IJ,KL) + 1/2 G(ij,kl) + 1/2 G(Ij,Kl) + 1/2 G(iJ,kL) */
101
dpd_buf4_init(&G1, CC_GAMMA, 0, 0, 0, 0, 0, 0, "GIjKl");
102
dpd_buf4_sort(&G1, CC_TMP0, qprs, 0, 0, "GjIKl");
103
dpd_buf4_init(&G2, CC_TMP0, 0, 0, 0, 0, 0, 0, "GjIKl");
104
dpd_buf4_sort(&G2, CC_TMP0, pqsr, 0, 0, "GiJkL");
106
dpd_buf4_init(&G2, CC_TMP0, 0, 0, 0, 0, 0, 0, "GiJkL");
107
dpd_buf4_axpy(&G2, &G1, 1.0);
109
dpd_buf4_init(&G2, CC_GAMMA, 0, 0, 0, 2, 2, 0, "GIJKL");
110
dpd_buf4_axpy(&G2, &G1, 1.0);
112
dpd_buf4_init(&G2, CC_GAMMA, 0, 0, 0, 2, 2, 0, "Gijkl");
113
dpd_buf4_axpy(&G2, &G1, 1.0);
115
dpd_buf4_scm(&G1, 0.5);
117
if(!params.aobasis) {
118
dpd_buf4_init(&A, CC_AINTS, 0, 0, 0, 0, 0, 0, "A <ij|kl>");
119
two_energy = dpd_buf4_dot(&A, &G1);
121
fprintf(outfile, "\tIJKL energy = %20.15f\n", two_energy);
122
total_two_energy += two_energy;
127
/* G(IJ,KA) <-- G(IJ,KA) + G(ij,ka) + G(Ij,Ka) + G(iJ,kA) */
128
dpd_buf4_init(&G1, CC_GAMMA, 0, 0, 10, 0, 10, 0, "GIjKa");
129
dpd_buf4_init(&G2, CC_GAMMA, 0, 0, 10, 0, 10, 0, "GiJkA");
130
dpd_buf4_axpy(&G2, &G1, 1.0);
132
dpd_buf4_init(&G2, CC_GAMMA, 0, 0, 10, 2, 10, 0, "GIJKA");
133
dpd_buf4_axpy(&G2, &G1, 1.0);
135
dpd_buf4_init(&G2, CC_GAMMA, 0, 0, 10, 2, 10, 0, "Gijka");
136
dpd_buf4_axpy(&G2, &G1, 1.0);
139
if(!params.aobasis) {
140
dpd_buf4_init(&E, CC_EINTS, 0, 0, 10, 0, 10, 0, "E <ij|ka>");
141
two_energy = dpd_buf4_dot(&E, &G1);
143
fprintf(outfile, "\tIJKA energy = %20.15f\n", 2*two_energy);
144
total_two_energy += 2*two_energy;
149
/* G(Ij,Ab) <-- G(Ij,Ab) - G(Ib,jA) */
150
dpd_buf4_init(&G1, CC_GAMMA, 0, 0, 5, 0, 5, 0, "GIjAb");
151
dpd_buf4_init(&G2, CC_GAMMA, 0, 10, 10, 10, 10, 0, "GIbjA");
152
dpd_buf4_sort(&G2, CC_TMP0, prsq, 0, 5, "GIbjA (Ij,Ab)");
154
dpd_buf4_init(&G2, CC_TMP0, 0, 0, 5, 0, 5, 0, "GIbjA (Ij,Ab)");
155
dpd_buf4_axpy(&G2, &G1, -1.0);
159
/* G(IJ,AB) <-- G(IJ,AB) - G(IB,JA) */
160
dpd_buf4_init(&G1, CC_GAMMA, 0, 0, 5, 2, 7, 0, "GIJAB");
161
dpd_buf4_copy(&G1, CC_TMP0, "G(IJ,AB)");
163
dpd_buf4_init(&G1, CC_TMP0, 0, 0, 5, 0, 5, 0, "G(IJ,AB)");
164
dpd_buf4_init(&G2, CC_GAMMA, 0, 10, 10, 10, 10, 0, "GIBJA");
165
dpd_buf4_sort(&G2, CC_TMP0, prsq, 0, 5, "GIBJA (IJ,AB)");
167
dpd_buf4_init(&G2, CC_TMP0, 0, 0, 5, 0, 5, 0, "GIBJA (IJ,AB)");
168
dpd_buf4_axpy(&G2, &G1, -1.0);
172
/* G(ij,ab) <-- G(ij,ab) - G(ib,ja) */
173
dpd_buf4_init(&G1, CC_GAMMA, 0, 0, 5, 2, 7, 0, "Gijab");
174
dpd_buf4_copy(&G1, CC_TMP0, "G(ij,ab)");
176
dpd_buf4_init(&G1, CC_TMP0, 0, 0, 5, 0, 5, 0, "G(ij,ab)");
177
dpd_buf4_init(&G2, CC_GAMMA, 0, 10, 10, 10, 10, 0, "Gibja");
178
dpd_buf4_sort(&G2, CC_TMP0, prsq, 0, 5, "Gibja (ij,ab)");
180
dpd_buf4_init(&G2, CC_TMP0, 0, 0, 5, 0, 5, 0, "Gibja (ij,ab)");
181
dpd_buf4_axpy(&G2, &G1, -1.0);
185
/* G(IJ,AB) <-- G(IJ,AB) + G(ij,ab) + G(Ij,Ab) + G(iJ,aB) */
186
dpd_buf4_init(&G1, CC_GAMMA, 0, 0, 5, 0, 5, 0, "GIjAb");
187
dpd_buf4_sort(&G1, CC_TMP0, qpsr, 0, 5, "G(jI,bA)");
188
dpd_buf4_init(&G2, CC_TMP0, 0, 0, 5, 0, 5, 0, "G(jI,bA)");
189
dpd_buf4_axpy(&G2, &G1, 1.0);
191
dpd_buf4_init(&G2, CC_TMP0, 0, 0, 5, 0, 5, 0, "G(IJ,AB)");
192
dpd_buf4_axpy(&G2, &G1, 1.0);
194
dpd_buf4_init(&G2, CC_TMP0, 0, 0, 5, 0, 5, 0, "G(ij,ab)");
195
dpd_buf4_axpy(&G2, &G1, 1.0);
198
if(!params.aobasis) {
199
dpd_buf4_init(&DInts, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
200
two_energy = dpd_buf4_dot(&DInts, &G1);
201
dpd_buf4_close(&DInts);
202
fprintf(outfile, "\tIJAB energy = %20.15f\n", two_energy);
203
total_two_energy += two_energy;
208
/* G(IB,JA) <-- G(IB,JA) + G(ib,ja) + G(Ib,Ja) + G(iB,jA) */
209
dpd_buf4_init(&G1, CC_GAMMA, 0, 10, 10, 10, 10, 0, "GIBJA");
210
dpd_buf4_init(&G2, CC_GAMMA, 0, 10, 10, 10, 10, 0, "Gibja");
211
dpd_buf4_axpy(&G2, &G1, 1.0);
213
dpd_buf4_init(&G2, CC_GAMMA, 0, 10, 10, 10, 10, 0, "GIbJa");
214
dpd_buf4_axpy(&G2, &G1, 1.0);
216
dpd_buf4_init(&G2, CC_GAMMA, 0, 10, 10, 10, 10, 0, "GiBjA");
217
dpd_buf4_axpy(&G2, &G1, 1.0);
220
if(!params.aobasis) {
221
dpd_buf4_init(&C, CC_CINTS, 0, 10, 10, 10, 10, 0, "C <ia|jb>");
222
two_energy = dpd_buf4_dot(&C, &G1);
224
fprintf(outfile, "\tIBJA energy = %20.15f\n", two_energy);
225
total_two_energy += two_energy;
230
/* G(CI,AB) <-- G(CI,AB) + G(ci,ab) + G(Ci,Ab) + G(cI,aB) */
231
dpd_buf4_init(&G1, CC_GAMMA, 0, 11, 5, 11, 5, 0, "GCiAb");
232
dpd_buf4_init(&G2, CC_GAMMA, 0, 11, 5, 11, 5, 0, "GcIaB");
233
dpd_buf4_axpy(&G2, &G1, 1.0);
235
dpd_buf4_init(&G2, CC_GAMMA, 0, 11, 5, 11, 7, 0, "GCIAB");
236
dpd_buf4_axpy(&G2, &G1, 1.0);
238
dpd_buf4_init(&G2, CC_GAMMA, 0, 11, 5, 11, 7, 0, "Gciab");
239
dpd_buf4_axpy(&G2, &G1, 1.0);
242
if(!params.aobasis) {
243
dpd_buf4_init(&FInts, CC_FINTS, 0, 10, 5, 10, 5, 0, "F <ia|bc>");
244
dpd_buf4_sort(&FInts, CC_TMP0, qprs, 11, 5, "F(CI,BA)");
245
dpd_buf4_close(&FInts);
246
dpd_buf4_init(&FInts, CC_TMP0, 0, 11, 5, 11, 5, 0, "F(CI,BA)");
247
dpd_buf4_sort(&FInts, CC_TMP1, pqsr, 11, 5, "F(CI,AB)");
248
dpd_buf4_close(&FInts);
249
dpd_buf4_init(&FInts, CC_TMP1, 0, 11, 5, 11, 5, 0, "F(CI,AB)");
250
two_energy = dpd_buf4_dot(&FInts, &G1);
251
dpd_buf4_close(&FInts);
252
fprintf(outfile, "\tCIAB energy = %20.15f\n", 2*two_energy);
253
total_two_energy += 2*two_energy;
258
/* G(Ab,Cd) <-- 1/2 G(AB,CD) + 1/2 G(ab,cd) + G(Ab,Cd) */
259
dpd_buf4_init(&G1, CC_GAMMA, 0, 5, 5, 5, 5, 0, "GAbCd");
260
dpd_buf4_sort(&G1, CC_TMP0, qprs, 5, 5, "G(bA,Cd)");
261
dpd_buf4_init(&G2, CC_TMP0, 0, 5, 5, 5, 5, 0, "G(bA,Cd)");
262
dpd_buf4_sort(&G2, CC_TMP1, pqsr, 5, 5, "G(aB,cD)");
264
dpd_buf4_init(&G2, CC_TMP1, 0, 5, 5, 5, 5, 0, "G(aB,cD)");
265
dpd_buf4_axpy(&G2, &G1, 1.0);
267
dpd_buf4_init(&G2, CC_GAMMA, 0, 5, 5, 7, 7, 0, "GABCD");
268
dpd_buf4_axpy(&G2, &G1, 1.0);
270
dpd_buf4_init(&G2, CC_GAMMA, 0, 5, 5, 7, 7, 0, "Gabcd");
271
dpd_buf4_axpy(&G2, &G1, 1.0);
273
dpd_buf4_scm(&G1, 0.5);
275
if(!params.aobasis) {
276
dpd_buf4_init(&B, CC_BINTS, 0, 5, 5, 5, 5, 0, "B <ab|cd>");
277
two_energy = dpd_buf4_dot(&B, &G1);
279
fprintf(outfile, "\tABCD energy = %20.15f\n", two_energy);
280
total_two_energy += two_energy;
285
if(!params.aobasis) {
286
fprintf(outfile, "\tTotal two-electron energy = %20.15f\n", total_two_energy);
288
fprintf(outfile, "\tCCSD correlation energy = %20.15f\n",
289
one_energy + total_two_energy);
290
fprintf(outfile, "\tTotal CCSD energy = %20.15f\n",
291
one_energy + total_two_energy + moinfo.eref);
294
fprintf(outfile, "\tTotal EOM CCSD correlation energy = %20.15f\n",
295
one_energy + total_two_energy);
296
fprintf(outfile, "\tCCSD correlation + EOM excitation energy = %20.15f\n",
297
moinfo.ecc + params.cceom_energy);
298
fprintf(outfile, "\tTotal EOM CCSD energy = %20.15f\n",
299
one_energy + total_two_energy + moinfo.eref);