libcommute/pycommute: A quantum operator algebra domain-specific language and exact diagonalization toolkit

I present libcommute, a C++11/14/17 template library that implements a domain-specific language for easy manipulating of polynomial operators used in the quantum many-body theory, as well as a software development toolkit for exact diagonalization codes. The library is written with expressiveness, extensibility and performance in mind. It features simple syntax for commonly used abstractions and algorithms, is well documented and covered by unit tests. libcommute is supplemented with Python 3 bindings called pycommute. They are useful for solving small scale diagonalization problems, rapid prototyping and wrapping of high performance libcommute-based computational cores in Python.


Motivation and significance
The last three decades have seen a boom in the field of computational quantum many-body theory. It has been driven by the limited efficiency of analytical approaches on the one hand, and by the ever growing availability of computer resources on the other. Exact diagonalization (ED) [1,2] has become one of the most important and commonly used families of numerical methods in the field. In a nutshell, an ED calculation amounts to reducing a quantum system of interest to a finite effective model, which retains some essential features of the original system, and whose Hamiltonian is amenable to numerical diagonalization. Computed energy levels and eigenstates are then combined to form expectation values of physical observables and correlation functions. In addition to being a stand-alone solver, ED is an integral part of more sophisticated methods, for instance the Hybridization Expansion Quantum Monte Carlo [3] and master equation solvers for non-equilibrium dynamics [4]. In many cases, ED can provide physical insights without revealing the full eigensystem of the problem. It is often sufficient to perform partial diagonalization within a sector defined by a set of quantum numbers, or to apply an iterative method [5,6,7] to compute only the ground and a few excited states.
In this paper, I present libcommute [8], a C++ library that aims at establishing a reliable and highly extensible framework for developing ED codes. Development of libcommute has been pursuing the following goals.
• Providing a lightweight, header-only library without external dependencies written in the well-established C++11 dialect (some advanced features require the more recent C++17).
• Employing C++11 meta-programming techniques to ensure expressiveness, ease of use, flexibility and a high degree of extensibility.
• Orientation towards high performance computing and solution of very large scale problems, especially in the context of iterative diagonalization methods.
• Robustness and good maintainability achieved via comprehensive unit testing (built on top of the Catch 2 framework [9]).
• Documentation covering all basic concepts of the library with detailed class and function descriptions and a set of usage examples.
It is worth noting that libcommute does not provide any container types to store matrices, or algorithms to find their eigensystems. These tasks are left to the plentiful and well optimized C++ linear algebra libraries, such as Eigen [10] and Blaze [11]. Instead, libcommute is focused on providing a convenient way to specify the relevant objects in the language of quantum operators. libcommute is accompanied by a set of Python 3 bindings called pycommute. These are meant to further simplify development of ED codes in one the following ways.
• pycommute can be directly used to solve moderately sized problems when performance is not a concern.
• pycommute is useful for rapid prototyping of libcommute-based codes.
• pycommute code may serve as a Python input layer for a high performance computational core written in C++ and wrapped in Python by means of pybind11 [12].
libcommute borrows and greatly extends upon some ideas from portions of the TRIQS library [13], which was co-written by the author of the present paper. The ±h.c. syntax feature was borrowed from the excellent TBTK library [14] written and maintained by Kristofer Björnson. Another related work is the C++ port of the SeQuant library [15], which features some facilities similar to those of libcommute's Domain-Specific Language (see Sec.2.2.1 for details). Other open-source ED construction toolkits for C++ include EDLib [16], pomerol [17], parts of the ALPS library [18], Qbasis [19] and the 'lattice-symmetries' package [20]. The most notable Python projects in this category are QuSpin [21,22], PYED [23], pybinding [24] and KWANT [25] (the latter two are specialized for diagonalization of non-interacting tightbinding models).

Software Architecture
libcommute is a C++11/14/17 template library that includes two major parts.
• A Domain-Specific Language (DSL) designed to easily construct and manipulate polynomial expressions built out of quantum-mechanical operators used in the quantum many-body theory (fermion and boson creation and annihilation operators, spin ladder operators). The most commonly used instances of such expressions are many-body Hamiltonians and operators representing physical observables.
• A representation of the polynomial expressions in a form of linear operators, which enables their action on state vectors in finite-dimensional Hilbert spaces. This feature, together with a set of supplementary tools, provides a basis for writing highly performant ED codes without loss of flexibility.
pycommute [26] is a Python package that uses pybind11 [12] to wrap a subset of libcommute's features. It also provides an extra module for easy construction of some model Hamiltonians widely used in the theory of quantum many-body systems, quantum optics and the theory of spin lattices. pycommute is available from the Python Package Index (PyPI) [27] and as part of a public Docker image [28].

Software Functionalities 2.2.1. The Domain-Specific Language
libcommute's DSL revolves around the notion of a polynomial expression (the expression class template).
Expressions are assembled from non-commuting algebra generators g α (such as fermion creation/annihilation operators c † i,j,k,... , c i,j,k,... ) and numerical coefficients C αβ... using a syntax that closely resembles mathematical notation. They support basic arithmetic operations +, −, · with simplification of their results, Hermitian conjugation, and the ±h.c. shorthand (plus/minus Hermitian conjugate). The operator-valued expressions are internally stored in the explicit form (1) rather than in a matrix representation, which means they do not suffer from the exponential dimensionality explosion, and their allowed degree is practically unlimited. Individual terms and generators in (1) are accessible via constant STL-like iterators.
Algebra generators carry statically typed lists of indices i, j, k, . . . which can represent integer lattice site coordinates, quasi-momentum components, spin or orbital string labels, or indices of a user-defined type. With C++17 available, it is also possible to construct expressions from generators with types of indices defined at runtime via the specially provided index type dynamic_indices.
libcommute supports building expressions out of fermionic ladder operators c, c † , bosonic operators a, a † (generators of a canonical (anti)commutation relation algebra), and arbitrary (half-)integer spin operators S z , S ± . Furthermore, generators of new algebras can be defined by deriving from the abstract base generator and implementing a few methods that describe (anti)commutation relations and simplification rules for the new algebra 1 . In general, libcommute's DSL can work with algebras, whose generators g α obey the commutation relations with real constants c, F αβ and f γ αβ . Such algebraic structures include the Lie and Clifford algebras.
Another customization point for the DSL is the type of the numerical coefficients C αβ... . Beside the common double precision real and complex numbers, one can use a custom numeric-like type. This feature can prove useful when working with time-dependent Hamiltonians. Coefficients of the corresponding expressions can be number-like objects storing interpolators or power series w.r.t. the time variable. Specializing the scalar_traits structure for the new coefficient type will teach libcommute how to use it in an expression. Mathematically speaking, any type whose values form a ring with operations + and · can serve as a coefficient type (in fact, a ring without a multiplicative identity would suffice).

Linear operators and Exact Diagonalization tools
The second major part of libcommute, which makes it practically useful for writing exact diagonalization solvers, is the linear operator framework. Instances of the loperator class represent action of polynomial expressions (1) on state vectors in finite dimensional Hilbert spaces. Such spaces are described by hilbert_space objects and constructed as ordered direct products of elementary spaces H i , Each elementary space H i corresponds to a single quantum degree of freedom (DOF) and can be one of the following.
• An elementary space corresponding to a DOF from a user-defined algebra.
The respective C++ types are classes implementing the elementary_space interface.
hilbert_space's API allows to explicitly build H out of elementary_space_* objects corresponding to the factors H i . Another option is to delegate this task to the function make_hilbert_space(), which automatically analyzes a polynomial expression and builds a minimal space required to accommodate the corresponding operator.
Instances of loperator are callable C++ objects that can be applied to a state vector object -a container or view type modeling a special StateVector concept. libcommute directly supports std::vector and some vector-like types of the Eigen 3 library [10] as state vector objects, including the Eigen::Map views. Types provided by other matrix-vector algebra libraries can readily be made loperator-compatible by overloading a few free functions for them. Requirements of the StateVector concept are rather lax and can be met even by such intricate data types as GPU memory-based [29,30,31] and distributed arrays [32,33,34,35].
Following the same design principle as expression, loperator does not internally store any matrices (neither dense nor sparse) and describes action on the state vectors only as a set of element selection and transformation rules. This solution is memory efficient and fits particularly well with widely used iterative methods of partial diagonalization, such as power iteration [5] and the Lanczos algorithm [6,36,7].
In addition to the basic building blocks of ED solvers -Hilbert spaces, state vectors and linear operatorslibcommute ships a handful of supplemental tools. sparse_state_vector is a lightweight sparse array type that can be used to implement iterative ED solvers with repeated elimination of negligible quantum amplitudes [37]. mapped_basis_view is a view of an underlying state vector container. It translates basis states of a Hilbert space according to a predefined map and is useful when only a specific sector of a model (block of its Hamiltonian) is to be diagonalized. A more specialized and less memory consuming version of such views is called n_fermion_(multi)sector_view and works for the N-particle sectors of models with multiple fermionic DOF. Finally, the space_partition utility class reveals invariant subspaces of Hamiltonians using the automatic space partition algorithm described in Section 4 of [38]. The algorithm reorders basis states in such a way that the Hamiltonian is brought to a block-diagonal form with the smallest possible block sizes. It generates input basis state maps for mapped_basis_view and, therefore, makes it possible to diagonalize individual sectors without knowing integrals of motion of the studied model.

pycommute
The Python package pycommute features three modules.
1. expression wraps the polynomial expression C++ types with real (class ExpressionR) and complex coefficients (ExpressionC), their arithmetics and functions to create generators of fermionic, bosonic and spin algebras. Support for user-defined algebras is planned for a future release. 2. loperator contains wrappers of Hilbert space types and linear operators acting on one-dimensional NumPy [39] arrays. It also exposes some of the state vector view types, the automatic space partitioning utilities, and a few overloads of the function make_matrix(), which returns a matrix representation of a linear operator in a form of a twodimensional NumPy array. 3. models is a set of factory functions for expressions of some widely used model Hamiltonians. These include tight-binding models, fermion and boson interaction terms, Ising, Heisenberg and Dzyaloshinskii-Moriya spin coupling terms, and a few spin-boson type models.
It is possible to pass expressions, linear operators and other pycommute objects between Python scripts and C++ functions defined in a pybind11module. This makes coding up interface layers of hybrid C++/Python ED codes a hassle-free task.

Illustrative Examples
This section briefly highlights a few libcommute/pycommute usage examples. Complete code listings of these and further examples with in-depth explanations are to be found on the respective documentation web-sites.

libcommute examples
The first example [40] concerns an integrable quantum system, a spin-1/2 Heisenberg chain. libcommute's DSL is used to verify some explicit expressions for higher charges of the model derived in [41].
The second example [42] demonstrates how to partially diagonalize a onedimensional Hubbard-Holstein model [43,44] using libcommute and Eigen 3 [10]. The model describes behavior of strongly correlated electrons on a lattice coupled to localized phonons. Diagonalization is performed within a sector with a fixed number of electrons.

pycommute examples
Both pycommute examples mentioned below are simple Python scripts that make use of functions from the models module to readily construct Hamiltonians, and of NumPy [39] to diagonalize them.
• [45]: Sector-wise diagonalization of an interaction term of electrons in an atomic d-shell parametrized by Slater integrals F k [46,47].
• [48]: Spectrum of a generalized Jaynes-Cummings (Tavis-Cummings) model [49] describing dynamics of two qubits coupled via an electromagnetic mode in a cavity.

Impact
The main anticipated impact of the presented software is a drastic reduction of development time for Exact Diagonalization codes written in C++ and (optionally) using Python scripting as their input/output layer. libcommute's domain specific language provides a very convenient and flexible way to deal with quantum-mechanical operators, while the loperator module abstracts out a lot of small and hard-to-keep-track-of details. Availability of an open source, well tested, well documented and highly reusable framework -such as libcommute -can make building ED codes a much less tedious and error-prone task. Since libcommute has been and will continue to be developed with speed and memory consumption in mind, its users can enjoy a decent performance level without investing much time into profiling and fine tuning.
The development time can be further reduced by first trying out new ideas on the level of pycommute-based Python scripts. As libcommute and pycommute share the same set of basic concepts, translating the scripts into HPC-ready C++ code is fairly easy.
Last but not least, pycommute allows one to define Hamiltonians and operators of physical observables and pass them to various libcommute/pycommutebased solvers in a uniform manner. It can, therefore, establish a universal input data format for such computational programs.
libcommute/pycommute is still a very young project. At the moment, it provides a foundation for the recent 2.0 release of the Pomerol ED library [17], as well as for a few privately developed solvers. Nonetheless, wider adoption of the framework by the community would certainly simplify development of modern and complex ED solvers, while improving reliability and reproducibility of their results.

Conclusions
To conclude, I have presented libcommute/pycommute, a combination of a C++11 template library and its Python bindings, whose primary goal is rapid development of complex Exact Diagonalization solvers for models of quantum many-body theory. Despite being a young project, it already has the potential to provide a solid foundation for future high performance scientific software, which is both easy to write and to use. Future development of libcommute will focus on expanding the arsenal of the ED construction tools. Newer versions of pycommute are going to expose more of libcommute's functionality, in particular, make it possible to define new operator algebras in the Python code.

Conflict of Interest
No conflict of interest exists: I wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.