CI  Coverage  License 

GaussianBasis offers highlevel utilities for molecular integral computations.
Current features include:
 Basis set parsing (
gbs
format)  Standard basis set files from BSE
 Oneelectron integral (1e)
 Twoelectron twocenter integral (2e2c)
 Twoelectrons threecenter integral (2e3c)
 Twoelectrons fourcenter integral (2e4c)
 Gradients (currrently under construction  watch out!)
Integral computations use by default the integral library libcint via libcint_jll.jl. A simple Juliawritten 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("sto3g", """
H 0.00 0.00 0.00
H 0.76 0.00 0.00""")
sto3g 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  
kinetic 
Kinetic integral  
nuclear 
Nuclear attraction integral  
ERI_2e4c 
Electron repulsion integral  returns a full rank4 tensor!  
sparseERI_2e4c 
Electron repulsion integral  returns nonzero elements along with a index tuple  
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) 

ERI_2e2c 
Electron repulsion integral over two centers 
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 (c_{n}), 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.7071067812⋅Y₁₋₁⋅r¹⋅exp(5.0⋅r²)
+ 0.7071067812⋅Y₁₋₁⋅r¹⋅exp(1.2⋅r²)
χ₁₀ = 0.7071067812⋅Y₁₀⋅r¹⋅exp(5.0⋅r²)
+ 0.7071067812⋅Y₁₀⋅r¹⋅exp(1.2⋅r²)
χ₁₁ = 0.7071067812⋅Y₁₁⋅r¹⋅exp(5.0⋅r²)
+ 0.7071067812⋅Y₁₁⋅r¹⋅exp(1.2⋅r²)
We can now check the fields (attributes):
julia> bf.l
1
julia> bf.coef
2element SVector{2, Float64} with indices SOneTo(2):
0.7071067811865475
0.7071067811865475
julia> bf.exp
2element SVector{2, Float64} with indices SOneTo(2):
5.0
1.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 theAtom
structure from Molecules.jl.
atoms = GaussianBasis.parse_string("""
H 0.00 0.00 0.00
H 0.76 0.00 0.00""")
BasisSet("sto3g", 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"
)
2element 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> bset.name == "UnequalHydrogens"
true
julia> bset.basis == shells
true
julia> bset.atoms == h2
true
Integrals over different basis sets
Functions such as ERI_2e3c
require two basis set as arguments. Looking at the corresponding equation
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("sto3g", """
H 0.00 0.00 0.00
H 0.76 0.00 0.00""")
julia> b2 = BasisSet("321g", """
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 elementwise
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] =
0.2845189435761272
julia> kinetic(b1, 1,2)
1×1 Matrix{Float64}:
0.2252049038643092
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