(a,b) | |||||
A1 | A2 | B1 | B2 | ||
A1 | X | 0 | 0 | 0 | |
A2 | 0 | X | 0 | 0 | |
(i,j) | |||||
B1 | 0 | 0 | X | 0 | |
B2 | 0 | 0 | 0 | X |
int dpd_init(int dpd_num, int nirreps, long int memory, int cachetype, int *cachefiles, int **cachelist, struct dpd_file4_cache_entry *priority, int num_subspaces, int *orbspi, int *orbsym, ...);
For four-index quantities, the dpd_init() function pre-computes a number of orbital lookup arrays about all pairwise combinations of the given orbital subspaces, including information regarding possible permutational symmetry or antisymmetry among the orbitals. libdpd then assigns a particular "pair number" to each possible combination of indices from the orbital subspace data provided to dpd_init(). For example, given only occupied and virtual orbital subspaces (as would be the case for RHF and ROHF reference wave functions), libdpd will automatically construct the following twelve possible index pairings:
Pair # | Left Subspace | Right Subspace | Permutational Symmetry | Index Packing |
0 | occupied | occupied | None | All p,q |
1 | occupied | occupied | Symmetric | p>q |
2 | occupied | occupied | Antisymmetric | p>q |
3 | occupied | occupied | Symmetric | p>=q |
4 | occupied | occupied | Antisymmetric | p>=q |
5 | virtual | virtual | None | All p,q |
6 | virtual | virtual | Symmetric | p>q |
7 | virtual | virtual | Antisymmetric | p>q |
8 | virtual | virtual | Symmetric | p>=q |
9 | virtual | virtual | Antisymmetric | p>=q |
10 | occupied | virtual | None | All p,q |
11 | virtual | occupied | None | All p,q |
Given the two-electron integral group <ij||ab> for example, we may wish to store these integrals in a matrix with the compound row index (i,j) and compound column index (a,b), as described earlier. Furthermore, since these integrals have perumtational antiysmmetry between indices i and j and between indices a and b, we may wish to store the integrals in a manner which avoids redundancy. To do this, we must choose from the above table the appropriate pair number for the row and column compound indices which reflects the desired permutational antisymmetry and index packing characteristics. For the current example, we would choose pair #2 for the row index and pair #7 for the column index (so that the i=j and a=b terms, which are zero, would be omitted from storage). On the other hand, if we wished to store the "normal" Dirac notation integrals <ia|jk>, which contain three occupied indices and do not have permutational symmetry, we would choose pair #10 for the row index and pair #0 for the column index.
For UHF references, for example, where one must use four orbital subspaces (alpha occupied, alpha virtual, beta occupied, and bet virtual) libdpd will automatically construct the following 32 possible index pairings:
Pair # | Left Subspace | Right Subspace | Permutational Symmetry | Index Packing (Shorthand) |
0 | alpha occupied | alpha occupied | None | I,J |
1 | alpha occupied | alpha occupied | Symmetric | I>J |
2 | alpha occupied | alpha occupied | Antisymmetric | I>J |
3 | alpha occupied | alpha occupied | Symmetric | I>=J |
4 | alpha occupied | alpha occupied | Antisymmetric | I>=J |
5 | alpha virtual | alpha virtual | None | A,B |
6 | alpha virtual | alpha virtual | Symmetric | A>B |
7 | alpha virtual | alpha virtual | Antisymmetric | A>B |
8 | alpha virtual | alpha virtual | Symmetric | A>=B |
9 | alpha virtual | alpha virtual | Antisymmetric | A>=B |
10 | beta occupied | beta occupied | None | i,j |
11 | beta occupied | beta occupied | Symmetric | i>j |
12 | beta occupied | beta occupied | Antisymmetric | i>j |
13 | beta occupied | beta occupied | Symmetric | i>=j |
14 | beta occupied | beta occupied | Antisymmetric | i>=j |
15 | beta virtual | beta virtual | None | a,b |
16 | beta virtual | beta virtual | Symmetric | a>b |
17 | beta virtual | beta virtual | Antisymmetric | a>b |
18 | beta virtual | beta virtual | Symmetric | a>=b |
19 | beta virtual | beta virtual | Antisymmetric | a>=b |
20 | alpha occupied | alpha virtual | None | I,A |
21 | alpha virtual | alpha occupied | None | A,I |
22 | alpha occupied | beta occupied | None | I,j |
23 | beta occupied | alpha occupied | None | j,I |
24 | alpha occupied | beta virtual | None | I,a |
25 | beta virtual | alpha occupied | None | a,I |
26 | alpha virtual | beta occupied | None | A,i |
27 | beta occupied | alpha virtual | None | i,A |
28 | alpha virtual | beta virtual | None | A,b |
29 | beta virtual | alpha virtual | None | b,A |
30 | beta occupied | beta virtual | None | i,a |
31 | beta virtual | beta occupied | None | a,i |
When the program is finished with all libdpd tasks, the global memory used by the library is deallocated using the int dpd_close(int dpdnum) function, which takes the DPD id as its argument.
typedef struct { int nirreps; /* No. of irreps */ int pqnum; /* Pair number for the row indices */ int rsnum; /* Pair number for the column indices */ int *rowtot; /* Row dimension for each submatrix */ int *coltot; /* Column dimension for each submatrix */ int **rowidx; /* Row index lookup array */ int **colidx; /* Column index lookup array */ int ***roworb; /* Row index -> orbital index lookup array */ int ***colorb; /* Column index -> orbital index lookup array */ int *ppi; /* Number of p indices per irrep */ int *qpi; /* Number of q indices per irrep */ int *rpi; /* Number of r indices per irrep */ int *spi; /* Number of s indices per irrep */ int *poff; /* Orbital offset for p */ int *qoff; /* Orbital offset for q */ int *roff; /* Orbital offset for r */ int *soff; /* Orbital offset for s */ int *psym; /* Orbital symmetry for index p */ int *qsym; /* Orbital symmetry for index q */ int *rsym; /* Orbital symmetry for index r */ int *ssym; /* Orbital symmetry for index s */ int perm_pq; /* Can p and q be permuted? (Boolean) */ int perm_rs; /* Can r and s be permuted? (Boolean) */ int peq; /* Can p and q be equal? (Boolean) */ int res; /* Can r and s be equal? (Boolean) */ } dpdparams4;It should be noted that the rowidx and colidx arrays require absolute orbital indices (within the appropriate subspace, of course), and rowborb and colorb arrays return absolute orbital indices. The offset arrays, poff, qoff, etc. may be used to compute the absolute indices from irrep-relative indices using a function of the form pabs = poff[irrep] + prel, where pabs is the absolute index and prel is the index within irrep.
For two-index quantities, the corresponding storage structure is:
typedef struct { int nirreps; /* No. of irreps */ int pnum; /* Orbital subspace for p */ int qnum; /* Orbital subspace for q */ int *rowtot; /* Row dimension for each submatrix */ int *coltot; /* Column dimension for each submatrix */ int *rowidx; /* Row index lookup array */ int *colidx; /* Column index lookup array */ int **roworb; /* Row index -> orbital index lookup array */ int **colorb; /* Column index -> orbital index lookup array */ int *ppi; /* Number of p indices per irrep */ int *qpi; /* Number of q indices per irrep */ int *poff; /* Orbital offset for p */ int *qoff; /* Orbital offset for q */ int *psym; /* Orbital symmetry for index p */ int *qsym; /* Orbital symmetry for index q */ } dpdparams2;The library interacts with data stored on disk via the dpdfile4 and dpdfile2 structures, which have dpdparams4 and dpdparams2 structures, respectively, as members:
typedef struct { char label[PSIO_KEYLEN]; /* Label needed by the I/O routines */ int filenum; /* The PSI unit number */ int my_irrep; /* Total irrep of this quantity */ psio_address *lfiles; /* The disk address of each irrep of data */ dpdparams4 *params; /* The current parameter structure */ int incore; /* Is this file4 already in memory? */ double ***matrix; /* Data */ } dpdfile4;The one-electron (two-index) counterpart of dpdfile4 is:
typedef struct { char label[PSIO_KEYLEN]; /* Label needed by the I/O routines */ int filenum; /* The PSI unit number */ int my_irrep /* The total irrep of this quantity */ psio_address *lfiles; /* The disk address of each irrep of data */ dpdparams2 *params; /* The current parameter structure */ int incore; /* Is this file2 already in memory? */ double ***matrix; /* Data */ } dpdfile2;As described above, one may need to store a given quantity in memory in a manner different from that which is used on disk. The library distinguishes between on-disk and in-memory storage using the dpdbuf4 structure (for four-index quantities only), which has a corresponding dpdfile4 as a member:
typedef struct { int anti; /* Boolean for on-the-fly antisymmetrization */ dpdparams4 *params; /* Parameters for in-memory storage */ dpdfile4 file; /* Structure for on-disk storage */ dpdshift4 shift; /* Structure for left or right index shifting */ double ***matrix; /* Data */ } dpdbuf4;The dpdbuf4 struct is the most common way by which data are transferred between disk and memory. The buffers are initialized using the dpd_buf4_init() function, which takes the following arguments (see dpd.h for the proper syntax):
int dpd_buf4_init(dpdbuf4 *Buf, int inputfile, int pqnum, int rsnum, int file_pqnum, int file_rsnum, int anti, char *label);
dpd_buf4_init(Buffer, PSIFILE, 0, 2, 5, 0, 5, 1, "<ij|ab> integrals");
These argument indicate that while the integrals are stored as (pqnum,rsnum) = (0, 5) on disk, they will be stored in memory as (2, 5), with the row indices packed. In addition, the data will be antisymmetrized as it is read into memory. Note, however, that the libdpd routines do not make sure that your intialization request makes sense (apart from some pre-processor defined debugging in a few high-level functions), and it is quite possible to erroneously pack indices for non-(anti)symmetric quantities.
When the program is finished with a given buffer, the int dpd_buf4_close(dpdbuf4 *Buf) function is called to deallocate the associated memory.
int dpd_contract444(dpdbuf4 *X, dpdbuf4 *Y, dpdbuf4 *Z, int target_X, int target_Y, double alpha, double beta): This function contracts two four-index buffers, X and Y, into a target four-index buffer, Z, using the general formula, alpha * X(pq,mn) * Y(mn,rs) = beta * Z(pq,rs). The current version of this function requires that the target (external) indices must both be contained in the bra (row) or ket (column) of X and Y. The value of target_X indicates that the target indices of X are contained in its bra (0) or its ket (1); target_Y is defined similarly.
For example, we may use dpd_contract444() to evaluate the following contraction found in the T2 amplitude equations from coupled cluster theory:
(A summation over the repeated indices c and d is implied.) Assuming that the <ab|cd> integrals exist on disk in a (pqnum,rsnum) = (5,5) format and the right-hand-side T amplitudes exist on disk in a (2,7) format (with i>j and c>d), the following code block will evaluate this contraction (ignoring the complications of spin factors):
#include <dpd.h> dpdbuf4 T2, T2new, Ints; dpd_buf4_init(&T2, T2FILE, 0, 2, 7, 2, 7, 0, "T2(ij,ab)"); dpd_buf4_init(&T2new, T2newFile, 0, 2, 7, 2, 7, 0, "T2new(ij,ab)"); dpd_buf4_init(&Ints, INTSFILE, 0, 7, 7, 5, 5, 1, "<ab|cd> integrals"); dpd_contract444(&T2, &Ints, &T2new, 0, 1, 1.0, 1.0); dpd_buf4_close(&Ints); dpd_buf4_close(&T2new); dpd_buf4_close(&T2);
The dpd_contract444() function will allocate memory for a single symmetry block of each of the two factors (T2 and Ints) as well as for the target (T2new). It then reads the data from disk for each buffer (including antisymmetrization and index packing for Ints), and then carries our the matrix multiplication. The current symmetry block of the product T2new is then written out to disk and the memory for all three buffers is deallocated. This procedure is repeated for all irreps.
int dpd_contract424(dpdbuf4 *X, dpdfile2 *Y, dpdbuf4 *Z, int sum_X, int sum_Y, int trans_Z, double alpha, double beta): This function evaluates contractions of the form alpha * X(pq,rm) * Y(m,s) = beta * Z(pq,rs). The summation index on X and Y is identified using sum_X and sum_Y, respectively, and the function automatically takes care of any sorting needed to shift the indices into appropriate ordering for matrix multiplication. In addition, the target, Z may be transposed using the trans_Z flag.
int dpd_contract244(dpdfile2 *X, dpdbuf4 *Y, dpdbuf4 *Z, int sum_X, int sum_Y, int trans_Z, double alpha, double beta): This function evaluates contractions of the form alpha * X(p,m) * Y(mq,rs) = beta * Z(pq,rs). The summation index on X and Y is identified using sum_X and sum_Y, respectively, and the function automatically takes care of any sorting needed to shift the indices into appropriate ordering for matrix multiplication. In addition, the target, Z may be transposed using the trans_Z flag.
int dpd_contract222(dpdfile2 *X, dpdfile2 *Y, dpdfile2 *Z, int target_X, int target_Y, double alpha, double beta): This function evaluates contractions of the form alpha * X(p,m) * Y(m,q) = beta * Z(p,q). The target indices on X and Y are identified using target_X and target_Y, respectively.
int dpd_contract422(dpdbuf4 *X, dpdfile2 *Y, dpdfile2 *Z, int trans_Y, int trans_Z, double alpha, double beta): This function evaluates contractions of the form alpha * X(pq,mn) * Y(m,n) = beta * Z(p,q). The bra indices on X must be the target indices and both indices of Y must be summed (of course). The target, Z, may be transposed using the trans_Z flag.
int dpd_contract442(dpdbuf4 *X, dpdbuf4 *Y, dpdfile2 *Z, int target_X, int target_Y, double alpha, double beta): This function evaluates contractions of the form alpha * X(pm,no) * Y(mn,oq) = beta * Z(p,q). The target indices on X and Y are identified using target_X and target_Y, respectively, and these factors are automatically sorted to the appropriate ordering for matrix multiplication to form Z.
int dpd_dot23(dpdfile2 *T, dpdbuf4 *I, dpdfile2 *Z, int transt, int transz, double alpha, double beta): This function evaluates contractions of the form alpha * T(q,r) * I(pq,rs) = beta * Z(p,s), where the summation indices correspond to indices 2 and 3 on the four-index factor. The transt and transz flags may be used to transpose the indices of the factor T or the target Z, as desired. Similar functions exist for other index patterns, including dpd_dot13(), dpd_dot14(), and dpd_dot24().
double dpd_buf4_dot(dpdbuf4 *BufA, dpdbuf4 *BufB): This function evaluates a dot product between two four-index buffers. That is, the function evaluates the sum of the products of coresponding elements of the two buffers.
int dpd_buf4_dirprd(dpdbuf4 *BufA, dpdbuf4 *BufB): This function evaluates a direct product between two four-index buffers and places the output in BufB. That is, the function evaluates each of the products of coresponding elements of the two buffers and places the result in the second buffer.
int dpd_buf4_axpy(dpdbuf4 *BufX, dpdbuf4 *BufY, double alpha): This function evaluates the standard "axpy" function, alpha * X(pq,rs) + Y(pq,rs) = Y(pq,rs).
enum indices {pqrs, pqsr, prqs, prsq, psqr, psrq, qprs, qpsr, qrps, qrsp, qspr, qsrp, rqps, rqsp, rpqs, rpsq, rsqp, rspq, sqrp, sqpr, srqp, srpq, spqr, sprq};
int dpd_buf4_copy(dpdbuf4 *InBuf, int outfilenum, char *label): This function copies a given buffer onto a new disk location identified by outfilenum and label.