Programmer's Manual for LIBDPD:
A Library for Including Spatial Symmetry in Quantum Chemistry Programs

T. Daniel Crawford
Version: June 18, 2001
crawdad@vt.edu


I. Introduction

In many-body theories such as coupled cluster or MBPT, one finds many complicated algebraic expressions involving products of multi-index quantities such as one- and two-electron integrals, cluster amplitudes, density matrices, etc. Efficient evaluation of these products can require a great deal of effort, particularly if one wishes to use Abelian point-group symmetry to reduce the number of terms that must be computed and stored. The direct-product decomposition library, LIBDPD, is designed to assist the programmer with this problem by providing (1) a symmetry-blocked, matrix-based storage scheme for all two- and four-index quantities, (2) a set of functions for evaluating various types of products among them, and (3) a set of utilities for sorting them to different index orderings. The library is currently used in the set of coupled cluster and perturbation theory energy and analytic gradient codes under development in the PSI package of quantum chemical programs. This manual describes the design of the library and provides a set of examples for its use. The header file dpd.h provides proper declarations for all structures and functions used in the library.

II. Fundamental Symmetry Concepts

In the current version of libdpd, I assume that all multi-index quantities are totally symmetric, i.e., the direct product of the irreducible representations (irreps) associated with the component orbital indices corresponds to the totally symmetric irrep of the given point group. (This assumption is in the process of being removed, however.) Furthermore, each orbital subset of interest (e.g., occupied or virtual orbitals) must be grouped by symmetry. This allows one to organize the given quantity in a symmetry-blocked matrix for which only the symmetry-allowed elements are actually stored on disk or in memory. As an example, consider the set of two-electron integrals (in antisymmetrized, Dirac notation) needed for the second-order MBPT energy expression, <ij||ab>, where i and j (a and b) represent occupied (virtual) orbitals. In C2v symmetry, a given integral is zero unless the direct product of the four irreps corresponding to indices, i, j, a, and b is A1, that is

i × j × a × b = A1.

Since the direct product of any irrep with itself also gives the totally symmetric irrep, we may write

i × j = a × b

If we organize the occupied and virtual orbitals into symmetry groups, we may represent the storage of the complete list of <ij||ab> integrals as a block-diagonal supermatrix with the form

(a,b)
A1A2B1B2
A1 X 0 0 0
A2 0 X 0 0
(i,j)
B1 0 0 X 0
B2 0 0 0 X

Here "X" and "0" indicate nonzero and zero submatrices, respectively, and the notation (p,q) is used to indicate the compound row or column indices. The compound indices may be computed based on a knowledge of the permutational symmetry (or antisymmetry) of the index pairs; if the given integrals are antisymmetric with respect to permutation of i and j, for example, we may compute (i,j) from the individual indices keeping i>j. Since we may partition the indices into any grouping we desire, we may select any grouping of the indices to use as row and column indices. For example, since

i × j × a = b

we could also construct a different supermatrix for these same integrals, using (i,j,a) as the compound row index and b as the column index. These same concepts apply to any multi-index quantity.

III. Library Initialization

Before the library may be used, it must first be initialized via the dpd_init() function. This function requires the following information (see dpd.h for proper syntax):

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 SubspaceRight SubspacePermutational SymmetryIndex Packing
0occupiedoccupiedNoneAll p,q
1occupiedoccupiedSymmetricp>q
2occupiedoccupiedAntisymmetricp>q
3occupiedoccupiedSymmetricp>=q
4occupiedoccupiedAntisymmetricp>=q
5virtualvirtualNoneAll p,q
6virtualvirtualSymmetricp>q
7virtualvirtualAntisymmetricp>q
8virtualvirtualSymmetricp>=q
9virtualvirtualAntisymmetricp>=q
10occupiedvirtualNoneAll p,q
11virtualoccupiedNoneAll 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 SubspaceRight SubspacePermutational SymmetryIndex Packing (Shorthand)
0alpha occupiedalpha occupiedNoneI,J
1alpha occupiedalpha occupiedSymmetricI>J
2alpha occupiedalpha occupiedAntisymmetricI>J
3alpha occupiedalpha occupiedSymmetricI>=J
4alpha occupiedalpha occupiedAntisymmetricI>=J
5alpha virtualalpha virtualNoneA,B
6alpha virtualalpha virtualSymmetricA>B
7alpha virtualalpha virtualAntisymmetricA>B
8alpha virtualalpha virtualSymmetricA>=B
9alpha virtualalpha virtualAntisymmetricA>=B
10beta occupiedbeta occupiedNonei,j
11beta occupiedbeta occupiedSymmetrici>j
12beta occupiedbeta occupiedAntisymmetrici>j
13beta occupiedbeta occupiedSymmetrici>=j
14beta occupiedbeta occupiedAntisymmetrici>=j
15beta virtualbeta virtualNonea,b
16beta virtualbeta virtualSymmetrica>b
17beta virtualbeta virtualAntisymmetrica>b
18beta virtualbeta virtualSymmetrica>=b
19beta virtualbeta virtualAntisymmetrica>=b
20alpha occupiedalpha virtualNoneI,A
21alpha virtualalpha occupiedNoneA,I
22alpha occupiedbeta occupiedNoneI,j
23beta occupiedalpha occupiedNonej,I
24alpha occupiedbeta virtualNoneI,a
25beta virtualalpha occupiedNonea,I
26alpha virtualbeta occupiedNoneA,i
27beta occupiedalpha virtualNonei,A
28alpha virtualbeta virtualNoneA,b
29beta virtualalpha virtualNoneb,A
30beta occupiedbeta virtualNonei,a
31beta virtualbeta occupiedNonea,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.

IV. Basic Data Structures and Storage Hierarchy

Since one may need to store a given quantity in memory in a manner different from that which is used on disk (e.g., one may need to pack or unpack indices), the library distinguishes between on-disk and in-memory storage of four-index quantities. For example, one might choose to store <ij|ab> integrals on disk using an (i,j)x(a,b) matrix with all values of i, j, a, and b (pair numbers 0 and 5 from the previous section). However, these can be automatically antisymmetrized as they are read into memory and the resulting <ij||ab> stored in a packed form with i>j and a>b. This section described the basic storage hierarchy that automates this process. For a general four-index quantity with dummy indices p, q, r, and s, the storage of the data (whether on disk or in memory) is defined by the dpdparams4 structure. Although the programmer should almost always be able to avoid direct interaction with the dpdparams4 structure, the library isn't perfectly object-oriented and so access may be necessary in some cases. The structure contains the following information:

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);

For example, consider again the Dirac-notation integrals <ij|ab> (which are not antisymmetrized) which are already stored on disk. If we wish to read the antisymmetrized form of these integrals into memory but to also pack the i and j indices (i.e. keep i>j), we might call dpd_buf4_init() in the following manner:

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.

V. Contraction Evaluation

The library provides a number of high-level functions for evaluating a variety of products among two- and four-index buffers. This section outlines these functions and provides simple examples of their use.

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:

tijab = tijcd <ab||cd>

(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).

VI. Sorting Utilities

There also exist a number of utilities for sorting two- and four-index quantities to various arrangements. int dpd_buf4_sort(dpdbuf4 *InBuf, int outfilenum, enum indices index, int pqnum, int rsnum, char *label): This function takes a given four-index InBuf and reorders indices according to the index variable, placing the result in outfilenum with pair values pqnum and rsnum. The argument index is an enumerated type with possible values:

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.


T. Daniel Crawford  / crawdad@vt.edu
Last modified: Mon Jun 18 17:33:08 2001