2
#include <libdpd/dpd.h>
6
/* WijmbL2(): Computes the contributions of the Wmnie HBAR matrix
7
** elements to the Lambda double de-excitation amplitude equations.
8
** These contributions are given in spin orbitals as:
10
** L_ij^ab = - P(ab) L_m^a Wijmb
12
** where Wijmb = Wmnie is defined as:
14
** Wmnie = <mn||ie> + t_i^f <mn||fe>
16
** [cf. Gauss and Stanton, JCP 103, 3561-3577 (1995)]
18
** All four spin cases are stored in (mn,ei) ordering (see Wmnie.c in
19
** the CCHBAR code). The three L_ij^ab spin cases are computed as:
21
** L(IJ,AB) <-- - L(M,A) W(IJ,BM) + L(M,B) W(IJ,AM) (one unique
23
** L(ij,ab) <-- - L(m,a) W(ij,bm) + L(m,b) W(ij,am) (one unique
25
** L(Ij,Ab) <-- - L(M,A) W(Ij,bM) - L(m,b) W(jI,Am)
30
void WijmbL2(int L_irr)
33
dpdbuf4 L2, newLijab, newLIJAB, newLIjAb;
34
dpdbuf4 W, WMNIE, Wmnie, WMnIe, WmNiE;
35
dpdbuf4 X1, X2, Z, Z1, Z2;
37
/* RHS += -P(ab) Lma * Wijmb */
38
if(params.ref == 0 || params.ref == 1) { /** RHF/ROHF **/
40
dpd_file2_init(&LIA, CC_OEI, L_irr, 0, 1, "LIA");
41
dpd_file2_init(&Lia, CC_OEI, L_irr, 0, 1, "Lia");
43
dpd_buf4_init(&WMNIE, CC_HBAR, 0, 2, 11, 2, 11, 0, "WMNIE");
44
dpd_buf4_init(&X1, CC_TMP1, L_irr, 2, 5, 2, 5, 0, "X(2,5) 1");
45
dpd_contract424(&WMNIE, &LIA, &X1, 3, 0, 0, 1.0, 0.0);
46
dpd_buf4_close(&WMNIE);
47
dpd_buf4_sort(&X1, CC_TMP1, pqsr, 2, 5, "X(2,5) 2");
48
dpd_buf4_init(&X2, CC_TMP1, L_irr, 2, 5, 2, 5, 0, "X(2,5) 2");
49
dpd_buf4_axpy(&X2, &X1, -1.0);
51
dpd_buf4_init(&newLIJAB, CC_LAMBDA, L_irr, 2, 5, 2, 7, 0, "New LIJAB");
52
dpd_buf4_axpy(&X1, &newLIJAB, 1.0);
53
dpd_buf4_close(&newLIJAB);
56
dpd_buf4_init(&Wmnie, CC_HBAR, 0, 2, 11, 2, 11, 0, "Wmnie");
57
dpd_buf4_init(&X1, CC_TMP1, L_irr, 2, 5, 2, 5, 0, "X(2,5) 1");
58
dpd_contract424(&Wmnie, &Lia, &X1, 3, 0, 0, 1.0, 0.0);
59
dpd_buf4_close(&Wmnie);
60
dpd_buf4_sort(&X1, CC_TMP1, pqsr, 2, 5, "X(2,5) 2");
61
dpd_buf4_init(&X2, CC_TMP1, L_irr, 2, 5, 2, 5, 0, "X(2,5) 2");
62
dpd_buf4_axpy(&X2, &X1, -1.0);
64
dpd_buf4_init(&newLijab, CC_LAMBDA, L_irr, 2, 5, 2, 7, 0, "New Lijab");
65
dpd_buf4_axpy(&X1, &newLijab, 1.0);
66
dpd_buf4_close(&newLijab);
68
dpd_buf4_init(&newLIjAb, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "New LIjAb");
70
dpd_buf4_init(&WMnIe, CC_HBAR, 0, 0, 11, 0, 11, 0, "WMnIe");
71
dpd_buf4_sort(&WMnIe, CC_TMP0, pqsr, 0, 10, "WMnIe (Mn,Ie)");
72
dpd_buf4_close(&WMnIe);
74
dpd_buf4_init(&WMnIe, CC_TMP0, 0, 0, 10, 0, 10, 0, "WMnIe (Mn,Ie)");
75
dpd_contract244(&LIA, &WMnIe, &newLIjAb, 0, 2, 1, -1.0, 1.0);
76
dpd_buf4_close(&WMnIe);
78
dpd_buf4_init(&WmNiE, CC_HBAR, 0, 0, 11, 0, 11, 0, "WmNiE");
79
dpd_buf4_sort(&WmNiE, CC_TMP0, qprs, 0, 11, "WmNiE (Nm,Ei)");
80
dpd_buf4_close(&WmNiE);
82
/* W(Nm,Ei) * L(i,b) --> L(Nm,Eb) */
83
dpd_buf4_init(&WmNiE, CC_TMP0, 0, 0, 11, 0, 11, 0, "WmNiE (Nm,Ei)");
84
dpd_contract424(&WmNiE, &Lia, &newLIjAb, 3, 0, 0, -1.0, 1.0);
85
dpd_buf4_close(&WmNiE);
87
dpd_buf4_close(&newLIjAb);
89
dpd_file2_close(&Lia);
90
dpd_file2_close(&LIA);
92
else if(params.ref == 2) { /** UHF **/
94
dpd_file2_init(&LIA, CC_OEI, L_irr, 0, 1, "LIA");
95
dpd_file2_init(&Lia, CC_OEI, L_irr, 2, 3, "Lia");
97
/** W(IJ,AM) L(M,B) --> Z(IJ,AB) **/
98
dpd_buf4_init(&Z, CC_TMP2, L_irr, 2, 5, 2, 5, 0, "Z'(IJ,AB)");
99
dpd_buf4_init(&W, CC_HBAR, 0, 2, 21, 2, 21, 0, "WMNIE");
100
dpd_contract424(&W, &LIA, &Z, 3, 0, 0, 1, 0);
102
/** Z(IJ,AB) --> Z(IJ,BA) **/
103
dpd_buf4_sort(&Z, CC_TMP2, pqsr, 2, 5, "Z'(IJ,BA)");
105
/** Z(IJ,AB) = Z(IJ,AB) - Z(IJ,BA) **/
106
dpd_buf4_init(&Z1, CC_TMP2, L_irr, 2, 5, 2, 5, 0, "Z'(IJ,AB)");
107
dpd_buf4_init(&Z2, CC_TMP2, L_irr, 2, 5, 2, 5, 0, "Z'(IJ,BA)");
108
dpd_buf4_axpy(&Z2, &Z1, -1);
110
/** Z(IJ,AB) --> New L(IJ,AB) **/
111
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 2, 5, 2, 7, 0, "New LIJAB");
112
dpd_buf4_axpy(&Z1, &L2, 1);
117
/** W(ij,am) L(m,b) --> Z(ij,ab) **/
118
dpd_buf4_init(&Z, CC_TMP2, L_irr, 12, 15, 12, 15, 0, "Z'(ij,ab)");
119
dpd_buf4_init(&W, CC_HBAR, 0, 12, 31, 12, 31, 0, "Wmnie");
120
dpd_contract424(&W, &Lia, &Z, 3, 0, 0, 1, 0);
122
/** Z(ij,ab) --> Z(ij,ba) **/
123
dpd_buf4_sort(&Z, CC_TMP2, pqsr, 12, 15, "Z'(ij,ba)");
125
/** Z(ij,ab) = Z(ij,ab) - Z(ij,ba) **/
126
dpd_buf4_init(&Z1, CC_TMP2, L_irr, 12, 15, 12, 15, 0, "Z'(ij,ab)");
127
dpd_buf4_init(&Z2, CC_TMP2, L_irr, 12, 15, 12, 15, 0, "Z'(ij,ba)");
128
dpd_buf4_axpy(&Z2, &Z1, -1);
130
/** Z(ij,ab) --> New L(ij,ab) **/
131
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 12, 15, 12, 17, 0, "New Lijab");
132
dpd_buf4_axpy(&Z1, &L2, 1);
137
/** Z(jI,Ab) = W(jI,Am) L(m,b) **/
138
dpd_buf4_init(&Z, CC_TMP2, L_irr, 23, 28, 23, 28, 0, "Z(jI,Ab)");
139
dpd_buf4_init(&W, CC_HBAR, 0, 23, 26, 23, 26, 0, "WmNiE");
140
dpd_contract424(&W, &Lia, &Z, 3, 0, 0, -1, 0);
142
/** Z(jI,Ab) --> New L(Ij,Ab) **/
143
dpd_buf4_sort_axpy(&Z, CC_LAMBDA, qprs, 22, 28, "New LIjAb", 1);
146
/** Z(Ij,bA) = W(Ij,bM) L(M,A) **/
147
dpd_buf4_init(&Z, CC_TMP2, L_irr, 22, 29, 22, 29, 0, "Z(Ij,bA)");
148
dpd_buf4_init(&W, CC_HBAR, 0, 22, 25, 22, 25, 0, "WMnIe");
149
dpd_contract424(&W, &LIA, &Z, 3, 0, 0, -1, 0);
151
/** Z(Ij,bA) --> New L(Ij,Ab) **/
152
dpd_buf4_sort_axpy(&Z, CC_LAMBDA, pqsr, 22, 28, "New LIjAb", 1);
155
dpd_file2_close(&Lia);
156
dpd_file2_close(&LIA);
added
removed
src/bin/cclambda/WijmbL2.c
