3
#include <libipv1/ip_lib.h>
4
#include <libciomr/libciomr.h>
5
#include <libiwl/iwl.h>
6
#include <libchkpt/chkpt.h>
11
#include "yoshimine.h"
13
#define INDEX(i,j) ((i>j) ? (ioff[(i)]+(j)) : (ioff[(j)]+(i)))
14
#define MAX0(a,b) (((a)>(b)) ? (a) : (b))
16
double fzc_energy_uhf(int nbfso, int *sosym, double *Pa, double *Pb,
17
double *Hca, double *Hcb,double *H,
18
int *first_so, int *ioff);
20
/* transform_two_uhf(): Carry out the transformation of the SO-basis
21
** two-electron integrals into the AA, BB, and AB MO-basis classes.
22
** In these comments, I use p,q,r,s to denote SO-basis indices,
23
** I,J,K,L to denote alpha-MO indices, and i,j,k,l to denote beta-MO
26
** (1) Yoshimine presort the SO-basis integrals into a supermatrix,
27
** with p>=q, r>=s, but with all pq and rs (just like the usual
28
** presort for integral transformations in transform_two.c).
30
** (2) Half-transform the integrals into the alpha-MO basis to
31
** form a (pq,IJ) matrix and Yoshimine sort the result into the
32
** transposed form, (IJ,pq).
34
** (3) Complete the second half-transformation for the AA and AB
35
** integral lists: (IJ,KL) and (IJ,kl).
37
** (4) Half-transform the SO-basis integrals into the beta-MO basis to
38
** form a (pq,ij) matrix and Yoshimine sort the result into the
39
** transposed form, (ij,pq).
41
** (5) Complete the second half-transformation for the BB integral
48
void transform_two_uhf(void)
51
long int maxcor, maxcord;
52
int max_buckets, first_tmp_file;
56
int *src_first, *src_last, *src_orbspi, *src_orbsym, src_orbs, src_ntri;
57
int *dst_first, *dst_last, *dst_orbspi, *dst_orbsym, dst_orbs, dst_ntri;
58
int i, j, ifirst, ilast, jfirst, jlast, ij;
59
int k, l, ksym, lsym, kfirst, klast, lfirst, llast, kl, klsym;
60
int p, q, psym, qsym, pfirst, plast, qfirst, qlast, pq, pqsym;
61
int r, s, rsym, ssym, rfirst, rlast, sfirst, slast, rs;
62
int P, Q, R, S, K, L, ktr, ltr;
63
struct yoshimine YBuffP;
64
struct yoshimine YBuffJ;
68
double *P_block, *J_block;
69
double **A, **B, ***C;
70
int A_cols, B_cols, *C_colspi;
71
int *reorder_alpha, *reorder_beta;
73
maxcor = params.maxcor;
74
maxcord = params.maxcord;
75
max_buckets = params.max_buckets;
76
first_tmp_file = params.first_tmp_file;
77
tolerance = params.tolerance;
78
print_lvl = params.print_lvl;
80
nirreps = moinfo.nirreps;
81
reorder_alpha = moinfo.order_alpha;
82
reorder_beta = moinfo.order_beta;
84
src_first = moinfo.first_so;
85
src_last = moinfo.last_so;
86
src_orbspi = moinfo.sopi;
87
src_orbsym = moinfo.sosym;
88
src_orbs = moinfo.nso;
89
src_ntri = src_orbs * (src_orbs+1)/2;
91
dst_first = (params.fzc && !params.do_all_tei) ? moinfo.fstact : moinfo.first;
92
dst_last = (!params.do_all_tei) ? moinfo.lstact : moinfo.last;
93
dst_orbspi = (!params.do_all_tei) ? moinfo.active : moinfo.orbspi;
94
dst_orbsym = moinfo.orbsym;
95
dst_orbs = (!params.do_all_tei) ? (moinfo.nmo - moinfo.nfzv - moinfo.nfzc) : moinfo.nmo;
96
dst_ntri = moinfo.nmo * (moinfo.nmo+1)/2;
97
C_colspi = dst_orbspi;
99
fzc_offset = (params.fzc && !params.do_all_tei) ? moinfo.nfzc : 0;
101
/** Presort the two-electron integrals **/
103
if (params.print_lvl) {
104
fprintf(outfile, "\n\tPre-sorting two-electron integrals...\n\n");
108
yosh_init(&YBuffP, src_ntri, src_ntri, maxcor, maxcord,
109
max_buckets, first_tmp_file, tolerance, outfile);
112
fprintf(outfile, "\tPresort Yoshimine Parameters");
113
yosh_print(&YBuffP, outfile);
114
fprintf(outfile, "\n");
118
yosh_init_buckets(&YBuffP);
120
yosh_rdtwo_uhf(&YBuffP, params.src_tei_file,
121
params.delete_src_tei, src_orbspi, nirreps, ioff, 0,
122
params.fzc && moinfo.nfzc, moinfo.fzc_density_alpha,
123
moinfo.fzc_density_beta, moinfo.fzc_operator_alpha,
124
moinfo.fzc_operator_beta, 1, (print_lvl > 5), outfile);
126
yosh_close_buckets(&YBuffP, 0);
128
yosh_sort(&YBuffP, params.presort_file, 0, ioff, NULL, src_orbs, src_ntri,
129
0, 1, 0, 0, 1, (print_lvl > 5), outfile);
133
/** Pre-sort complete **/
135
/** Compute the frozen core energy **/
137
if (params.fzc && moinfo.nfzc) {
139
fzc_energy_uhf(moinfo.nso, moinfo.sosym, moinfo.fzc_density_alpha,
140
moinfo.fzc_density_beta, moinfo.fzc_operator_alpha,
141
moinfo.fzc_operator_beta, moinfo.oe_ints,
142
moinfo.first_so, ioff);
143
free(moinfo.fzc_density_alpha);
144
free(moinfo.fzc_density_beta);
147
/* Write efzc to chkpt file */
148
chkpt_init(PSIO_OPEN_OLD);
149
chkpt_wt_efzc(moinfo.efzc);
152
if (params.print_lvl)
153
fprintf(outfile, "\n\tFrozen core energy = %20.15lf\n", moinfo.efzc);
155
P_block = init_array(src_ntri);
156
J_block = init_array(MAX0(src_ntri,dst_ntri));
158
A_cols = MAX0(src_orbs, dst_orbs);
160
A = block_matrix(A_cols,A_cols);
161
B = block_matrix(src_orbs, dst_orbs);
163
/** AA/AB First-half transformation **/
165
fprintf(outfile, "\n\tBeginning AA/AB two-electron transform...\n");
167
C = moinfo.evects_alpha;
169
iwl_buf_init(&PBuff, params.presort_file, tolerance, 1, 1);
170
yosh_init(&YBuffJ, dst_ntri, src_ntri, maxcor, maxcord, max_buckets,
171
first_tmp_file, tolerance, outfile);
173
yosh_init_buckets(&YBuffJ);
176
fprintf(outfile, "\tHalf-transform Yoshimine Parameters");
177
yosh_print(&YBuffJ, outfile);
178
fprintf(outfile, "\n");
182
for (psym=0; psym < nirreps; psym++) {
183
pfirst = src_first[psym];
184
plast = src_last[psym];
185
for (p=pfirst; p <= plast; p++) {
186
for (qsym=0; qsym <= psym; qsym++) {
187
qfirst = src_first[qsym];
188
qlast = src_last[qsym];
190
for (q=qfirst; (q<=qlast) && (q <= p); q++) {
193
zero_arr(P_block,src_ntri);
194
iwl_buf_rd(&PBuff, pq, P_block, ioff, ioff, 0, (print_lvl>4),
197
for (rsym=0; rsym < nirreps; rsym++) {
198
rfirst = src_first[rsym];
199
rlast = src_last[rsym];
200
kfirst = dst_first[rsym];
201
klast = dst_last[rsym];
204
sfirst = src_first[ssym];
205
slast = src_last[ssym];
206
lfirst = dst_first[ssym];
207
llast = dst_last[ssym];
209
for(r=rfirst,R=0; r <= rlast; r++,R++) {
210
for(s=sfirst,S=0; (s <= slast) && (s <= r);
213
A[R][S] = P_block[rs];
214
A[S][R] = P_block[rs];
219
for (r=rfirst,R=0; r <= rlast; r++,R++) {
220
for (s=sfirst,S=0; s <= slast; s++,S++) {
222
A[R][S] = P_block[rs];
227
if (C_colspi[ssym] > 0)
228
C_DGEMM('n', 'n', src_orbspi[rsym],
229
dst_orbspi[ssym], src_orbspi[ssym], 1.0,
230
A[0], A_cols, C[ssym][0], C_colspi[ssym], 0.0,
232
if (C_colspi[rsym] > 0)
233
C_DGEMM('t', 'n', dst_orbspi[rsym],
234
dst_orbspi[ssym], src_orbspi[rsym], 1.0,
235
C[rsym][0], C_colspi[rsym], B[0], B_cols, 0.0,
238
zero_arr(J_block, dst_ntri);
239
for (k=kfirst,K=0; k <= klast; k++,K++) {
240
for (l=lfirst,L=0; (l <= llast) && (l <= k); l++,L++) {
242
J_block[kl] = A[K][L];
246
yosh_wrt_arr(&YBuffJ, p, q, pq, pqsym, J_block,
247
moinfo.nmo, ioff, dst_orbsym, dst_first, dst_last,
248
1, (print_lvl > 4), outfile);
258
iwl_buf_close(&PBuff, 1);
260
if (params.print_lvl)
261
fprintf(outfile, "\tSorting AA/AB half-transformed integrals...\n");
264
yosh_close_buckets(&YBuffJ, 0);
265
yosh_sort(&YBuffJ, params.jfile, 0, ioff, NULL, src_orbs, src_ntri, 0, 1, 0, 0,
266
1, (print_lvl > 5), outfile);
270
fprintf(outfile, "\tFinished AA/AB half-transformation...\n");
274
/** AA Second-half transformation **/
276
C = moinfo.evects_alpha;
278
iwl_buf_init(&JBuff, params.jfile, tolerance, 1, 1);
279
iwl_buf_init(&MBuff, params.aa_mfile, tolerance, 0, 0);
281
for (ksym=0; ksym < nirreps; ksym++) {
282
kfirst = dst_first[ksym];
283
klast = dst_last[ksym];
284
for (k=kfirst; k <= klast; k++) {
285
for (lsym=0; lsym <= ksym; lsym++) {
286
lfirst = dst_first[lsym];
287
llast = dst_last[lsym];
289
for (l=lfirst; (l <= llast) && (l <= k); l++) {
291
if (!params.backtr) {
292
ktr = reorder_alpha[k] - fzc_offset;
293
ltr = reorder_alpha[l] - fzc_offset;
296
zero_arr(J_block, src_ntri);
297
iwl_buf_rd(&JBuff, kl, J_block, ioff, ioff, 0, 0, outfile);
299
for (psym=0; psym < nirreps; psym++) {
300
pfirst = src_first[psym];
301
plast = src_last[psym];
302
ifirst = dst_first[psym];
303
ilast = dst_last[psym];
306
qfirst = src_first[qsym];
307
qlast = src_last[qsym];
308
jfirst = dst_first[qsym];
309
jlast = dst_last[qsym];
310
for (p=pfirst,P=0; p <= plast; p++,P++) {
311
for (q=qfirst,Q=0; q <= qlast; q++,Q++) {
313
A[P][Q] = J_block[pq];
318
if (C_colspi[qsym] > 0)
319
C_DGEMM('n', 'n', src_orbspi[psym],
320
dst_orbspi[qsym], src_orbspi[qsym], 1.0,
321
A[0], A_cols, C[qsym][0], C_colspi[qsym], 0.0,
323
if (C_colspi[psym] > 0)
324
C_DGEMM('t', 'n', dst_orbspi[psym],
325
dst_orbspi[qsym], src_orbspi[psym], 1.0,
326
C[psym][0], C_colspi[psym], B[0], B_cols, 0.0,
329
iwl_buf_wrt_mat(&MBuff, ktr, ltr, A,
330
ifirst, ilast, jfirst, jlast,
331
reorder_alpha, fzc_offset, params.print_te_ints,
340
iwl_buf_close(&JBuff, 1);
342
iwl_buf_flush(&MBuff, 1);
343
iwl_buf_close(&MBuff, 1);
346
fprintf(outfile, "\n\tAA Transformation finished.\n");
347
fprintf(outfile, "\tTwo-electron AA integrals written to file%d.\n",params.aa_mfile);
351
/** AB Second-half transformation **/
353
C = moinfo.evects_beta;
355
iwl_buf_init(&JBuff, params.jfile, tolerance, 1, 1);
356
iwl_buf_init(&MBuff, params.ab_mfile, tolerance, 0, 0);
358
for (ksym=0; ksym < nirreps; ksym++) {
359
kfirst = dst_first[ksym];
360
klast = dst_last[ksym];
361
for (k=kfirst; k <= klast; k++) {
362
for (lsym=0; lsym <= ksym; lsym++) {
363
lfirst = dst_first[lsym];
364
llast = dst_last[lsym];
366
for (l=lfirst; (l <= llast) && (l <= k); l++) {
368
if (!params.backtr) {
369
ktr = reorder_alpha[k] - fzc_offset;
370
ltr = reorder_alpha[l] - fzc_offset;
373
zero_arr(J_block, src_ntri);
374
iwl_buf_rd(&JBuff, kl, J_block, ioff, ioff, 0, 0, outfile);
376
for (psym=0; psym < nirreps; psym++) {
377
pfirst = src_first[psym];
378
plast = src_last[psym];
379
ifirst = dst_first[psym];
380
ilast = dst_last[psym];
383
qfirst = src_first[qsym];
384
qlast = src_last[qsym];
385
jfirst = dst_first[qsym];
386
jlast = dst_last[qsym];
387
for (p=pfirst,P=0; p <= plast; p++,P++) {
388
for (q=qfirst,Q=0; q <= qlast; q++,Q++) {
390
A[P][Q] = J_block[pq];
395
if (C_colspi[qsym] > 0)
396
C_DGEMM('n', 'n', src_orbspi[psym],
397
dst_orbspi[qsym], src_orbspi[qsym], 1.0,
398
A[0], A_cols, C[qsym][0], C_colspi[qsym], 0.0,
400
if (C_colspi[psym] > 0)
401
C_DGEMM('t', 'n', dst_orbspi[psym],
402
dst_orbspi[qsym], src_orbspi[psym], 1.0,
403
C[psym][0], C_colspi[psym], B[0], B_cols, 0.0,
406
iwl_buf_wrt_mat2(&MBuff, ktr, ltr, A,
407
ifirst, ilast, jfirst, jlast,
408
reorder_beta, fzc_offset, params.print_te_ints,
417
iwl_buf_close(&JBuff, 0);
419
iwl_buf_flush(&MBuff, 1);
420
iwl_buf_close(&MBuff, 1);
423
fprintf(outfile, "\n\tAB Transformation finished.\n");
424
fprintf(outfile, "\tTwo-electron AB integrals written to file%d.\n",params.ab_mfile);
428
/** BB First-half transformation **/
430
fprintf(outfile, "\n\tBeginning BB two-electron transform...\n");
432
C = moinfo.evects_beta;
434
iwl_buf_init(&PBuff, params.presort_file, tolerance, 1, 1);
435
yosh_init(&YBuffJ, dst_ntri, src_ntri, maxcor, maxcord, max_buckets,
436
first_tmp_file, tolerance, outfile);
438
yosh_init_buckets(&YBuffJ);
441
fprintf(outfile, "\tHalf-transform Yoshimine Parameters");
442
yosh_print(&YBuffJ, outfile);
443
fprintf(outfile, "\n");
447
for (psym=0; psym < nirreps; psym++) {
448
pfirst = src_first[psym];
449
plast = src_last[psym];
450
for (p=pfirst; p <= plast; p++) {
451
for (qsym=0; qsym <= psym; qsym++) {
452
qfirst = src_first[qsym];
453
qlast = src_last[qsym];
455
for (q=qfirst; (q<=qlast) && (q <= p); q++) {
458
zero_arr(P_block,src_ntri);
459
iwl_buf_rd(&PBuff, pq, P_block, ioff, ioff, 0, (print_lvl>4),
462
for (rsym=0; rsym < nirreps; rsym++) {
463
rfirst = src_first[rsym];
464
rlast = src_last[rsym];
465
kfirst = dst_first[rsym];
466
klast = dst_last[rsym];
469
sfirst = src_first[ssym];
470
slast = src_last[ssym];
471
lfirst = dst_first[ssym];
472
llast = dst_last[ssym];
474
for(r=rfirst,R=0; r <= rlast; r++,R++) {
475
for(s=sfirst,S=0; (s <= slast) && (s <= r);
478
A[R][S] = P_block[rs];
479
A[S][R] = P_block[rs];
484
for (r=rfirst,R=0; r <= rlast; r++,R++) {
485
for (s=sfirst,S=0; s <= slast; s++,S++) {
487
A[R][S] = P_block[rs];
492
if (C_colspi[ssym] > 0)
493
C_DGEMM('n', 'n', src_orbspi[rsym],
494
dst_orbspi[ssym], src_orbspi[ssym], 1.0,
495
A[0], A_cols, C[ssym][0], C_colspi[ssym], 0.0,
497
if (C_colspi[rsym] > 0)
498
C_DGEMM('t', 'n', dst_orbspi[rsym],
499
dst_orbspi[ssym], src_orbspi[rsym], 1.0,
500
C[rsym][0], C_colspi[rsym], B[0], B_cols, 0.0,
503
zero_arr(J_block, dst_ntri);
504
for (k=kfirst,K=0; k <= klast; k++,K++) {
505
for (l=lfirst,L=0; (l <= llast) && (l <= k); l++,L++) {
507
J_block[kl] = A[K][L];
511
yosh_wrt_arr(&YBuffJ, p, q, pq, pqsym, J_block,
512
moinfo.nmo, ioff, dst_orbsym, dst_first, dst_last,
513
1, (print_lvl > 4), outfile);
523
iwl_buf_close(&PBuff, params.keep_presort);
525
if (params.print_lvl)
526
fprintf(outfile, "\tSorting BB half-transformed integrals...\n");
529
yosh_close_buckets(&YBuffJ, 0);
530
yosh_sort(&YBuffJ, params.jfile, 0, ioff, NULL, src_orbs, src_ntri, 0, 1, 0, 0,
531
1, (print_lvl > 5), outfile);
535
fprintf(outfile, "\tFinished BB half-transformation...\n");
539
/** BB Second-half transformation **/
541
C = moinfo.evects_beta;
543
iwl_buf_init(&JBuff, params.jfile, tolerance, 1, 1);
544
iwl_buf_init(&MBuff, params.bb_mfile, tolerance, 0, 0);
546
for (ksym=0; ksym < nirreps; ksym++) {
547
kfirst = dst_first[ksym];
548
klast = dst_last[ksym];
549
for (k=kfirst; k <= klast; k++) {
550
for (lsym=0; lsym <= ksym; lsym++) {
551
lfirst = dst_first[lsym];
552
llast = dst_last[lsym];
554
for (l=lfirst; (l <= llast) && (l <= k); l++) {
556
if (!params.backtr) {
557
ktr = reorder_beta[k] - fzc_offset;
558
ltr = reorder_beta[l] - fzc_offset;
561
zero_arr(J_block, src_ntri);
562
iwl_buf_rd(&JBuff, kl, J_block, ioff, ioff, 0, 0, outfile);
564
for (psym=0; psym < nirreps; psym++) {
565
pfirst = src_first[psym];
566
plast = src_last[psym];
567
ifirst = dst_first[psym];
568
ilast = dst_last[psym];
571
qfirst = src_first[qsym];
572
qlast = src_last[qsym];
573
jfirst = dst_first[qsym];
574
jlast = dst_last[qsym];
575
for (p=pfirst,P=0; p <= plast; p++,P++) {
576
for (q=qfirst,Q=0; q <= qlast; q++,Q++) {
578
A[P][Q] = J_block[pq];
583
if (C_colspi[qsym] > 0)
584
C_DGEMM('n', 'n', src_orbspi[psym],
585
dst_orbspi[qsym], src_orbspi[qsym], 1.0,
586
A[0], A_cols, C[qsym][0], C_colspi[qsym], 0.0,
588
if (C_colspi[psym] > 0)
589
C_DGEMM('t', 'n', dst_orbspi[psym],
590
dst_orbspi[qsym], src_orbspi[psym], 1.0,
591
C[psym][0], C_colspi[psym], B[0], B_cols, 0.0,
594
iwl_buf_wrt_mat(&MBuff, ktr, ltr, A,
595
ifirst, ilast, jfirst, jlast,
596
reorder_beta, fzc_offset, params.print_te_ints,
605
iwl_buf_close(&JBuff, 0);
607
iwl_buf_flush(&MBuff, 1);
608
iwl_buf_close(&MBuff, 1);
611
fprintf(outfile, "\n\tBB Transformation finished.\n");
612
fprintf(outfile, "\tTwo-electron BB integrals written to file%d.\n",params.bb_mfile);
added
removed
src/bin/transqt/transform_two_uhf.c
