CheMPS2: A free open-source spin-adapted implementation of the density matrix renormalization group for ab initio quantum chemistry

Abstract

The density matrix renormalization group (DMRG) has become an indispensable numerical tool to find exact eigenstates of finite-size quantum systems with strong correlation. In the fields of condensed matter, nuclear structure and molecular electronic structure, it has significantly extended the system sizes that can be handled compared to full configuration interaction, without losing numerical accuracy. For quantum chemistry (QC), the most efficient implementations of DMRG require the incorporation of particle number, spin and point group symmetries in the underlying matrix product state (MPS) ansatz, as well as the use of so-called complementary operators. The symmetries introduce a sparse block structure in the MPS ansatz and in the intermediary contracted tensors. If a symmetry is non-abelian, the Wigner–Eckart theorem allows to factorize a tensor into a Clebsch–Gordan coefficient and a reduced tensor. In addition, the fermion signs have to be carefully tracked. Because of these challenges, implementing DMRG efficiently for QC is not straightforward. Efficient and freely available implementations are therefore highly desired. In this work we present CheMPS2, our free open-source spin-adapted implementation of DMRG for ab initio QC. Around CheMPS2, we have implemented the augmented Hessian Newton–Raphson complete active space self-consistent field method, with exact Hessian. The bond dissociation curves of the 12 lowest states of the carbon dimer were obtained at the DMRG(28 orbitals, 12 electrons, DSU(2)=2500)/cc-pVDZ level of theory. The contribution of 1s core correlation to the X1Σg+ bond dissociation curve of the carbon dimer was estimated by comparing energies at the DMRG(36o, 12e, DSU(2)=2500)/cc-pCVDZ and DMRG-SCF(34o, 8e, DSU(2)=2500)/cc-pCVDZ levels of theory. Program summaryProgram title: CheMPS2Catalogue identifier: AESE_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AESE_v1_0.htmlProgram obtainable from: CPC Program Library, Queen’s University, Belfast, N. IrelandLicensing provisions: GNU General Public License, version 2No. of lines in distributed program, including test data, etc.: 131472No. of bytes in distributed program, including test data, etc.: 1645700Distribution format: tar.gzProgramming language: C++.Computer: x86-64.Operating system: Scientific Linux 6.0.RAM: 10 MB–64 GBClassification: 16.1.External routines: Basic Linear Algebra Subprograms (BLAS), Linear Algebra Package (LAPACK), GNU Scientific Library (GSL), and Hierarchical Data Format Release 5 (HDF5)Nature of problem:The many-body Hilbert space grows exponentially with the number of single-particle states. Exact diagonalization solvers can therefore only handle small systems, of up to 18 electrons in 18 orbitals. Interesting active spaces are often significantly larger.Solution method:The density matrix renormalization group allows the extension of the size of active spaces, for which numerically exact solutions can be found, to about 40 electrons in 40 orbitals. In addition, it provides a rigorous variational upper bound to energies, as it has an underlying wavefunction ansatz, the matrix product state.Restrictions:Our implementation of the density matrix renormalization group is spin-adapted. This means that targeted eigenstates in the active space are exact eigenstates of the total electronic spin operator. Hamiltonians which break this symmetry (a magnetic field term for example) cannot be handled by our code. As electron repulsion integrals in Gaussian basis sets have eightfold permutation symmetry, we have used this property in our code.Unusual features:The nature of the matrix product state ansatz allows for exact spin coupling. In CheMPS2, the total electronic spin is imposed (not just the spin projection), in addition to the particle-number and abelian point-group symmetries.Running time:The running time depends on the size of the targeted active space, the number of desired eigenstates, their symmetry, the density of states, the individual orbital symmetries, the orbital ordering, the desired level of convergence, and the chosen convergence scheme. To converge a single point of one of the dissociation curves of the carbon dimer (D∞h→D2h symmetry) in the cc-pVDZ basis (28 orbitals; their ordering is described in Section 5.3) with 2500 reduced renormalized basis states (see the convergence scheme in Section 5.4; the variational energy then lies 0.1mEh above the fully converged result) takes about 8 h on a single node with a dual-socket octa-core Intel Xeon Sandy Bridge (E5-2670) (16 cores at 2.6 GHz), and requires 6 GB of RAM.