One-Electron Integrals

Overlap

For the atomic orbitals $\chi_\mu, \chi_\nu$, the overlap integral is calculated as:

\[S_{\mu\nu} = \int \chi_\mu(\mathbf{r})\, \chi_\nu(\mathbf{r})\, d\mathbf{r}\]

A full-dense overlap matrix can be calculated directly from a basis set

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
              """)
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

Alternatively, an in-place (mutating) version can be used:

julia> out = zeros(bset.nbas, bset.nbas)
julia> overlap!(out, bset)
julia> out == S
true

A per-shell option is also available:

julia> S43 = overlap(bset, 4, 3) # 4 and 3 are shell indexes
1×3 Matrix{Float64}:
 0.0572785  0.0447509  0.0

For maximum efficiency, you may need to use the most primitive call, mutating and using a pre-allocated output

julia> out = zeros(num_basis(bset[4]), num_basis(bset[3]))
julia> overlap!(out, bset, 4, 3) 
julia> out == S43
true

The overlap can also be computed using two basis set. This will calculate $\langle \chi_\mu | \chi_\nu \rangle$ for $\mu \in A$ and $\nu \in B$, where $A$ and $B$ are two different basis set. This can be useful when projecting results from one basis onto another.

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

julia> overlap(bset, bset2)
7×13 Matrix{Float64}:
 0.998058    0.255215   0.0        0.0        0.0       0.207535   …  0.0      3.07065e-6  0.00788467   3.07065e-6   0.00788467
 0.205853    0.853358   0.0        0.0        0.0       0.980775      0.0      0.0105854   0.112372     0.0105854    0.112372
 0.0         0.0        0.886749   0.0        0.0       0.0           0.0      0.0161306   0.0936416   -0.0161306   -0.0936416
 0.0         0.0        0.0        0.886749   0.0       0.0           0.0      0.0126026   0.0731609    0.0126026    0.0731609
 0.0         0.0        0.0        0.0        0.886749  0.0           0.83337  0.0         0.0          0.0          0.0
 0.00354245  0.021765   0.0137339  0.0107301  0.0       0.0800824  …  0.0      0.914077    0.899458     0.00124358   0.0171688
 0.00354245  0.021765  -0.0137339  0.0107301  0.0       0.0800824     0.0      0.00124358  0.0171688    0.914077     0.899458

GaussianBasis.overlap — Function

overlap(BS::BasisSet) -> Matrix{Float64}
overlap(BS::BasisSet, i, j) -> Matrix{Float64}
overlap(BS1::BasisSet, BS2::BasisSet) -> Matrix{Float64}
overlap(BS1::BasisSet, BS2::BasisSet, i, j) -> Matrix{Float64}

Compute the AO overlap matrix S_{μν} = ⟨χ_μ|χ_ν⟩.

Methods

overlap(BS): full, dense, symmetric nbas × nbas matrix for BS.
overlap(BS, i, j): just the (Ni,Nj) block for shells i and j of BS (shell indices, not AO indices).
overlap(BS1, BS2): mixed-basis matrix, μ running over BS1 and ν over BS2 – a dense BS1.nbas × BS2.nbas matrix, not generally symmetric. Useful e.g. for projecting a density or orbitals from one basis onto another.
overlap(BS1, BS2, i, j): mixed-basis block between shell i of BS1 and shell j of BS2.

For repeated calls (e.g. in a hot loop or a geometry scan), see overlap!, which writes into a preallocated array instead of allocating.

source

GaussianBasis.overlap! — Function

overlap!(out, BS::BasisSet, i, j)
overlap!(out, BS::BasisSet)
overlap!(out, BS1::BasisSet, BS2::BasisSet, i, j)
overlap!(out, BS1::BasisSet, BS2::BasisSet)

Mutating counterpart of overlap: writes into the caller-supplied out instead of allocating. This shell-pair, single-basis form is the primitive every other method builds on.

Methods

overlap!(out, BS, i, j): out must be (Ni,Nj), the block for shells i/j of BS (shell indices, not AO indices).
overlap!(out, BS): out must be a dense nbas × nbas matrix.
overlap!(out, BS1, BS2, i, j): out must be (Ni,Nj), the block for shell i of BS1 and shell j of BS2.
overlap!(out, BS1, BS2): out must be BS1.nbas × BS2.nbas.

source

Kinetic Energy

For the atomic orbitals $\chi_\mu, \chi_\nu$, the kinetic energy integral is calculated as:

\[T_{\mu\nu} = -\frac{1}{2}\int \chi_\mu(\mathbf{r})\, \nabla^2 \chi_\nu(\mathbf{r})\, d\mathbf{r}\]

A full-dense kinetic energy matrix can be calculated directly from a basis set

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
              """)
julia> T = kinetic(bset)
7×7 Matrix{Float64}:
 29.0032      -0.168011    0.0          0.0         0.0      -0.00160233  -0.00160233
 -0.168011     0.808128    0.0          0.0         0.0      -0.0167939   -0.0167939
  0.0          0.0         2.52873      0.0         0.0      -0.00612697   0.00612697
  0.0          0.0         0.0          2.52873     0.0      -0.00478691  -0.00478691
  0.0          0.0         0.0          0.0         2.52873   0.0          0.0
 -0.00160233  -0.0167939  -0.00612697  -0.00478691  0.0       0.760032    -0.00447755
 -0.00160233  -0.0167939   0.00612697  -0.00478691  0.0      -0.00447755   0.760032

Alternatively, a in-place (mutating) version can be used:

julia> out = zeros(bset.nbas, bset.nbas)
julia> kinetic!(out, bset)
julia> out == T
true

A per-shell option is also available:

julia> T43 = kinetic(bset, 4, 3) # 4 and 3 are shell indexes
1×3 Matrix{Float64}:
 -0.00612697  -0.00478691  0.0

For maximum efficiency, you may need to use the most primitive call, mutating and using a pre-allocated output

julia> out = zeros(num_basis(bset[4]), num_basis(bset[3]))
julia> kinetic!(out, bset, 4, 3) 
julia> out == T43
true

The kinetic energy can also be computed using two basis set. This will calculate $\langle \chi_\mu | -\frac{1}{2}\nabla^2 | \chi_\nu \rangle$ for $\mu \in A$ and $\nu \in B$, where $A$ and $B$ are two different basis set. This can be useful when projecting results from one basis onto another.

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

julia> kinetic(bset, bset2)
7×13 Matrix{Float64}:
 29.5083      -2.60308      0.0         …  -4.54054e-5  -0.00358451
 -0.234446     1.16236      0.0            -0.0150019   -0.0147396
  0.0          0.0          3.37591         0.0166165   -0.00997474
  0.0          0.0          0.0            -0.0129822    0.00779312
  0.0          0.0          0.0             0.0          0.0
 -0.00139057  -0.00795767  -0.00295226  …  -0.0013665   -0.00713319
 -0.00139057  -0.00795767   0.00295226      1.03401      0.314867

GaussianBasis.kinetic — Function

kinetic(BS::BasisSet) -> Matrix{Float64}
kinetic(BS::BasisSet, i, j) -> Matrix{Float64}
kinetic(BS1::BasisSet, BS2::BasisSet) -> Matrix{Float64}
kinetic(BS1::BasisSet, BS2::BasisSet, i, j) -> Matrix{Float64}

Compute the AO kinetic energy matrix T.

Methods

kinetic(BS): full, dense, symmetric nbas × nbas matrix for BS.
kinetic(BS, i, j): just the (Ni,Nj) block for shells i and j of BS (shell indices, not AO indices).
kinetic(BS1, BS2): mixed-basis matrix, BS1.nbas × BS2.nbas. See overlap for the general mixed-basis convention.
kinetic(BS1, BS2, i, j): mixed-basis block between shell i of BS1 and shell j of BS2.

For repeated calls, see kinetic!, which writes into a preallocated array instead of allocating.

source

GaussianBasis.kinetic! — Function

kinetic!(out, BS::BasisSet, i, j)
kinetic!(out, BS::BasisSet)
kinetic!(out, BS1::BasisSet, BS2::BasisSet, i, j)
kinetic!(out, BS1::BasisSet, BS2::BasisSet)

Mutating counterpart of kinetic: writes into the caller-supplied out instead of allocating. This shell-pair, single-basis form is the primitive every other method builds on.

Methods

kinetic!(out, BS, i, j): out must be (Ni,Nj), the block for shells i/j of BS (shell indices, not AO indices).
kinetic!(out, BS): out must be a dense nbas × nbas matrix.
kinetic!(out, BS1, BS2, i, j): out must be (Ni,Nj), the block for shell i of BS1 and shell j of BS2.
kinetic!(out, BS1, BS2): out must be BS1.nbas × BS2.nbas.

source

Nuclear-electron Attraction

For the atomic orbitals $\chi_\mu, \chi_\nu$, the nuclear attraction integral, summed over every nucleus $C$ in the molecule with charge $Z_C$ at position $\mathbf{R}_C$, is calculated as:

\[V_{\mu\nu} = -\sum_C Z_C \int \chi_\mu(\mathbf{r}) \frac{1}{|\mathbf{r}-\mathbf{R}_C|} \chi_\nu(\mathbf{r})\, d\mathbf{r}\]

A full-dense nuclear attraction matrix can be calculated directly from a basis set

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
              """)
julia> V = nuclear(bset)
7×7 Matrix{Float64}:
 -61.1276      -7.3036      0.0       -0.00405309   0.0      -0.129311   -0.129311
  -7.3036      -9.56106     0.0       -0.0511636    0.0      -0.390434   -0.390434
   0.0          0.0        -9.49705    0.0          0.0      -0.289694    0.289694
  -0.00405309  -0.0511636   0.0       -9.48994      0.0      -0.229019   -0.229019
   0.0          0.0         0.0        0.0         -9.47879   0.0         0.0
  -0.129311    -0.390434   -0.289694  -0.229019     0.0      -3.42421    -0.0375252
  -0.129311    -0.390434    0.289694  -0.229019     0.0      -0.0375252  -3.42421

Alternatively, a in-place (mutating) version can be used:

julia> out = zeros(bset.nbas, bset.nbas)
julia> nuclear!(out, bset)
julia> out == V
true

A per-shell option is also available:

julia> V43 = nuclear(bset, 4, 3) # 4 and 3 are shell indexes
1×3 Matrix{Float64}:
 -0.289694  -0.229019  0.0

For maximum efficiency, you may need to use the most primitive call, mutating and using a pre-allocated output

julia> out = zeros(num_basis(bset[4]), num_basis(bset[3]))
julia> nuclear!(out, bset, 4, 3) 
julia> out == V43
true

The nuclear attraction can also be computed using two basis set (the nuclei are always taken from the first BasisSet). This will calculate $\langle \chi_\mu | \sum_C -Z_C/|\mathbf{r}-\mathbf{R}_C| | \chi_\nu \rangle$ for $\mu \in A$ and $\nu \in B$, where $A$ and $B$ are two different basis set. This can be useful when projecting results from one basis onto another.

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

julia> nuclear(bset, bset2)
7×13 Matrix{Float64}:
 -61.8862      -5.33044      0.0       …  -0.247703   -6.51069e-5  -0.247703
  -6.94425     -9.79189      0.0          -0.706771   -0.0420387   -0.706771
   0.0          0.0        -10.433        -0.499966    0.0630876    0.499966
  -0.00309527  -0.0329122    0.0          -0.395625   -0.0493949   -0.395625
   0.0          0.0          0.0           0.0         0.0          0.0
  -0.122874    -0.21832     -0.114557  …  -2.78587    -0.00443586  -0.0649447
  -0.122874    -0.21832      0.114557     -0.0649447  -3.40551     -2.78587

GaussianBasis.nuclear — Function

nuclear(BS::BasisSet) -> Matrix{Float64}
nuclear(BS::BasisSet, i, j) -> Matrix{Float64}
nuclear(BS1::BasisSet, BS2::BasisSet) -> Matrix{Float64}
nuclear(BS1::BasisSet, BS2::BasisSet, i, j) -> Matrix{Float64}

Compute the AO nuclear attraction matrix V (potential from every nucleus in BS.atoms, summed).

Methods

nuclear(BS): full, dense, symmetric nbas × nbas matrix for BS.
nuclear(BS, i, j): just the (Ni,Nj) block for shells i and j of BS (shell indices, not AO indices).
nuclear(BS1, BS2): mixed-basis matrix (nuclei taken from BS1.atoms), BS1.nbas × BS2.nbas. See overlap for the general mixed-basis convention.
nuclear(BS1, BS2, i, j): mixed-basis block between shell i of BS1 and shell j of BS2, using the nuclei of BS1.

For repeated calls, see nuclear!, which writes into a preallocated array instead of allocating.

source

GaussianBasis.nuclear! — Function

nuclear!(out, BS::BasisSet, i, j)
nuclear!(out, BS::BasisSet)
nuclear!(out, BS1::BasisSet, BS2::BasisSet, i, j)
nuclear!(out, BS1::BasisSet, BS2::BasisSet)

Mutating counterpart of nuclear: writes into the caller-supplied out instead of allocating. This shell-pair, single-basis form is the primitive every other method builds on.

Methods

nuclear!(out, BS, i, j): out must be (Ni,Nj), the block for shells i/j of BS (shell indices, not AO indices).
nuclear!(out, BS): out must be a dense nbas × nbas matrix.
nuclear!(out, BS1, BS2, i, j): out must be (Ni,Nj), the block for shell i of BS1 and shell j of BS2.
nuclear!(out, BS1, BS2): out must be BS1.nbas × BS2.nbas.

source

Multipole Integrals

Theory

Multipole integrals are AO integrals of Cartesian powers of the position operator $\mathbf{r} = (x,y,z)$ about the origin. The dipole, quadrupole, octupole, and hexadecapole integrals are

\[\langle \mu | r_a | \nu \rangle, \quad \langle \mu | r_a r_b | \nu \rangle, \quad \langle \mu | r_a r_b r_c | \nu \rangle, \quad \langle \mu | r_a r_b r_c r_d | \nu \rangle \qquad a,b,c,d \in \{x,y,z\}\]

respectively, each component indexed by which Cartesian direction(s) the $r$ factor(s) refer to. They are the building blocks for molecular electric multipole moments (e.g. contracting the dipole integrals with a density matrix, plus the nuclear point-charge contribution, gives the molecular dipole moment) and for response properties such as polarizabilities.

All multipole integrals are evaluated about the coordinate origin – shift the molecule's geometry first if a different reference point (e.g. the center of mass) is needed.

Usage

GaussianBasis.dipole — Function

dipole(BS::BasisSet) -> Array{Float64,3}
dipole(BS::BasisSet, i, j) -> Array{Float64,3}

Compute the AO electric dipole integral tensor, in atomic units (bohr) about the origin. Combine with a density matrix and nuclear charges to get a molecular dipole moment.

Methods

dipole(BS): full tensor for BS. Returns a dense nbas × nbas × 3 array, the trailing axis indexing the x,y,z Cartesian components.
dipole(BS, i, j): just the block for shells i and j of BS (shell indices, not AO indices). Returns an (Ni,Nj,3) array.

For repeated calls, see dipole!, which writes into a preallocated array instead of allocating.

source

GaussianBasis.dipole! — Function

dipole!(out, BS::BasisSet, i, j)
dipole!(out, BS::BasisSet)

Mutating counterpart of dipole: writes into the caller-supplied out instead of allocating. This shell-pair form is the primitive the full-tensor form builds on.

Methods

dipole!(out, BS, i, j): out must be an (Ni,Nj,3) array, the block for shells i/j of BS (shell indices, not AO indices), the trailing axis indexing the x,y,z Cartesian components.
dipole!(out, BS): out must be a dense nbas × nbas × 3 array.

source

GaussianBasis.quadrupole — Function

quadrupole(BS::BasisSet) -> Array{Float64,4}
quadrupole(BS::BasisSet, i, j) -> Array{Float64,4}

Compute the AO electric quadrupole integral tensor, in atomic units about the origin; see dipole for the general convention.

Methods

quadrupole(BS): full tensor for BS. Returns a dense nbas × nbas × 3 × 3 array.
quadrupole(BS, i, j): just the block for shells i and j of BS (shell indices, not AO indices). Returns an (Ni,Nj,3,3) array.

For repeated calls, see quadrupole!, which writes into a preallocated array instead of allocating.

source

GaussianBasis.quadrupole! — Function

quadrupole!(out, BS::BasisSet, i, j)
quadrupole!(out, BS::BasisSet)

Mutating counterpart of quadrupole: writes into the caller-supplied out instead of allocating. This shell-pair form is the primitive the full-tensor form builds on.

Methods

quadrupole!(out, BS, i, j): out must be an (Ni,Nj,3,3) array, the block for shells i/j of BS (shell indices, not AO indices).
quadrupole!(out, BS): out must be a dense nbas × nbas × 3 × 3 array.

source

GaussianBasis.octupole — Function

octupole(BS::BasisSet) -> Array{Float64,5}
octupole(BS::BasisSet, i, j) -> Array{Float64,5}

Compute the AO electric octupole integral tensor, in atomic units about the origin; see dipole for the general convention.

Methods

octupole(BS): full tensor for BS. Returns a dense nbas × nbas × 3 × 3 × 3 array.
octupole(BS, i, j): just the block for shells i and j of BS (shell indices, not AO indices). Returns an (Ni,Nj,3,3,3) array.

For repeated calls, see octupole!, which writes into a preallocated array instead of allocating.

source

GaussianBasis.octupole! — Function

octupole!(out, BS::BasisSet, i, j)
octupole!(out, BS::BasisSet)

Mutating counterpart of octupole: writes into the caller-supplied out instead of allocating. This shell-pair form is the primitive the full-tensor form builds on.

Methods

octupole!(out, BS, i, j): out must be an (Ni,Nj,3,3,3) array, the block for shells i/j of BS (shell indices, not AO indices).
octupole!(out, BS): out must be a dense nbas × nbas × 3 × 3 × 3 array.

source

GaussianBasis.hexadecapole — Function

hexadecapole(BS::BasisSet) -> Array{Float64,6}
hexadecapole(BS::BasisSet, i, j) -> Array{Float64,6}

Compute the AO electric hexadecapole integral tensor, in atomic units about the origin; see dipole for the general convention.

Methods

hexadecapole(BS): full tensor for BS. Returns a dense nbas × nbas × 3 × 3 × 3 × 3 array.
hexadecapole(BS, i, j): just the block for shells i and j of BS (shell indices, not AO indices). Returns an (Ni,Nj,3,3,3,3) array.

For repeated calls, see hexadecapole!, which writes into a preallocated array instead of allocating.

source

GaussianBasis.hexadecapole! — Function

hexadecapole!(out, BS::BasisSet, i, j)
hexadecapole!(out, BS::BasisSet)

Mutating counterpart of hexadecapole: writes into the caller-supplied out instead of allocating. This shell-pair form is the primitive the full-tensor form builds on.

Methods

hexadecapole!(out, BS, i, j): out must be an (Ni,Nj,3,3,3,3) array, the block for shells i/j of BS (shell indices, not AO indices).
hexadecapole!(out, BS): out must be a dense nbas × nbas × 3 × 3 × 3 × 3 array.

source

Example

julia> using GaussianBasis, Molecules

julia> water = Molecules.parse_string("""
       O        0.000000000000     -0.143225816552      0.000000000000
       H        1.638036840407      1.136548822547     -0.000000000000
       H       -1.638036840407      1.136548822547     -0.000000000000
       """);

julia> bset = BasisSet("sto-3g", water);

julia> D = dipole(bset);

julia> size(D)
(7, 7, 3)

julia> D[1, 1, :]   # ⟨1s_O|r|1s_O⟩ in atomic units (x, y, z)
3-element Vector{Float64}:
  0.0
 -0.2706575672723027
  0.0

octupole and hexadecapole follow the same convention, with 3 and 4 trailing Cartesian axes respectively.