SpinAdaptedSecondQuantization.jl

SpinAdaptedSecondQuantization.jl is a Julia package for doing symbolic second quantization mainly targeted at quantum chemistry methods, such as post Hartree-Fock methods.

To install the package open a Julia terminal, type a ] to open the package mode, then write add SpinAdaptedSecondQuantization to install.

(v1.8) pkg> add SpinAdaptedSecondQuantization

Julia version 1.8 or higher is required.

The package can then be loaded

julia> using SpinAdaptedSecondQuantization

Basic Structure

The core type of the package is the Expression type. This contains an array of Term types. A term is a product that can contain a combination of various elements. These are:

• A scalar multiplier
• An array of operators
• An array of tensors
• An array of Kronecker deltas
• An array of summation indices
• A dictionary of index constraints

Indices

Most types will contain indices. For convenience these are represented as integers, but they are purely symbolic. When constrained to a certain index space they will be printed with names according to the semi-standard names for general molecular orbital (MO) indices pqrstuv, abcdefg for virtual and ijklmno for occupied.

By default indices are unconstrained and will be printed as subscript indices ₁₂...

julia> E(1, 2)E_₁₂

Usually, however, one wants to constrain indices to a certain subspace of indices. By default the available spaces are GeneralOrbital, VirtualOrbital and OccupiedOrbital which can be used for example like:

julia> E(1, 2) * constrain(1 => VirtualOrbital, 2 => OccupiedOrbital)E_ai

The shorthand functions electron(...), virtual(...) and occupied(...) can also be used

julia> E(1, 2) * electron(1, 2)E_pq
julia> E(1, 2) * E(3, 4) * virtual(1, 3) * occupied(2, 4)E_ai E_bj

Operators

A small variety of operator types is supplied by default listed in the following table.

A few examples are:

julia> E(1, 2) * electron(1, 2)E_pq
julia> E(1, 2) * E(3, 4) * virtual(1, 3) * occupied(2, 4)E_ai E_bj
julia> e(1, 2, 3, 4) * electron(1, 2, 3, 4)E_pq E_rs
- δ_qr E_ps
julia> fermion(1, α) * occupied(1)a⁻_iα
julia> fermiondag(2, β) * virtual(2)a†_aβ
julia> bosondag() * boson()b† b⁻
julia> τ(1, 2) * electron(1, 2)T_pq

Tensors

The simplest tensor type is the RealTensor which has a name and an array of indices.

Examples:

julia> real_tensor("h", 1, 2) * electron(1, 2)h_pq
julia> real_tensor("g", 1, 2, 3, 4) * electron(1, 2, 3, 4)g_pqrs

Some other types of tensors are supported with various symmetries in the indices. Supported tensor types are listed below

rsym_tensor is typically used when full real orbital symmetry of integrals is required, such as 2-fold symmetry of hpq and 8-fold symmetry of gpqrs. psym_tensor is useful when doing coupled cluster theory, where symmetries within index pairs do not exist because the integrals are T1 transformed.

Examples of symmetric tensors:

julia> rsym_tensor("g", 4, 3, 1, 2) * electron(1, 2, 3, 4)g_pqrs
julia> psym_tensor("g", 4, 3, 1, 2) * electron(1, 2, 3, 4)g_pqsr

Kronecker Deltas

Kronecker deltas are used to constrain indices to be equal, otherwise the term would be zero. They can have two or more indices each.

Example:

julia> delta(1, 2)δ_₁₂
julia> δ(1, 2) # Unicode version equivalent to the one aboveδ_₁₂
julia> d1 = δ(1, 2)δ_₁₂
julia> d2 = δ(2, 3)δ_₂₃
julia> d1 * d2 # Compacts delta expression since all are equalδ_₁₂₃

Summation Indices

A term can represent a sum over all values of a specific (or many) indices. A sum can be constructed using the summation function, or the unicode aliases ∑ or Σ.

Example:

julia> a = summation(E(1, 2) * electron(1, 2), [1])∑_q(E_qp)
julia> b = ∑(E(1, 2) * electron(1, 2), 1:2) # Unicode∑_pq(E_pq)
julia> a * b # Automatically renames summation indices to not collide∑_qrs(E_qp E_rs)

Hartree-Fock energy expression

As a simple example of use of the package, here is how to derive the HF energy expression.

julia> h = ∑(real_tensor("h", 1, 2) * E(1, 2) * electron(1, 2), 1:2)∑_pq(h_pq E_pq)
julia> g = 1//2 * simplify(
           ∑(psym_tensor("g", 1:4...) * e(1:4...) * electron(1:4...), 1:4)
       )- 1/2 ∑_pqr(g_prrq E_pq)
+ 1/2 ∑_pqrs(g_pqrs E_pq E_rs)
julia> H = h + g + real_tensor("h_nuc")h_nuc
+ ∑_pq(h_pq E_pq)
- 1/2 ∑_pqr(g_prrq E_pq)
+ 1/2 ∑_pqrs(g_pqrs E_pq E_rs)
julia> E_hf = simplify_heavy(hf_expectation_value(H))h_nuc
+ 2 ∑_i(h_ii)
+ 2 ∑_ij(g_iijj)
- ∑_ij(g_ijji)