This package is a work in progress! It mainly exists as an internal experiment by the ITensor developers to test out the capabilities of ITensorMPS.jl for solving quantum chemistry problems and is not meant for external usage, except as a reference for those interested in trying out similar calculations. The design is very naïve, basically we call out to the Hartree-Fock functionality in PySCF using PythonCall.jl and construct a (symbolic) second quantized Hamiltonian in the resulting molecular orbital basis. We then convert that second quantized Hamiltonian to an MPO and run DMRG. One major outstanding issue we came across is that in standard Gaussian basis sets, the number of terms in the quantum chemistry Hamiltonian grows very quickly (generically, it scales as |
The main functionality of this package is outputting a second quantized quantum chemistry Hamiltonian in the molecular orbital basis, given a molecule and atomic orbital basis.
Under the hood, the package uses Hartree-Fock implemented in PySCF, which we call using PythonCall.jl, to obtain the molecular orbital basis and one-electron and two-electron integrals.
The main output is an OpSum
from ITensors.jl, which is a representation of the second quantized Hamiltonian. This can be converted into a variety of other formats, such as a matrix product operator (MPO) to run DMRG, quantum circuit, full matrix representation for exact diagonalization (ED) for full configuration interaction (FCI) calculations, etc.
julia> using Pkg
julia> Pkg.add(; url="https://github.com/mtfishman/ITensorChemistry.jl")
using ITensors, ITensorMPS
using ITensorChemistry
using Plots
function energy_at_bond(r)
# define molecule geometry
molecule = Molecule([("H", 0.0, 0.0, 0.0),
("H", r, 0.0, 0.0)])
# build electronic hamiltonian and solve HF
hf = molecular_orbital_hamiltonian(molecule; basis="sto-3g")
hamiltonian = hf.hamiltonian
hartree_fock_state = hf.hartree_fock_state
hartree_fock_energy = hf.hartree_fock_energy
# hilbert space
s = siteinds("Electron", 2; conserve_qns=true)
H = MPO(hamiltonian, s)
# initialize MPS to HF state
ψhf = MPS(s, hartree_fock_state)
# run dmrg
dmrg_kwargs = (;
nsweeps=10,
maxdim=[10,20,30,40,50,100],
cutoff=1e-8,
noise=[1e-6, 1e-7, 1e-8, 0.0],
)
dmrg_energy, _ = dmrg(H, ψhf; nsweeps=10, outputlevel=0)
return hartree_fock_energy, dmrg_energy
end
# bond distances
r⃗ = 0.3:0.03:3.0
energies = []
for r in r⃗
push!(energies, energy_at_bond(r))
end
using ITensors, ITensorMPS
using ITensorChemistry
# Nitrogen molecule
molecule = Molecule("N₂")
basis = "sto-3g"
@show molecule
hf = molecular_orbital_hamiltonian(molecule; basis);
hamiltonian = hf.hamiltonian;
hartree_fock_state = hf.hartree_fock_state
println("Number of orbitals = ", length(hartree_fock_state))
@show hamiltonian[end];
println("Number of fermionic operators = ", length(hamiltonian))
println("Hartree-Fock state |HF⟩ = |", prod(string.(hartree_fock_state)),"⟩")
qubit_hamiltonian = jordanwigner(hamiltonian);
qubit_state = jordanwigner(hartree_fock_state)
@show qubit_hamiltonian[end];
println("Number of qubit operators = ", length(qubit_hamiltonian))
println("Hartree-Fock state |HF⟩ = |", prod(string.(qubit_state)),"⟩")
# --------------------------------------------------------------------------
# molecule = Molecule
# Atom 1: N, r⃗ = (0.0, 0.0, 0.550296)
# Atom 2: N, r⃗ = (0.0, 0.0, -0.550296)
# Number of orbitals = 10
# Number of fermionic operators = 14181
# |HF⟩ = |4444444111⟩
# Number of qubit operators = 17005
# |HF⟩ = |22222222222222111111⟩