3
\brief Enter brief description of file here
6
** YOSHIMINE.C: Functions for the Yoshimine Sort Object
9
** Center for Computational Quantum Chemistry, UGA
12
** Made slightly more general to handle MP2 restricted transformations by
22
#include <libciomr/libciomr.h>
24
#include <libiwl/iwl.h>
27
#include "yoshimine.h"
30
namespace psi { namespace transqt {
32
extern struct MOInfo moinfo;
34
#define MIN0(a,b) (((a)<(b)) ? (a) : (b))
35
#define MAX0(a,b) (((a)>(b)) ? (a) : (b))
36
#define INDEX(i,j) ((i>j) ? (ioff[(i)]+(j)) : (ioff[(j)]+(i)))
39
** YOSH_INIT(): This function initializes a Yoshimine sort object.
40
** The data is contained in a structure YBuff.
43
** YBuff = pointer to struct which will hold data for the object
44
** bra_indices = the number of bra_index pairs for the two-electron
45
** integrals being sorted
46
** ket_indices = the number of ket_index pairs for the two-electron
47
** integrals being sorted
48
** maxcor = the core memory, in bytes
49
** maxcord = the core memory, in doubles
50
** max_buckets = the max number of buckets to use (may be limited due
51
** to the fact that each bucket requires a consecutively
52
** numbered binary temp file).
53
** first_tmp_file = the number of the first tmp file used in the
54
** Yoshimine sort (e.g. 80 for file80).
55
** cutoff = minimum value to be kept for any value during the sort
56
** outfile = the text output file
60
** Note: bra_indices and ket_indices replace nbstri in an attempt to somewhat
61
** generalize the sort for four-index quantities which index pairs may or
62
** may not be canonicalizable, e.g. integrals of (ov|ov) type, as may be
63
** found in MP2 energy calculations...the first and second (or third and
64
** fourth) indices are not necessarily interchangeable. Additionally,
65
** this modification has the added benefit that it will work when the
66
** number of left and right basis functions are not equal (e.g., when
67
** one has half-backtransformed integrals of the type (ao ao | mo mo)
68
** where nbfao != nbfmo.
70
void yosh_init(struct yoshimine *YBuff, unsigned bra_indices,
71
unsigned ket_indices, long maxcor,
72
long maxcord, const int max_buckets,
73
unsigned int first_tmp_file,
74
double cutoff, FILE *outfile)
76
unsigned long long int twoel_array_size; /*--- Although on 32-bit systems one can only allocate 2GB arrays
77
in 1 process space, one can store much bigger integrals files on disk ---*/
78
unsigned int nbuckets;
80
unsigned long int bytes_per_bucket, free_bytes_per_bucket;
82
YBuff->first_tmp_file = first_tmp_file;
83
twoel_array_size = bra_indices; twoel_array_size *= ket_indices;
84
YBuff->core_loads = (twoel_array_size - 1) / maxcord + 1 ;
85
nbuckets = YBuff->core_loads ;
86
YBuff->nbuckets = nbuckets;
87
YBuff->cutoff = cutoff;
88
YBuff->bra_indices = bra_indices;
89
YBuff->ket_indices = ket_indices;
90
if (nbuckets > max_buckets) {
91
fprintf(stderr, "(yosh_init): maximum number of buckets exceeded\n") ;
92
fprintf(outfile, "(yosh_init): maximum number of buckets exceeded\n") ;
93
fprintf(outfile, " wanted %d buckets\n", nbuckets) ;
95
exit(PSI_RETURN_FAILURE) ;
98
/* if the number of pq does not divide evenly among the buckets, then
99
* the last bucket will have the remainder of the pq's.
101
YBuff->pq_per_bucket = bra_indices / nbuckets ;
104
bytes_per_bucket = ((unsigned long int) (4*sizeof(int) + sizeof(double))) *
105
((unsigned long int) twoel_array_size) + (unsigned long int) (sizeof(struct iwlbuf)
106
+ IWL_INTS_PER_BUF * (4*sizeof(Label) + sizeof(Value)));
107
if (bytes_per_bucket > (unsigned long int) (maxcor/nbuckets))
108
bytes_per_bucket = (unsigned long int) (maxcor / nbuckets);
111
bytes_per_bucket = (unsigned long int) (maxcor / nbuckets);
113
free_bytes_per_bucket = bytes_per_bucket -
114
(unsigned long int) (sizeof(struct iwlbuf) + IWL_INTS_PER_BUF * (4*sizeof(Label) + sizeof(Value)));
115
YBuff->bucketsize = free_bytes_per_bucket / (4 * sizeof(int) +
117
YBuff->buckets = (struct bucket *) malloc(nbuckets * sizeof(struct bucket));
118
YBuff->bucket_for_pq = (int *) malloc (bra_indices * sizeof(int));
119
for (i=0,pq=0; i<nbuckets; i++) {
120
if (i != (nbuckets - 1)) {
121
for (j=0; j<YBuff->pq_per_bucket; j++)
122
YBuff->bucket_for_pq[pq++] = i;
124
else for (pq=pq; pq<bra_indices; pq++) YBuff->bucket_for_pq[pq] = i;
127
for (i=0,j=0; i<(nbuckets-1); i++) {
128
YBuff->buckets[i].in_bucket = 0;
129
YBuff->buckets[i].lo = j;
130
YBuff->buckets[i].hi = YBuff->buckets[i].lo + YBuff->pq_per_bucket - 1;
131
j += YBuff->pq_per_bucket;
133
YBuff->buckets[i].in_bucket = 0;
134
YBuff->buckets[i].lo = j;
135
YBuff->buckets[i].hi = bra_indices - 1;
138
/* YOSH_DONE(): Free allocated memory and reset all options for a
139
** given Yoshimine sorting structure.
141
void yosh_done(struct yoshimine *YBuff)
143
YBuff->core_loads = 0;
145
free(YBuff->bucket_for_pq);
146
YBuff->bucketsize = 0;
147
free(YBuff->buckets);
148
YBuff->first_tmp_file = 0;
149
YBuff->pq_per_bucket = 0;
150
YBuff->bra_indices = 0;
151
YBuff->ket_indices = 0;
155
void yosh_init_buckets(struct yoshimine *YBuff)
159
for (i=0; i<(YBuff->nbuckets); i++) {
160
YBuff->buckets[i].p = init_int_array(YBuff->bucketsize);
161
YBuff->buckets[i].q = init_int_array(YBuff->bucketsize);
162
YBuff->buckets[i].r = init_int_array(YBuff->bucketsize);
163
YBuff->buckets[i].s = init_int_array(YBuff->bucketsize);
164
YBuff->buckets[i].val = init_array(YBuff->bucketsize);
165
iwl_buf_init(&(YBuff->buckets[i].IWLBuf),YBuff->first_tmp_file+i,
171
** YOSH_PRINT(): Print out the Yoshimine structure
174
void yosh_print(struct yoshimine *YBuff, FILE *outfile)
176
fprintf(outfile, " Yoshimine structure:\n");
177
fprintf(outfile, "\tbra_indices = %10d\n", YBuff->bra_indices);
178
fprintf(outfile, "\tket_indices = %10d\n", YBuff->ket_indices);
179
fprintf(outfile, "\tbin size = %10d\n", YBuff->bucketsize) ;
180
fprintf(outfile, "\tbins = %10d\n", YBuff->nbuckets) ;
181
fprintf(outfile, "\tcore loads = %10d\n", YBuff->core_loads) ;
182
fprintf(outfile, "\tpq/bin = %10d\n", YBuff->pq_per_bucket) ;
183
fprintf(outfile, "\tcutoff = %10.2E\n", YBuff->cutoff) ;
188
** YOSH_CLOSE_BUCKETS(): Close the temporary binary files used for the
189
** Yoshimine sort (not the final output file).
192
** YBuff = pointer to yoshimine object
193
** erase = 1 to erase tmp files, else 0
197
** This function was formerly called yosh_close_tmp_files().
199
void yosh_close_buckets(struct yoshimine *YBuff, int erase)
203
for (i=0; i<YBuff->nbuckets; i++) { /* close but keep */
204
iwl_buf_close(&(YBuff->buckets[i].IWLBuf), !erase);
205
free(YBuff->buckets[i].p);
206
free(YBuff->buckets[i].q);
207
free(YBuff->buckets[i].r);
208
free(YBuff->buckets[i].s);
209
free(YBuff->buckets[i].val);
216
** YOSH_RDTWO() : Read two-electron integrals from
217
** file33 (in IWL format) and prepare them for Yoshimine sorting.
219
** Adapted from Ed Seidl's CSCF rdtwo.c routine
221
** Center for Computational Quantum Chemistry, UGA
224
** Modified February 1995 to use new Yoshimine data structure
225
** Modified March 1995 to use nsoff array instead of call to abs_orb()
228
** YBuff = Yoshimine object pointer
229
** itapERI = unit number for two el. file (33)
230
** num_so = array of number of symm orbs in each irrep (for reindex)
231
** nirreps = number of irreps
232
** ioff = standard lexical index array
233
** elbert = 1 for Elbert ordering, 0 for canonical ordering
234
** P = frozen core density matrix (lower triangle)
235
** Hc = frozen core operator (lower triangle)
236
** matrix = 1 for all rs for given pq, 0 otherwise
237
** (for matrix multiplication algorithm)
238
** del_tei_file = 1 to delete the tei file (33), 0 otherwise
239
** printflag = 1 for printing (for debugging only!) else 0
240
** outfile = file to print integrals to (if printflag is set)
242
void yosh_rdtwo(struct yoshimine *YBuff, int itapERI, int del_tei_file, int *num_so,
243
int nirreps, int *ioff, int elbert, int fzcflag, double *P, double *Hc,
244
int matrix, int printflag, FILE *outfile)
248
int ior, ism, jor, jsm;
249
int kor, ksm, lor, lsm;
250
int iabs, jabs, kabs, labs ;
254
struct bucket *bptr ;
256
int whichbucket, lastflag = 0, firstfile;
258
int a,b,c,d,ab,cd,ad,bc,dum,found=0;
259
int al[8], bl[8], cl[8], dl[8];
264
fprintf(outfile, "Yoshimine rdtwo routine entered\n");
265
fprintf(outfile, "Two-electron integrals from file%d:\n",itapERI);
268
firstfile = YBuff->first_tmp_file;
270
iwl_buf_init(&ERIIN,itapERI,0.0,1,1);
272
nsoff = init_int_array(nirreps);
274
for (i=1; i<nirreps; i++) {
275
nsoff[i] = nsoff[i-1] + num_so[i-1];
279
/* read a buffer full */
280
ilsti = ERIIN.lastbuf;
284
for (i=0; i < nbuf ; i++,tmp += 2) { /* do funky stuff to unpack ints */
285
iabs = abs(ERIIN.labels[fi]);
286
jabs = ERIIN.labels[fi+1];
287
kabs = ERIIN.labels[fi+2];
288
labs = ERIIN.labels[fi+3];
289
value = ERIIN.values[i];
292
/* calculate ijkl lexical index */
293
ij = ioff[iabs] + jabs;
294
kl = ioff[kabs] + labs;
295
ijkl = ioff[ij] + kl;
297
/* newly added March 1995: construct the frozen core operator */
303
ab = ioff[MAX0(a,b)] + MIN0(a,b);
304
cd = ioff[MAX0(c,d)] + MIN0(c,d);
305
bc = ioff[MAX0(b,c)] + MIN0(b,c);
306
ad = ioff[MAX0(a,d)] + MIN0(a,d);
307
Hc[cd] += 2.0 * P[ab] * value;
308
if (b >= c) Hc[bc] -= P[ad] * value;
314
if (!(a == al[0] && b == bl[0] && c == cl[0] && d == dl[0])) {
315
ab = ioff[MAX0(a,b)] + MIN0(a,b);
316
cd = ioff[MAX0(c,d)] + MIN0(c,d);
317
bc = ioff[MAX0(b,c)] + MIN0(b,c);
318
ad = ioff[MAX0(a,d)] + MIN0(a,d);
319
if (c >= d) Hc[cd] += 2.0 * P[ab] * value;
320
if (b >= c) Hc[bc] -= P[ad] * value;
327
for (dum=0, found=0; dum < 2 && !found; dum++) {
328
if (a==al[dum] && b==bl[dum] && c==cl[dum] && d==dl[dum]) found=1;
331
ab = ioff[MAX0(a,b)] + MIN0(a,b);
332
cd = ioff[MAX0(c,d)] + MIN0(c,d);
333
bc = ioff[MAX0(b,c)] + MIN0(b,c);
334
ad = ioff[MAX0(a,d)] + MIN0(a,d);
335
if (c >= d) Hc[cd] += 2.0 * P[ab] * value;
336
if (b >= c) Hc[bc] -= P[ad] * value;
343
for (dum=0, found=0; dum < 3 && !found; dum++) {
344
if(a==al[dum] && b==bl[dum] && c==cl[dum] && d==dl[dum]) found=1;
347
ab = ioff[MAX0(a,b)] + MIN0(a,b);
348
cd = ioff[MAX0(c,d)] + MIN0(c,d);
349
bc = ioff[MAX0(b,c)] + MIN0(b,c);
350
ad = ioff[MAX0(a,d)] + MIN0(a,d);
351
if (c >= d) Hc[cd] += 2.0 * P[ab] * value;
352
if (b >= c) Hc[bc] -= P[ad] * value;
359
for (dum=0, found=0; dum < 4 && !found; dum++) {
360
if(a==al[dum] && b==bl[dum] && c==cl[dum] && d==dl[dum]) found=1;
363
ab = ioff[MAX0(a,b)] + MIN0(a,b);
364
cd = ioff[MAX0(c,d)] + MIN0(c,d);
365
bc = ioff[MAX0(b,c)] + MIN0(b,c);
366
ad = ioff[MAX0(a,d)] + MIN0(a,d);
367
if (c >= d) Hc[cd] += 2.0 * P[ab] * value;
368
if (b >= c) Hc[bc] -= P[ad] * value;
375
for (dum=0, found=0; dum < 5 && !found; dum++) {
376
if (a==al[dum] && b==bl[dum] && c==cl[dum] && d==dl[dum]) found=1;
379
ab = ioff[MAX0(a,b)] + MIN0(a,b);
380
cd = ioff[MAX0(c,d)] + MIN0(c,d);
381
bc = ioff[MAX0(b,c)] + MIN0(b,c);
382
ad = ioff[MAX0(a,d)] + MIN0(a,d);
383
if (c >= d) Hc[cd] += 2.0 * P[ab] * value;
384
if (b >= c) Hc[bc] -= P[ad] * value;
391
for (dum=0, found=0; dum < 6 && !found; dum++) {
392
if (a==al[dum] && b==bl[dum] && c==cl[dum] && d==dl[dum]) found=1;
395
ab = ioff[MAX0(a,b)] + MIN0(a,b);
396
cd = ioff[MAX0(c,d)] + MIN0(c,d);
397
bc = ioff[MAX0(b,c)] + MIN0(b,c);
398
ad = ioff[MAX0(a,d)] + MIN0(a,d);
399
if (c >= d) Hc[cd] += 2.0 * P[ab] * value;
400
if (b >= c) Hc[bc] -= P[ad] * value;
407
for (dum=0, found=0; dum < 7 && !found; dum++) {
408
if(a==al[dum] && b==bl[dum] && c==cl[dum] && d==dl[dum]) found=1;
411
ab = ioff[MAX0(a,b)] + MIN0(a,b);
412
cd = ioff[MAX0(c,d)] + MIN0(c,d);
413
bc = ioff[MAX0(b,c)] + MIN0(b,c);
414
ad = ioff[MAX0(a,d)] + MIN0(a,d);
415
if (c >= d) Hc[cd] += 2.0 * P[ab] * value;
416
if (b >= c) Hc[bc] -= P[ad] * value;
418
} /* end construction of frozen core operator */
420
/* figure out what bucket to put it in, and do so
422
* Elbert wants us to sort by the lower index (kl)
423
* i.e. for us, ij > kl (guaranteed in 33), but for them kl > ij
427
if (elbert) whichbucket = YBuff->bucket_for_pq[kl] ;
428
else whichbucket = YBuff->bucket_for_pq[ij] ;
430
bptr = YBuff->buckets+whichbucket ;
431
tmpi = (bptr->in_bucket)++ ;
433
/* switch things around here for Elbert (k->p, l->q, i->r, j->s) */
435
bptr->p[tmpi] = kabs;
436
bptr->q[tmpi] = labs;
437
bptr->r[tmpi] = iabs;
438
bptr->s[tmpi] = jabs;
441
bptr->p[tmpi] = iabs;
442
bptr->q[tmpi] = jabs;
443
bptr->r[tmpi] = kabs;
444
bptr->s[tmpi] = labs;
447
bptr->val[tmpi] = value;
450
fprintf(outfile, "%4d %4d %4d %4d %4d %10.6lf\n",
451
iabs, jabs, kabs, labs, ijkl, value) ;
452
if ((tmpi+1) == YBuff->bucketsize) { /* need to flush bucket to disk */
453
flush_bucket(bptr, 0);
457
if(matrix && ij != kl) {
458
whichbucket = YBuff->bucket_for_pq[kl] ;
459
bptr = YBuff->buckets+whichbucket ;
460
tmpi = (bptr->in_bucket)++;
461
bptr->p[tmpi] = kabs;
462
bptr->q[tmpi] = labs;
463
bptr->r[tmpi] = iabs;
464
bptr->s[tmpi] = jabs;
465
bptr->val[tmpi] = value;
466
if ((tmpi+1) == YBuff->bucketsize) {
467
flush_bucket(bptr, 0);
474
iwl_buf_fetch(&ERIIN);
477
/* flush partially filled buckets */
478
/* Ok, after "matrix" was added above, we ran into the possibility of
479
* flushing TWO buffers with the lastflag set. This would be bad,
480
* because the second buffer would never be read. Therefore, I have
481
* always passed a lastflag of 0 to flush_bucket() in the code above,
482
* and now I flush all buckets here with lastflag set to 1. There
483
* is a small possibility that I will write a buffer of all zeroes.
484
* This should not actually cause a problem, the way the iwl buf reads
485
* currently work. Make sure to be careful if rewriting iwl routines!
487
for (i=0; i<YBuff->nbuckets; i++) {
488
flush_bucket((YBuff->buckets)+i, 1);
493
iwl_buf_close(&ERIIN, !del_tei_file);
497
** YOSH_RDTWO_UHF() : Read two-electron integrals from file33 (in IWL
498
** format) and prepare them for Yoshimine sorting.
501
** YBuff = Yoshimine object pointer
502
** itapERI = unit number for two el. file (33)
503
** num_so = array of number of symm orbs in each irrep (for reindex)
504
** nirreps = number of irreps
505
** ioff = standard lexical index array
506
** elbert = 1 for Elbert ordering, 0 for canonical ordering
507
** Pa = alpha frozen core density matrix (lower triangle)
508
** Pb = beta frozen core density matrix (lower triangle)
509
** Hca = alpha frozen core operator (lower triangle)
510
** Hcb = beta frozen core operator (lower triangle)
511
** matrix = 1 for all rs for given pq, 0 otherwise
512
** (for matrix multiplication algorithm)
513
** del_tei_file = 1 to delete the tei file (33), 0 otherwise
514
** printflag = 1 for printing (for debugging only!) else 0
515
** outfile = file to print integrals to (if printflag is set)
517
void yosh_rdtwo_uhf(struct yoshimine *YBuff, int itapERI, int del_tei_file, int *num_so,
518
int nirreps, int *ioff, int elbert, int fzcflag, double *Pa, double *Pb,
519
double *Hca, double *Hcb, int matrix, int printflag, FILE *outfile)
523
int ior, ism, jor, jsm;
524
int kor, ksm, lor, lsm;
525
int iabs, jabs, kabs, labs ;
529
struct bucket *bptr ;
531
int whichbucket, lastflag = 0, firstfile;
533
int a,b,c,d,ab,cd,ad,bc,dum,found=0;
534
int al[8], bl[8], cl[8], dl[8];
539
fprintf(outfile, "Yoshimine rdtwo routine entered\n");
540
fprintf(outfile, "Two-electron integrals from file%d:\n",itapERI);
543
firstfile = YBuff->first_tmp_file;
545
iwl_buf_init(&ERIIN,itapERI,0.0,1,1);
547
nsoff = init_int_array(nirreps);
549
for (i=1; i<nirreps; i++) {
550
nsoff[i] = nsoff[i-1] + num_so[i-1];
554
/* read a buffer full */
555
ilsti = ERIIN.lastbuf;
559
for (i=0; i < nbuf ; i++,tmp += 2) { /* do funky stuff to unpack ints */
560
iabs = abs(ERIIN.labels[fi]);
561
jabs = ERIIN.labels[fi+1];
562
kabs = ERIIN.labels[fi+2];
563
labs = ERIIN.labels[fi+3];
564
value = ERIIN.values[i];
567
/* calculate ijkl lexical index */
568
ij = ioff[iabs] + jabs;
569
kl = ioff[kabs] + labs;
570
ijkl = ioff[ij] + kl;
572
/* construct the UHF frozen core operator */
578
ab = ioff[MAX0(a,b)] + MIN0(a,b);
579
cd = ioff[MAX0(c,d)] + MIN0(c,d);
580
bc = ioff[MAX0(b,c)] + MIN0(b,c);
581
ad = ioff[MAX0(a,d)] + MIN0(a,d);
582
Hca[cd] += (Pa[ab] + Pb[ab]) * value;
583
Hcb[cd] += (Pa[ab] + Pb[ab]) * value;
585
Hca[bc] -= Pa[ad] * value;
586
Hcb[bc] -= Pb[ad] * value;
593
if (!(a == al[0] && b == bl[0] && c == cl[0] && d == dl[0])) {
594
ab = ioff[MAX0(a,b)] + MIN0(a,b);
595
cd = ioff[MAX0(c,d)] + MIN0(c,d);
596
bc = ioff[MAX0(b,c)] + MIN0(b,c);
597
ad = ioff[MAX0(a,d)] + MIN0(a,d);
599
Hca[cd] += (Pa[ab] + Pb[ab]) * value;
600
Hcb[cd] += (Pa[ab] + Pb[ab]) * value;
603
Hca[bc] -= Pa[ad] * value;
604
Hcb[bc] -= Pb[ad] * value;
612
for (dum=0, found=0; dum < 2 && !found; dum++) {
613
if (a==al[dum] && b==bl[dum] && c==cl[dum] && d==dl[dum]) found=1;
616
ab = ioff[MAX0(a,b)] + MIN0(a,b);
617
cd = ioff[MAX0(c,d)] + MIN0(c,d);
618
bc = ioff[MAX0(b,c)] + MIN0(b,c);
619
ad = ioff[MAX0(a,d)] + MIN0(a,d);
621
Hca[cd] += (Pa[ab] + Pb[ab]) * value;
622
Hcb[cd] += (Pa[ab] + Pb[ab]) * value;
625
Hca[bc] -= Pa[ad] * value;
626
Hcb[bc] -= Pb[ad] * value;
634
for (dum=0, found=0; dum < 3 && !found; dum++) {
635
if(a==al[dum] && b==bl[dum] && c==cl[dum] && d==dl[dum]) found=1;
638
ab = ioff[MAX0(a,b)] + MIN0(a,b);
639
cd = ioff[MAX0(c,d)] + MIN0(c,d);
640
bc = ioff[MAX0(b,c)] + MIN0(b,c);
641
ad = ioff[MAX0(a,d)] + MIN0(a,d);
643
Hca[cd] += (Pa[ab] + Pb[ab]) * value;
644
Hcb[cd] += (Pa[ab] + Pb[ab]) * value;
647
Hca[bc] -= Pa[ad] * value;
648
Hcb[bc] -= Pb[ad] * value;
656
for (dum=0, found=0; dum < 4 && !found; dum++) {
657
if(a==al[dum] && b==bl[dum] && c==cl[dum] && d==dl[dum]) found=1;
660
ab = ioff[MAX0(a,b)] + MIN0(a,b);
661
cd = ioff[MAX0(c,d)] + MIN0(c,d);
662
bc = ioff[MAX0(b,c)] + MIN0(b,c);
663
ad = ioff[MAX0(a,d)] + MIN0(a,d);
665
Hca[cd] += (Pa[ab] + Pb[ab]) * value;
666
Hcb[cd] += (Pa[ab] + Pb[ab]) * value;
669
Hca[bc] -= Pa[ad] * value;
670
Hcb[bc] -= Pb[ad] * value;
678
for (dum=0, found=0; dum < 5 && !found; dum++) {
679
if (a==al[dum] && b==bl[dum] && c==cl[dum] && d==dl[dum]) found=1;
682
ab = ioff[MAX0(a,b)] + MIN0(a,b);
683
cd = ioff[MAX0(c,d)] + MIN0(c,d);
684
bc = ioff[MAX0(b,c)] + MIN0(b,c);
685
ad = ioff[MAX0(a,d)] + MIN0(a,d);
687
Hca[cd] += (Pa[ab] + Pb[ab]) * value;
688
Hcb[cd] += (Pa[ab] + Pb[ab]) * value;
691
Hca[bc] -= Pa[ad] * value;
692
Hcb[bc] -= Pb[ad] * value;
700
for (dum=0, found=0; dum < 6 && !found; dum++) {
701
if (a==al[dum] && b==bl[dum] && c==cl[dum] && d==dl[dum]) found=1;
704
ab = ioff[MAX0(a,b)] + MIN0(a,b);
705
cd = ioff[MAX0(c,d)] + MIN0(c,d);
706
bc = ioff[MAX0(b,c)] + MIN0(b,c);
707
ad = ioff[MAX0(a,d)] + MIN0(a,d);
709
Hca[cd] += (Pa[ab] + Pb[ab]) * value;
710
Hcb[cd] += (Pa[ab] + Pb[ab]) * value;
713
Hca[bc] -= Pa[ad] * value;
714
Hcb[bc] -= Pb[ad] * value;
722
for (dum=0, found=0; dum < 7 && !found; dum++) {
723
if(a==al[dum] && b==bl[dum] && c==cl[dum] && d==dl[dum]) found=1;
726
ab = ioff[MAX0(a,b)] + MIN0(a,b);
727
cd = ioff[MAX0(c,d)] + MIN0(c,d);
728
bc = ioff[MAX0(b,c)] + MIN0(b,c);
729
ad = ioff[MAX0(a,d)] + MIN0(a,d);
731
Hca[cd] += (Pa[ab] + Pb[ab]) * value;
732
Hcb[cd] += (Pa[ab] + Pb[ab]) * value;
735
Hca[bc] -= Pa[ad] * value;
736
Hcb[bc] -= Pb[ad] * value;
739
} /* end construction of frozen core operator */
741
/* figure out what bucket to put it in, and do so
743
* Elbert wants us to sort by the lower index (kl)
744
* i.e. for us, ij > kl (guaranteed in 33), but for them kl > ij
748
if (elbert) whichbucket = YBuff->bucket_for_pq[kl] ;
749
else whichbucket = YBuff->bucket_for_pq[ij] ;
751
bptr = YBuff->buckets+whichbucket ;
752
tmpi = (bptr->in_bucket)++ ;
754
/* switch things around here for Elbert (k->p, l->q, i->r, j->s) */
756
bptr->p[tmpi] = kabs;
757
bptr->q[tmpi] = labs;
758
bptr->r[tmpi] = iabs;
759
bptr->s[tmpi] = jabs;
762
bptr->p[tmpi] = iabs;
763
bptr->q[tmpi] = jabs;
764
bptr->r[tmpi] = kabs;
765
bptr->s[tmpi] = labs;
768
bptr->val[tmpi] = value;
771
fprintf(outfile, "%4d %4d %4d %4d %4d %10.6lf\n",
772
iabs, jabs, kabs, labs, ijkl, value) ;
773
if ((tmpi+1) == YBuff->bucketsize) { /* need to flush bucket to disk */
774
flush_bucket(bptr, 0);
778
if(matrix && ij != kl) {
779
whichbucket = YBuff->bucket_for_pq[kl] ;
780
bptr = YBuff->buckets+whichbucket ;
781
tmpi = (bptr->in_bucket)++;
782
bptr->p[tmpi] = kabs;
783
bptr->q[tmpi] = labs;
784
bptr->r[tmpi] = iabs;
785
bptr->s[tmpi] = jabs;
786
bptr->val[tmpi] = value;
787
if ((tmpi+1) == YBuff->bucketsize) {
788
flush_bucket(bptr, 0);
795
iwl_buf_fetch(&ERIIN);
798
/* flush partially filled buckets */
799
/* Ok, after "matrix" was added above, we ran into the possibility of
800
* flushing TWO buffers with the lastflag set. This would be bad,
801
* because the second buffer would never be read. Therefore, I have
802
* always passed a lastflag of 0 to flush_bucket() in the code above,
803
* and now I flush all buckets here with lastflag set to 1. There
804
* is a small possibility that I will write a buffer of all zeroes.
805
* This should not actually cause a problem, the way the iwl buf reads
806
* currently work. Make sure to be careful if rewriting iwl routines!
808
for (i=0; i<YBuff->nbuckets; i++) {
809
flush_bucket((YBuff->buckets)+i, 1);
813
iwl_buf_close(&ERIIN, !del_tei_file);
817
** YOSH_RDTWO_BACKTR() : Read two-electron integrals from an IWL file and
818
** prepare them for Yoshimine sorting. We have removed support for
819
** Elbert loops and frozen core, since the former is not currently
820
** being used and the latter should not apply to backtransforms.
821
** The main issue is whether we need to symmetrize the twopdm, or
822
** whether it has already been symmetrized. We assume that it always
823
** has the symmetry (pq|rs) = (rs|pq), but it may not have the other
824
** left-pair and right-pair permutational symmetries (the CI twopdm
825
** for example does not have this symmetry naturally). The transform
826
** needs this symmetry (as does the derivative program), so we will
827
** enforce it here if the user specifies. We currently assume that
828
** the tpdm on disk has only unique (pq|rs) pairs, and so for our
829
** purposes we need to generate (rs|pq) just like we do when reading
830
** the AO integrals in the regular forwards transform.
832
** Based on the YOSH_RDTWO34() function
834
** Created August 1997
837
** YBuff = Yoshimine object pointer
838
** tei_file = unit number for two-electron integrals
839
** ioff = standard lexical index array
840
** symmetrize = symmetrize the incoming 2pdm
841
** add_ref_pt = Add the factors arising from a reference determinant
842
** (n.b. assumes lowest MO's occupied)
843
** del_tei_file = 1 to delete the tei file, 0 otherwise
844
** printflag = 1 for printing (for debugging only!) else 0
845
** outfile = file to print integrals to (if printflag is set)
847
void yosh_rdtwo_backtr(struct yoshimine *YBuff, int tei_file, int *ioff,
848
int symmetrize, int add_ref_pt, int del_tei_file,
849
int printflag, FILE *outfile)
852
int iabs, jabs, kabs, labs;
855
int whichbucket, lastbuf, idx;
862
fprintf(outfile, "Yoshimine rdtwo_backtr routine entered\n");
863
fprintf(outfile, "Two-particle density from file %d:\n", tei_file);
866
iwl_buf_init(&Inbuf, tei_file, YBuff->cutoff, 1, 0);
867
lblptr = Inbuf.labels;
868
valptr = Inbuf.values;
871
iwl_buf_fetch(&Inbuf);
872
lastbuf = Inbuf.lastbuf;
873
for (idx=4*Inbuf.idx; Inbuf.idx < Inbuf.inbuf; Inbuf.idx++) {
874
iabs = (int) lblptr[idx++];
875
jabs = (int) lblptr[idx++];
876
kabs = (int) lblptr[idx++];
877
labs = (int) lblptr[idx++];
879
iabs = moinfo.corr2pitz_nofzv[iabs];
880
jabs = moinfo.corr2pitz_nofzv[jabs];
881
kabs = moinfo.corr2pitz_nofzv[kabs];
882
labs = moinfo.corr2pitz_nofzv[labs];
884
value = valptr[Inbuf.idx];
886
if (iabs != jabs) value *= 0.5;
887
if (kabs != labs) value *= 0.5;
890
/* calculate ijkl lexical index (make no i>=j assumptions for now) */
891
ij = INDEX(iabs,jabs);
892
kl = INDEX(kabs,labs);
893
ijkl = INDEX(ij,kl); /* ijkl needed only in the print routine */
895
/* figure out what bucket to put it in, and do so */
897
whichbucket = YBuff->bucket_for_pq[ij];
898
bptr = YBuff->buckets+whichbucket;
899
tmpi = (bptr->in_bucket)++;
900
bptr->p[tmpi] = iabs;
901
bptr->q[tmpi] = jabs;
902
bptr->r[tmpi] = kabs;
903
bptr->s[tmpi] = labs;
904
bptr->val[tmpi] = value;
907
fprintf(outfile, "%4d %4d %4d %4d %4d %10.6lf\n",
908
iabs, jabs, kabs, labs, ijkl, value) ;
909
if ((tmpi+1) == YBuff->bucketsize) { /* need to flush bucket to disk */
910
flush_bucket(bptr, 0);
915
/* this generates (kl|ij) from (ij|kl) if necessary and puts it out */
916
/* anaglogous to the "matrix" option in yosh_rdtwo() */
917
if (iabs != kabs || jabs != labs) {
918
whichbucket = YBuff->bucket_for_pq[kl];
919
bptr = YBuff->buckets+whichbucket;
920
tmpi = (bptr->in_bucket)++;
921
bptr->p[tmpi] = kabs;
922
bptr->q[tmpi] = labs;
923
bptr->r[tmpi] = iabs;
924
bptr->s[tmpi] = jabs;
925
bptr->val[tmpi] = value;
927
fprintf(outfile, "%4d %4d %4d %4d %4d %10.6lf\n",
928
kabs, labs, iabs, jabs, ijkl, value) ;
929
if ((tmpi+1) == YBuff->bucketsize) {
930
flush_bucket(bptr, 0);
939
/* now add in contributions from the reference determinant if requested */
940
if (add_ref_pt) add_2pdm_ref_pt(YBuff, ioff, printflag, outfile);
942
/* flush partially filled buckets */
943
for (i=0; i<YBuff->nbuckets; i++) {
944
flush_bucket((YBuff->buckets)+i, 1);
947
iwl_buf_close(&Inbuf, !del_tei_file);
951
** YOSH_RDTWO_BACKTR_UHF() : Read two-particle density elements from
952
** an IWL file and prepare them for Yoshimine sorting. The sorted
953
** twopdm elements, G(pqrs), produced by this code have unique p-q and
954
** r-s combinations, but no pq-rs packing. However, the input
955
** twopdm's may not have this same structure, so the boolean arguments
956
** swap_bk and symm_pq are used to correct this problem. There are
957
** two circumstances to consider:
959
** (1) If the MO twopdm lacks pq-rs symmetry (e.g., the AB twopdm),
960
** its input file should include all rs for each pq. The swap_bk flag
961
** should be set to "0" in this case so that no "extra" G(rs,pq)
962
** components are written to the sorting buffers. If the input twopdm
963
** has pq-rs symmetry (e.g., the AA and BB twopdms), only unique pq-rs
964
** combinations should be included and the swap_bk flag should be set
967
** (2) If the MO density lacks p-q and r-s symmetry, the input file
968
** should include all combinations of p,q and r,s, and the symm_pq
969
** flag should be set to "1". If the MO density has p-q and r-s
970
** symmetry, then its input file should include only unique p,q and
971
** r,s combinations, and the symm_pq flag should be set to "0".
973
** Note that "intermediate" symmetry cases, where the MO twopdm has
974
** p-q symmetry but not r-s symmetry, for example, are not included
977
** Also, this code assumes the input indices are in QTS ordering and
978
** converts them automatically to Pitzer.
982
** Based on the YOSH_RDTWO_BACKTR() function above by
986
** YBuff = Yoshimine object pointer
987
** tei_file = unit number for two-electron integrals
988
** ioff = standard lexical index array
989
** swap_bk = sort both G(pq,rs) and G(rs,pq) combinations
990
** symm_pq = symmetrize both p,q and r,s combinations
991
** del_tei_file = 1 to delete the tei file, 0 otherwise
992
** printflag = 1 for printing (for debugging only!) else 0
993
** outfile = file to print integrals to (if printflag is set)
995
void yosh_rdtwo_backtr_uhf(const char *spin, struct yoshimine *YBuff, int tei_file, int *ioff,
996
int swap_bk, int symm_pq, int del_tei_file,
997
int printflag, FILE *outfile)
1000
int iabs, jabs, kabs, labs;
1002
struct bucket *bptr;
1003
int whichbucket, lastbuf, idx;
1005
struct iwlbuf Inbuf;
1008
int *iorder, *jorder, *korder, *lorder;
1011
fprintf(outfile, "Yoshimine rdtwo_backtr routine entered\n");
1012
fprintf(outfile, "Two-particle density from file %d:\n", tei_file);
1015
if(!strcmp(spin, "AA")) {
1016
iorder = moinfo.corr2pitz_nofzv_a;
1017
jorder = moinfo.corr2pitz_nofzv_a;
1018
korder = moinfo.corr2pitz_nofzv_a;
1019
lorder = moinfo.corr2pitz_nofzv_a;
1021
else if(!strcmp(spin, "BB")) {
1022
iorder = moinfo.corr2pitz_nofzv_b;
1023
jorder = moinfo.corr2pitz_nofzv_b;
1024
korder = moinfo.corr2pitz_nofzv_b;
1025
lorder = moinfo.corr2pitz_nofzv_b;
1027
else if(!strcmp(spin, "AB")) {
1028
iorder = moinfo.corr2pitz_nofzv_a;
1029
jorder = moinfo.corr2pitz_nofzv_a;
1030
korder = moinfo.corr2pitz_nofzv_b;
1031
lorder = moinfo.corr2pitz_nofzv_b;
1034
fprintf(outfile, "\n\tInvalid spin cases requested for backtransformation!\n");
1035
exit(PSI_RETURN_FAILURE);
1038
iwl_buf_init(&Inbuf, tei_file, YBuff->cutoff, 1, 0);
1039
lblptr = Inbuf.labels;
1040
valptr = Inbuf.values;
1043
iwl_buf_fetch(&Inbuf);
1044
lastbuf = Inbuf.lastbuf;
1045
for (idx=4*Inbuf.idx; Inbuf.idx < Inbuf.inbuf; Inbuf.idx++) {
1046
iabs = (int) lblptr[idx++];
1047
jabs = (int) lblptr[idx++];
1048
kabs = (int) lblptr[idx++];
1049
labs = (int) lblptr[idx++];
1051
iabs = iorder[iabs];
1052
jabs = jorder[jabs];
1053
kabs = korder[kabs];
1054
labs = lorder[labs];
1056
value = valptr[Inbuf.idx];
1059
if (iabs != jabs) value *= 0.5;
1060
if (kabs != labs) value *= 0.5;
1063
/* calculate ijkl lexical index (make no i>=j assumptions for now) */
1064
ij = INDEX(iabs,jabs);
1065
kl = INDEX(kabs,labs);
1066
ijkl = INDEX(ij,kl); /* ijkl needed only in the print routine */
1068
/* figure out what bucket to put it in, and do so */
1070
whichbucket = YBuff->bucket_for_pq[ij];
1071
bptr = YBuff->buckets+whichbucket;
1072
tmpi = (bptr->in_bucket)++;
1073
bptr->p[tmpi] = iabs;
1074
bptr->q[tmpi] = jabs;
1075
bptr->r[tmpi] = kabs;
1076
bptr->s[tmpi] = labs;
1077
bptr->val[tmpi] = value;
1080
fprintf(outfile, "%4d %4d %4d %4d %4d %10.6lf\n",
1081
iabs, jabs, kabs, labs, ijkl, value) ;
1082
if ((tmpi+1) == YBuff->bucketsize) { /* need to flush bucket to disk */
1083
flush_bucket(bptr, 0);
1084
bptr->in_bucket = 0;
1087
/* this generates (kl|ij) from (ij|kl) if necessary and puts it out */
1088
/* anaglogous to the "matrix" option in yosh_rdtwo() */
1089
if (swap_bk && (iabs != kabs || jabs != labs)) {
1090
whichbucket = YBuff->bucket_for_pq[kl];
1091
bptr = YBuff->buckets+whichbucket;
1092
tmpi = (bptr->in_bucket)++;
1093
bptr->p[tmpi] = kabs;
1094
bptr->q[tmpi] = labs;
1095
bptr->r[tmpi] = iabs;
1096
bptr->s[tmpi] = jabs;
1097
bptr->val[tmpi] = value;
1099
fprintf(outfile, "%4d %4d %4d %4d %4d %10.6lf\n",
1100
kabs, labs, iabs, jabs, ijkl, value) ;
1101
if ((tmpi+1) == YBuff->bucketsize) {
1102
flush_bucket(bptr, 0);
1103
bptr->in_bucket = 0;
1111
/* flush partially filled buckets */
1112
for (i=0; i<YBuff->nbuckets; i++) {
1113
flush_bucket((YBuff->buckets)+i, 1);
1116
iwl_buf_close(&Inbuf, !del_tei_file);
1122
** This function adds in the contributions to the two-particle density
1123
** matrix from the reference determinant, as might be required in MBPT
1124
** or CC theory. Assume the reference is made up of the lowest-lying
1125
** orbitals as specified by DOCC and SOCC arrays. Assume restricted
1126
** orbitals. Assume correlated ordering is docc, socc, virt...may
1127
** not always be true, but order array is already wiped out by this
1128
** point for backtransforms, would have to fetch it again.
1130
** David Sherrill, Feb 1998
1132
void add_2pdm_ref_pt(struct yoshimine *YBuff, int *ioff, int pflg,
1135
int i, j, ii, jj, ij;
1139
ndocc = moinfo.ndocc; nsocc = moinfo.nsocc;
1141
/* closed-shell part */
1142
for (i=0; i<ndocc; i++) {
1143
iabs = moinfo.corr2pitz_nofzv[i];
1144
ii = ioff[iabs] + iabs;
1146
for (j=0; j<i; j++) {
1147
jabs = moinfo.corr2pitz_nofzv[j];
1148
jj = ioff[jabs] + jabs;
1149
ij = ioff[iabs] + jabs;
1151
yosh_buff_put_val(YBuff,ioff,ii,iabs,iabs,jabs,jabs, 2.0,pflg,outfile);
1152
yosh_buff_put_val(YBuff,ioff,jj,jabs,jabs,iabs,iabs, 2.0,pflg,outfile);
1153
yosh_buff_put_val(YBuff,ioff,ij,iabs,jabs,jabs,iabs,-0.25,pflg,outfile);
1154
yosh_buff_put_val(YBuff,ioff,ij,jabs,iabs,iabs,jabs,-0.25,pflg,outfile);
1158
jabs = moinfo.corr2pitz_nofzv[j];
1159
jj = ioff[jabs] + jabs;
1160
ij = ioff[iabs] + jabs;
1161
yosh_buff_put_val(YBuff,ioff,ii,iabs,iabs,jabs,jabs, 1.0,pflg,outfile);
1166
/* open-shell part, if any */
1167
for (i=ndocc; i<ndocc+nsocc; i++) {
1168
iabs = moinfo.corr2pitz_nofzv[i];
1169
ii = ioff[iabs] + iabs;
1171
for (j=0; j<ndocc; j++) {
1172
jabs = moinfo.corr2pitz_nofzv[j];
1173
jj = ioff[jabs] + jabs;
1174
ij = ioff[iabs] + jabs;
1176
yosh_buff_put_val(YBuff,ioff,ii,iabs,iabs,jabs,jabs, 1.0,pflg,outfile);
1177
yosh_buff_put_val(YBuff,ioff,jj,jabs,jabs,iabs,iabs, 1.0,pflg,outfile);
1178
yosh_buff_put_val(YBuff,ioff,ij,iabs,jabs,jabs,iabs,-0.125,pflg,outfile);
1179
yosh_buff_put_val(YBuff,ioff,ij,jabs,iabs,iabs,jabs,-0.125,pflg,outfile);
1182
for (j=ndocc; j<i; j++) {
1183
jabs = moinfo.corr2pitz_nofzv[j];
1184
jj = ioff[jabs] + jabs;
1185
ij = ioff[iabs] + jabs;
1187
yosh_buff_put_val(YBuff,ioff,ii,iabs,iabs,jabs,jabs, 0.5,pflg,outfile);
1188
yosh_buff_put_val(YBuff,ioff,jj,jabs,jabs,iabs,iabs, 0.5,pflg,outfile);
1189
yosh_buff_put_val(YBuff,ioff,ij,iabs,jabs,jabs,iabs,-0.125,pflg,outfile);
1190
yosh_buff_put_val(YBuff,ioff,ij,jabs,iabs,iabs,jabs,-0.125,pflg,outfile);
1200
** YOSH_BUFF_PUT_VAL
1202
** This function puts a value and its associated indices to a yoshimine
1205
void yosh_buff_put_val(struct yoshimine *YBuff, int *ioff, int pq,
1206
int p, int q, int r, int s, double value, int prtflg,
1209
struct bucket *bptr;
1213
whichbucket = YBuff->bucket_for_pq[pq];
1214
bptr = YBuff->buckets+whichbucket;
1215
tmpi = (bptr->in_bucket)++;
1220
bptr->val[tmpi] = value;
1223
fprintf(outfile, "%4d %4d %4d %4d %10.6lf\n", p, q, r, s,
1226
if ((tmpi+1) == YBuff->bucketsize) { /* need to flush bucket to disk */
1227
flush_bucket(bptr, 0);
1228
bptr->in_bucket = 0;
1237
** YOSH_SORT(): Sort all the buckets in the Yoshimine sorting algorithm.
1238
** The call to sortbuf() will cause the intermediate files to be
1239
** deleted unless keep_bins is set to 1.
1242
** YBuff = pointer to Yoshimine object
1243
** out_tape = number for binary output file
1244
** keep_bins = keep the intermediate tmp files
1245
** ioff = the usual offset array unless no_pq_perm, in which
1246
** case it is appropriate for the left indices, i.e.,
1248
** ioff2 = the Elbert ioff2 array if elbert=true. If no_pq_perm,
1249
** then this is the usual ioff array, appropriate
1250
** for the right indices
1251
** nbfso = number of basis fns in symmetry orbitals
1252
** ket_indices = number of ket indices (usually ntri)
1253
** elbert = are inputs in Elbert ordering? p>=q, r>=s, rs>=pq
1254
** intermediate = 1 if sorting an intermediate in the transformation
1255
** (argument to sortbuf()). This implies that the full
1256
** set of rs are available for a given pq.
1257
** no_pq_perm = if p and q are not interchangeable (e.g., one is
1258
** occupied and one is virtual, as in MP2)
1259
** qdim = the number of possible values for index q
1260
** add = do additions of integrals during the sort?
1261
** print_lvl = verbosity level (how much to print)
1262
** outfile = text output file
1266
void yosh_sort(struct yoshimine *YBuff, int out_tape, int keep_bins,
1267
int *ioff, int *ioff2, int nbfso, int ket_indices, int elbert,
1268
int intermediate, int no_pq_perm, int qdim, int add, int print_lvl,
1273
struct iwlbuf inbuf, outbuf;
1275
/* may be slightly more than pq_per_bucket pq's in each bucket
1276
* if the pq's didn't divide evenly among the buckets. The remainder
1277
* will go to the last bucket.
1279
max_pq = YBuff->buckets[YBuff->core_loads-1].hi -
1280
YBuff->buckets[YBuff->core_loads-1].lo
1283
twoel_ints = init_array(max_pq * ket_indices) ;
1284
iwl_buf_init(&outbuf, out_tape, YBuff->cutoff, 0, 0);
1286
for (i=0; i<YBuff->core_loads-1; i++) {
1287
if (print_lvl > 1) fprintf(outfile, "Sorting bin %d\n", i+1);
1288
iwl_buf_init(&inbuf, YBuff->first_tmp_file+i, YBuff->cutoff, 1, 0);
1289
sortbuf(&inbuf, &outbuf, twoel_ints, (YBuff->buckets)[i].lo,
1290
(YBuff->buckets)[i].hi, ioff, ioff2, nbfso, elbert,
1291
intermediate, no_pq_perm, qdim, add, (print_lvl > 4), outfile);
1292
zero_arr(twoel_ints, max_pq * ket_indices);
1293
/* zero_arr(twoel_ints, YBuff->pq_per_bucket * YBuff->bra_indices); */
1294
iwl_buf_close(&inbuf, keep_bins);
1298
if (print_lvl > 1) fprintf(outfile, "Sorting bin %d\n", i+1) ;
1299
iwl_buf_init(&inbuf, YBuff->first_tmp_file+i, YBuff->cutoff, 1, 0);
1300
sortbuf(&inbuf, &outbuf, twoel_ints, (YBuff->buckets)[i].lo,
1301
(YBuff->buckets)[i].hi, ioff, ioff2, nbfso, elbert,
1302
intermediate, no_pq_perm, qdim, add, (print_lvl > 4), outfile);
1303
iwl_buf_close(&inbuf, keep_bins);
1305
if (print_lvl > 1) fprintf(outfile, "Done sorting.\n");
1307
iwl_buf_flush(&outbuf, 1);
1308
iwl_buf_close(&outbuf, 1);
1314
** YOSH_FREE(): Free up a Yoshimine object. Free any dynamically-allocated
1317
void yosh_free(struct yoshimine *YBuff)
1319
free(YBuff->buckets);
1324
** YOSH_FLUSH(): Flush any of the buckets and tag them as 'last buffer'
1325
** for each bucket tmp file.
1328
void yosh_flush(struct yoshimine *YBuff)
1332
for (i=0; i<YBuff->nbuckets; i++) {
1333
flush_bucket((YBuff->buckets)+i, 1);
1340
** FLUSH_BUCKET(): This function flushes a Yoshimine bucket to a temporary
1341
** binary file. Must be careful to call with lastbuf==1 when the
1342
** last buffer has been reached, so that the iwl buffer can set the
1343
** `last buffer' flag on output.
1345
void flush_bucket(struct bucket *bptr, int lastbuf)
1348
iwl_buf_wrt_arr(&(bptr->IWLBuf), bptr->val, bptr->p, bptr->q,
1349
bptr->r, bptr->s, bptr->in_bucket);
1350
iwl_buf_flush(&(bptr->IWLBuf), lastbuf);
1355
** YOSH_WRT_ARR(): Write an array to a Yoshimine object...useful for writing
1356
** intermediates to disk from a transformation program. Write only
1357
** nonzero elements so that less disk space is required, and so that
1358
** zeroes for i,j,k,l and value for the first element in a buffer can
1359
** indicate an empty buffer. (Eventually I should probably use
1360
** the second byte of the flags field in each buffer to indicate a
1361
** zero buffer, else use the second and third bytes to denote the
1362
** number of integrals in each buffer. File34 had a field like
1363
** this, but it assumes an integral number of integer words per
1364
** double, which might not be true in the long run. IWL is currently
1365
** free of this assumption).
1367
** The sortbuf() routine wants to have blocks containing all rs for
1368
** a given pq. If instead, we want all pq for a given rs, we
1369
** need to reverse the order of pq and rs in the tmp files. Set
1373
** YBuff = pointer to Yoshimine object
1374
** p = common p value for array
1375
** q = common q value for array
1376
** pq = compound index determined from p and q
1377
** pqsym = the direct product of the symmetries of p and q
1378
** arr = the array containing the data for a given pq
1379
** rmax = the maximum value of the third index
1380
** ioff = the standard offset array
1381
** orbsym = the irrep for each orbital
1382
** firsti = the first orbital for each irrep
1383
** lasti = last orbital for each irrep
1384
** sortby_rs = described above
1385
** printflag = verbosity level for printing
1386
** outfile = text output file
1389
void yosh_wrt_arr(struct yoshimine *YBuff, int p, int q, int pq, int pqsym,
1390
double *arr, int rmax, int *ioff, int *orbsym, int *firsti,
1391
int *lasti, int sortby_rs, int printflag, FILE *outfile)
1393
int r, s, rs, rsym, ssym, smax;
1396
int lastflag = 0, firstfile;
1397
struct bucket *bptr;
1401
fprintf(outfile, "\nyosh_wrt_arr called for p=%d,q=%d\n", p, q);
1404
firstfile = YBuff->first_tmp_file;
1406
for (r=0; r<rmax; r++) {
1408
ssym = pqsym ^ rsym;
1409
if (ssym > rsym) continue;
1411
smax = (rsym == ssym) ? r : lasti[ssym];
1413
for (s=firsti[ssym]; s<=smax; s++) {
1417
if (fabs(value) > YBuff->cutoff) {
1418
/* figure out which bucket the integral should go in */
1419
if (sortby_rs) whichbucket = YBuff->bucket_for_pq[rs];
1420
else whichbucket = YBuff->bucket_for_pq[pq];
1422
bptr = YBuff->buckets+whichbucket ;
1423
tmpi = (bptr->in_bucket)++ ;
1425
/* sortbuf wants to sort by the first index. If we want to
1426
* sort by rs, we need to map r->p, s->q, p->r, q->s
1441
bptr->val[tmpi] = value;
1444
fprintf(outfile, "%4d %4d %4d %4d %10.6lf\n",
1445
p, q, r, s, arr[rs]);
1447
/* if we need to flush bucket to disk */
1448
if ((tmpi+1) == YBuff->bucketsize) {
1449
flush_bucket(bptr, lastflag);
1450
bptr->in_bucket = 0 ;
1454
* What if this was the last buffer, and we didn't set the lastbuf
1455
* flag? Must write a buffer of zeroes using the yosh_flush()
1456
* routine. Must be careful to always have the last buffer
1457
* set the lastbuf flag (to avoid read error) and must also be
1458
* sure that an exactly-filled buffer is written out. Carelessness
1459
* might cause a failure to write the last buffer if it has
1460
* exactly 'bucketsize' elements. ---CDS
1463
} /* end if (fabs(arr[rs])) ... */
1465
} /* end loop over s */
1467
} /* end loop over r */
1469
if (printflag) fflush(outfile);
1475
** YOSH_WRT_ARR2(): Write an array to a Yoshimine object (nonzero
1476
** elements only). The values of the integrals are in array 'arr',
1477
** and the indices are input as follows: all integrals have the
1478
** common indices 'p' and 'q'. The remaining indices are input
1479
** in arrays 'rlist' and 'slist'. The number of integrals in the
1480
** current block is given by the parameter 'size'. As currently
1481
** written, the routine will write out the integrals in canonical
1482
** ordering p>=q, r>=s, pq>=rs. It might conceivably be useful
1483
** in the future to use a 'sortby_rs' flag as in yosh_wrt_arr().
1486
** YBuff = pointer to Yoshimine object
1487
** size = number of integrals to process in the array
1488
** arr = the array containing the data for a given pq
1489
** p = common p value for array
1490
** q = common q value for array
1491
** rlist = list of indices r
1492
** slist = list of indices s (no guarantee rlist[i] >= slist[i])
1493
** ioff = the standard offset array
1494
** printflag = verbosity level for printing
1495
** outfile = text output file
1498
void yosh_wrt_arr2(struct yoshimine *YBuff, int size, double *arr,
1499
int p, int q, int *rlist, int *slist, int *ioff, int printflag,
1503
int i1, j1, k1, l1, ij1, kl1;
1505
int i2, j2, k2, l2, ij2, kl2;
1508
int lastflag = 0, firstfile;
1509
struct bucket *bptr;
1512
firstfile = YBuff->first_tmp_file;
1516
ij1 = ioff[i1] + j1;
1518
for (x=0; x<size; x++) {
1524
if (fabs(value) < YBuff->cutoff) continue;
1526
k1 = MAX0(ktmp,ltmp);
1527
l1 = MIN0(ktmp,ltmp);
1528
kl1 = ioff[k1] + l1;
1547
whichbucket = YBuff->bucket_for_pq[ij2];
1548
bptr = YBuff->buckets+whichbucket ;
1549
tmpi = (bptr->in_bucket)++ ;
1556
bptr->val[tmpi] = value;
1558
if (printflag) fprintf(outfile, "%4d %4d %4d %4d %10.6lf\n",
1559
i2, j2, k2, l2, value);
1561
/* if we need to flush bucket to disk */
1562
if ((tmpi+1) == YBuff->bucketsize) {
1563
flush_bucket(bptr, lastflag);
1564
bptr->in_bucket = 0;
1567
} /* end loop over x */
1572
** YOSH_WRT_ARR_MP2(): Write an array to a Yoshimine object...useful for
1573
** writing intermediates to disk from a transformation program. Write
1574
** only nonzero elements so that less disk space is required, and so that
1575
** zeroes for i,j,k,l and value for the first element in a buffer can
1576
** indicate an empty buffer.
1578
** The sortbuf() routine wants to have blocks containing all rs for
1579
** a given pq. If instead, we want all pq for a given rs, we
1580
** need to reverse the order of pq and rs in the tmp files. Set
1583
** This routine has been altered to work _only_ for MP2 restricted
1584
** sorts. It expects the r-index to correlate to occupied orbitals and
1585
** the s-index to virtual orbitals. This routine may not, in general
1586
** be used with other orbital distributions.
1590
** YBuff = pointer to Yoshimine object
1591
** p = common p value for array
1592
** q = common q value for array
1593
** pq = compound index determined from p and q
1594
** pqsym = the direct product of the symmetries of p and q
1595
** arr = the array containing the data for a given pq
1596
** rsym = the irrep of the r-indices
1597
** firstr = the first r-index for each irrep
1598
** lastr = the last r-index for each irrep
1599
** firsts = the first s-index for each irrep
1600
** lasts = the last s-index for each irrep
1601
** sortby_rs = described above
1602
** ndocc = the number of doubly-occupied orbitals
1603
** nvirt = the number of virtual orbitals
1604
** occ = the Pitzer -> QTS ordering array for the occupied orbitals
1605
** vir = the Pitzer -> QTS ordering array for the virtual orbitals
1606
** ioff3 = offset array for ndocc*nvirt arrays
1607
** printflag = verbosity level for printing
1608
** outfile = text output file
1611
void yosh_wrt_arr_mp2(struct yoshimine *YBuff, int p, int q, int pq,
1612
int pqsym, double **arr, int rsym, int *firstr,
1613
int *lastr, int *firsts, int *lasts, int sortby_rs,
1614
int ndocc, int nvirt, int *occ, int *vir, int *ioff3,
1615
int printflag, FILE *outfile)
1622
int lastflag = 0, firstfile;
1623
struct bucket *bptr;
1626
firstfile = YBuff->first_tmp_file;
1627
ssym = pqsym ^ rsym;
1628
for (r=firstr[rsym], R=0; r <= lastr[rsym]; r++, R++) {
1630
for (s=firsts[ssym], S=0; s <=lasts[ssym]; s++, S++) {
1632
rs = ioff3[rnew] + snew;
1635
if (fabs(value) > YBuff->cutoff) {
1636
/* figure out which bucket the integral should go in */
1637
if (sortby_rs) whichbucket = YBuff->bucket_for_pq[rs];
1638
else whichbucket = YBuff->bucket_for_pq[pq];
1640
bptr = YBuff->buckets+whichbucket ;
1641
tmpi = (bptr->in_bucket)++ ;
1643
/* sortbuf wants to sort by the first index. If we want to
1644
* sort by rs, we need to map r->p, s->q, p->r, q->s
1645
* We also write out the QTS ordered indices for r and s, and
1646
* not the Pitzer indices.
1649
bptr->p[tmpi] = rnew;
1650
bptr->q[tmpi] = snew;
1657
bptr->r[tmpi] = rnew;
1658
bptr->s[tmpi] = snew;
1661
bptr->val[tmpi] = value;
1664
fprintf(outfile, "%4d %4d %4d %4d %10.6lf\n",
1667
/* if we need to flush bucket to disk */
1668
if ((tmpi+1) == YBuff->bucketsize) {
1669
flush_bucket(bptr, lastflag);
1670
bptr->in_bucket = 0;
1673
} /* end if (fabs(arr[rs])) ... */
1675
} /* end loop over s */
1677
} /* end loop over r */
1684
** YOSH_WRT_ARR_MP2R12A(): Write an array to a Yoshimine object...useful for
1685
** writing intermediates to disk from a transformation program. Write
1686
** only nonzero elements so that less disk space is required, and so that
1687
** zeroes for i,j,k,l and value for the first element in a buffer can
1688
** indicate an empty buffer.
1690
** The sortbuf() routine wants to have blocks containing all rs for
1691
** a given pq. If instead, we want all pq for a given rs, we
1692
** need to reverse the order of pq and rs in the tmp files. Set
1695
** This routine has been altered to work _only_ for MP2R12A restricted
1696
** sorts. It expects the r-index to correlate to occupied orbitals.
1697
** This routine may not, in general
1698
** be used with other orbital distributions.
1702
** YBuff = pointer to Yoshimine object
1703
** p = common p value for array
1704
** q = common q value for array
1705
** pq = compound index determined from p and q
1706
** pqsym = the direct product of the symmetries of p and q
1707
** arr = the array containing the data for a given pq
1708
** rsym = the irrep of the r-indices
1709
** firstr = the first r-index for each irrep
1710
** lastr = the last r-index for each irrep
1711
** firsts = the first s-index for each irrep
1712
** lasts = the last s-index for each irrep
1713
** sortby_rs = described above
1714
** occ = the Pitzer -> QTS ordering array for the occupied orbitals
1715
** ioff3 = offset array for ndocc*nmo arrays
1716
** printflag = verbosity level for printing
1717
** outfile = text output file
1720
void yosh_wrt_arr_mp2r12a(struct yoshimine *YBuff, int p, int q, int pq,
1721
int pqsym, double **arr, int rsym, int *firstr,
1722
int *lastr, int *firsts, int *lasts, int sortby_rs,
1723
int *occ, int *ioff3,
1724
int printflag, FILE *outfile)
1731
int lastflag = 0, firstfile;
1732
struct bucket *bptr;
1735
firstfile = YBuff->first_tmp_file;
1736
ssym = pqsym ^ rsym;
1737
for (r=firstr[rsym], R=0; r <= lastr[rsym]; r++, R++) {
1739
for (s=firsts[ssym], S=0; s <=lasts[ssym]; s++, S++) {
1741
rs = ioff3[rnew] + snew;
1744
if (fabs(value) > YBuff->cutoff) {
1745
/* figure out which bucket the integral should go in */
1746
if (sortby_rs) whichbucket = YBuff->bucket_for_pq[rs];
1747
else whichbucket = YBuff->bucket_for_pq[pq];
1749
bptr = YBuff->buckets+whichbucket ;
1750
tmpi = (bptr->in_bucket)++ ;
1752
/* sortbuf wants to sort by the first index. If we want to
1753
* sort by rs, we need to map r->p, s->q, p->r, q->s
1754
* We also write out the QTS ordered index for r, and
1755
* not the Pitzer index.
1758
bptr->p[tmpi] = rnew;
1759
bptr->q[tmpi] = snew;
1766
bptr->r[tmpi] = rnew;
1767
bptr->s[tmpi] = snew;
1770
bptr->val[tmpi] = value;
1773
fprintf(outfile, "%4d %4d %4d %4d %10.6lf\n",
1776
/* if we need to flush bucket to disk */
1777
if ((tmpi+1) == YBuff->bucketsize) {
1778
flush_bucket(bptr, lastflag);
1779
bptr->in_bucket = 0;
1782
} /* end if (fabs(arr[rs])) ... */
1784
} /* end loop over s */
1786
} /* end loop over r */
1789
}} // end namespace psi::transqt
added
removed
src/bin/transqt/yoshimine.cc
