11
13
dpdbuf4 Z, Z1, Z2, I, Y2, L2, W;
16
int nirreps, Gbm, Gef, Gjf, Ge, Gf, Gj, bm, ef, jf;
17
int Gfe, Gb, Gm, b, m, B, M, e, f, fe, nrows, ncols, nlinks;
15
20
sprintf(lbl, "Z_%s_%1s_MI", pert_y, cart_y);
16
21
dpd_file2_init(&z1, CC_TMP0, irrep_y, 0, 0, lbl);
130
135
sprintf(lbl, "X_%s_%1s_IA (%5.3f)", pert_x, cart_x, omega_x);
131
136
dpd_file2_init(&X1, CC_OEI, irrep_x, 0, 1, lbl);
133
sprintf(lbl, "Z_%s_%1s_MbEj (Mb,jE)", pert_x, cart_x);
134
dpd_buf4_init(&Z1, CC_TMP0, irrep_x, 10, 10, 10, 10, 0, lbl);
135
dpd_buf4_init(&W, CC_HBAR, 0, 10, 5, 10, 5, 0, "WAmEf");
136
dpd_contract244(&X1, &W, &Z1, 1, 2, 1, 1, 0);
138
sprintf(lbl, "Z_%s_%1s_MbEj (bM,jE)", pert_x, cart_x);
139
dpd_buf4_init(&Z1, CC_TMP0, irrep_x, 11, 10, 11, 10, 0, lbl);
140
dpd_buf4_init(&W, CC_HBAR, 0, 11, 5, 11, 5, 0, "WAmEf");
141
/* dpd_contract244(&X1, &W, &Z1, 1, 2, 1, 1, 0); */
142
/* add out-of-core algorithm here */
143
/* Z(bm,jf) <-- X1(j,e) * W(bm,ef) */
144
dpd_file2_mat_init(&X1);
145
dpd_file2_mat_rd(&X1);
146
for(Gbm=0; Gbm < moinfo.nirreps; Gbm++) {
147
Gef = Gbm; Gjf = Gbm^irrep_x;
149
dpd_buf4_mat_irrep_init(&Z1, Gbm);
150
dpd_buf4_mat_irrep_row_init(&W, Gbm);
151
for(bm=0; bm < W.params->rowtot[Gbm]; bm++) {
152
dpd_buf4_mat_irrep_row_rd(&W, Gbm, bm);
154
for(Gj=0; Gj < moinfo.nirreps; Gj++) {
155
Gf = Gjf^Gj; Ge = Gef^Gf;
157
ef = W.col_offset[Gbm][Ge];
158
jf = Z1.col_offset[Gbm][Gj];
160
nrows = moinfo.occpi[Gj];
161
ncols = moinfo.virtpi[Gf];
162
nlinks = moinfo.virtpi[Ge];
164
if(nrows && ncols && nlinks)
165
C_DGEMM('n','n', nrows, ncols, nlinks, 1.0, &(X1.matrix[Gj][0][0]), nlinks,
166
&(W.matrix[Gbm][0][ef]), ncols, 0.0, &(Z1.matrix[Gbm][bm][jf]), ncols);
169
dpd_buf4_mat_irrep_row_close(&W, Gbm);
170
dpd_buf4_mat_irrep_wrt(&Z1, Gbm);
171
dpd_buf4_mat_irrep_close(&Z1, Gbm);
174
dpd_file2_mat_close(&X1);
175
/* out-of-core algorithm done */
137
176
dpd_buf4_close(&W);
138
178
sprintf(lbl, "Z_%s_%1s_MbEj (ME,jb)", pert_x, cart_x);
139
dpd_buf4_sort(&Z1, CC_TMP0, psrq, 10, 10, lbl);
179
dpd_buf4_sort(&Z1, CC_TMP0, qsrp, 10, 10, lbl);
140
180
dpd_buf4_close(&Z1);
142
182
sprintf(lbl, "Z_%s_%1s_WMbEj (bM,Ej)", pert_x, cart_x);
143
183
dpd_buf4_init(&Z1, CC_TMP0, irrep_x, 11, 11, 11, 11, 0, lbl);
144
dpd_buf4_init(&W, CC_HBAR, 0, 0, 11, 0, 11, 0, "WMnIe");
184
dpd_buf4_init(&W, CC_HBAR, 0, 0, 11, 0, 11, 0, "WMnIe (Mn,eI)");
145
185
dpd_contract244(&X1, &W, &Z1, 0, 0, 0, -1, 0);
146
186
dpd_buf4_close(&W);
147
187
sprintf(lbl, "Z_%s_%1s_MbEj (ME,jb)", pert_x, cart_x);
151
191
sprintf(lbl, "Z_%s_%1s_MbeJ (Mb,eJ)", pert_x, cart_x);
152
192
dpd_buf4_init(&Z1, CC_TMP0, irrep_x, 10, 11, 10, 11, 0, lbl);
153
dpd_buf4_init(&W, CC_HBAR, 0, 0, 11, 0, 11, 0, "WMnIe");
193
dpd_buf4_init(&W, CC_HBAR, 0, 0, 11, 0, 11, 0, "WMnIe (Mn,eI)");
154
194
dpd_contract424(&W, &X1, &Z1, 1, 0, 1, 1, 0);
155
195
dpd_buf4_close(&W);
156
dpd_buf4_init(&W, CC_HBAR, 0, 10, 5, 10, 5, 0, "WAmEf");
157
dpd_contract424(&W, &X1, &Z1, 3, 1, 0, -1, 1);
159
196
sprintf(lbl, "Z_%s_%1s_MbeJ (Me,Jb)", pert_x, cart_x);
160
197
dpd_buf4_sort(&Z1, CC_TMP0, prsq, 10, 10, lbl);
161
198
dpd_buf4_close(&Z1);
200
sprintf(lbl, "Z_%s_%1s_MbeJ (bM,eJ)", pert_x, cart_x);
201
dpd_buf4_init(&Z1, CC_TMP0, irrep_x, 11, 11, 11, 11, 0, lbl);
202
dpd_buf4_init(&W, CC_HBAR, 0, 11, 5, 11, 5, 0, "WAmEf");
203
dpd_contract424(&W, &X1, &Z1, 3, 1, 0, -1, 0);
205
sprintf(lbl, "Z_%s_%1s_MbeJ (Me,Jb)", pert_x, cart_x);
206
dpd_buf4_sort_axpy(&Z1, CC_TMP0, qrsp, 10, 10, lbl, 1);
163
209
dpd_file2_close(&X1);
165
211
sprintf(lbl, "Z_%s_%1s_MbEj (ME,jb)", pert_x, cart_x);
231
277
sprintf(lbl, "Z_%s_%1s_mj" , pert_x, cart_x);
232
278
dpd_file2_init(&z, CC_TMP0, irrep_x, 0, 0, lbl);
233
dpd_buf4_init(&W, CC_HBAR, 0, 0, 11, 0, 11, 0, "2WMnIe - WnMIe");
279
dpd_buf4_init(&W, CC_HBAR, 0, 0, 11, 0, 11, 0, "2WMnIe - WnMIe (Mn,eI)");
234
280
dpd_dot23(&X1, &W, &z, 0, 0, 1, 0);
235
281
dpd_buf4_close(&W);
236
282
dpd_buf4_init(&Z1, CC_TMP0, 0, 0, 5, 0, 5, 0, "Z(Ij,Ab)");
243
289
dpd_file2_close(&z);
247
sprintf(lbl, "Z_%s_%1s_ae", pert_x, cart_x);
292
sprintf(lbl, "Z_%s_%1s_AE", pert_x, cart_x);
248
293
dpd_file2_init(&z, CC_TMP0, irrep_x, 1, 1, lbl);
250
dpd_buf4_init(&W, CC_HBAR, 0, 11, 5, 11, 5, 0, "WAmEf 2(Am,Ef) - (Am,fE)");
251
dpd_dot24(&X1, &W, &z, 0, 0, 1, 0);
295
/* dpd_buf4_init(&W, CC_HBAR, 0, 11, 5, 11, 5, 0, "WAmEf 2(Am,Ef) - (Am,fE)"); */
296
/* dpd_dot24(&X1, &W, &z, 0, 0, 1, 0); */
297
dpd_buf4_init(&W, CC_HBAR, 0, 11, 5, 11, 5, 0, "WAmEf");
298
dpd_file2_scm(&z, 0);
299
dpd_file2_mat_init(&z);
300
dpd_file2_mat_init(&X1);
301
dpd_file2_mat_rd(&X1);
302
for(Gbm=0; Gbm < moinfo.nirreps; Gbm++) {
303
Gfe = Gbm; /* W is totally symmetric */
304
dpd_buf4_mat_irrep_row_init(&W, Gbm);
305
X = init_array(W.params->coltot[Gfe]);
306
for(bm=0; bm < W.params->rowtot[Gbm]; bm++) {
307
dpd_buf4_mat_irrep_row_rd(&W, Gbm, bm);
308
b = W.params->roworb[Gbm][bm][0];
309
m = W.params->roworb[Gbm][bm][1];
310
Gb = W.params->psym[b];
314
B = b - moinfo.vir_off[Gb];
315
M = m - moinfo.occ_off[Gm];
316
zero_arr(X, W.params->coltot[Gfe]);
317
for(fe=0; fe < W.params->coltot[Gfe]; fe++) {
318
f = W.params->colorb[Gfe][fe][0];
319
e = W.params->colorb[Gfe][fe][1];
320
ef = W.params->colidx[e][f];
321
X[fe] = 2.0 * W.matrix[Gbm][0][fe] - W.matrix[Gbm][0][ef];
323
nrows = moinfo.virtpi[Gf];
324
ncols = moinfo.virtpi[Ge];
326
C_DGEMV('n',nrows,ncols,1,&X[W.col_offset[Gfe][Gf]],ncols,
327
X1.matrix[Gm][M],1,1,z.matrix[Gb][B],1);
331
dpd_buf4_mat_irrep_row_close(&W, Gbm);
333
dpd_file2_mat_close(&X1);
334
dpd_file2_mat_wrt(&z);
335
dpd_file2_mat_close(&z);
252
336
dpd_buf4_close(&W);
253
338
dpd_buf4_init(&Z1, CC_TMP0, 0, 0, 5, 0, 5, 0, "Z(Ij,Ab) temp");
254
339
dpd_contract424(&Y2, &z, &Z1, 3, 1, 0, 1, 0);
255
340
dpd_buf4_init(&Z, CC_TMP0, 0, 0, 5, 0, 5, 0, "Z(Ij,Ab) Final");
267
352
dpd_file2_init(&X1, CC_OEI, irrep_x, 0, 1, lbl);
268
353
sprintf(lbl, "Z_%s_%1s_MnjI", pert_x, cart_x);
269
354
dpd_buf4_init(&Z1, CC_TMP0, irrep_x, 0, 0, 0, 0, 0, lbl);
270
dpd_buf4_init(&W, CC_HBAR, 0, 0, 11, 0, 11, 0, "WMnIe");
355
dpd_buf4_init(&W, CC_HBAR, 0, 0, 11, 0, 11, 0, "WMnIe (Mn,eI)");
271
356
dpd_contract244(&X1, &W, &Z1, 1, 2, 1, 1, 0);
272
357
dpd_buf4_close(&W);
273
358
sprintf(lbl, "Z_%s_%1s_MnIj", pert_x, cart_x);
286
371
dpd_buf4_close(&Z1);
287
372
dpd_buf4_close(&Z);
289
sprintf(lbl, "Z_%s_%1s_AmIj (mA,Ij)", pert_y, cart_y);
290
dpd_buf4_init(&Z1, CC_TMP0, irrep_y, 10, 0, 10, 0, 0, lbl);
374
sprintf(lbl, "Z_%s_%1s_AmIj (Am,Ij)", pert_y, cart_y);
375
dpd_buf4_init(&Z1, CC_TMP0, irrep_y, 11, 0, 11, 0, 0, lbl);
291
376
sprintf(lbl, "X_%s_%1s_IjAb (%5.3f)", pert_y, cart_y, omega_y);
292
377
dpd_buf4_init(&Y2, CC_LR, irrep_y, 0, 5, 0, 5, 0, lbl);
293
dpd_buf4_init(&W, CC_HBAR, 0, 10, 5, 10, 5, 0, "WAmEf");
378
dpd_buf4_init(&W, CC_HBAR, 0, 11, 5, 11, 5, 0, "WAmEf");
294
379
dpd_contract444(&W, &Y2, &Z1, 0, 0, 1, 0);
295
380
dpd_buf4_close(&W);
296
381
dpd_buf4_close(&Y2);
298
dpd_buf4_init(&Z, CC_TMP0, 0, 0, 5, 0, 5, 0, "Z(Ij,bA)");
383
dpd_buf4_init(&Z, CC_TMP0, 0, 0, 5, 0, 5, 0, "Z(Ij,Ab) Final");
299
384
sprintf(lbl, "X_%s_%1s_IA (%5.3f)", pert_x, cart_x, omega_x);
300
385
dpd_file2_init(&X1, CC_OEI, irrep_x, 0, 1, lbl);
301
dpd_contract244(&X1, &Z1, &Z, 0, 0, 1, -2, 0);
386
dpd_contract424(&Z1, &X1, &Z, 1, 0, 0, -2, 1);
302
387
dpd_file2_close(&X1);
303
388
dpd_buf4_close(&Z1);
304
dpd_buf4_sort_axpy(&Z, CC_TMP0, pqsr, 0, 5, "Z(Ij,Ab) Final", 1);
305
389
dpd_buf4_close(&Z);
307
391
dpd_buf4_init(&Z, CC_TMP0, 0, 0, 5, 0, 5, 0, "Z(Ij,Ab) Final");