GaussianBasis.jl

GaussianBasis.jl evaluates molecular integrals (and their nuclear-coordinate derivatives) over contracted Gaussian-type atomic orbitals, using libcint as its default backend. It is the integral engine behind Fermi.jl.

Installation

using Pkg
Pkg.add("GaussianBasis")

Quickstart

A BasisSet pairs a set of atoms with a set of contracted Gaussian basis functions read from a standard basis set library (e.g. sto-3g, cc-pvdz):

julia> using GaussianBasis

julia> bset = BasisSet("sto-3g", """
       O        0.000000000000     -0.143225816552      0.000000000000
       H        1.638036840407      1.136548822547     -0.000000000000
       H       -1.638036840407      1.136548822547     -0.000000000000
       """);

sto-3g Basis Set
Type: Spherical   Backend: Libcint

Number of shells: 5
Number of basis:  7

O: 1s 2s 1p
H: 1s
H: 1s

Geometries are Cartesian coordinates in Angstrom. From a BasisSet, you can compute:

julia> S = overlap(bset);
7×7 Matrix{Float64}:
 1.0         0.236704    0.0        0.0        0.0  0.00410862   0.00410862
 0.236704    1.0         0.0        0.0        0.0  0.0644883    0.0644883
 0.0         0.0         1.0        0.0        0.0  0.0572785   -0.0572785
 0.0         0.0         0.0        1.0        0.0  0.0447509    0.0447509
 0.0         0.0         0.0        0.0        1.0  0.0          0.0
 0.00410862  0.0644883   0.0572785  0.0447509  0.0  1.0          0.0100209
 0.00410862  0.0644883  -0.0572785  0.0447509  0.0  0.0100209    1.0

Electron Repulsion Integrals (ERI) are returned as full $N^4$ arrays.

julia> ERI_2e4c(bset) |> size
(7, 7, 7, 7)

To accommodate memory sensitive operations, such as integral direct (on-the-fly) methods, integrals can be computed in shell batches of size $(2l+1)^4$, where $l$ is the angular momentum of the shell.

The formula above is only valid for Spherical basis set. For more information check CartesianShell.

julia> ERI_2e4c(bset, 4,4,4,3) # Note that those are shell indexes, not basis.
1×1×1×3 Array{Float64, 4}:
[:, :, 1, 1] =
 0.026307122405698283

[:, :, 1, 2] =
 0.02055337661033345

[:, :, 1, 3] =
 0.0

For an even more allocation friendly option, check out ERI_2e4c!

Basis Set

A basis set is a collection of basis functions that, in turn, represent one-electron wave functions for an electron bound to a nucleus. A BasisSet object can be conveniently created using a string of molecular coordinates and the name of a standard basis set. The available basis sets can be found in lib/.

julia> BasisSet("cc-pvdz", """
              F        0.000000000000     -0.143225816552      0.000000000000
              H        1.638036840407      1.136548822547     -0.000000000000
""")
cc-pvdz Basis Set
Type: Spherical   Backend: Libcint

Number of shells: 9
Number of basis:  19

F: 1s 2s 3s 1p 2p 1d 
H: 1s 2s 1p

You may add new basis sets by simply creating the .gbs file inside lib/. These files can be found in the Basis Set Exchange. GaussianBasis.jl uses the Psi4 format.

In GaussianBasis.jl, a BasisSet object is the fundamental object necessary to compute integrals. In simple terms, a basis set is simply a container of Shells (see SphericalShell and CartesianShell), which represent a group of similar wave functions with the same $l$ value but different $m_l$ values.

julia> for b in bset
           println(b)
       end
S shell on Fluorine (1 basis function, 9 primitives)
S shell on Fluorine (1 basis function, 9 primitives)
S shell on Fluorine (1 basis function, 1 primitive)
P shell on Fluorine (3 basis functions, 4 primitives)
P shell on Fluorine (3 basis functions, 1 primitive)
D shell on Fluorine (5 basis functions, 1 primitive)
S shell on Hydrogen (1 basis function, 4 primitives)
S shell on Hydrogen (1 basis function, 1 primitive)
P shell on Hydrogen (3 basis functions, 1 primitive)

There is no individual basis function object. That is because integrals are calculated in shell batches. It is important to pay attention to the difference between number of basis and number of shells. Likewise, the index of a shell and the index of a basis are fundamentally different.

GaussianBasis.BasisSet — Type
BasisSet

Object holding a set of ShellFunction objects associated with an array of atoms.

Fields

NameTypeDescription
atomsVector{Atom}An array of Molecules.Atom objects
nameStringString holding the basis set name
shellsVector{Vector{ShellFunction}}An Array of arrays with ShellFunction
natomsInt32Number of atoms in the BasisSet
nbasInt32Number of basis functions.
nshellsInt32Number of shells, i.e. ShellFunction objects
lc_atomsArray{Int32,1}Integer array mapping data to libcint
lc_basArray{Int32,1}Integer array mapping data to libcint
lc_envArray{Float64,1}Float64 array mapping data to libcint

Example

Build a basis set from default options

julia> water = """ 
 O        1.2091536548      1.7664118189     -0.0171613972
 H        2.1984800075      1.7977100627      0.0121161719
 H        0.9197881882      2.4580185570      0.6297938832"""

julia> bset = BasisSet("sto-3g", water)
sto-3g Basis Set
Type: Spherical   Backend: Libcint

Number of shells: 5
Number of basis:  7

O: 1s 2s 1p 
H: 1s 
H: 1s

The BasisSet object can be accessed as vector, returning the n-th shell.

julia> bset[1] # Show the first shell
S shell on Oxygen at position [1.2091536548, 1.7664118189, -0.0171613972] Å
Contains 1 basis function built from 3 primitive gaussians

χ₀₀  =    0.8486970052⋅Y₀₀⋅exp(-5.033151319⋅r²)
     +    1.1352008076⋅Y₀₀⋅exp(-1.169596125⋅r²)
     +    0.8567529838⋅Y₀₀⋅exp(-0.38038896⋅r²)

You can also create your own crazy mix! Let us create one S and one P basis functions for H

julia> H1 = GaussianBasis.Atom(:H, [0.0, 0.0, 0.0]);
julia> H2 = GaussianBasis.Atom(:H, [0.0, 0.0, 0.7]);
julia> s = ShellFunction(0, [0.5215367271], [0.122], H1)
julia> p = ShellFunction(1, [1.9584045349], [0.727], H2)

The basis set is constructed with an array of atoms (Vector{Atom}) and a corresponding array of Vector{ShellFunction} holding all basis functions for that particular atom. In this example, we consider an unequal treatment for the two atoms in the H₂ molecule.

julia> shells = [s, p]  # One s function on the first hydrogen
                        # One p function on the second hydrogen
2-element Vector{SphericalShell{Molecules.Atom{Float64, Float64}}}:
 S shell on Hydrogen (1 basis function, 1 primitive)
 P shell on Hydrogen (3 basis functions, 1 primitive)
julia> BasisSet("UnequalHydrogens", shells)
UnequalHydrogens Basis Set
Type: Spherical   Backend: Libcint

Number of shells: 2
Number of basis:  4

H: 1s 
H: 1p

Note that the number of shells and the number of basis functions are different. Integrals are computed over shells, while the number of basis functions is the actual number of functions in the basis set and dictates the size of your integral arrays.

source
GaussianBasis.SphericalShell — Type
 SphericalShell

Object representing a shell of basis functions in Spherical coordinates. The shells types are s, p, d, ... corresponding to angular momentum numbers l = 0, 1, 2, ... For each shell, there are 2l+1 basis functions corresponding to the magnetic quantum number m = -l, ..., l.

Fields

NameDescription
lAngular momentum number
coefArray holding expansion coefficients for the primitive functions
expArray holding exponents for primitive functions

Examples

julia> atom = GaussianBasis.Atom(1, 1.008, [0.0, 0.0, 0.0])
julia> p_shell = SphericalShell(1, [1/√2, 1/√2], [5.0, 1.2], atom)
P shell on Hydrogen at position [0.0, 0.0, 0.0] Å
Contains 3 basis functions 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²)
source
GaussianBasis.CartesianShell — Type
 CartesianShell

Object representing a shell of basis functions in Cartesian coordinates. The shells types are s, p, d, ... corresponding to angular momentum numbers l = 0, 1, 2, ... Note that the number of basis functions in a Cartesian shell is given by (l+1)(l+2)/2, which is larger than the number of basis functions in a Spherical shell (2l+1) for l > 1. That is, m_l values are not well defined in Cartesian shells.

Fields

NameDescription
lPseudo-angular momentum number l = a+b+c for xᵃyᵇzᶜ
coefArray holding expansion coefficients for the primitive functions
expArray holding exponents for primitive functions

Examples

julia> atom = GaussianBasis.Atom(1, 1.008, [0.0, 0.0, 0.0])
julia> d_shell = CartesianShell(2, [1/√2], [5.0], atom)
D shell on Hydrogen at position [0.0, 0.0, 0.0] Å
Contains 6 basis functions built from 1 primitive gaussian

χ(x²) =    0.7071067812⋅x²⋅exp(-5.0⋅r²)

χ(xy) =    0.7071067812⋅xy⋅exp(-5.0⋅r²)

χ(xz) =    0.7071067812⋅xz⋅exp(-5.0⋅r²)

χ(y²) =    0.7071067812⋅y²⋅exp(-5.0⋅r²)

χ(yz) =    0.7071067812⋅yz⋅exp(-5.0⋅r²)

χ(z²) =    0.7071067812⋅z²⋅exp(-5.0⋅r²)
source
GaussianBasis.atomic_orbital_amplitude — Function
atomic_orbital_amplitude(B::BasisSet, i::Integer, r::AbstractVector) -> Real
atomic_orbital_amplitude(B::BasisSet, i::Integer, r::AbstractArray) -> AbstractArray

Returns the amplitude of the ith basis function at a set of real space positions. r can be either a 3-vector containing a single position or a 3 × … array of positions for evaluating many points at once more efficiently.

Examples

julia> bset = BasisSet("sto-3g", "H 0 0 0");
julia> atomic_orbital_amplitude(bset, 1, [0.0,0.0,0.0])
0.6282468778403579
julia> atomic_orbital_amplitude(bset, 1, [0.1;0.2;0.3;;-0.1;0.3;-0.2])
2-element Vector{Float64}:
 0.49840190793869554
 0.49840190793869554
source
GaussianBasis.Libcint — Module
Libcint

Minimal wrap around the integral library libcint. This module exposes libcint functions to the Julia interface.

source
GaussianBasis.on_atom_flags — Function
on_atom_flags(BS::BasisSet, iA::Int, shells...)

Whether each of shells (any number of shell indices) sits on atom iA, as an NTuple{N,Bool}.

source
GaussianBasis.merge_basis — Function
merge_basis(BS1::BasisSet, BS2::BasisSet) -> BasisSet

The concatenation of BS1 and BS2 into one BasisSet – regular shells first, auxiliary second – as libcint's 3-center kernels require, since they resolve all three shell indices against a single basis. An auxiliary shell k of BS2 is therefore addressed as k + BS1.nshells in the result.

Depends only on BS1/BS2, never on atoms or shells, so routines that take a Bmerged keyword should be handed one built once outside the loop.

source

Contributing

Contributions in the form of bug reports, feature requests, and pull requests are all welcome.

Reporting Issues and Requesting Features

If you run into a bug, are looking for a feature, or simply need some help with GaussianBasis.jl, feel free to open an issue on GitHub. You can also reach out directly to Gustavo Aroeira.

Development Setup

Clone the repository and instantiate its environment:

julia --project=. -e 'using Pkg; Pkg.instantiate()'

Running Tests

From the package root:

julia --project=. -e 'using Pkg; Pkg.test()'

or, from a Julia REPL started with --project=.:

julia> ]
pkg> test

Building the Docs

The docs live under docs/ and are built with Documenter.jl. To build and preview them locally:

julia --project=docs -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
julia --project=docs docs/make.jl

This generates docs/build/; open docs/build/index.html in a browser to preview your changes.

Submitting a Pull Request

  • Branch off main.
  • Make sure Pkg.test() passes locally – CI runs the same suite on Julia 1.9 and the latest 1.x release, on Linux and macOS.
  • Keep PRs minimal, do not bundle several unrelated changes into a single PR.
  • There's no enforced code formatter for this repo, so just match the style of the surrounding code.
  • Use of AI is allowed, but you must have reviewed the entire code and assured its quality.
  • Your new code must be tested. Make sure to wire your test suite into test/runtests.jl.
  • At the minimum, you should provide docstrings for your user-facing functions.