Package to handle integrals over Gaussian-type atomic orbitals.
Author FermiQC
22 Stars
Updated Last
7 Months Ago
Started In
August 2021

CI Coverage License

GaussianBasis offers high-level utilities for molecular integral computations.

Current features include:

  • Basis set parsing (gbs format)
  • Standard basis set files from BSE
  • One-electron integral (1e)
  • Two-electron two-center integral (2e2c)
  • Two-electrons three-center integral (2e3c)
  • Two-electrons four-center integral (2e4c)
  • Gradients (currrently under construction - watch out!)

Integral computations use by default the integral library libcint via libcint_jll.jl. A simple Julia-written integral module Acsint.jl is also available, but it is significantly slower than the libcint.

Basic Usage

The simplest way to use the code is by first creating a BasisSet object. For example

julia> bset = BasisSet("sto-3g", """
              H        0.00      0.00     0.00                 
              H        0.76      0.00     0.00""")
sto-3g Basis Set
Type: Spherical   Backend: Libcint
Number of shells: 2
Number of basis:  2

H: 1s 
H: 1s

Next, call the desired integral function with the BasisSet object as the argument. Let's take the overlap function as an example:

julia> overlap(bset)
2×2 Matrix{Float64}:
 1.0       0.646804
 0.646804  1.0
Function Description Formula
overlap Overlap between two basis functions S
kinetic Kinetic integral T
nuclear Nuclear attraction integral V
ERI_2e4c Electron repulsion integral - returns a full rank-4 tensor! ERI
sparseERI_2e4c Electron repulsion integral - returns non-zero elements along with a index tuple sERI
ERI_2e3c Electron repulsion integral over three centers. Note: this function requires another basis set as the second argument (that is the auxiliary basis set in Density Fitting). It must be called as ERI_2c3c(bset, aux) 3cERI
ERI_2e2c Electron repulsion integral over two centers 2cERI

Advanced Usage

Basis Functions

BasisFunction object is the central data type within this package. Here, BasisFunction is an abstract type with two concrete structures: SphericalShell and CartesianShell. By default SphericalShell is created. In general a spherical basis function is


where the sum goes over primitive functions. A BasisFunction object contains the data to reproduce the mathematical object, i.e. the angular momentum number (l), expansion coefficients (cn), and exponential factors (ξn). We can create a basis function by passing these arguments orderly:

julia> using StaticArrays
julia> atom = GaussianBasis.Atom(8, 16.0, [1.0, 0.0, 0.0])
julia> bf = BasisFunction(1, SVector(1/√2, 1/√2), SVector(5.0, 1.2), atom)
P shell with 3 basis built from 2 primitive gaussians

χ₁₋₁ =    0.7071067812Y₁₋₁exp(-5.0r²)
     +    0.7071067812Y₁₋₁exp(-1.2r²)

χ₁₀  =    0.7071067812Y₁₀exp(-5.0r²)
     +    0.7071067812Y₁₀exp(-1.2r²)

χ₁₁  =    0.7071067812Y₁₁exp(-5.0r²)
     +    0.7071067812Y₁₁exp(-1.2r²)

We can now check the fields (attributes):

julia> bf.l

julia> bf.coef
2-element SVector{2, Float64} with indices SOneTo(2):

julia> bf.exp
2-element SVector{2, Float64} with indices SOneTo(2):

Note that exp and coef are expected to be SVector from StaticArrays.

Basis Set

The BasisSet object is the main ingredient for integrals. It can be created in a number of ways:

  • The highest level approach takes two strings as arguments, one for the basis set name and another for the XYZ file. See Basic Usage.

  • You can pass your vector of Atom structures instead of an XYZ string as the second argument. GaussianBasis uses the Atom structure from Molecules.jl.

atoms = GaussianBasis.parse_string("""
            H        0.00      0.00     0.00                 
            H        0.76      0.00     0.00""")
BasisSet("sto-3g", atoms)
  • Finally, instead of searching into GaussianBasis/lib for a basis set file matching the desired name, you can construct your own from scratch. We further discuss this approach below.

Basis sets are mainly composed of two arrays: a vector of atoms and a vector of basis functions objects. We can construct both manually for maximum flexibility:

julia> h2 = GaussianBasis.parse_string(
  "H 0.0 0.0 0.0
   H 0.0 0.0 0.7"
2-element Vector{Atom{Int16, Float64}}:
Atom{Int16, Float64}(1, 1.008, [0.0, 0.0, 0.0])
Atom{Int16, Float64}(1, 1.008, [0.0, 0.0, 0.7])

Next, we create a vector of basis functions.

julia> shells = [BasisFunction(0, SVector(0.5215367271), SVector(0.122), h2[1]),
BasisFunction(0, SVector(0.5215367271), SVector(0.122), h2[2]),
BasisFunction(1, SVector(1.9584045349), SVector(0.727), h2[2])];

Finally, we create the basis set object. Note that, you got to make sure your procedure is consistent. The atoms used to construct the basis set object must be in the atom vector, otherwise unexpected results may arise.

julia> bset = BasisSet("UnequalHydrogens", h2, shells)
UnequalHydrogens Basis Set
Type: Spherical{Molecules.Atom, 1, Float64}   Backend: Libcint
Number of shells: 3
Number of basis:  5

H: 1s 
H: 1s 1p

The most import fields here are:

julia> == "UnequalHydrogens"
julia> bset.basis == shells 
julia> bset.atoms == h2

Integrals over different basis sets

Functions such as ERI_2e3c require two basis set as arguments. Looking at the corresponding equation 3cERI we see two basis set: Χ and P. If your first basis set has 2 basis functions and the second has 4, your output array is a 2x2x4 tensor. For example

julia> b1 = BasisSet("sto-3g", """
              H        0.00      0.00     0.00                 
              H        0.76      0.00     0.00""")
julia> b2 = BasisSet("3-21g", """
              H        0.00      0.00     0.00                 
              H        0.76      0.00     0.00""")
julia> ERI_2e3c(b1,b2)
2×2×4 Array{Float64, 3}:
[:, :, 1] =
 3.26737  1.85666
 1.85666  2.44615

[:, :, 2] =
 6.18932  3.83049
 3.83049  5.60161

[:, :, 3] =
 2.44615  1.85666
 1.85666  3.26737

[:, :, 4] =
 5.60161  3.83049
 3.83049  6.18932

One electron integrals can also be employed with different basis set.

julia> overlap(b1, b2)
2×4 Matrix{Float64}:
 0.914077  0.899458  0.473201  0.708339
 0.473201  0.708339  0.914077  0.899458

julia> kinetic(b1, b2)
2×4 Matrix{Float64}:
 1.03401  0.314867  0.20091  0.203163
 0.20091  0.203163  1.03401  0.314867

This can be useful when working with projections from one basis set onto another.

Computing integrals element-wise

For all integrals, you can get the full array by using the general syntax integral(basisset) (e.g. overlap(bset) or ERI_2e4c(bset)). Alternatively, you can specify a shell combination for which the integral must be computed

julia> ERI_2e4c(b1, 1,2,2,1)
1×1×1×1 Array{Float64, 4}:
[:, :, 1, 1] =

julia> kinetic(b1, 1,2)
1×1 Matrix{Float64}:

Mutating versions of the functions are also available

julia> S = zeros(2,2);
julia> overlap!(S, b1)
julia> S
2×2 Matrix{Float64}:
 1.0       0.646804
 0.646804  1.0

Used By Packages