7
#include <libciomr/libciomr.h>
10
#define MIN0(a,b) (((a)<(b)) ? (a) : (b))
11
#define MAX0(a,b) (((a)>(b)) ? (a) : (b))
18
** Function reads a file of two-electron integrals into
19
** core and writes them back out again in canonical order. Used in
20
** Yoshimine sorts where we have a file containing all rs for a few
21
** values of pq, but the ints are not in canonical order. At the
22
** very least, we need to sort to make sure that all (pq|rs) for a
23
** given pq are grouped together, since the transformation program
24
** wants to work with all rs for a single pq value at one time.
25
** We may or may not use the restriction pq >= rs (not used if
26
** intermediate = 1, which is how this routine is always called
29
** One interesting issue here is that the intermediate array ('ints')
30
** must be big enough to hold the integrals in the current buffer, but
31
** we don't generally want it to be much larger than necessary! Thus
32
** we calculate an 'offset' which is the canonical index of the first
33
** integral in the buffer, and we use this so that the first integral
34
** in the buffer is stored in ints[0]. What's different for the
35
** Elbert ordering is that we have all rs for a given pq but rs >= pq!
36
** If we want our canonical indices to be consecutive WITHIN THE
37
** CURRENT BUFFER, we MUST use upper triangle ordering rather than
38
** lower triangle! That's what ioff2 is used for. Obviously the
39
** offset must also be calculated with ioff2 for the Elbert order.
40
** Formula for ioff2: ioff2[0] = 0; ioff2[i] = ioff2[i-1] + n - i;
41
** Note that this is not the case when this routine is used for MP2
42
** sorts. The definitions of the ioff arrays can be confusing, so
43
** care should be taken when using this routine.
45
** \param Inbuf = IWL buffer for input
46
** \param Outbuf = IWL buffer for output
47
** \param ints = array to hold integrals in
48
** \param fpq = first pq for this tape
49
** \param lpq = last pq for this tape
50
** \param ioff = offset array for the left indices
51
** \param ioff2 = offset array for Elbert sorts or for the right
52
** indices when no_pq_perm=1
53
** \param nbfso = number of basis functions in SO's
54
** \param lastsort = 1 if this is the last intape, 0 otherwise
55
** \param elbert = integrals obey rs >= pq. Use ioff2 to get offset.
56
** \param intermediate= 1 if sorting a intermediate in the transformation
57
** which is indexed as X[ij][kl] where ij runs from
58
** fpq to lpq and kl runs from 0 to nbstri
59
** \param no_pq_perm = don't use permutational symmetry to swap p and q
60
** (appropriate for MP2 where one is occ and one is virt)
61
** \param qdim = dimensions for the q index...nvirt for MP2
62
** \param add = add contributions to the same integral during sort
63
** \param printflg = 1 for printing, 0 otherwise
64
** \param outfile = output file pointer
68
** Revised 6/27/96 by CDS for new IWL format
69
** N.B. Now need to iwl_flush the output buffer...not done in here!!
72
void sortbuf(struct iwlbuf *Inbuf, struct iwlbuf *Outbuf,
73
double *ints, int fpq, int lpq, int *ioff, int *ioff2,
74
int nbfso, int elbert, int intermediate, int no_pq_perm,
75
int qdim, int add, int printflg, FILE *outfile)
78
Value *valptr; /* array of integral values */
79
Label *lblptr; /* array of integral labels */
80
int idx; /* index for curr integral (0..ints_per_buf) */
81
int lastbuf; /* last buffer flag */
82
int p, q, qmax, qmin, r, rmin, rmax, s, smin, smax, pq, rs;
83
long int pqrs, offset;
84
int first_p, first_q, first_pq, last_p, last_q;
88
fprintf(outfile, "\nsortbuf for pq=%d to %d\n", fpq, lpq);
91
if (no_pq_perm && !intermediate) {
92
fprintf(outfile,"(sortbuf): illegal parameter combination.\n");
93
fprintf(stderr, "(sortbuf): illegal parameter combination.\n");
96
nbstri = nbfso * (nbfso + 1) / 2;
98
/* figure out ranges on things */
99
/* I believe this section works fine, even with different ioff arrays */
101
while (fpq >= ioff[i] && i < BIGNUM) i++;
103
fprintf(outfile, "(sortbuf): parameter error\n") ;
106
first_p = i-1 ; first_q = fpq - ioff[i-1];
107
first_pq = ioff[first_p] + first_q;
108
if (first_pq != fpq) {
109
fprintf(outfile, "(sortbuf): fpq != first_pq.\n");
110
fprintf(stderr, "(sortbuf): fpq != first_pq.\n");
114
if (elbert) offset = ioff2[first_pq] + first_pq ;
115
else offset = ioff[first_pq] ;
120
while (lpq >= ioff[i] && i < BIGNUM) i++ ;
122
fprintf(outfile, "(sortbuf): parameter error\n") ;
125
last_p = i-1 ; last_q = lpq - ioff[i-1] ;
128
lblptr = Inbuf->labels;
129
valptr = Inbuf->values;
131
/* read a buffer at a time until we're done */
134
iwl_buf_fetch(Inbuf);
135
lastbuf = Inbuf->lastbuf;
136
for (idx=4*Inbuf->idx; Inbuf->idx<Inbuf->inbuf; Inbuf->idx++) {
137
p = (int) lblptr[idx++];
138
q = (int) lblptr[idx++];
139
r = (int) lblptr[idx++];
140
s = (int) lblptr[idx++];
142
/* if (no_pq_perm) ioff is the appropriate offset array for the left
143
indices (ioff[p] = nvirt * p for MP2); ioff2 is then the usual
144
ioff offset array, used for the right indices */
147
rs = ioff2[MAX0(r,s)] + MIN0(r,s);
150
pq = ioff[MAX0(p,q)] + MIN0(p,q);
151
rs = ioff[MAX0(r,s)] + MIN0(r,s);
156
pqrs = ioff2[pq] + rs;
158
pqrs = ioff[MAX0(pq,rs)];
163
pqrs = (pq - first_pq);
168
if (printflg && ints[pqrs-offset] != 0.0)
169
fprintf(outfile, "Adding %10.6lf to el %d %d %d %d = %10.6lf\n",
170
valptr[Inbuf->idx], p, q, r, s, ints[pqrs-offset]);
172
if (add) ints[pqrs-offset] += valptr[Inbuf->idx];
173
else ints[pqrs-offset] += valptr[Inbuf->idx];
176
fprintf(outfile, "<%d %d %d %d | %d %d [%ld] = %10.6lf\n",
177
p, q, r, s, pq, rs, pqrs, ints[pqrs-offset]) ;
181
/* now write them out again, in order */
182
lblptr = Outbuf->labels;
183
valptr = Outbuf->values;
187
for (p=first_p; p<=last_p; p++) {
188
qmax = (p==last_p) ? last_q : p ;
189
qmin = (p==first_p) ? first_q : 0 ;
191
qmax = (p==last_p) ? last_q : (qdim-1);
192
qmin = (p==first_p) ? first_q: 0;
194
for (q=qmin; q<=qmax; q++) {
195
pq = ioff[p] + q ; /* This should be fine even with MP2 */
198
rmin = (elbert) ? p : 0 ;
199
rmax = (elbert) ? nbfso : p+1 ;
201
else { /* This should be fine with MP2, also */
206
for (r=rmin; r<rmax; r++) {
211
smin = (p==r) ? q : 0 ;
214
smax = (p==r) ? (q+1) : (r+1) ;
218
else { /* This should be fine with MP2, also */
223
for (s=smin; s < smax; s++) {
224
if(no_pq_perm) rs = ioff2[r] + s;
225
else rs = ioff[r] + s;
227
/* Again, this should be fine with MP2 */
228
if (elbert) pqrs = ioff2[pq] + rs ;
229
else if (intermediate) {
230
pqrs = (pq - first_pq);
234
else pqrs = ioff[pq] + rs ;
236
if (fabs(ints[pqrs-offset]) > Outbuf->cutoff) {
242
valptr[Outbuf->idx] = ints[pqrs-offset];
244
fprintf(outfile, ">%d %d %d %d | %d %d [%ld] = %10.6lf\n",
245
p, q, r, s, pq, rs, pqrs, ints[pqrs-offset]) ;
248
if (Outbuf->idx == Outbuf->ints_per_buf) {
250
Outbuf->inbuf = Outbuf->idx;
added
removed
src/lib/libiwl/sortbuf.c
