library for evaluating polarization integrals according to Schwerdtfeger et al. (1988) (Reference [CPP]).

Polarization Integrals

The polarization integrals are defined as

$$\langle \text{CGTO}_i \vert \frac{x^{m_x} y^{m_y} z^{m_z}}{r^k} \left(1 - e^{-\alpha r^2} \right)^q \vert \text{CGTO}_j \rangle$$

between unnormalized primitive Cartesian Gaussian functions

CGTO_i(x, y, z) = (x − x_i)^n_xi(y − y_i)^n_yi(z − z_i)^n_ziexp (−β_i(r⃗−r⃗_i)²)

and

CGTO_j(x, y, z) = (x − x_j)^n_xj(y − y_j)^n_yj(z − z_j)^n_zjexp (−β_j(r⃗−r⃗_j)²)

for k > 2. The power of the cutoff function q has to satisfy

q ≥ κ(k) − κ(m_x) − κ(m_y) − κ(m_z) − 1

where

$$\begin{aligned} \kappa = \begin{cases} \frac{n}{2} \quad \text{, if n is even } \\\ \frac{n+1}{2} \quad \text{, if n is odd } \end{cases} \end{aligned}$$

Requirements

• python3
• pybind11

For parallel evaluation of integrals on a GPU you need in addition

• CUDA Toolkit (nvcc)
• a CUDA supported GPU device

Getting Started

The package is installed by running

$ pip install -e .

Tests

To check the proper functioning of the code, it is recommended to run a set of tests. In order to compute the numerical integrals needed for comparison one first has to install the python package becke from https://github.com/humeniuka/becke_multicenter_integration . Then the test suite is run with

$ cd tests
$ python -m unittest

Example

The following code evaluates the matrix elements of the operator r^− 4 between an unnormalized s-orbital at the origin and a shell of unnnormalized p-orbitals at the point r⃗_j = (0.0, 0.0, 1.0).

First a PolarizationIntegral is created, which needs to know the centers, angular momenta and radial exponents of the two shells, as well as the polarization operator and the cutoff function:

from polarization_integrals import PolarizationIntegral

# s- and p-shell
li, lj = 0, 1
# Op = r^{-4}
k, mx,my,mz = 4, 0,0,0
# cutoff function
alpha = 50.0
q = 4

I = PolarizationIntegral(0.0, 0.0, 0.0,  li,  0.5,  
                         0.0, 0.0, 1.0,  lj,  0.5,
                         k, mx,my,mz,
                         alpha, q)

Then the integrals can be evaluated for all combinations of unnormalized primitives from each shell.

# <s|Op|px>
print( I.compute_pair(0,0,0,  1,0,0) )
# <s|Op|py>
print( I.compute_pair(0,0,0,  0,1,0) )
# <s|Op|pz>
print( I.compute_pair(0,0,0,  0,0,1) )

GPU Support

Polarization integrals can be calculated in parallel on a GPU which supports CUDA. The kernels with python binding are located in the folder src_gpu/. The kernels and python wrapper are compiled with

$ cd src_gpu
$ make

The correctness of the GPU integrals should be verified by comparison with the CPU implementation by running a set of tests with

$ cd tests_gpu
$ python -m unittest

References

[CPP] P. Schwerdtfeger, H. Silberbach,

'Multicenter integrals over long-range operators using Cartesian Gaussian functions', Phys. Rev. A 37, 2834 https://doi.org/10.1103/PhysRevA.37.2834

[CPP-Erratum] Phys. Rev. A 42, 665

https://doi.org/10.1103/PhysRevA.42.665

[CPP-Erratum2] Phys. Rev. A 103, 069901

https://doi.org/10.1103/PhysRevA.103.069901

[library] A. Humeniuk, W. Glover,

'Efficient CPU and GPU implementations of multicenter integrals over long-range operators using Cartesian Gaussian functions', submitted

Files

The directory tree below depicts the structure of the source code package:

│   setup.py                             script for installing python package
│   README.rst                           instructions on installation and usage
│   LICENSE.txt                          MIT license
│
└───polarization_integrals               python package
│   │   __init__.py
│
└───src
│   │   export.cc                      - python wrapper around CPU implementation
│   │   Faddeeva.hh                    - error and Dawson function for x >= 6.0
│   │   Faddeeva.cc
│   │   polarization.h                 - definition of PolarizationIntegral class
│   │   polarization.cc                  CPU implementation of polarization integrals
│   │   radial_overlap.h               - analytical overlap integrals between squares
│   │   radial_overlap.cc                of radial Gaussian basis functions
│   │   Makefile                       - compile python module for CPU integrals
│   │
│   
└───src_gpu
│   │   export.cc                      - python wrapper around GPU implementation
│   │   Dawson.cu                      - implementation of Dawson function for x >= 6.0
│   │   Dawson_real.cu                 - automatically generated from Dawson.cu
│   │   polarization.h                 - definition of Primitive, PrimitivePair and PolarizationIntegral classes
│   │   polarization.cu                - GPU implementation of polarization integrals
│   │   double_to_real_cast.sh         - convenience script for replacing numerical constants in source code
│   │   Makefile                       - compile python module for GPU integrals with single and double precision
│   │   code_generator.py              - Parts of the C++ code were generated automatically using this script.
│   │   README.txt                     - instructions on profiling GPU implementation
│   │   test.cu                        - example for using GPU code
│   │   run_tests.sh                   - run and benchmark GPU code in single and double precision
│
└───tests
│   │   example.py                     - simple example demonstrating how to use the library
│   │   make_table_I.py                - produce table I in the article
│   │   polarization_ints_numerical.py - numerical integrals using Becke's scheme
│   │   polarization_ints_reference.py - pure python implementation of some integrals
│   │   upper_bounds.py                - numerical integrals of right-hand side of Cauchy-Schwarz inequality
│   │   radial_overlap.py              - numerical overlaps using Becke's scheme
│   │   test_###.py                    - different tests
│   │   ...
│
└───tests_gpu
│   │   errors.py                      - comparison of numerical errors between CPU and GPU implementations
│   │   test_integrals_gpu_sp.py       - test cases in single precision
│   │   test_integrals_gpu_dp.py       - test cases double precision