3
\brief Enter brief description of file here
6
#include <libdpd/dpd.h>
13
namespace psi { namespace cchbar {
15
/* Wamef_build(): Computes all contributions to the Wamef HBAR matrix
16
** elements, whose spin-orbital definition is:
18
** Wamef = <am||ef> - t_n^a <nm||ef>
20
** (cf. Gauss and Stanton, JCP 103, 3561-3577 (1995).)
22
** The storage and naming convention for each spin case are
25
** Spin Case Storage Name
26
** ---------- --------- --------
27
** WAMEF (MA,E>F) "WAMEF"
28
** Wamef (ma,e>f) "Wamef"
29
** WAmEf (mA,Ef) "WAmEf"
30
** WaMeF (Ma,eF) "WaMeF"
31
** -----------------------------------
34
** RHF Cases: Note that only the WAmEf spin case is required, and
35
** we store it AS WRITTEN, (Am,Ef).
38
** For all spin cases, we now use only the following
39
** WAMEF (AM,E>F) "WAMEF"
40
** Wamef (am,e>f) "Wamef"
41
** WAmEf (Am,Ef) "WAmEf"
42
** WaMeF (aM,eF) "WaMeF"
45
** For CC3, these are computed by ccenergy in file CC_HET1
49
void Wamef_build(void) {
50
dpdbuf4 Wamef, WAMEF, WAmEf, WaMeF, W;
53
int h, Ga, Gn, Gm, A, a, row, nrows, ncols;
57
dpd_buf4_init(&F, CC_FINTS, 0, 10, 5, 10, 5, 0, "F <ia|bc>");
58
dpd_buf4_sort(&F, CC_HBAR, qpsr, 11, 5, "WAmEf");
61
dpd_buf4_init(&W, CC_HBAR, 0, 11, 5, 11, 5, 0, "WAmEf");
62
dpd_file2_init(&tIA, CC_OEI, 0, 0, 1, "tIA");
63
dpd_file2_mat_init(&tIA);
64
dpd_file2_mat_rd(&tIA);
65
dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
67
/* dpd_contract244(&tIA, &D, &W, 0, 0, 0, -1, 1); */
68
/** OOC code below added 05/04/05, -TDC **/
69
for(h=0; h < moinfo.nirreps; h++) { /* h = Gam = Gnm = Gef */
71
dpd_buf4_mat_irrep_init(&D, h);
72
dpd_buf4_mat_irrep_rd(&D, h);
75
for(Ga=0; Ga < moinfo.nirreps; Ga++) {
77
Gn = Ga; /* T1 is totally symmetric */
79
W.matrix[h] = dpd_block_matrix(moinfo.occpi[Gm], W.params->coltot[h]);
81
for(A=0; A < moinfo.virtpi[Ga]; A++) {
82
a = moinfo.vir_off[Ga] + A;
84
dpd_buf4_mat_irrep_rd_block(&W, h, W.row_offset[h][a], moinfo.occpi[Gm]);
86
nrows = moinfo.occpi[Gn];
87
ncols = moinfo.occpi[Gm] * W.params->coltot[h];
90
C_DGEMV('t',nrows,ncols,-1.0,&(D.matrix[h][row][0]),ncols,&(tIA.matrix[Gn][0][A]),
91
moinfo.virtpi[Ga],1.0, W.matrix[h][0],1);
94
dpd_buf4_mat_irrep_wrt_block(&W, h, W.row_offset[h][a], moinfo.occpi[Gm]);
97
row += moinfo.occpi[Gn] * moinfo.occpi[Gm];
99
dpd_free_block(W.matrix[h], moinfo.occpi[Gm], W.params->coltot[h]);
102
dpd_buf4_mat_irrep_close(&D, h);
106
dpd_file2_mat_close(&tIA);
107
dpd_file2_close(&tIA);
111
else if(params.ref == 1) { /** ROHF **/
113
dpd_file2_init(&tIA, CC_OEI, 0, 0, 1, "tIA");
114
dpd_file2_init(&tia, CC_OEI, 0, 0, 1, "tia");
116
/* <AM||EF> --> W(AM,E>F) */
117
/* <am||ef> --> W(am,e>f) */
118
dpd_buf4_init(&F, CC_FINTS, 0, 11, 7, 11, 5, 1, "F <ai|bc>");
119
dpd_buf4_copy(&F, CC_HBAR, "WAMEF");
120
dpd_buf4_copy(&F, CC_HBAR, "Wamef");
123
/* T(N,A) <NM||EF> --> W(AM,E>F) */
124
dpd_buf4_init(&WAMEF, CC_HBAR, 0, 11, 7, 11, 7, 0, "WAMEF");
125
dpd_buf4_init(&D, CC_DINTS, 0, 0, 7, 0, 7, 0, "D <ij||ab> (ij,a>b)");
126
dpd_contract244(&tIA, &D, &WAMEF, 0, 0, 0, -1, 1);
128
dpd_buf4_close(&WAMEF);
130
/* T(n,a) <nm||ef> --> W(am,e>f) */
131
dpd_buf4_init(&Wamef, CC_HBAR, 0, 11, 7, 11, 7, 0, "Wamef");
132
dpd_buf4_init(&D, CC_DINTS, 0, 0, 7, 0, 7, 0, "D <ij||ab> (ij,a>b)");
133
dpd_contract244(&tia, &D, &Wamef, 0, 0, 0, -1, 1);
135
dpd_buf4_close(&Wamef);
137
/* <Am|Ef> --> W(Am,Ef) */
138
dpd_buf4_init(&F, CC_FINTS, 0, 11, 5, 11, 5, 0, "F <ai|bc>");
139
dpd_buf4_copy(&F, CC_HBAR, "WAmEf");
142
/* <aM|eF> --> W(aM,eF) */
143
dpd_buf4_init(&F, CC_FINTS, 0, 11, 5, 11, 5, 0, "F <ai|bc>");
144
dpd_buf4_copy(&F, CC_HBAR, "WaMeF");
147
/* T(N,A) <Nm|Ef> --> W(Am,Ef) */
148
dpd_buf4_init(&WAmEf, CC_HBAR, 0, 11, 5, 11, 5, 0, "WAmEf");
149
dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
150
dpd_contract244(&tIA, &D, &WAmEf, 0, 0, 0, -1, 1);
152
dpd_buf4_close(&WAmEf);
154
/* T(n,a) <nM|eF> --> W(aM,eF) */
155
dpd_buf4_init(&WaMeF, CC_HBAR, 0, 11, 5, 11, 5, 0, "WaMeF");
156
dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
157
dpd_contract244(&tia, &D, &WaMeF, 0, 0, 0, -1, 1);
159
dpd_buf4_close(&WaMeF);
161
dpd_file2_close(&tIA);
162
dpd_file2_close(&tia);
165
else if(params.ref == 2) { /** UHF **/
167
dpd_file2_init(&tIA, CC_OEI, 0, 0, 1, "tIA");
168
dpd_file2_init(&tia, CC_OEI, 0, 2, 3, "tia");
170
/* <AM||EF> --> W(AM,E>F) */
171
dpd_buf4_init(&F, CC_FINTS, 0, 21, 7, 21, 5, 1, "F <AI|BC>");
172
dpd_buf4_copy(&F, CC_HBAR, "WAMEF");
175
/* <am||ef> --> W(am,e>f) */
176
dpd_buf4_init(&F, CC_FINTS, 0, 31, 17, 31, 15, 1, "F <ai|bc>");
177
dpd_buf4_copy(&F, CC_HBAR, "Wamef");
180
/* T(N,A) <NM||EF> --> W(AM,E>F) */
181
dpd_buf4_init(&WAMEF, CC_HBAR, 0, 21, 7, 21, 7, 0, "WAMEF");
182
dpd_buf4_init(&D, CC_DINTS, 0, 0, 7, 0, 7, 0, "D <IJ||AB> (IJ,A>B)");
183
dpd_contract244(&tIA, &D, &WAMEF, 0, 0, 0, -1, 1);
185
dpd_buf4_close(&WAMEF);
187
/* T(n,a) <nm||ef> --> W(am,e>f) */
188
dpd_buf4_init(&Wamef, CC_HBAR, 0, 31, 17, 31, 17, 0, "Wamef");
189
dpd_buf4_init(&D, CC_DINTS, 0, 10, 17, 10, 17, 0, "D <ij||ab> (ij,a>b)");
190
dpd_contract244(&tia, &D, &Wamef, 0, 0, 0, -1, 1);
192
dpd_buf4_close(&Wamef);
194
/* <Am|Ef> --> W(Am,Ef) */
195
dpd_buf4_init(&F, CC_FINTS, 0, 26, 28, 26, 28, 0, "F <Ai|Bc>");
196
dpd_buf4_copy(&F, CC_HBAR, "WAmEf");
199
/* <aM|eF> --> W(aM,eF) */
200
dpd_buf4_init(&F, CC_FINTS, 0, 25, 29, 25, 29, 0, "F <aI|bC>");
201
dpd_buf4_copy(&F, CC_HBAR, "WaMeF");
204
/* T(N,A) <Nm|Ef> --> W(Am,Ef) */
205
dpd_buf4_init(&WAmEf, CC_HBAR, 0, 26, 28, 26, 28, 0, "WAmEf");
206
dpd_buf4_init(&D, CC_DINTS, 0, 22, 28, 22, 28, 0, "D <Ij|Ab>");
207
dpd_contract244(&tIA, &D, &WAmEf, 0, 0, 0, -1, 1);
209
dpd_buf4_close(&WAmEf);
211
/* T(n,a) <nM|eF> --> W(aM,eF) */
212
dpd_buf4_init(&WaMeF, CC_HBAR, 0, 25, 29, 25, 29, 0, "WaMeF");
213
dpd_buf4_init(&D, CC_DINTS, 0, 23, 29, 23, 29, 0, "D <iJ|aB>");
214
dpd_contract244(&tia, &D, &WaMeF, 0, 0, 0, -1, 1);
216
dpd_buf4_close(&WaMeF);
218
dpd_file2_close(&tIA);
219
dpd_file2_close(&tia);
225
}} // namespace psi::cchbar