This project contains a Julia implementation of the DMRG (Density Matrix Renormalization Group) algorithm and some miscellaneous Julia experiments. The experiments were used during the development of the DMRG code, and to explore the structure the molecular hamiltonian matrix.
The DMRG implementation exactly follows the algorithm presented in Schollwöck's paper, The density-matrix renormalization group in the age of matrix product states:
https://arxiv.org/abs/1008.3477
The code is organized into several directories:
-
dmrg: contains all the DMRG code -
experiments: contains various Julia experiments -
itensor: ITensor examples -
lapack: contains an example of how to call a ScaLAPACK function from Julia(qrcp.jl). The example specifically calls the ScaLAPACK functiondlaqps(). -
operator_mpo: contains an MPO implementation that follows Snajberk's paper (https://edoc.ub.uni-muenchen.de/21974/1/Snajberg_Philipp.pdf). This was initially developed to generate MPOs for use with DMRG. However, ITensor's MPO generation mechanism is vastly superior. The tests contain examples of using both MPO implementations, however this code doesn't have much value. It's definitely recommended to generate MPOs using ITensor.
This code has been tested with Julia Version 1.6.0 on macOS:
Julia Version 1.6.0
Commit f9720dc2eb (2021-03-24 12:55 UTC)
Platform Info:
OS: macOS (x86_64-apple-darwin19.6.0)
CPU: Intel(R) Core(TM) i7-4960HQ CPU @ 2.60GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-11.0.1 (ORCJIT, haswell)
Julia packages have notoriously long load times. This code uses some fairly heavy packages and so, not surprisingly, the initial startup time is very long. There are some ways to reduce this initial load time.
For ITensor specifically, one can use its compile() function:
using ITensors
ITensors.compile()
For more details see section 2.7 of https://arxiv.org/pdf/2007.14822.pdf.
A more general solution is to use PackageCompiler. This can speed up loading of any package. The only penalty is that one needs to rerun the compilation to upgrade the package to a new version.
See https://julialang.github.io/PackageCompiler.jl/dev for details.
To install the Julia packages required by the DMRG code and its associated tests, execute the following commands from the Julia REPL:
import Pkg
Pkg.add("TensorOperations")
Pkg.add("KrylovKit")
Pkg.add("ITensors")
The Snajberk-based MPO code requires some additional packages:
import Pkg
Pkg.add("SymPy")
Pkg.add("Parameters")
To directly compute the molecular hamiltonian matrix and
compute its eigenvalue decomposition, install Arpack.
import Pkg
Pkg.add("Arpack")
The following code runs DMRG on a Heisenberg spin chain containing 2 spins (sites):
include("src/dmrg/dmrg.jl")
mpo = MPO(create_heisenberg_mpo(2))
mps = MPS(randn(2, 2))
energy = compute_ground_state!(mpo, mps)
The test test_dmrg.jl contains some good examples of DMRG execution on
a few different types of systems. From the root directory of the project run
julia
include("test/dmrg/test_dmrg.jl")
Execution time for this test is pretty long. It's partially due to package loading. The bulk of the time is spent computing the ground state of a water molecule.
The Julia version of ITensor (https://github.com/ITensor/ITensors.jl) provides a well optimized implementation of the DMRG algorithm.
The file itensor_utils.jl provides some convenience methods for
the execution of DMRG with ITensor. Note also that the examples and
tests in this project use ITensor to generate MPOs. For example,
h,v = read_electron_integral_tensors("data/h2/h2.ezfio.FCIDUMP")
energy,state = molecular_dmrg(h, v)
nuclear_repulsion_energy = 0.71510433908108118
total_energy = energy + nuclear_repulsion_energy
@test isapprox(total_energy, -1.1372838344894254, atol=0.00000001)
As noted above, this project contains code (OperatorMPO) that generates MPOs
using the Snajberk-based construction. However, this code is effectively worthless.
The ITensor AutoMPO method is vastly superior.
The file fcicump.jl provides a function to read one-electron and two-electron
tensor from an FCIDUMP file:
read_electron_integral_tensors(fcidump_filename::String)
This code was specifically written to read FCIDUMP files generated by Quantum Package (https://quantumpackage.github.io/qp2).
Note that the two-electron is permuted at the end of the function:
two_electron_integral_tensor = permutedims(two_electron_integral_tensor, [3,2,1,4])
It's not clear why this permutation is necessary - and it may be possible that it's actually not necessary. However, quantum chemistry calculations using the non-permuted two-electron tensor seem to be incorrect. Additionally, the datasets at Quantaggle (https://github.com/qulacs/Quantaggle_dataset/tree/master/datasets/Small_Molecules_1) don't require this permutation.
One logical explanation is that the difference is due to the ordering of the orbitals. However, if that's the case, then the ground state energy should be the same regardless of the permutation of orbitals. So the problem seems to lie elsewhere.
Note that the hydrogen and water molecule calculations provided in this project use the permuted two-electron tensor. And the results are correct - at least as far as the author can tell.
In any case, there's some weirdness somewhere. So if you use this function for molecules other than hydrogen and water, be sure to consult with a chemist regarding this little mystery.
As noted above Quantaggle provides quantum chemistry data:
https://github.com/qulacs/Quantaggle_dataset/tree/master/datasets/Small_Molecules_1
The data is conveniently stored in HDF5 format. The code below reads the
one-electron and two-electron tensors (and full configuration interaction Energy)
for a water molecule. First install the HDF5 package:
import Pkg
Pkg.add("HDF5")
And then:
using HDF5
io = h5open("data/h2o/H2O_sto-3g_singlet_0.96_104.5deg_0.96.hdf5", "r")
data = read(io)
one_electron_tensor = data["one_body_integrals"]
two_electron_tensor = data["two_body_integrals"]
fci_energy = data["fci_energy"]
As also noted above, there's some confusion - for the author - regarding the correct representation of the two-electron tensor. In the case of Quantaggle, it's not clear if the two-electron tensor is stored in row-major or column-major format. And, given the symmetries of this tensor, it's hard to figure out which data storage format is being used.
In any case, as previously stated, there's some weirdness somewhere. So if you use this data for molecules other than hydrogen and water, be sure to consult with a chemist regarding this other little mystery.
The file molecular_hamiltonian_matrix.jl contains functions to directly
generate the molecular hamiltonian matrix for any (not too large) molecule.
The matrix is computed by directly evaluating the second quantization
expression given arbitrary one-electron and two-electron tensors.
The hamiltonian matrix is stored in sparse format, and it's eigenvalue
decomposition can be readily computed using Arpack.
For example, after installing Arpack:
import Arpack
include("src/experiments/molecules/molecular_hamiltonian_matrix.jl")
h,v = read_electron_integral_tensors("data/h2/h2.ezfio.FCIDUMP")
hamiltonian = molecular_hamiltonian_matrix(h, v)
# add nuclear repulsion part to hamiltonian
m_qp_full = m_qp + (I(16) * 0.71510433908108118)
eigenvalues,eigenvectors = Arpack.eigs(hamiltonian)
@test isapprox(eigenvalues[1], -1.1372838344894252, atol=0.0000001)
Note also that the computation of the hamiltonian matrix is optionally multithreaded. To exploit this, simply launch Julia with several threads. For example,
julia --threads 4
The file data/h2/h2_matrix_and_mpo.jl contains the hamiltonian for a
hydrogen molecule, both in matrix format and MPO format.
This data is generated as follows:
h,v = read_electron_integral_tensors("data/h2/h2.ezfio.FCIDUMP")
h2_matrix = molecular_hamiltonian_matrix(h, v)
sites = ITensors.siteinds("Electron", 2)
one_electron_itensor_mpo = create_one_electron_mpo(sites, h)
two_electron_itensor_mpo = create_two_electron_mpo(sites, v)
one_electron_mpo = MPO(one_electron_itensor_mpo)
two_electron_mpo = MPO(two_electron_itensor_mpo)
h2_mpo = one_electron_mpo + two_electron_mpo