Skip to Main Content

Skip Nav Destination

Article navigation

Rapid Communication| August 06 2014

Communication: A reduced scaling J-engine based reformulation of SOS-MP2 using graphics processing units

S. A. Maurer;

S. A. Maurer

1Chair of Theoretical Chemistry, Department of Chemistry,

University of Munich (LMU)

, Butenandtstr. 7, D-81377 München,

Germany

2Center for Integrated Protein Science (CIPSM) at the Department of Chemistry,

University of Munich (LMU)

, Butenandtstr. 5–13, D-81377 München,

Germany

Search for other works by this author on:

PubMed

Google Scholar

J. Kussmann;

J. Kussmann

1Chair of Theoretical Chemistry, Department of Chemistry,

University of Munich (LMU)

, Butenandtstr. 7, D-81377 München,

Germany

2Center for Integrated Protein Science (CIPSM) at the Department of Chemistry,

University of Munich (LMU)

, Butenandtstr. 5–13, D-81377 München,

Germany

Search for other works by this author on:

PubMed

Google Scholar

C. Ochsenfeld

C. Ochsenfeld ^a)

1Chair of Theoretical Chemistry, Department of Chemistry,

University of Munich (LMU)

, Butenandtstr. 7, D-81377 München,

Germany

2Center for Integrated Protein Science (CIPSM) at the Department of Chemistry,

University of Munich (LMU)

, Butenandtstr. 5–13, D-81377 München,

Germany

Search for other works by this author on:

PubMed

Google Scholar

Author & Article Information

Electronic mail: [email protected]

J. Chem. Phys. 141, 051106 (2014)

https://doi.org/10.1063/1.4891797

We present a low-prefactor, cubically scaling scaled-opposite-spin second-order Møller-Plesset perturbation theory (SOS-MP2) method which is highly suitable for massively parallel architectures like graphics processing units (GPU). The scaling is reduced from

${\cal O\mbox{(N^{5})}}$

O (N^{5})

to

${\cal O\mbox{(N^{3})}}$

O (N^{3})

by a reformulation of the MP2-expression in the atomic orbital basis via Laplace transformation and the resolution-of-the-identity (RI) approximation of the integrals in combination with efficient sparse algebra for the 3-center integral transformation. In contrast to previous works that employ GPUs for post Hartree-Fock calculations, we do not simply employ GPU-based linear algebra libraries to accelerate the conventional algorithm. Instead, our reformulation allows to replace the rate-determining contraction step with a modified J-engine algorithm, that has been proven to be highly efficient on GPUs. Thus, our SOS-MP2 scheme enables us to treat large molecular systems in an accurate and efficient manner on a single GPU-server.

I. INTRODUCTION

Second-order Møller-Plesset perturbation theory (MP2)¹ is an important method in the field of ab initio quantum chemistry, however, the

${\cal O\mbox{(N^{5})}}$

O (N^{5})

scaling with molecular size N hampers its applicability to large systems. Therefore, many efficient low- or linear-scaling methods have been developed over the last decades, e.g., see Refs. 2–8. Despite its success, reducing the prefactors is a central issue which is the focus of our present work.

Our method is based on the Laplace scaled-opposite-spin MP2 approach within the resolution-of-the-identity (RI) approximation (SOS-RI-MP2) of Jung et al.⁹ where the essential idea is a reversed order of summation compared to a conventional RI-MP2 calculation, i.e., the molecular orbital (MO) indices are contracted first before the summations over auxiliary functions are carried out. In this work, we use an atomic orbital-based reformulation which leads to an asymptotic cubic scaling, that is one order of magnitude smaller than in the original MO-based SOS-RI-MP2 approach, and which allows for efficient evaluation on graphics processing units (GPUs). In contrast to the asymptotically cubic scaling local SOS-RI-MP2 variant of Jung et al.,¹⁰ our method does not rely on a local metric or orbital localization, but uses local atomic orbitals (AOs).

In the last years several methods focusing on using GPUs for post Hartree-Fock (HF) methods have been published^11–15 that employ GPUs only for linear algebra operations. In this work, however, our reformulation for SOS-RI-MP2 allows to replace the computationally demanding contraction step with a modified J-engine algorithm, which has already been proven to be highly efficient on GPUs within self-consistent field (SCF) calculations.^16,17 This integral-direct approach also enables a complete avoidance of disk-IO.

II. THEORY

The conventional RI-MP2 expression for the opposite-spin MP2 term reads

\begin{eqnarray}&&\hspace*{-7pt} E_{RI-MP2}^{OS}\nonumber\\&&\hspace*{-7pt}\quad =\! -\! \sum _{i j a b} \sum _{R S R^{\prime } S^{\prime }}\!\! \frac{(ia|R) [ {\bf J}^{-1} ]_{R S} (S|jb) (ia|R^{\prime }) [ {\bf J}^{-1} ]_{R^{\prime } S^{\prime }} (S^{\prime }|jb)}{\epsilon _{a} + \epsilon _{b} - \epsilon _{i} - \epsilon _{j}},\nonumber\\\end{eqnarray}

\begin{matrix} E_{R I - M P 2}^{O S} \\ = - \sum_{i j a b} \sum_{R S R^{'} S^{'}} \frac{(i a | R) {[J^{- 1}]}_{R S} (S | j b) (i a | R^{'}) {[J^{- 1}]}_{R^{'} S^{'}} (S^{'} | j b)}{ε_{a} + ε_{b} - ε_{i} - ε_{j}}, \end{matrix}

(1)

where the matrix J contains the two-center two-electron integrals in the auxiliary basis (R, S). Following the Laplace AO-MP2 approach^2–4 this expression can be transformed into the AO basis

\begin{eqnarray}E_{RI-AO-MP2}^{OS} &=& - \sum _{\alpha } \sum _{\genfrac{}{}{0.0pt}{}{\mu \nu \lambda \sigma }{\mu ^{\prime }\nu ^{\prime }\lambda ^{\prime }\sigma ^{\prime }}} \sum _{R S R^{\prime } S^{\prime }} \underline{P}_{\mu \mu ^{\prime }} \overline{P}_{\nu \nu ^{\prime }} \underline{P}_{\lambda \lambda ^{\prime }} \overline{P}_{\sigma \sigma ^{\prime }}\nonumber\\&&\times\,(\mu \nu |R) [ {\bf J}^{-1} ]_{R S} (S|\lambda \sigma ) (\mu ^{\prime }\nu ^{\prime }|R^{\prime }) [ {\bf J}^{-1} ]_{R^{\prime } S^{\prime }}\nonumber\\&&\times\, (S^{\prime }|\lambda ^{\prime }\sigma ^{\prime }).\end{eqnarray}

\begin{matrix} E_{R I - A O - M P 2}^{O S} & = & - \sum_{α} \sum_{\binom{μ ν λ σ}{μ^{'} ν^{'} λ^{'} σ^{'}}} \sum_{R S R^{'} S^{'}} {\underset{̲}{P}}_{μ μ^{'}} {\bar{P}}_{ν ν^{'}} {\underset{̲}{P}}_{λ λ^{'}} {\bar{P}}_{σ σ^{'}} \\ \times (μ ν | R) {[J^{- 1}]}_{R S} (S | λ σ) (μ^{'} ν^{'} | R^{'}) {[J^{- 1}]}_{R^{'} S^{'}} \\ \times (S^{'} | λ^{'} σ^{'}) . \end{matrix}

(2)

The sum over α corresponds to the quadrature sum of the Laplace transform of the energy denominator where usually 5-6 points have been found to provide good accuracy.^4,18,19 The pseudo-densities

$\underline{\bf P}$

\underset{̲}{P}

and

$\overline{\bf P}$

\bar{P}

are defined as

\begin{equation}\begin{split} \underline{P}_{\mu \nu } &= (\omega _{\alpha })^{\frac{1}{4}} \sum _{i}^{\text{occ}} c_{\mu i} e^{\epsilon _{i} t_{\alpha }} c_{\nu i}, \\\overline{P}_{\mu \nu } &= (\omega _{\alpha })^{\frac{1}{4}} \sum _{a}^{\text{virt}} c_{\mu a} e^{-\epsilon _{a} t_{\alpha }} c_{\nu a}. \end{split}\end{equation}

\begin{matrix} {\underset{̲}{P}}_{μ ν} & = {(ω_{α})}^{\frac{1}{4}} \sum_{i}^{occ} c_{μ i} e^{ε_{i} t_{α}} c_{ν i}, \\ {\bar{P}}_{μ ν} & = {(ω_{α})}^{\frac{1}{4}} \sum_{a}^{virt} c_{μ a} e^{- ε_{a} t_{α}} c_{ν a} . \end{matrix}

(3)

For an efficient evaluation, we define intermediates

\begin{eqnarray}Z_{R S} &=& \sum _{\mu \nu \mu ^{\prime }\nu ^{\prime }} (R|\mu ^{\prime }\nu ^{\prime }) \underline{P}_{\mu \mu ^{\prime }} \overline{P}_{\nu \nu ^{\prime }} (\mu \nu |S)\nonumber\\&=& \sum _{\mu \nu } (R|\underline{\mu }\overline{\nu }) (\mu \nu |S),\end{eqnarray}

\begin{matrix} Z_{R S} & = & \sum_{μ ν μ^{'} ν^{'}} (R | μ^{'} ν^{'}) {\underset{̲}{P}}_{μ μ^{'}} {\bar{P}}_{ν ν^{'}} (μ ν | S) \\ = & \sum_{μ ν} (R | \underset{̲}{μ} \bar{ν}) (μ ν | S), \end{matrix}

(4)

which allow to calculate the OS-MP2 energy via

\begin{equation}{\bf \tilde{Z}} = {\bf Z J^{-1}}\end{equation}

\tilde{Z} = Z J^{- 1}

(5)

and

\begin{equation}E_{RI-AO-MP2}^{OS} = - \sum _{\alpha } \sum _{R S} \tilde{Z}_{R S} \tilde{Z}_{S R}.\end{equation}

E_{R I - A O - M P 2}^{O S} = - \sum_{α} \sum_{R S} {\tilde{Z}}_{R S} {\tilde{Z}}_{S R} .

(6)

The crucial step is the calculation of the intermediate matrix Z which can be calculated with cubic scaling cost as shown below. It should be noted that the matrix Z differs from the matrices X and Y used in the methods of Jung et al. (see Eqs. (12) and (19) in Ref. 10), where the inverse, or inverse square root, respectively, of the matrix J is included.

The formation of the Z matrix requires transformed three-center integrals

$(R|\underline{\mu }\overline{\nu })$

(R | \underset{̲}{μ} \bar{ν})

which we calculate in an asymptotically quadratic step from the untransformed integrals (R|μν) using sparse matrix multiplications in the highly efficient block-compressed sparse row (BCSR) format²⁰ (see Ref. 21 for details of our BCSR implementation). For an efficient transformation, it is furthermore advantageous to perform a Cholesky decomposition of the occupied pseudo-density,

\begin{equation}\underline{\bf P} = \underline{\bf L} \underline{\bf L} ^T,\end{equation}

\underset{̲}{P} = \underset{̲}{L} {\underset{̲}{L}}^{T},

(7)

which provides a coefficient matrix

$\underline{\bf L}$

\underset{̲}{L}

for local occupied pseudo-MOs¹⁹ that can be used as an intermediate basis. The transformation starts with a sparse square matrix of integrals for a given auxiliary index M^(R) with

$M^{(R)}_{\mu \nu } = (R|\mu \nu )$

M_{μ ν}^{(R)} = (R | μ ν)

where the sparsity reflects the linear-scaling number of significant basis function products. The transformation to the transformed integral block

\begin{equation}\underline{\overline{\bf M}}^{(R)} = \underline{\bf P} {\bf M}^{(R)} \overline{\bf P}\end{equation}

{\underset{̲}{\bar{M}}}^{(R)} = \underset{̲}{P} M^{(R)} \bar{P}

(8)

with

$\underline{\overline{M}}^{(R)} = (R|\underline{\mu }\overline{\nu })$

{\underset{̲}{\bar{M}}}^{(R)} = (R | \underset{̲}{μ} \bar{ν})

is then performed in three consecutive multiplications,

\begin{equation}\underline{\overline{\bf M}}^{(R)} = \underline{\bf L} ((\underline{\bf L} ^T {\bf M}^{(R)} ) \overline{\bf P} ).\end{equation}

{\underset{̲}{\bar{M}}}^{(R)} = \underset{̲}{L} (({\underset{̲}{L}}^{T} M^{(R)}) \bar{P}) .

(9)

It should be noted that the number of occupied Cholesky pseudo-MOs is equal or, in some cases, slightly smaller than the number of occupied MOs (see Ref. 19 for an in-depth discussion) and the dimension of the matrices is therefore largely reduced by the first transformation with the Cholesky coefficients. In any of the three multiplications, two matrix dimensions correspond to the number of atomic orbitals and one corresponds to the number of occupied pseudo-MOs. Further savings in the final transformation step are possible, if one exploits the restriction of the subsequent contraction Eq. (4) to significant basis function products: The final multiplication in Eq. (9) can be restricted to matrix blocks which contain (⁠

$\underline{\mu }, \overline{\nu }$

\underset{̲}{μ}, \bar{ν}

⁠) pairs that match significant basis function pairs. All matrices have a sparse structure and the matrix multiplications therefore scale asymptotically linearly with system size for a given auxiliary index which leads to an asymptotically quadratic scaling of the total transformation step.

An outline of the steps in our method is given in Table I and timings on linear alkanes and DNA systems are given in Tables II and III. Steps (1)–(4), (7), and (8) have negligible cost (<10% in all calculations, 2% for the largest DNA system) and the effective scaling for these steps is at most cubically with system size. The cubic scaling steps (3), (4), and (7) are performed using highly optimized linear algebra routines. The significant steps are the transformation of the three-center integrals, step (5), and the contraction with the untransformed integrals in step (6) which is performed on GPU – details of our algorithm are presented in Sec. III. The scaling of the transformation step (5) is asymptotically quadratic and this behavior is observed as the general trend for large linear alkanes but superimposed with considerable fluctuations due to the blocking in the BCSR format with the chosen block size of around 100 × 100. At the present stage of development, we do not adjust the block structure for the rectangular Cholesky matrices, which leads to some variable overhead in the multiplications which can be significant for smaller systems. For the DNA systems, the transformation matrices are less sparse and the quadratic scaling regime is not reached but an effective cubic scaling is observed.

Table I.

Outline of the algorithm for SOS-RI-AO-MP2. Steps (4)–(7) are repeated for every Laplace point.

(1)

Calculation of (R|μν)

$\mathcal {O}(\mathrm{N^2})$

O (N^{2})

(2)

Calculation of J_RS = (R|S)

$\mathcal {O}(\mathrm{N^2})$

O (N^{2})

(3)

Calculation of

$\bf J^{-1}$

J^{- 1}

$\mathcal {O}(\mathrm{N^3})$

O (N^{3})

(4)

Calculation of pseudo-densities

$\mathcal {O}(\mathrm{N^3})$

O (N^{3})

(5)

Transformation of (R|μν) to

$(R|\underline{\mu }\overline{\nu })$

(R | \underset{̲}{μ} \bar{ν})

$\mathcal {O}(\mathrm{N^2})$

O (N^{2})

(6)

Contraction

$\sum _{\mu \nu } (R|\underline{\mu }\overline{\nu })(\mu \nu |S)$

\sum_{μ ν} (R | \underset{̲}{μ} \bar{ν}) (μ ν | S)

(on GPU)

$\mathcal {O}(\mathrm{N^3})$

O (N^{3})

(7)

Multiplication Z J⁻¹

$\mathcal {O}(\mathrm{N^3})$

O (N^{3})

(8)

Contraction

$\sum _{R S} \tilde{Z}_{R S} \tilde{Z}_{S R}$

\sum_{R S} {\tilde{Z}}_{R S} {\tilde{Z}}_{S R}

$\mathcal {O}(\mathrm{N^2})$

O (N^{2})

Table II.

Wall times and scaling behavior with respect to the number of basis functions for SOS-RI-AO-MP2 calculations on linear alkanes in a def2-SVP and def2-TZVP basis on a computing node with two Intel Xeon E5-2620 processors (12 CPU cores) and four Nvidia GeForce GTX Titan GPUs. Timings and scaling exponents [

${\cal O\mbox{(N^{x})}}$

O (N^{x})

] are given for the whole calculation as well as selected steps of the algorithm defined in Table I.

def2-SVP	# Basis	Total		Transformation (5)		Contraction (6) on GPU		Steps (1)–(4), (7), (8)
system	functions	time [s]	${\cal O\mbox{(N^{x})}}$ $O (N^{x})$	time [s]	${\cal O\mbox{(N^{x})}}$ $O (N^{x})$	time [s]	${\cal O\mbox{(N^{x})}}$ $O (N^{x})$	time [s]	${\cal O\mbox{(N^{x})}}\,$ $O (N^{x})$
C₂₀H₄₂	490	98	…	51	…	39	…	8	…
C₄₀H₈₂	970	575	2.59	374	2.92	171	2.16	30	1.94
C₈₀H₁₆₂	1930	2248	1.98	1051	1.50	1048	2.64	149	2.33
C₁₆₀H₃₂₂	3850	14 134	2.66	6206	2.57	7127	2.78	801	2.44
C₃₂₀H₆₄₂	7690	95 765	2.77	33 881	2.45	56 358	2.99	5526	2.79
def2-TZVP	# Basis	Total		Transformation (5)		Contraction (6) on GPU		Steps (1)–(4), (7), (8)
system	functions	time [s]	${\cal O\mbox{(N^{x})}}$ $O (N^{x})$	time [s]	${\cal O\mbox{(N^{x})}}$ $O (N^{x})$	time [s]	${\cal O\mbox{(N^{x})}}$ $O (N^{x})$	time [s]	${\cal O\mbox{(N^{x})\,}}$ $O (N^{x})$
C₂₀H₄₂	872	406	…	261	…	120	…	25	…
C₄₀H₈₂	1732	1831	2.19	932	1.85	784	2.74	115	2.22
C₈₀H₁₆₂	3452	11 800	2.70	5790	2.65	5445	2.81	565	2.31

def2-SVP	# Basis	Total		Transformation (5)		Contraction (6) on GPU		Steps (1)–(4), (7), (8)
system	functions	time [s]	${\cal O\mbox{(N^{x})}}$ $O (N^{x})$	time [s]	${\cal O\mbox{(N^{x})}}$ $O (N^{x})$	time [s]	${\cal O\mbox{(N^{x})}}$ $O (N^{x})$	time [s]	${\cal O\mbox{(N^{x})}}\,$ $O (N^{x})$
C₂₀H₄₂	490	98	…	51	…	39	…	8	…
C₄₀H₈₂	970	575	2.59	374	2.92	171	2.16	30	1.94
C₈₀H₁₆₂	1930	2248	1.98	1051	1.50	1048	2.64	149	2.33
C₁₆₀H₃₂₂	3850	14 134	2.66	6206	2.57	7127	2.78	801	2.44
C₃₂₀H₆₄₂	7690	95 765	2.77	33 881	2.45	56 358	2.99	5526	2.79
def2-TZVP	# Basis	Total		Transformation (5)		Contraction (6) on GPU		Steps (1)–(4), (7), (8)
system	functions	time [s]	${\cal O\mbox{(N^{x})}}$ $O (N^{x})$	time [s]	${\cal O\mbox{(N^{x})}}$ $O (N^{x})$	time [s]	${\cal O\mbox{(N^{x})}}$ $O (N^{x})$	time [s]	${\cal O\mbox{(N^{x})\,}}$ $O (N^{x})$
C₂₀H₄₂	872	406	…	261	…	120	…	25	…
C₄₀H₈₂	1732	1831	2.19	932	1.85	784	2.74	115	2.22
C₈₀H₁₆₂	3452	11 800	2.70	5790	2.65	5445	2.81	565	2.31

Table III.

Wall times and scaling behavior with respect to the number of basis functions for SOS-RI-AO-MP2 calculations on DNA systems in a def2-SVP basis on a computing node with two Intel Xeon E5-2620 processors (12 CPU cores) and four Nvidia GeForce GTX Titan GPUs. Timings and scaling exponents [

${\cal O\mbox{(N^{x})}}$

O (N^{x})

] are given for the whole calculation as well as selected steps of the algorithm defined in Table I.

				Transformation		Contraction		Steps
	# Basis	Total		(5)		(6) on GPU		(1)–(4), (7), (8)		Reference^a
System	functions	time [s]	${\cal O\mbox{(N^{x})}}$ $O (N^{x})$	time [s]	${\cal O\mbox{(N^{x})}}$ $O (N^{x})$	time [s]	${\cal O\mbox{(N^{x})}}$ $O (N^{x})$	time [s]	${\cal O\mbox{(N^{x})}}$ $O (N^{x})$	time [s]	${\cal O\mbox{(N^{x})}}$ $O (N^{x})$
DNA₁	625	238	…	143	…	82	…	13	…	438	…
DNA₂	1332	2096	2.88	1279	2.90	721	2.87	96	2.64	8290	3.88
DNA₄	2746	19 601	3.09	13 289	3.24	5750	2.87	562	2.44	146 000	3.96
DNA₈	5574	185 659	3.18	126 060	3.18	55 882	3.21	3717	2.67	…	…

				Transformation		Contraction		Steps
	# Basis	Total		(5)		(6) on GPU		(1)–(4), (7), (8)		Reference^a
System	functions	time [s]	${\cal O\mbox{(N^{x})}}$ $O (N^{x})$	time [s]	${\cal O\mbox{(N^{x})}}$ $O (N^{x})$	time [s]	${\cal O\mbox{(N^{x})}}$ $O (N^{x})$	time [s]	${\cal O\mbox{(N^{x})}}$ $O (N^{x})$	time [s]	${\cal O\mbox{(N^{x})}}$ $O (N^{x})$
DNA₁	625	238	…	143	…	82	…	13	…	438	…
DNA₂	1332	2096	2.88	1279	2.90	721	2.87	96	2.64	8290	3.88
DNA₄	2746	19 601	3.09	13 289	3.24	5750	2.87	562	2.44	146 000	3.96
DNA₈	5574	185 659	3.18	126 060	3.18	55 882	3.21	3717	2.67	…	…

^a

Reference calculation with conventional SOS-RI-MP2 (1 core@Intel Xeon E5-2620, only serial version available).

It should be noted that the scaling of the computational cost given in Table I is based on increased sparsity in large molecular systems with a fixed choice of basis set. More generally, one can also express the scaling behavior in terms of the number of atom-centers (N_c) and number of AOs per center (N_ao/c). The number of auxilliary AOs per center can usually be given in terms of N_ao/c multiplied by a constant c_aux, i.e., as c_aux × N_ao/c. Thus, we obtain for the two rate-determining steps (5) and (6) a scaling behavior of c_occ × c_aux ×

$\mathcal {O}(\mathrm{N_{c}^{2}N_{ao/c}^{3}})$

O (N_{c}^{2} N_{ao / c}^{3})

and

$c_{aux}^{2}\times\mathcal {O}(\mathrm{N_{c}^{3}N_{ao/c}^{4}})$

c_{a u x}^{2} \times O (N_{c}^{3} N_{ao / c}^{4})

⁠, respectively, where either

$\rm N_c$

N_{c}

or

$\rm N_{ao/c}$

N_{ao / c}

is constant (scaling with basis-set size or system size, respectively). The first expression with the prefactor c_occ comes from the use of the Cholesky factor

$\underline{\mathbf {L}}$

\underset{̲}{L}

in the sparse multiplications of step (5). These expressions show the scaling for a fixed system size with increasing basis set size, i.e.,

$\rm N_{ao/c}^{3}$

N_{ao / c}^{3}

for step (5) and

$\rm N_{ao/c}^{4}$

N_{ao / c}^{4}

for step (6). Nevertheless, the asymptotic scaling with respect to the system size (i.e.,

$\rm N_{ao/c}$

N_{ao / c}

is constant) is a quadratic (5) or cubic (6) scaling behavior independent of the choice of the basis set.

In order to provide some measure of comparison for the efficiency of our algorithm, we performed SOS-RI-MP2 calculations with a development version of Q-Chem²² for the first three DNA-fragments (1 core, Intel Xeon [email protected]), showing a speed-up of approx. 7.5 for DNA₄. Here, it has to be stressed that a fair and balanced comparison of two different algorithms aiming for different computing architectures is very difficult. For larger systems, however, the speed-up of the presented algorithm will be far larger due to the less favorable

${\cal O\mbox{(N^{4})}}$

O (N^{4})

of the conventional SOS-RI-MP2 algorithm.

In our present implementation, the untransformed three-center integrals are saved to disk once and read in the transformation step for every Laplace point. An implementation without any I/O can be easily realized by performing steps (3), (4), and (5) in batches of R. For each batch of this auxiliary index, all three-center integrals can be kept in memory. The transformation can be performed in memory and in step (5) the untransformed integrals are calculated on-the-fly and directly contracted so that all operations can be performed without any hard disk access. It also has to be stressed that the current implementation is strictly sequential with respect to the single steps given in Table I, i.e., the GPU-kernels block the further execution of the code in order to provide meaningful timings for the single steps of the algorithm. However, the overall computational time can be significantly reduced by executing the GPU-kernels in a non-blocking fashion, so that the contraction step (6) and the transformation step (5) can be executed simultaneously.

III. DETAILS: GPU-BASED ALGORITHM

For an efficient evaluation of the contraction step (6), a modification of the J-engine^23,24 based GPU-algorithm for the evaluation of the Coulomb matrix in SCF calculations¹⁷ is employed. The algorithm is based on the 1 thread/1 primitive integral (1T1PI) approach and the shell-pair ordering as proposed by Ufimtsev and Martínez.^16,25 To be consistent with this work with respect to the choice of bra and ket distributions and to make the analogy to the Coulomb matrix construction most clear, we rewrite Eq. (4) (for real functions) as

\begin{equation}\begin{split} Z_{R S} &= \sum _{\mu \nu } (S|\mu \nu ) \underline{\overline{M}}^{(R)}_{\mu \nu } , \end{split}\end{equation}

\begin{matrix} Z_{R S} & = \sum_{μ ν} (S | μ ν) {\underset{̲}{\bar{M}}}_{μ ν}^{(R)}, \end{matrix}

(10)

with

$\underline{\overline{M}}^{(R)}_{\mu \nu } = (R|\underline{\mu }\overline{\nu })$

{\underset{̲}{\bar{M}}}_{μ ν}^{(R)} = (R | \underset{̲}{μ} \bar{ν})

(cf. Eq. ). In contrast to the regular Coulomb matrix construction, we execute a Coulomb-type contraction for any auxiliary index R to calculate one row of Z using the matrix of transformed three-center integrals

$\underline{\overline{\bf M}}^{(R)}$

{\underset{̲}{\bar{M}}}^{(R)}

in place of the density matrix.

In the regular algorithm to determine the Coulomb matrix for SCF calculations, the ket shell-pairs are contracted with a single density matrix only. Here, however, N_aux density-like matrices have to be contracted, so that the use of the conventional algorithm would show the same large prefactor and some adjustments of the algorithm were required.

Crucial to the performance of the algorithm is the use of shared memory to hold the intermediate values J_[pq], however, this memory is strictly limited on GPUs. Thus, we decided to retain the two-dimensional thread-block setup but we extend the grid by a further dimension to represent the offset to the corresponding indices of the current R-batch. Therefore, we strongly reduce the communication between host and GPU-devices. It should be noted that in this case the bra-data does not represent shell-pairs but only the shells of the auxilliary basis.

The complete algorithm is also depicted in Fig. 1. Note that the elements of

$(R|\underline{\mu }\overline{\nu })$

(R | \underset{̲}{μ} \bar{ν})

are already transformed to the cartesian basis and non-axial normalization coefficients have been incorporated. Furthermore, the elements are sorted according to the l-quantum number combination of

$\underline{\mu }\overline{\nu }$

\underset{̲}{μ} \bar{ν}

⁠, so that the batches can be read from disk in a coalesced fashion.

FIG. 1.

FIG. 1. Scheme of the algorithm for the contraction step Eq. (10). If nothing else is indicated, the steps are executed on CPU.

View large Download slide

Scheme of the algorithm for the contraction step Eq. (10). If nothing else is indicated, the steps are executed on CPU.

IV. CONCLUSION

We present a reformulation of the SOS-RI-MP2 algorithm which not only reduces the computational effort to

${\cal O\mbox{(N^{3})}}$

O (N^{3})

⁠, but also allows for an efficient evaluation of the rate-determining contraction step on GPUs by employing a modified J-engine algorithm. As a – to our knowledge – first reformulation of a post-HF method which is highly suitable for GPUs beyond the mere use of GPU-accelerated linear algebra libraries, our ansatz may also be applicable to other post-HF methods in order to enable an efficient use of massively parallel processors.

ACKNOWLEDGMENTS

C.O. acknowledges financial support by the “Deutsche Forschungsgemeinschaft” (DFG) for the project Oc35/4-1. Further funding was provided by the DFG via the initiatives SFB 749 “Dynamik und Intermediate molekularer Transformationen” (project C7) and the DFG cluster of excellence EXC 114 “Center for Integrated Protein Science Munich” (CIPSM).

REFERENCES

1.

C.

Møller

and

M. S.

Plesset

,

Phys. Rev.

46

,

618

(

1934

).

https://doi.org/10.1103/PhysRev.46.618

2.

J.

Almlöf

,

Chem. Phys. Lett.

181

,

319

(

1991

).

https://doi.org/10.1016/0009-2614(91)80078-C

3.

M.

Häser

and

J.

Almlöf

,

J. Chem. Phys.

96

,

489

(

1992

).

https://doi.org/10.1063/1.462485

4.

M.

Häser

,

Theor. Chim. Acta

87

,

147

(

1993

).

https://doi.org/10.1007/BF01113535

5.

P. Y.

Ayala

and

G. E.

Scuseria

,

J. Chem. Phys.

110

,

3660

(

1999

).

https://doi.org/10.1063/1.478256

6.

S.

Saebø

and

P.

Pulay

,

J. Chem. Phys.

115

,

3975

(

2001

).

https://doi.org/10.1063/1.1389291

7.

H.-J.

Werner

,

F. R.

Manby

, and

P. J.

Knowles

,

J. Chem. Phys.

118

,

8149

(

2003

).

https://doi.org/10.1063/1.1564816

8.

B.

Doser

,

D. S.

Lambrecht

,

J.

Kussmann

, and

C.

Ochsenfeld

,

J. Chem. Phys.

130

,

064107

(

2009

).

https://doi.org/10.1063/1.3072903

9.

Y.

Jung

,

R. C.

Lochan

,

A. D.

Dutoi

, and

M.

Head-Gordon

,

J. Chem. Phys.

121

,

9793

(

2004

).

https://doi.org/10.1063/1.1809602

10.

Y.

Jung

,

Y.

Shao

, and

M.

Head-Gordon

,

J. Comput. Chem.

28

,

1953

(

2007

).

https://doi.org/10.1002/jcc.20590

11.

R.

Olivares-Amaya

,

M. A.

Watson

,

R. G.

Edgar

,

L.

Vogt

,

Y.

Shao

, and

A.

Aspuru-Guzik

,

J. Chem. Theory Comput.

6

,

135

(

2010

).

https://doi.org/10.1021/ct900543q

12.

E. A.

DePrince

III and

J. R.

Hammond

,

J. Chem. Theory Comput.

7

,

1287

(

2011

).

https://doi.org/10.1021/ct100584w

13.

W.

Ma

,

S.

Krishnamoorthy

,

O.

Villa

, and

K.

Kowalski

,

J. Chem. Theory Comput.

7

,

1316

(

2011

).

https://doi.org/10.1021/ct1007247

14.

K.

Bhaskaran-Nair

,

W.

Ma

,

S.

Krishnamoorthy

,

O.

Villa

,

H. J. J.

van Dam

,

E.

Aprà

, and

K.

Kowalski

,

J. Chem. Theory Comput.

9

,

1949

(

2013

).

https://doi.org/10.1021/ct301130u

15.

M.

Del Ben

,

J.

Hutter

, and

J.

VandeVondele

,

J. Chem. Theory Comput.

9

,

2654

(

2013

).

https://doi.org/10.1021/ct4002202

16.

I. S.

Ufimtsev

and

T. J.

Martinez

,

J. Chem. Theory Comput.

5

,

1004

(

2009

).

https://doi.org/10.1021/ct800526s

17.

J.

Kussmann

and

C.

Ochsenfeld

,

J. Chem. Phys.

138

,

134114

(

2013

).

https://doi.org/10.1063/1.4796441

18.

S. A.

Maurer

,

D. S.

Lambrecht

,

J.

Kussmann

, and

C.

Ochsenfeld

,

J. Chem. Phys.

138

,

014101

(

2013

).

https://doi.org/10.1063/1.4770502

19.

S.

Maurer

,

L.

Clin

, and

C.

Ochsenfeld

,

J. Chem. Phys.

140

,

224112

(

2014

).

https://doi.org/10.1063/1.4881144

20.

F. G.

Gustavson

,

ACM Trans. Math. Softw.

4

,

250

(

1978

).

https://doi.org/10.1145/355791.355796

21.

J.

Kussmann

and

C.

Ochsenfeld

,

J. Chem. Phys.

127

,

054103

(

2007

).

https://doi.org/10.1063/1.2749509

22.

Development version of the Q-Chem program package, http://www.q-chem.com.

23.

C. A.

White

and

M.

Head-Gordon

,

J. Chem. Phys.

104

,

2620

(

1996

).

https://doi.org/10.1063/1.470986

24.

G. R.

Ahmadi

and

J.

Almlöf

,

Chem. Phys. Lett.

246

,

364

(

1995

).

https://doi.org/10.1016/0009-2614(95)01127-4

25.

I. S.

Ufimtsev

and

T. J.

Martínez

,

J. Chem. Theory Comput.

4

,

222

(

2008

).

https://doi.org/10.1021/ct700268q

© 2014 AIP Publishing LLC.

2014

AIP Publishing LLC