wrHopHopupHopdnQuick Start Code on GitHub
"},{"location":"#overview","title":"Overview","text":"XDiag is a library for performing Exact Diagonalizations of quantum many-body systems. Key features include optimized combinatorical algorithms for navigating Hilbert spaces, iterative linear algebra algorithms, shared and distributed memory parallelization. It consist of two packages:
Please support our work by citing XDiag and the implemented algorithms if it is used in your published research.
@article{Wietek2018,\n title = {Sublattice coding algorithm and distributed memory parallelization for large-scale exact diagonalizations of quantum many-body systems},\n author = {Wietek, Alexander and L\\\"auchli, Andreas M.},\n journal = {Phys. Rev. E},\n volume = {98},\n issue = {3},\n pages = {033309},\n numpages = {10},\n year = {2018},\n month = {Sep},\n publisher = {American Physical Society},\n doi = {10.1103/PhysRevE.98.033309},\n url = {https://link.aps.org/doi/10.1103/PhysRevE.98.033309}\n}\nSupporting material for lecture held at Quant24 master's school at MPI PKS. Consists of a Jupyter notebook and a sample lattice file describing the \\(N=12\\) site triangular lattice Heisenberg model:
This notebook uses the Julia verision of XDiag and covers the basic functionality:
Hello World
Prints out a greeting containing information on the version of the code.
source
Groundstate energy
Computes the ground state energy of a simple Heisenberg spin \\(S=1/2\\) chain
source
\\(t\\)-\\(J\\) time evolution
Computes the time evolution of a state in the \\(t\\)-\\(J\\) model with distributed parallelization
source
Normal XDiag
Template CMakeLists.txt which can be used to compile applications with the normal XDiag library.
source
Distributed XDiag
Template CMakeLists.txt which can be used to compile applications with the distributed XDiag library.
source
Enter the package mode using ] in the Julia REPL and type:
add XDiag\nThat's it!
"},{"location":"installation/#c-compilation","title":"C++ Compilation","text":"Using XDiag with C++ is a two-step process. First the xdiag library needs to be compiled and installed. Therafter, application codes are compiled in a second step. Here we explain how to compile the library.
g++, clang, or Intel's icpx)Download the source code using git
cd /path/to/where/xdiag/should/be\ngit clone https://github.com/awietek/xdiag.git\nCompile the default library
cd xdiag\ncmake -S . -B build\ncmake --build build\ncmake --install build\ninstall. Compile the distributed library
To use the distributed computing features of xdiag, the distributed library has to be built which requires MPI.
cd xdiag\ncmake -S . -B build -D XDIAG_DISTRIBUTED=On\ncmake --build build\ncmake --install build\nInfo
It might be necessary to explicitly define MPI compiler, e.g. mpicxx like this
cmake -S . -B build -D XDIAG_DISTRIBUTED=On -D CMAKE_CXX_COMPILER=mpicxx\nParallel compilation To speed up the compilation process, the build step can be performed in parallel using the -j flag
cmake --build build -j\nListing compile options
The available compilation options can be displayed using
cmake -L .\nChoosing a certain compiler
The compiler (e.g. icpx) can be specified using
cmake -S . -B build -D CMAKE_CXX_COMPILER=icpx\nWarning
If the xdiag library is compiled with a certain compiler, it is advisable to also compile the application codes with the same compiler.
Setting the install path
In the installation step, the install directory can be set in the following way
cmake --install build --prefix /my/install/prefix\nDisabling HDF5/OpenMP
To disable support for HDF5 or OpenMP support, use
cmake -S . -B build -D XDIAG_DISABLE_OPENMP=On -D XDIAG_DISABLE_HDF5=On\nBuilding and running tests
To compile and run the testing programs, use
cmake -S . -B build -D BUILD_TESTING=On\ncmake --build build\nbuild/tests/tests\nBuilding the Julia wrapper locally
First, get the path to the CxxWrap package of julia. To do so, enter the Julia REPL,
julia\nusing CxxWrap\nCxxWrap.prefix_path()\n/path/to/libcxxwrap-julia-prefix. This is then used to configure the cmake compilation. cmake -S . -B build -D XDIAG_JULIA_WRAPPER=On -D CMAKE_PREFIX_PATH=/path/to/libcxxwrap-julia-prefix\ncmake --build build\ncmake --install build\nlibxdiagjl.so, (or the corresponding library format on non-Linux systems). The source files for the documentation can be found in the directory docs. The documentation is built using Material for MKDocs. To work on it locally, it can be served using
mkdocs serve\nfrom the xdiag root source directory. A local build of the documentation can then be accessed in a webbrowser at the adress
127.0.0.1:8000\nLet us set up our first program using the xdiag library.
using XDiag\nsay_hello()\n#include <xdiag/all.hpp>\n\nusing namespace xdiag;\n\nint main() try {\n say_hello();\n} catch (Error e) {\n error_trace(e);\n}\nThe function say_hello() prints out a welcome message, which also contains information which exact XDiag version is used. In Julia this is all there is to it.
For the C++ code we need to create two files to compile the program. The first is the actual C++ code. What is maybe a bit unfamiliar is the try / catch block. XDiag implements a traceback mechanism for runtime errors, which is activated by this idiom. While not stricly necessary here, it is a good practice to make use of this.
Now that the application program is written, we next need to set up the compilation instructions using CMake. To do so we create a second file called CMakeLists.txt in the same directory.
cmake_minimum_required(VERSION 3.19)\n\nproject(\n hello_world\n)\n\nfind_package(xdiag REQUIRED HINTS \"/path/to/xdiag/install\")\nadd_executable(main main.cpp)\ntarget_link_libraries(main PRIVATE xdiag::xdiag)\nYou should replace \"/path/to/xdiag/install\" with the appropriate directory where your XDiag library is installed after compilation. This exact CMakeLists.txt file can be used to compile any XDiag application.
Info
For using the distributed XDiag library the last line of the above CMakeLists.txt should be changed to
target_link_libraries(main PUBLIC xdiag::xdiag_distributed)\nWe then compile the application code,
cmake -S . -B build\ncmake --build build\nand finally run our first xdiag application.
./build/main\nWe compute the ground state energy of the \\(S=1/2\\) Heisenberg chain on a periodic chain lattice in one dimension. The Hamiltonian is given by
\\[ H = J\\sum_{\\langle i,j \\rangle} \\mathbf{S}_i \\cdot \\mathbf{S}_j\\]where \\(\\mathbf{S}_i = (S_i^x, S_i^y, S_i^z)\\) are the spin \\(S=1/2\\) operators and \\(\\langle i,j \\rangle\\) denotes summation over nearest-meighbor sites \\(i\\) and \\(j\\).
The following code, sets up the Hilbert space, defines the Hamiltonian and finally calls an iterative eigenvalue solver to compute the ground state energy.
JuliaC++using XDiag\n\nlet \n N = 16\n nup = N \u00f7 2\n block = Spinhalf(N, nup)\n\n # Define the nearest-neighbor Heisenberg model\n ops = OpSum()\n for i in 1:N\n ops += Op(\"HB\", \"J\", [i-1, i % N])\n end\n ops[\"J\"] = 1.0\n\n set_verbosity(2) # set verbosity for monitoring progress\n e0 = eigval0(ops, block) # compute ground state energy\n\n println(\"Ground state energy: $e0\")\nend\n#include <xdiag/all.hpp>\n\nusing namespace xdiag;\n\nint main() try {\n int N = 16;\n int nup = N / 2;\n Spinhalf block(N, nup);\n\n // Define the nearest-neighbor Heisenberg model\n OpSum ops;\n for (int i = 0; i < N; ++i) {\n ops += Op(\"HB\", \"J\", {i, (i + 1) % N});\n }\n ops[\"J\"] = 1.0;\n\n set_verbosity(2); // set verbosity for monitoring progress\n double e0 = eigval0(ops, block); // compute ground state energy\n\n Log(\"Ground state energy: {:.12f}\", e0);\n\n} catch (Error e) {\n error_trace(e);\n}\nSep. 9, 2024
Introduced 1-indexing everywhere in Julia version
Aug. 27, 2024
Lanczos routines and multicolumn States
Aug. 16, 2024
Small patch release providing small utility functions
Aug. 15, 2024
Basic functionality for three Hilbert space types, Spinhalf, tJ, and Electron, has been implemented. Features are:
A step-by-step guide to using XDiag
Installation Language Julia installation How to install the XDiag in Julia C++ compilation How to compile the XDiag C++ library First steps Language Writing code How to write an application using XDiag Compilation How to compile an application using XDiag in C++ Hilbert spaces How to define a Hilbert space Operators How to define an operator Ground state How compute a ground state (energy) Measurements How to measure an observable on a state Input/Output Language TOML How to read input from a TOML file"},{"location":"documentation/","title":"Documentation","text":"XDiag uses the C++ library Armadillo as a linear algebra backend. Documentation for linear algebra operations can, therefore, be found in the Armadillo Documentation.
"},{"location":"documentation/#algebra","title":"Algebra","text":"Name Description Language matrix Creates the full matrix representation of an operator on a block apply Applies an operator to a state \\(\\vert \\phi \\rangle = O \\vert \\psi\\rangle\\) dot Computes the dot product between two states inner Computes an expectation value \\(\\langle \\psi \\vert O \\vert \\psi \\rangle\\) norm Computes the 2-norm of a state norm1 Computes the 1-norm of a state norminf Computes the \\(\\infty\\)-norm of a state"},{"location":"documentation/#algorithms","title":"Algorithms","text":"Diagonalization
Name Description Language eigval0 Computes the lowest lying eigenvalue of an operator eig0 Computes the lowest lying eigenvalue and eigenvector of an operator eigvals_lanczos Performs an iterative eigenvalue calculation using the Lanczos algorithm eigs_lanczos Performs an iterative eigenvalue calculation building eigenvectors using the Lanczos algorithmTime evolution
Name Description Language time_evolve Performs a real-time evolution \\(e^{ -iHt} \\vert \\psi \\rangle\\) of a State with a given Hermitian operator \\(H\\) imaginary_time_evolve Performs a imaginary-time evolution \\(e^{ -\\tau H}\\vert\\psi\\rangle\\) of a State with a given Hermitian operator \\(H\\) evolve_lanczos Computes the exponential $e^{z H}\\vert\\psi\\rangle $ of a Hermitian operator times a State for a real or complex \\(z\\) using the Lanczos algorithm time_evolve_expokit Performs a real-time evolution \\(e^{ -iHt} \\vert \\psi \\rangle\\) using a highly accurate Lanczos algorithm"},{"location":"documentation/#blocks","title":"Blocks","text":"Shared memory
Name Description Language Spinhalf Block of a spin \\(S=1/2\\) type Hilbert space tJ Block of a \\(t-J\\) type Hilbert space Electron Block of a Electron type Hilbert spaceDistributed memory
Name Description Language SpinhalfDistributed Block of a spin \\(S=1/2\\) type Hilbert space (distributed computing) tJDistributed Block of a \\(t-J\\) type Hilbert space (distributed computing)"},{"location":"documentation/#building","title":"Building","text":"Name Description Language Compilation Advanced settings for compilation of the C++ library using CMake Documentation How to build and work on the documentation locally Julia Wrapper How to build and develop the Julia wrapper locally"},{"location":"documentation/#input-output","title":"Input / Output","text":"TOML
Name Description Language FileToml A file handler for TOML files read_opsum reads an OpSum from a TOML file read_permutation_group reads a PermutationGroup from a TOML file read_representation reads a Representation from a TOML fileHDF5
Name Description Language FileH5 A file handler for hdf5 files"},{"location":"documentation/#operators","title":"Operators","text":"Operator types A summary of all the operator types defined in XDiag Name Description Language Op A linear operator acting on the Hilbert space OpSum Sum of couplings times operators hc Returns the hermitian conjugate of an Op or OpSum. symmetrize Symmetrizes an operator with a PermutationGroup or Representation Scalar A scalar number which can be either real or complex Coupling Describes a coupling associated with an operator, either string or scalar"},{"location":"documentation/#states","title":"States","text":"Name Description Language State A generic state describing a quantum wave function ProductState A product state of local configurations RandomState A random state with normal distributed coefficients fill Fill a state with a given model state product_state Creates a filled product state random_state Create a filled random state with normal distributed coefficients zero_state Create a filled state with all zero entries zero Set all coefficients of a given state to zero"},{"location":"documentation/#symmetries","title":"Symmetries","text":"Name Description Language Permutation Permutations of indices or lattice sites PermutationGroup A group of permutations Representation A (1D) irreducible representation of a finite group"},{"location":"documentation/#utilities","title":"Utilities","text":"Name Description Language set_verbosity Sets how much information is printed during computations say_hello Prints a nice welcome message with version number print_version Prints the plain version number Logging Controling what is written to standard output Timing Measurng wall time straightforwardly XDIAG_SHOW Macro for printing debugging information"},{"location":"documentation/algebra/algebra/","title":"Algebra routines","text":"Several basic algebraic operations for states and operators.
Sources algebra.hpp, algebra.cpp
"},{"location":"documentation/algebra/algebra/#dot","title":"dot","text":"Computes the dot product \\(\\langle v \\vert w \\rangle\\) between two states \\(\\vert v \\rangle\\) and \\(\\vert w \\rangle\\). In C++, please use the dotC function if one of the two states is expected to be complex.
C++Juliadouble dot(State const &v, State const &w);\ncomplex dotC(State const &v, State const &w);\ndot(v::State, w::State)\nComputes the expectation value \\(\\langle v | O |v \\rangle\\) of an operator \\(O\\) and a state \\(|v\\rangle\\). The operator can either be an Op or an OpSum object. In C++, please use the innerC function if either the operator or the state are complex.
C++Juliadouble inner(Op const &op, State const &v);\ndouble inner(OpSum const &ops, State const &v);\ncomplex innerC(Op const &op, State const &v);\ncomplex innerC(OpSum const &ops, State const &v);\ninner(op::Op, v::State)\ninner(ops::OpSum, v::State)\nComputes the 2-norm \\(\\parallel |v \\rangle \\parallel_2\\) of a state \\(|v \\rangle\\) defined as
\\[ \\parallel |v \\rangle \\parallel_2 = \\sum_n |\\langle n | v \\rangle |^2, \\]where \\(\\{ |n\\rangle \\}\\) denotes an orthonormal basis of the block.
C++Juliadouble norm(State const &v);\nnorm(state::State)\nComputes the 1-norm \\(\\parallel |v \\rangle \\parallel_1\\) of a state \\(|v \\rangle\\) defined as
\\[ \\parallel |v \\rangle \\parallel_1 = \\sum_n |\\langle n | v \\rangle |, \\]where \\(\\{ |n\\rangle \\}\\) denotes an orthonormal basis of the block.
C++Juliadouble norm1(State const &v);\nnorm1(state::State)\nComputes the \\(\\infty\\)-norm \\(\\parallel |v \\rangle \\parallel_\\infty\\) of a state \\(|v \\rangle\\) defined as
\\[ \\parallel |v \\rangle \\parallel_\\infty = \\max_n |\\langle n | v \\rangle |, \\]where \\(\\{ |n\\rangle \\}\\) denotes an orthonormal basis of the block.
C++Juliadouble norminf(State const &v);\nnorminf(state::State)\nint N = 8;\nauto block = Spinhalf(N, N / 2);\nauto ops = OpSum();\nfor (int i=0; i<N; ++i) {\n ops += Op(\"SdotS\", {i, (i+1)%N});\n}\nauto [e0, psi] = eig0(ops, block);\n\nXDIAG_SHOW(norm(psi));\nXDIAG_SHOW(norm1(psi));\nXDIAG_SHOW(norminf(psi));\n\nXDIAG_SHOW(dot(psi, psi));\nXDIAG_SHOW(e0);\nXDIAG_SHOW(inner(ops, psi));\n\nauto phi = random_state(block);\nXDIAG_SHOW(phi.vector());\nXDIAG_SHOW(psi.vector());\nXDIAG_SHOW((psi + 2.0*phi).vector());\nXDIAG_SHOW((psi*complex(0,3.0) + phi/2.0).vectorC());\nlet \n N = 8\n block = Spinhalf(N, N \u00f7 2)\n ops = OpSum()\n for i in 1:N\n ops += Op(\"HB\", 1.0, [i, mod1(i+1, N)])\n end\n e0, psi = eig0(ops, block);\n\n @show norm(psi)\n @show norm1(psi)\n @show norminf(psi)\n\n @show dot(psi, psi)\n @show e0, inner(ops, psi)\n\n phi = rand(block)\n display(vector(phi))\n display(vector(psi))\n display(vector(psi + 2.0*phi))\n display(vector(psi*3.0im + phi/2.0))\nend\nApplies an operator given as an Op or OpSum to a State \\(\\vert w \\rangle = \\mathcal{O} \\vert v\\rangle\\).
Sources apply.hpp, apply.cpp
"},{"location":"documentation/algebra/apply/#definition","title":"Definition","text":"An operator \\(\\mathcal{O}\\) can be applied to a state \\(\\vert v\\rangle\\) in two ways:
Only the input state on which the operator acts is defined is given. The block of the output state is calculated and eventually created automatically.
C++JuliaState apply(Op const &op, State const &v);\nState apply(OpSum const &ops, State const &v);\napply(op::Op, v::State)\napply(ops::OpSum, v::State)\nThe output state is also handed as an argument which is overwritten. The compatibility of quantum numbers is checked. This way the output block is not created automatically and, thus, can be used to save computation time if the output block appears repeatedly in the computation.
C++Juliavoid apply(Op const &op, State const &v, State &w);\nvoid apply(OpSum const &ops, State const &v, State &w);\napply(op::Op, v::State, w::State)\napply(ops::OpSum, v::State, w::State)\nint N = 8;\nauto block = Spinhalf(N, N / 2);\nauto ops = OpSum();\nfor (int i=0; i<N; ++i){\n ops += Op(\"SdotS\", {i, (i+1)%N});\n}\nauto [e0, psi] = eig0(ops, block);\nauto phi = apply(Op(\"S+\", 2), psi);\nXDIAG_SHOW(inner(ops, psi));\nXDIAG_SHOW(inner(ops, phi));\nlet \n N = 8\n block = Spinhalf(N, N \u00f7 2)\n ops = OpSum()\n for i in 1:N\n ops += Op(\"HB\", 1.0, [i, mod1(i+1, N)])\n end\n e0, psi = eig0(ops, block);\n\n blockp = Spinhalf(N, N \u00f7 2 + 1)\n phi = zeros(blockp)\n apply(Op(\"S+\", 1.0, 2), psi, phi)\n @show inner(ops, psi)\n @show inner(ops, phi)\nend\nCreates a numerical matrix with real (matrix) or complex (matrixC) coefficients given an Op or OpSum on a certain block.
Sources matrix.hpp, matrix.cpp
"},{"location":"documentation/algebra/matrix/#definition","title":"Definition","text":"A matrix can be created in two ways:
Only the input block on which the operator is defined is given. The output block is calculated and eventually created automatically.
C++Juliatemplate <class block_t>\narma::mat matrix(Op const &op, block_t const &block);\n\ntemplate <class block_t>\narma::mat matrix(OpSum const &ops, block_t const &block);\n\ntemplate <class block_t>\narma::cx_mat matrixC(Op const &op, block_t const &block);\n\ntemplate <typename block_t>\narma::cx_mat matrixC(OpSum const &ops, block_t const &block);\nmatrix(op::Op, block::Block)\nmatrix(ops::OpSum, block::Block)\nThe output block is also handed as an argument. The compatibility of quantum numbers is checked. This way the output block is not created automatically and, thus, can be used to save computation time if the output block appears repeatedly in the computation.
C++Juliatemplate <class block_t>\narma::mat matrix(Op const &op, block_t const &block_in, \n block_t const &block_out);\n\ntemplate <class block_t>\narma::mat matrix(OpSum const &ops, block_t const &block_in, \n block_t const &block_out);\n\ntemplate <class block_t>\narma::cx_mat matrixC(Op const &op, block_t const &block_in, \n block_t const &block_out);\n\ntemplate <typename block_t>\narma::cx_mat matrixC(OpSum const &ops, block_t const &block_in, \n block_t const &block_out);\nmatrix(op::Op, block_in::Block, block_out::Block)\nmatrix(ops::OpSum, block_in::Block, block_out::Block)\nComment: In Julia, depending on whether a real/complex matrix is generated also a real/complex matrix is returned. The C++ version has to return a fixed type. If a real matrix is desired, use the function matrix. If a complex matrix is desired, use the function matrixC.
// Creates matrix H_{k=2} in Eq (18.23) of https://link.springer.com/content/pdf/10.1007/978-3-540-74686-7_18.pdf\nint N = 4;\nint nup = 3;\nint ndn = 2;\n\n// Define a Hubbard chain model\nauto ops = OpSum();\nfor (int i=0; i< N; ++i){\n ops += \"T\" * Op(\"Hop\", {i, (i+1) % N});\n}\nops+= \"U\" * Op(\"HubbardU\");\nops[\"T\"] = 1.0;\nops[\"U\"] = 5.0;\n\n// Create the a permutation group\nauto p1 = Permutation({0, 1, 2, 3});\nauto p2 = Permutation({1, 2, 3, 0});\nauto p3 = Permutation({2, 3, 0, 1});\nauto p4 = Permutation({3, 0, 1, 2});\nauto group = PermutationGroup({p1, p2, p3, p4});\nauto irrep = Representation(group, arma::vec{1.0, -1.0, 1.0, -1.0});\nauto block = Electron(N, nup, ndn, irrep);\nauto H = matrix(ops, block);\nH.print();\nlet\n # Creates matrix H_{k=2} in Eq (18.23) of https://link.springer.com/content/pdf/10.1007/978-3-540-74686-7_18.pdf\n N = 4\n nup = 3\n ndn = 2\n\n # Define a Hubbard chain model\n ops = OpSum()\n for i in 1:N\n ops += Op(\"Hop\", \"T\", [i, mod1(i+1, N)])\n end\n ops[\"T\"] = 1.0;\n ops[\"U\"] = 5.0;\n\n # Create the a permutation group\n p1 = Permutation([1, 2, 3, 4])\n p2 = Permutation([2, 3, 4, 1])\n p3 = Permutation([3, 4, 1, 2])\n p4 = Permutation([4, 1, 2, 3])\n group = PermutationGroup([p1, p2, p3, p4])\n irrep = Representation([1, -1, 1, -1])\n block = Electron(N, nup, ndn, group, irrep)\n\n H = matrix(ops, block)\n display(H)\nend\nComputes the groud state energy and the ground state of a Hermitian operator on a block by using an iterative Lanczos algorithm. This function is a shortcut for the eigs_lanczos function. We refer to eigs_lanczos for further details on the algorithm and the convergence criterion.
Sources sparse_diag.hpp, sparse_diag.cpp
"},{"location":"documentation/algorithms/eig0/#definition","title":"Definition","text":"C++Juliastd::tuple<double, State> eig0(OpSum const &ops, Block const &block,\n double precision = 1e-12, int64_t max_iterations = 1000,\n bool force_complex = false, int64_t random_seed = 42);\nfunction eig0(\n ops::OpSum,\n block::Block;\n precision::Real = 1e-12,\n maxiter::Int64 = 1000,\n force_complex::Bool = false,\n seed::Int64 = 42,\n)\nops State groundstate"},{"location":"documentation/algorithms/eig0/#usage-example","title":"Usage Example","text":"C++Julia int N = 8;\nint nup = N / 2;\nauto block = Spinhalf(N, nup);\n\n// Define the nearest-neighbor Heisenberg model\nauto ops = OpSum();\nfor (int i=0; i<N; ++i) {\n ops += \"J\" * Op(\"SdotS\", {i, (i+1) % N});\n}\nops[\"J\"] = 1.0;\nauto [e0, gs] = eig0(ops, block);\nlet \n N = 8\n nup = N \u00f7 2\n block = Spinhalf(N, nup)\n\n # Define the nearest-neighbor Heisenberg model\n ops = OpSum()\n for i in 1:N\n ops += Op(\"HB\", \"J\", [i, mod1(i+1, N)])\n end\n ops[\"J\"] = 1.0;\n\n e0, gs = eig0(ops, block);\nend\nPerforms an iterative eigenvalue calculation building eigenvectors using the Lanczos algorithm. Returns the tridiagonal matrix, eigenvalues, number of iterations and the stopping criterion. The Lanczos interations are performed twice, where at the second run the eigenvectors are built.
Sources eigs_lanczos.hpp, eigs_lanczos.cpp
"},{"location":"documentation/algorithms/eigs_lanczos/#definition","title":"Definition","text":"The Lanczos algorithm can be run in two distinct ways:
A random intial state \\(|\\psi_0\\rangle = |r\\rangle\\) with normal distributed entries is used.
C++Juliaeigs_lanczos_result_t\neigs_lanczos(OpSum const &ops, Block const &block, int64_t neigvals = 1,\n double precision = 1e-12, int64_t max_iterations = 1000,\n bool force_complex = false, double deflation_tol = 1e-7,\n int64_t random_seed = 42);\nfunction eigs_lanczos(\n ops::OpSum,\n block::Block;\n neigvals::Int64 = 1,\n precision::Float64 = 1e-12,\n max_iterations::Int64 = 1000,\n force_complex::Bool = false,\n deflation_tol::Float64 = 1e-7,\n random_seed::Int64 = 42,\n)\nThe initial state \\(|\\psi_0\\rangle\\) is explicitly specified
C++eigs_lanczos_result_t \neigs_lanczos(OpSum const &ops, State &state0, int64_t neigvals = 1,\n double precision = 1e-12, int64_t max_iterations = 1000,\n bool force_complex = false, double deflation_tol = 1e-7);\nA struct with the following entries
Entry Description alphas diagonal elements of the tridiagonal matrix betas off-diagonal elements of the tridiagonal matrix eigenvalues the computed Ritz eigenvalues of the tridiagonal matrix eigenvectors State of shape $D \\times $neigvals holding all low-lying eigenvalues up to neigvals niterations number of iterations performed criterion string denoting the reason why the algorithm stopped"},{"location":"documentation/algorithms/eigs_lanczos/#convergence-criterion","title":"Convergence criterion","text":"The algorithm terminates if the \\(k\\)-th (\\(k\\) is the argument neigvals) approximate eigenvalue changes only by a fraction smaller than \\(\\epsilon\\) (\\(k\\) is the argument precision), i.e.
Here, \\(\\tilde{e}_k^{(n)}\\) denotes the Lanczos approximation to the \\(k\\)-th eigenvalue after \\(n\\) iterations.
"},{"location":"documentation/algorithms/eigs_lanczos/#usage-example","title":"Usage Example","text":"C++int N = 8;\nint nup = N / 2;\nauto block = Spinhalf(N, nup);\n\n// Define the nearest-neighbor Heisenberg model\nauto ops = OpSum();\nfor (int i=0; i<N; ++i) {\n ops += \"J\" * Op(\"SdotS\", {i, (i+1) % N});\n}\nops[\"J\"] = 1.0;\n\n// With random intial state\nauto res = eigs_lanczos(ops, block);\nXDIAG_SHOW(res.alphas);\nXDIAG_SHOW(res.betas);\nXDIAG_SHOW(res.eigenvalues);\nXDIAG_SHOW(res.eigenvectors);\nComputes the groud state energy of a Hermitian operator on a block by using an iterative Lanczos algorithm. This function is a shortcut for the eigvals_lanczos function. We refer to eigvals_lanczos for further details on the algorithm and the convergence criterion.
Sources sparse_diag.hpp, sparse_diag.cpp
"},{"location":"documentation/algorithms/eigval0/#definition","title":"Definition","text":"C++Juliadouble eigval0(OpSum const &ops, Block const &block, double precision = 1e-12,\n int64_t max_iterations = 1000, bool force_complex = false,\n int64_t random_seed = 42);\nfunction eigval0(\n ops::OpSum,\n block::Block;\n precision::Real = 1e-12,\n maxiter::Integer = 1000,\n force_complex::Bool = false,\n seed::Integer = 42,\n)\nops"},{"location":"documentation/algorithms/eigval0/#usage-example","title":"Usage Example","text":"C++Julia int N = 8;\nint nup = N / 2;\nauto block = Spinhalf(N, nup);\n\n// Define the nearest-neighbor Heisenberg model\nauto ops = OpSum();\nfor (int i=0; i<N; ++i) {\n ops += \"J\" * Op(\"SdotS\", {i, (i+1) % N});\n}\nops[\"J\"] = 1.0;\ndouble e0 = eigval0(ops, block);\nlet \n N = 8\n nup = N \u00f7 2\n block = Spinhalf(N, nup)\n\n # Define the nearest-neighbor Heisenberg model\n ops = OpSum()\n for i in 1:N\n ops += Op(\"HB\", \"J\", [i, mod1(i+1, N)])\n end\n ops[\"J\"] = 1.0\n\n e0 = eigval0(ops, block);\nend\nPerforms an iterative eigenvalue calculation using the Lanczos algorithm. Returns the tridiagonal matrix, eigenvalues, number of iterations and the stopping criterion.
Sources eigvals_lanczos.hpp, eigvals_lanczos.cpp
"},{"location":"documentation/algorithms/eigvals_lanczos/#definition","title":"Definition","text":"The Lanczos algorithm can be run in thre distinct ways:
A random intial state \\(|\\psi_0\\rangle = |r\\rangle\\) with normal distributed entries is used.
C++Juliaeigvals_lanczos_result_t\neigvals_lanczos(OpSum const &ops, Block const &block, int64_t neigvals = 1,\n double precision = 1e-12, int64_t max_iterations = 1000,\n bool force_complex = false, double deflation_tol = 1e-7,\n int64_t random_seed = 42);\nfunction eigvals_lanczos(\n ops::OpSum,\n block::Block;\n neigvals::Int64 = 1,\n precision::Float64 = 1e-12,\n max_iterations::Int64 = 1000,\n force_complex::Bool = false,\n deflation_tol::Float64 = 1e-7,\n random_seed::Int64 = 42,\n)\nThe initial state \\(|\\psi_0\\rangle\\) is explicitly specified.
C++eigvals_lanczos_result_t \neigvals_lanczos(OpSum const &ops, State psi0, int64_t neigvals = 1,\n double precision = 1e-12, int64_t max_iterations = 1000,\n bool force_complex = false, double deflation_tol = 1e-7);\nNotice this version copies the initial state, which requires memory but keeps the orginal state intact.
The initial state \\(|\\psi_0\\rangle\\) is explicitly specified and overwritten in the process. This version can save memory, but the initial state \\(|\\psi_0\\rangle\\) cannot be used later.
C++eigvals_lanczos_result_t \neigvals_lanczos_inplace(OpSum const &ops, State &psi0, int64_t neigvals = 1,\n double precision = 1e-12, int64_t max_iterations = 1000,\n bool force_complex = false, double deflation_tol = 1e-7);\nA struct of type eigvals_lanczos_result_t with the following entries.
The algorithm terminates if the \\(k\\)-th (\\(k\\) is the argument neigvals) approximate eigenvalue changes only by a fraction smaller than \\(\\epsilon\\) (\\(k\\) is the argument precision), i.e.
Here, \\(\\tilde{e}_k^{(n)}\\) denotes the Lanczos approximation to the \\(k\\)-th eigenvalue after \\(n\\) iterations.
"},{"location":"documentation/algorithms/eigvals_lanczos/#usage-example","title":"Usage Example","text":"C++int N = 8;\nint nup = N / 2;\nauto block = Spinhalf(N, nup);\n\n// Define the nearest-neighbor Heisenberg model\nauto ops = OpSum();\nfor (int i=0; i<N; ++i) {\n ops += \"J\" * Op(\"SdotS\", {i, (i+1) % N});\n}\nops[\"J\"] = 1.0;\n\n// With random intial state\nauto res = eigvals_lanczos(ops, block);\nXDIAG_SHOW(res.alphas);\nXDIAG_SHOW(res.betas);\nXDIAG_SHOW(res.eigenvalues);\n\n// With specific initial state\nauto psi0 = product_state(block, {\"Up\", \"Dn\", \"Up\", \"Dn\", \"Up\", \"Dn\", \"Up\", \"Dn\"});\nauto res2 = eigvals_lanczos(ops, psi0);\nXDIAG_SHOW(res.alphas);\nXDIAG_SHOW(res.betas);\nXDIAG_SHOW(res.eigenvalues);\nComputes the exponential of a Hermitian operator \\(H\\) with an arbitrary real or complex prefactor \\(z\\) applied to a State \\(\\vert \\psi_0\\rangle\\),
\\[\\vert \\psi(z) \\rangle = e^{z(H - \\delta)} \\vert \\psi_0\\rangle.\\]Here, \\(\\delta\\) denotes a real number shifting the spectrum of \\(H\\). The algorithm implemented is described in the following publication.
On Krylov Subspace Approximations to the Matrix Exponential Operator Marlis Hochbruck and Christian Lubich SIAM Journal on Numerical Analysis, Vol. 34, Iss. 5 (1997) DOI: 10.1137/S0036142995280572
Sources evolve_lanczos.hpp, evolve_lanczos.cpp
"},{"location":"documentation/algorithms/evolve_lanczos/#definition","title":"Definition","text":"The method is provided in two variants:
Returning a new state while the input state remains untouched. This variant is safe to use and simple to code.
C++evolve_lanczos_result_t\nevolve_lanczos(OpSum const &H, State psi, double z, double precision = 1e-12,\n double shift = 0., bool normalize = false,\n int64_t max_iterations = 1000, double deflation_tol = 1e-7);\n\nevolve_lanczos_result_t\nevolve_lanczos(OpSum const &H, State psi, complex z, double precision = 1e-12,\n double shift = 0., bool normalize = false,\n int64_t max_iterations = 1000, double deflation_tol = 1e-7);\nAn inplace variant evolve_lanczos_inplace, where the input state is overwritten and contains the time evolved state upon exit. This version is more memory efficient than evolve_lanczos.
evolve_lanczos_inplace_result_t\nevolve_lanczos_inplace(OpSum const &H, State &psi, double z, \n double precision = 1e-12, double shift = 0.,\n bool normalize = false, int64_t max_iterations = 1000, \n double deflation_tol = 1e-7);\n\nevolve_lanczos_inplace_result_t\nevolve_lanczos_inplace(OpSum const &H, State &psi, complex z, \n double precision = 1e-12, double shift = 0.,\n bool normalize = false, int64_t max_iterations = 1000, \n double deflation_tol = 1e-7);\nThe parameter shift can be used to turn all eigenvalues of the matrix \\(H - \\delta \\;\\textrm{Id}\\) positive whenever \\(\\delta < E_0\\), where \\(E_0\\) denotes the ground state energy of \\(H\\).
A struct with the following entries
Entry Description alphas diagonal elements of the Lanczos tridiagonal matrix betas off-diagonal elements of the Lanczos tridiagonal matrix eigenvalues the computed Ritz eigenvalues of the tridiagonal matrix niterations number of iterations performed criterion string denoting the reason why the Lanczosalgorithm stopped state time-evolved State \\(\\vert \\psi(t)\\rangle\\) (not defined for inplace variant)"},{"location":"documentation/algorithms/evolve_lanczos/#convergence-criterion","title":"Convergence criterion","text":"The algorithm is estimating the following error,
\\[ \\varepsilon = \\parallel \\vert \\tilde{\\psi}(t)\\rangle - e^{z(H - \\delta)} \\vert \\psi_0\\rangle \\parallel_2, \\]where \\(\\vert \\tilde{\\psi}(t) \\rangle\\) denotes the approximation computed during the algorithm. As the exact solution is not available this error is estimated using the method described by Algorithm 2 in
Expokit: A Software Package for Computing Matrix Exponentials Roger B. Sidje ACM Trans. Math. Softw., 24(1):130-156, 1998. (1998) DOI: 10.1145/285861.285868
"},{"location":"documentation/algorithms/evolve_lanczos/#usage-example","title":"Usage Example","text":"C++int N = 8;\nint nup = N / 2;\nauto block = Spinhalf(N, nup);\n\n// Define the nearest-neighbor Heisenberg model\nauto ops = OpSum();\nfor (int i=0; i<N; ++i) {\n ops += Op(\"SdotS\", {i, (i+1) % N});\n}\n\n// Compute ground state energy\ndouble e0 = eigval0(ops, block);\n\nauto psi0 = product_state(block, {\"Up\", \"Dn\", \"Up\", \"Dn\", \"Up\", \"Dn\", \"Up\", \"Dn\"});\ndouble time = 1.0;\ndouble precision = 1e-12;\nauto res = evolve_lanczos(ops, psi0, time, precision, e0, true, 500);\nXDIAG_SHOW(res.alphas);\nXDIAG_SHOW(res.betas);\nComputes the imaginary-time evolution,
\\[\\vert \\psi(\\tau) \\rangle = e^{-(H - \\delta) \\tau} \\vert \\psi_0\\rangle,\\]of a State \\(\\vert \\psi_0 \\rangle\\) and a Hermitian operator \\(H\\) using an iterative algorithm. \\(\\delta\\) here denotes a real number which can be chosen as the ground state energy \\(\\delta=E_0\\) of \\(H\\).
Sources imaginary_time_evolve.hpp, imaginary_time_evolve.cpp
"},{"location":"documentation/algorithms/imaginary_time_evolve/#definition","title":"Definition","text":"The method is provided in two variants:
Returning a new state while the input state remains untouched. This variant is safe to use and simple to code.
C++State imaginary_time_evolve(OpSum const &H, State psi0, double time,\n double precision = 1e-12, double shift = 0.);\nAn inplace variant imaginary_time_evolve_inplace, where the input state is overwritten and contains the time evolved state upon exit. This version is more memory efficient than imaginary_time_evolve.
void imaginary_time_evolve_inplace(OpSum const &H, State &psi0, double time,\n double precision = 1e-12, shift = 0.);\nThe routine calls the subroutine evolve_lanczos implementing a Lanczos algorithm to perform the evolution. This routine can also be called explicitly if more control is desired. Please also confer to the page evolve_lanczos for further details on the specifics of the algorithm. The parameter shift can be used to turn all eigenvalues of the matrix \\(H - \\delta \\;\\textrm{Id}\\) positive whenever \\(\\delta < E_0\\), where \\(E_0\\) denotes the ground state energy of \\(H\\).
int N = 8;\nint nup = N / 2;\nauto block = Spinhalf(N, nup);\n\n// Define the nearest-neighbor Heisenberg model\nauto ops = OpSum();\nfor (int i=0; i<N; ++i) {\n ops += Op(\"SdotS\", {i, (i+1) % N});\n}\n\n// Compute ground state energy\ndouble e0 = eigval0(ops, block);\n\nauto psi0 = product_state(block, {\"Up\", \"Dn\", \"Up\", \"Dn\", \"Up\", \"Dn\", \"Up\", \"Dn\"});\ndouble time = 1.0;\ndouble precision = 1e-12;\nauto psi = imaginary_time_evolve(ops, psi0, time, precision, e0);\nimaginary_time_evolve_inplace(ops, psi0, time, precision, e0);\nXDIAG_SHOW(isapprox(psi0, psi));\nComputes the real-time evolution,
\\[\\vert \\psi(t) \\rangle = e^{-iHt} \\vert \\psi_0\\rangle,\\]of a State \\(\\vert \\psi_0 \\rangle\\) and a Hermitian operator \\(H\\) using an iterative algorithm.
Sources time_evolve.hpp, time_evolve.cpp
"},{"location":"documentation/algorithms/time_evolve/#definition","title":"Definition","text":"The method is provided in two variants:
Returning a new state while the input state remains untouched. This variant is safe to use and simple to code.
C++State time_evolve(OpSum const &H, State psi0, double time,\n double precision = 1e-12,\n std::string algorithm = \"lanczos\");\nAn inplace variant time_evolve_inplace, where the input state is overwritten and contains the time evolved state upon exit. This version is more memory efficient than time_evolve.
void time_evolve_inplace(OpSum const &H, State &psi0, double time,\n double precision = 1e-12,\n std::string algorithm = \"lanczos\");\nlanczos or expokit lanczos The algorithm parameter decised which backend is run. If lanczos is chosen, the evolve_lanczos routine is called with the standard arguments. Alternatively, expokit chooses the time_evolve_expokit routine. For a detailed documentation of the algorithms we refer to the evolve_lanczos and time_evolve_expokit pages. Broadly speaking, the expokit can yield higher precision states at arbitrarily long times at the cost of increased memory and computing time. In practice, we recommend analysing the effect of the precision parameters on the time evolution series obtained in both cases.
int N = 8;\nint nup = N / 2;\nauto block = Spinhalf(N, nup);\n\n// Define the nearest-neighbor Heisenberg model\nauto ops = OpSum();\nfor (int i=0; i<N; ++i) {\n ops += Op(\"SdotS\", {i, (i+1) % N});\n}\n\nauto psi0 = product_state(block, {\"Up\", \"Dn\", \"Up\", \"Dn\", \"Up\", \"Dn\", \"Up\", \"Dn\"});\ndouble time = 1.0;\nauto psi = time_evolve(ops, psi0, time);\ntime_evolve_inplace(ops, psi0, time);\nXDIAG_SHOW(isapprox(psi0, psi));\nComputes the real-time evolution,
\\[\\vert \\psi(t) \\rangle = e^{-iHt} \\vert \\psi_0\\rangle,\\]of a State \\(\\vert \\psi_0 \\rangle\\) and a Hermitian operator \\(H\\) using the iterative algorithm implemented by Expokit
Expokit: A Software Package for Computing Matrix Exponentials Roger B. Sidje ACM Trans. Math. Softw., 24(1):130-156, 1998. (1998) DOI: 10.1145/285861.285868
The algorithm features automatic stepsize control and computes approximate solutions with high precision according to our tests. Yet, the evolve_lanczos implementation is currently faster and more memory efficient.
Sources time_evolve_expokit.hpp, time_evolve_expokit.cpp
"},{"location":"documentation/algorithms/time_evolve_expokit/#definition","title":"Definition","text":"The method is provided in two variants:
Returning a new state while the input state remains untouched. This variant is safe to use and simple to code.
C++time_evolve_expokit_return_t time_evolve_expokit(\n OpSum const &ops, State state, double time, double precision = 1e-12,\n int64_t m = 30, double anorm = 0., int64_t nnorm = 2);\nAn inplace variant time_evolve_expokit_inplace, where the input state is overwritten and contains the time evolved state upon exit. This version is more memory efficient than time_evolve_expokit.
time_evolve_expokit_inplace_return_t time_evolve_expokit_inplace(\n OpSum const &ops, State &state, double time, double precision = 1e-12,\n int64_t m = 30, double anorm = 0., int64_t nnorm = 2);\nA struct with the following entries
Entry Description error the computed error estimate during evolution hump the \"hump\" as defined in Expokit 10.1145/285861.285868 state time-evolved State \\(\\vert \\psi(t)\\rangle\\) (not defined for inplace variant)"},{"location":"documentation/algorithms/time_evolve_expokit/#convergence-criterion","title":"Convergence criterion","text":"The algorithm is estimating the following error,
\\[ \\varepsilon = \\parallel \\vert \\tilde{\\psi}(t)\\rangle - e^{z(H - \\delta)} \\vert \\psi_0\\rangle \\parallel_2, \\]where \\(\\vert \\tilde{\\psi}(t) \\rangle\\) denotes the approximation computed during the algorithm. As the exact solution is not available this error is estimated using the method described by Algorithm 2 in
Expokit: A Software Package for Computing Matrix Exponentials Roger B. Sidje ACM Trans. Math. Softw., 24(1):130-156, 1998. (1998) DOI: 10.1145/285861.285868
"},{"location":"documentation/algorithms/time_evolve_expokit/#usage-example","title":"Usage Example","text":"C++int N = 8;\nint nup = N / 2;\nauto block = Spinhalf(N, nup);\n\n// Define the nearest-neighbor Heisenberg model\nauto ops = OpSum();\nfor (int i=0; i<N; ++i) {\n ops += Op(\"SdotS\", {i, (i+1) % N});\n}\n\nauto psi0 = product_state(block, {\"Up\", \"Dn\", \"Up\", \"Dn\", \"Up\", \"Dn\", \"Up\", \"Dn\"});\ndouble time = 1.0;\nauto res1 = time_evolve_expokit(ops, psi0, time, 20);\nauto res2 = time_evolve_expokit_inplace(ops, psi0, time, 20);\nXDIAG_SHOW(isapprox(psi0, res1.state));\nXDIAG_SHOW(res1.error);\nXDIAG_SHOW(res1.hump);\nA block in an Electron (fermions with \\(\\uparrow, \\downarrow\\) spin) Hilbert space.
Sources electron.hpp, electron.cpp
"},{"location":"documentation/blocks/electron/#constructors","title":"Constructors","text":"C++JuliaElectron(int64_t nsites, std::string backend = \"auto\");\nElectron(int64_t nsites, int64_t nup, int64_t ndn, std::string backend = \"auto\");\nElectron(int64_t nsites, Representation irrep, std::string backend = \"auto\");\nElectron(int64_t nsites, int64_t nup, int64_t ndn, Representation irrep, std::string backend = \"auto\");\nElectron(nsites::Integer, backend::String)\nElectron(nsites::Integer, nup::Integer, ndn::Integer, backend::String)\nElectron(nsites::Integer, irrep::Representation, backend::String)\nElectron(nsites::Integer, nup::Integer, ndn::Integer, irrep::Representation, backend::String)\nauto The parameter backend chooses how the block is coded internally. By using the default parameter auto the backend is chosen automatically. Alternatives are 32bit, 64bit
An Electron block can be iterated over, where at each iteration a ProductState representing the corresponding basis state is returned.
C++Juliaauto block = Electron(4, 2, 2);\nfor (auto pstate : block) {\n Log(\"{} {}\", to_string(pstate), index(block, pstate));\n}\nblock = Electron(4, 2, 2)\nfor pstate in block\n @show pstate, index(block, pstate) \nend\nReturns the index of a given ProductState in the basis of the tJ block.
C++Juliaint64_t index(tJ const &block, ProductState const &pstate);\nindex(block::tJ, pstate::ProductState)\n1-indexing
In the C++ version, the index count starts from \"0\" whereas in Julia the index count starts from \"1\".
"},{"location":"documentation/blocks/electron/#nsites","title":"nsites","text":"Returns the number of sites of the block.
C++Juliaint64_t nsites(tJ const &block);\nnsites(block::tJ)\nReturns the size of the block, i.e. its dimension.
C++Juliaint64_t size(tJ const &block) const;\nsize(block::tJ)\nReturns the dimension of the block, same as \"size\" for non-distributed blocks.
C++Juliaint64_t dim(tJ const &block) const;\ndim(block::tJ)\nReturns whether the block can be used with real arithmetic. Complex arithmetic is needed when a Representation is genuinely complex.
C++Juliabool isreal(tJ const &block);\nisreal(block::tJ)\nint N = 4;\nint nup = 2;\nint ndn = 1;\n\n// without number conservation\nauto block = Electron(N);\nXDIAG_SHOW(block);\n\n// with number conservation\nauto block_np = Electron(N, nup, ndn);\nXDIAG_SHOW(block_np);\n\n// with symmetries, without number conservation\nauto p1 = Permutation({0, 1, 2, 3});\nauto p2 = Permutation({1, 2, 3, 0});\nauto p3 = Permutation({2, 3, 0, 1});\nauto p4 = Permutation({3, 0, 1, 2});\nauto group = PermutationGroup({p1, p2, p3, p4});\nauto irrep = Representation(group, arma::vec{1, -1, 1, -1});\nauto block_sym = Electron(N, irrep);\nXDIAG_SHOW(block_sym);\n\n// with symmetries and number conservation\nauto block_sym_np = Electron(N, nup, ndn, irrep);\nXDIAG_SHOW(block_sym_np);\nXDIAG_SHOW(block_sym_np.nsites());\nXDIAG_SHOW(block_sym_np.size());\n\n// Iteration\nfor (auto pstate : block_sym_np) {\n Log(\"{} {}\", to_string(pstate), index(block_sym_np, pstate));\n}\nN = 4\nnup = 2\nndn = 1\n\n# without number conservation\nblock = Electron(N)\n@show block\n\n# with number conservation\nblock_np = Electron(N, nup, ndn)\n@show block_np\n\n# with symmetries, without number conservation\np1 = Permutation([1, 2, 3, 4])\np2 = Permutation([2, 3, 4, 1])\np3 = Permutation([3, 4, 1, 2])\np4 = Permutation([4, 1, 2, 3])\ngroup = PermutationGroup([p1, p2, p3, p4])\nrep = Representation([1, -1, 1, -1])\nblock_sym = Electron(N, group, rep)\n@show block_sym\n\n# with symmetries and number conservation\nblock_sym_np = Electron(N, nup, ndn, group, rep)\n@show block_sym_np\n\n@show nsites(block_sym_np)\n@show size(block_sym_np)\n\n# Iteration\nfor pstate in block_sym_np\n @show pstate, index(block_sym_np, pstate)\nend\n@show permutation_group(block_sym_np)\n@show irrep(block_sym_np)\nA block in a spin \\(S=1/2\\) Hilbert space.
Sources spinhalf.hpp, spinhalf.cpp
"},{"location":"documentation/blocks/spinhalf/#constructors","title":"Constructors","text":"C++JuliaSpinhalf(int64_t nsites, std::string backend = \"auto\");\nSpinhalf(int64_t nsites, int64_t nup, std::string backend = \"auto\");\nSpinhalf(int64_t nsites, Representation const &irrep, std::string backend = \"auto\");\nSpinhalf(int64_t nsites, int64_t nup, Representation const &irrep, std::string backend = \"auto\");\nSpinhalf(nsites::Integer, backend::String)\nSpinhalf(nsites::Integer, nup::Integer, backend::String)\nSpinhalf(nsites::Integer, irrep::Representation, backend::String)\nSpinhalf(nsites::Integer, nup::Integer, irrep::Representation, backend::String)\nauto The parameter backend chooses how the block is coded internally. By using the default parameter auto the backend is chosen automatically. Alternatives are 32bit, 64bit, 1sublattice, 2sublattice, 3sublattice, 4sublattice, and 5sublattice. The backends xsublattice implement the sublattice coding algorithm described in Wietek, L\u00e4uchli, Phys. Rev. E 98, 033309 (2018). The sublattice coding algorithms impose certain constraints on the symmetries used, as described in the reference.
An Spinhalf block can be iterated over, where at each iteration a ProductState representing the corresponding basis state is returned.
C++Juliaauto block = Spinhalf(4, 2);\nfor (auto pstate : block) {\n Log(\"{} {}\", to_string(pstate), block.index(pstate));\n}\nblock = Spinhalf(4, 2)\nfor pstate in block\n @show pstate, index(block, pstate) \nend\nReturns the index of a given ProductState in the basis of the Spinhalf block.
C++Juliaint64_t index(Spinhalf const &block, ProductState const &pstate);\nindex(block::Spinhalf, pstate::ProductState)\n1-indexing
In the C++ version, the index count starts from \"0\" whereas in Julia the index count starts from \"1\".
"},{"location":"documentation/blocks/spinhalf/#nsites","title":"nsites","text":"Returns the number of sites of the block.
C++Juliaint64_t nsites(Spinhalf const &block);\nnsites(block::Spinhalf)\nReturns the size of the block, i.e. its dimension.
C++Juliaint64_t size(Spinhalf const &block) const;\nsize(block::Spinhalf)\nReturns the dimension of the block, same as \"size\" for non-distributed blocks.
C++Juliaint64_t dim(Spinhalf const &block) const;\ndim(block::Spinhalf)\nReturns whether the block can be used with real arithmetic. Complex arithmetic is needed when a Representation is genuinely complex.
C++Juliabool isreal(Spinhalf const &block);\nisreal(block::Spinhalf)\nint N = 4;\nint nup = 2;\n\n// without Sz conservation\nauto block = Spinhalf(N);\nXDIAG_SHOW(block);\n\n// with Sz conservation\nauto block_sz = Spinhalf(N, nup);\nXDIAG_SHOW(block_sz);\n\n// with symmetries, without Sz\nPermutation p1 = {0, 1, 2, 3};\nPermutation p2 = {1, 2, 3, 0};\nPermutation p3 = {2, 3, 0, 1};\nPermutation p4 = {3, 0, 1, 2};\nauto group = PermutationGroup({p1, p2, p3, p4});\nauto irrep = Representation(group, arma::vec({1, -1, 1, -1}));\nauto block_sym = Spinhalf(N, irrep);\nXDIAG_SHOW(block_sym);\n\n// with symmetries and Sz\nauto block_sym_sz = Spinhalf(N, nup, irrep);\nXDIAG_SHOW(block_sym_sz);\n\nXDIAG_SHOW(block_sym_sz.nsites());\nXDIAG_SHOW(block_sym_sz.size());\n\n// Iteration\nfor (auto pstate : block_sym_sz) {\n Log(\"{} {}\", to_string(pstate), index(block_sym_sz, pstate));\n}\nN = 4\nnup = 2\n\n# without Sz conservation\nblock = Spinhalf(N)\n@show block\n\n\n# with Sz conservation\nblock_sz = Spinhalf(N, nup)\n@show block_sz\n\n# with symmetries, without Sz\np1 = Permutation([1, 2, 3, 4])\np2 = Permutation([2, 3, 4, 1])\np3 = Permutation([3, 4, 1, 2])\np4 = Permutation([4, 1, 2, 3])\ngroup = PermutationGroup([p1, p2, p3, p4])\nrep = Representation([1, -1, 1, -1])\nblock_sym = Spinhalf(N, group, rep)\n@show block_sym\n\n# with symmetries and Sz\nblock_sym_sz = Spinhalf(N, nup, group, rep)\n@show block_sym_sz\n\n@show nsites(block_sym_sz)\n@show size(block_sym_sz)\n\n# Iteration\nfor pstate in block_sym_sz\n @show pstate, index(block_sym_sz, pstate)\nend\n@show permutation_group(block_sym_sz)\n@show irrep(block_sym_sz)\nA block in a spin \\(S=1/2\\) Hilbert space with distributed computing capabilities.
Sources spinhalf_distributed.hpp, spinhalf_distributed.cpp
"},{"location":"documentation/blocks/spinhalf_distributed/#constructors","title":"Constructors","text":"C++SpinhalfDistributed(int64_t nsites, int64_t nup, std::string backend = \"auto\");\nauto The parameter backend chooses how the block is coded internally. By using the default parameter auto the backend is chosen automatically. Alternatives are 32bit, 64bit.
An SpinhalfDistributed block can be iterated over, where at each iteration a ProductState representing the corresponding basis state is returned.
C++auto block = SpinhalfDistributed(4, 2);\nfor (auto pstate : block) {\n Log(\"{} {}\", to_string(pstate), block.index(pstate));\n}\nReturns the index of a given ProductState in the basis of the SpinhalfDistributed block.
C++int64_t index(SpinhalfDistributed const &block, ProductState const &pstate);\nReturns the number of sites of the block.
C++int64_t nsites(SpinhalfDistributed const &block);\nReturns the size of the block on a local process.
C++int64_t size(SpinhalfDistributed const &block) const;\nReturns the dimension of the block, i.e. the sum of all sizes across all processes.
C++int64_t dim(SpinhalfDistributed const &block) const;\nReturns whether the block can be used with real arithmetic. Complex arithmetic is needed when a Representation is genuinely complex.
C++bool isreal(SpinhalfDistributed const &block);\nA block in a \\(t-J\\) type Hilbert space, i.e. fermions with \\(\\uparrow, \\downarrow\\) spin excluding doubly occupied sites.
Sources tj.hpp, tj.cpp
"},{"location":"documentation/blocks/tJ/#constructors","title":"Constructors","text":"C++JuliatJ(int64_t nsites, int64_t nup, int64_t ndn, std::string backend = \"auto\");\ntJ(int64_t nsites, int64_t nup, int64_t ndn, Representation const &irrep, std::string backend = \"auto\");\ntJ(nsites::Integer, nup::Integer, ndn::Integer, backend::String)\ntJ(nsites::Integer, nup::Integer, ndn::Integer, irrep::Representation, backend::String)\nauto The parameter backend chooses how the block is coded internally. By using the default parameter auto the backend is chosen automatically. Alternatives are 32bit, 64bit.
An tJ block can be iterated over, where at each iteration a ProductState representing the corresponding basis state is returned.
C++Juliaauto block = tJ(4, 2, 1);\nfor (auto pstate : block) {\n Log(\"{} {}\", to_string(pstate), index(block, pstate));\n}\nblock = tJ(4, 2, 1)\nfor pstate in block\n @show pstate, index(block, pstate) \nend\nReturns the index of a given ProductState in the basis of the tJ block.
C++Juliaint64_t index(tJ const &block, ProductState const &pstate);\nindex(block::tJ, pstate::ProductState)\n1-indexing
In the C++ version, the index count starts from \"0\" whereas in Julia the index count starts from \"1\".
"},{"location":"documentation/blocks/tJ/#nsites","title":"nsites","text":"Returns the number of sites of the block.
C++Juliaint64_t nsites(tJ const &block);\nnsites(block::tJ)\nReturns the size of the block, i.e. its dimension.
C++Juliaint64_t size(tJ const &block) const;\nsize(block::tJ)\nReturns the dimension of the block, same as \"size\" for non-distributed blocks.
C++Juliaint64_t dim(tJ const &block) const;\ndim(block::tJ)\nReturns whether the block can be used with real arithmetic. Complex arithmetic is needed when a Representation is genuinely complex.
C++Juliabool isreal(tJ const &block);\nisreal(block::tJ)\nint N = 4;\nint nup = 2;\nint ndn = 1;\n\n// without permutation symmetries\nauto block = tJ(N, nup, ndn);\nXDIAG_SHOW(block);\n\n// with permutation symmetries\nauto p1 = Permutation({0, 1, 2, 3});\nauto p2 = Permutation({1, 2, 3, 0});\nauto p3 = Permutation({2, 3, 0, 1});\nauto p4 = Permutation({3, 0, 1, 2});\nauto group = PermutationGroup({p1, p2, p3, p4});\nauto irrep = Representation(group, arma::vec{1, -1, 1, -1});\nauto block_sym = tJ(N, nup, ndn, irrep);\nXDIAG_SHOW(block_sym);\nXDIAG_SHOW(block_sym.nsites());\nXDIAG_SHOW(block_sym.size());\n\n// Iteration\nfor (auto pstate : block_sym) {\n Log(\"{} {}\", to_string(pstate), index(block_sym, pstate));\n}\nN = 4\nnup = 2\nndn = 1\n\n# without permutation symmetries\nblock = tJ(N, nup, ndn)\n@show block\n\n# with permutation symmetries\np1 = Permutation([1, 2, 3, 4])\np2 = Permutation([2, 3, 4, 1])\np3 = Permutation([3, 4, 1, 2])\np4 = Permutation([4, 1, 2, 3])\ngroup = PermutationGroup([p1, p2, p3, p4])\nrep = Representation([1, -1, 1, -1])\nblock_sym = tJ(N, nup, ndn, group, rep)\n@show block_sym\n\n@show nsites(block_sym)\n@show size(block_sym)\n\n# Iteration\nfor pstate in block_sym\n @show pstate, index(block_sym, pstate)\nend\n@show permutation_group(block_sym)\n@show irrep(block_sym)\nA block in a \\(t-J\\) type Hilbert space, i.e. fermions with \\(\\uparrow, \\downarrow\\) spin excluding doubly occupied sites with distributed computing capabilities.
Sources tj_distributed.hpp, tj_distributed.cpp
"},{"location":"documentation/blocks/tJ_distributed/#constructors","title":"Constructors","text":"C++tJDistributed(int64_t nsites, int64_t nup, int64_t ndn, std::string backend = \"auto\");\nauto The parameter backend chooses how the block is coded internally. By using the default parameter auto the backend is chosen automatically. Alternatives are 32bit, 64bit.
An tJDistributed block can be iterated over, where at each iteration a ProductState representing the corresponding basis state is returned.
C++auto block = tJDistributed(4, 2, 1);\nfor (auto pstate : block) {\n Log(\"{} {}\", to_string(pstate), index(block, pstate));\n}\nReturns the index of a given ProductState in the basis of the tJDistributed block.
C++int64_t index(tJDistributed const &block, ProductState const &pstate);\nReturns the number of sites of the block.
C++int64_t nsites(tJDistributed const &block);\nReturns the size of the block on a local process.
C++int64_t size(tJDistributed const &block) const;\nReturns the dimension of the block, i.e. the sum of all sizes across all processes.
C++int64_t dim(tJDistributed const &block) const;\nReturns whether the block can be used with real arithmetic. Complex arithmetic is needed when a Representation is genuinely complex.
C++bool isreal(tJDistributed const &block);\nDownload the source code using git
cd /path/to/where/xdiag/should/be\ngit clone https://github.com/awietek/xdiag.git\nCompile the default library
cd xdiag\ncmake -S . -B build\ncmake --build build\ncmake --install build\ninstall. Compile the distributed library
To use the distributed computing features of xdiag, the distributed library has to be built which requires MPI.
cd xdiag\ncmake -S . -B build -D XDIAG_DISTRIBUTED=On\ncmake --build build\ncmake --install build\nInfo
It might be necessary to explicitly define MPI compiler, e.g. mpicxx like this
cmake -S . -B build -D XDIAG_DISTRIBUTED=On -D CMAKE_CXX_COMPILER=mpicxx\nParallel compilation To speed up the compilation process, the build step can be performed in parallel using the -j flag
cmake --build build -j\nListing compile options
The available compilation options can be displayed using
cmake -L .\nShared / Static libraries
XDiag can be built either as a static or shared library. By default, a static library is built. To build a shared library, use the option XDIAG_SHARED_LIBS, e.g.
cmake -S . -B build -D XDIAG_SHARED_LIBS=On\nChoosing a certain compiler
The compiler (e.g. icpx) can be specified using
cmake -S . -B build -D CMAKE_CXX_COMPILER=icpx\nWarning
If the xdiag library is compiled with a certain compiler, it is advisable to also compile the application codes with the same compiler.
Setting the install path
In the installation step, the install directory can be set in the following way
cmake --install build --prefix /my/install/prefix\nDisabling HDF5/OpenMP
To disable support for HDF5 or OpenMP support, use
cmake -S . -B build -D XDIAG_DISABLE_OPENMP=On -D XDIAG_DISABLE_HDF5=On\nBuilding and running tests
To compile and run the testing programs, use
cmake -S . -B build -D BUILD_TESTING=On\ncmake --build build\nbuild/tests/tests\nmissing documentation
"},{"location":"documentation/compilation/building_documentation/","title":"Building Documentation","text":"The source files for the documentation can be found in the directory docs. The documentation is built using Material for MKDocs. To install it please run the following commands,
pip install mkdocs\npip install mkdocs-material\nTo work on the documentation locally, it can be served using
mkdocs serve\nfrom the xdiag root source directory. A local build of the documentation can then be accessed in a webbrowser at the adress
127.0.0.1:8000\nIn order to develop and extend the julia wrapper, one should work locally and build a local version of the xdiag Julia binaries. First, get the path to the CxxWrap package of julia. To do so, enter the Julia REPL,
julia\nusing CxxWrap\nCxxWrap.prefix_path()\n/path/to/libcxxwrap-julia-prefix. This is then used to configure the cmake compilation. cmake -S . -B build -D XDIAG_JULIA_WRAPPER=On -D CMAKE_PREFIX_PATH=/path/to/libcxxwrap-julia-prefix\ncmake --build build\ncmake --install build\nlibxdiagjl.so, (or the corresponding library format on non-Linux systems)."},{"location":"documentation/io/file_h5/","title":"FileH5","text":"A file handler for hdf5 files. The proper tool to write results of XDiag simulations to disk. Only provided for the C++ version as Julia already provides good tools handling hdf5 files.
Sources file_h5.hpp, file_h5.cpp
"},{"location":"documentation/io/file_h5/#constructors","title":"Constructors","text":"C++FileH5(std::string filename, std::string iomode = \"r\");\nw There are four possible values of iomode
r: read-only modew: secure write mode, new file is created if it does not existw!: forced write mode, existing file will be overwrittena: append mode, new datasets in an existing file can be createdReturns a handler to a value to be read or written from or to the hdf5 file.
C++hdf5::FileH5Handler operator[](std::string key);\nstd::string filename = XDIAG_DIRECTORY \"/misc/data/hdf5/write.h5\";\nauto fl = FileH5(filename, \"w!\");\n\n// Write output to the hdf5 file\nfl[\"val\"] = 12;\nfl[\"test/to\"] = 22;\nfl[\"test/to2/group\"] = 32;\nfl[\"test/to3/group2/asdf\"] = 42;\n\nauto mat = arma::cx_mat(3, 5, arma::fill::randn);\nfl[\"a/b/c/mat\"] = mat;\nA file handler for TOML files. Most funtionality is only provided for the C++ version as Julia already provides good tools handling TOML files.
Sources file_toml.hpp, file_toml.cpp
"},{"location":"documentation/io/file_toml/#constructors","title":"Constructors","text":"C++JuliaFileToml(const char *filename);\nFileToml(std::string filename);\nFileToml(filename::String)\nReturns whether or not the TOML file has a certain key defined.
C++bool defined(FileToml const &fl, std::string key);\nReturns a handler to a value to be read from the TOML file.
C++io::FileTomlHandler operator[](std::string key);\nauto fl = FileToml(XDIAG_DIRECTORY \"/misc/data/toml/input.toml\");\nXDIAG_SHOW(defined(fl, \"N\"));\n\nint N = fl[\"N\"].as<int>();\nint nup = fl[\"nup\"].as<int>();\ndouble J1 = fl[\"J1\"].as<double>();\ndouble J2 = fl[\"J2\"].as<double>();\n\nauto block = Spinhalf(N, nup);\nauto H = OpSum();\nfor (int i=0; i<N; ++i){\n H += J1 * Op(\"SdotS\", {i, (i+1)%N});\n H += J2 * Op(\"SdotS\", {i, (i+2)%N});\n}\ndouble e0 = eigval0(H, block);\nXDIAG_SHOW(e0);\nReads an OpSum object from a TOML file.
Sources read.hpp, read.cpp
"},{"location":"documentation/io/read_opsum/#definition","title":"Definition","text":"C++JuliaOpSum read_opsum(FileToml file, std::string tag);\nread_opsum(file::FileToml, tag::String);\nAn OpSum can be defined in a TOML file as a simple list. The entries of the list are themselves also lists, which contain two or more entries:
Here is typical example of a OpSum specification in a TOML file:
Interactions = [\n ['J1', 'SdotS', 0, 1],\n ['J1', 'SdotS', 1, 2],\n ['J1', 'SdotS', 2, 0],\n ['Jchi', 'ScalarChirality', 0, 1, 2],\n]\nstd::string file = XDIAG_DIRECTORY \"/misc/data/triangular.9.hop.sublattices.tsl.toml\";\nauto fl = FileToml(file);\nauto ops = read_opsum(fl, \"Interactions\");\nXDIAG_SHOW(ops);\n\nReads an PermutationGroup object from a TOML file.
Sources read.hpp, read.cpp
"},{"location":"documentation/io/read_permutation_group/#definition","title":"Definition","text":"C++JuliaOpSum read_permutation_group(FileToml file, std::string tag);\nread_permutation_group(file::FileToml, tag::String);\nA PermutationGroup can be defined in a TOML file as an integer matrix, where the rows are the integers of the permutation.
A typical specification of a \\(C_4\\) PermutationGroup is shown here:
Symmetries = [\n [0, 1, 2, 3],\n [1, 2, 3, 0],\n [2, 3, 0, 1],\n [3, 0, 1, 2]\n]\n1-indexing in Julia / 0-indexing in C++
To enumerate the sites of a Permutation, we start counting at 1 in Julia and 0 in C++.
"},{"location":"documentation/io/read_permutation_group/#usage-example","title":"Usage Example","text":"C++Juliastd::string file = XDIAG_DIRECTORY \"/misc/data/triangular.9.hop.sublattices.tsl.toml\";\nauto fl = FileToml(file);\nauto group = read_permutation_group(fl, \"Symmetries\");\nXDIAG_SHOW(group);\n\nReads an Representation object from a TOML file.
Sources read.hpp, read.cpp
"},{"location":"documentation/io/read_representation/#definition","title":"Definition","text":"C++JuliaOpSum read_representation(FileToml file, std::string irrep_tag, \n std::string group_tag = \"Symmetries\");\nread_representation(file::FileToml, irrep_tag::String, \n group_tag::String = \"Symmetries\");\nA Representation can be defined in a TOML file by specifying up two or three things:
characters of the representation allowed_symmetries of the representation, i.e. a list of the number of symmetries used in the irrep. By default all symmetries of the group are used in the representation.A typical specification of several Representations is shown here:
Symmetries = [\n [0, 1, 2, 3],\n [1, 2, 3, 0],\n [2, 3, 0, 1],\n [3, 0, 1, 2]\n ]\n\n # real irrep, momentum 0\n [k_0]\n characters=[1.0, 1.0, 1.0, 1.0]\n\n # complex irrep, momentum pi/2\n [k_pi2]\n characters=[\n [1.0, 0.0], \n [0.0, 1.0], \n [-1.0, 0.0],\n [0.0, -1.0]\n ]\n\n # real irrep, momentum pi\n [k_pi]\n characters=[1.0, -1.0, 1.0, -1.0]\n\n # real irrep on subgroup\n [k_pi2_half]\n characters=[1.0, -1.0]\n allowed_symmetries=[0, 2]\n1-indexing in Julia / 0-indexing in C++
To enumerate the sites of a Permutation, we start counting at 1 in Julia and 0 in C++.
"},{"location":"documentation/io/read_representation/#usage-example","title":"Usage Example","text":"The example reads the representation defined in the irreps.toml file, whose contents are shown up in the section Data format.
std::string file = XDIAG_DIRECTORY \"/misc/data/irreps.toml\";\nauto fl = FileToml(file);\n\nauto k_0 = read_representation(fl, \"k_0\");\nXDIAG_SHOW(k_0);\nXDIAG_SHOW(isreal(k_0));\n\nauto k_pi2 = read_representation(fl, \"k_pi2\");\nXDIAG_SHOW(k_pi2);\nXDIAG_SHOW(isreal(k_pi2));\n\nauto k_pi = read_representation(fl, \"k_pi\");\nXDIAG_SHOW(k_pi);\nXDIAG_SHOW(isreal(k_pi));\n\nauto k_pi2_half = read_representation(fl, \"k_pi2_half\");\nXDIAG_SHOW(k_pi2_half);\nXDIAG_SHOW(isreal(k_pi2_half));\n\nDescribes the coupling of an operator. A coupling can either be a string or a Scalar, which is either a real or complex double precision floating point number.
Sources coupling.hpp, coupling.cpp
"},{"location":"documentation/operators/coupling/#constructors","title":"Constructors","text":"C++Coupling(std::string value);\nCoupling(double value);\nCoupling(complex value);\nCoupling(Scalar value););\nReturns whether or not the coupling is a Scalar, i.e. a real or complex number.
C++bool isscalar(Coupling const& c)\nReturns whether or not the coupling is a string.
C++bool isstring(Coupling const& c)\nReturns the Scalar if the Coupling holds a scalar. Otherwise throws an xdiag::Error
Scalar scalar(Coupling const& c)\nReturns the string if the Coupling holds a string. Otherwise throws an xdiag::Error
std::string string(Coupling const& c)\nReturns the hermitian conjugate \\(\\mathcal{O}^\\dagger\\) of an operator \\(\\mathcal{O}\\) represented by an Op or OpSum object. Please note the details when conjugating complex couplings, outlined in OpSum # Complex couplings.
Sources hc.hpp, hc.cpp
"},{"location":"documentation/operators/hc/#definition","title":"Definition","text":"C++Op hc(Op const &op)\nOpSum hc(OpSum const &ops)\nauto cdagup = Op(\"Cdagup\", 0);\nauto sdots = Op(\"SdotS\", {0, 1});\nauto hop = (1.0 + 1.0i) * Op(\"Hop\", {0, 1});\nXDIAG_SHOW(cdagup == hc(cdagup));\nXDIAG_SHOW(sdots == hc(sdots));\nXDIAG_SHOW(hop == hc(hop));\nObject describing a single linear operator acting on a Hilbert space.
Sources op.hpp, op.cpp
Every operator is defined by up to three paramaters:
The type of the operator. This is a string argument which determines what kind of operator is represented. A detailed overview of the available types can be found at Operator types
The sites of the operator. This defines which physical sites (or orbitals) the operator acts upon. While most operator types require sites, there are also operator types (e.g. HubbardU) which do not need to define sites.
For special interactions, it can be necessary to additionally specify a numerical matrix, which can either be real or complex. An example is the operator type Matrix defining generic spin interactions.
Op(std::string type);\nOp(std::string type, int64_t site);\nOp(std::string type, std::vector<int64_t> const &sites);\nOp(std::string type, int64_t site, arma::mat const &matrix);\nOp(std::string type, std::vector<int64_t> const &sites, arma::mat const &matrix);\nOp(std::string type, int64_t site, arma::cx_mat const &matrix);\nOp(std::string type, std::vector<int64_t> const &sites, arma::cx_mat const &matrix);\nOp(type::String)\nOp(type::String, site::Int64)\nOp(type::String, sites::Vector{Int64})\nOp(type::String, site::Int64, matrix::Matrix{Float64})\nOp(type::String, sites::Vector{Int64}, matrix::Matrix{Float64})\nOp(type::String, site::Int64, matrix::Matrix{ComplexF64})\nOp(type::String, sites::Vector{Int64}, matrix::Matrix{ComplexF64})\n1-indexing in Julia / 0-indexing in C++
To enumerate the sites of an Op, we start counting at 1 in Julia and 0 in C++.
"},{"location":"documentation/operators/op/#methods","title":"Methods","text":""},{"location":"documentation/operators/op/#isreal","title":"isreal","text":"Returns whether an Op is a real operator.
C++Juliabool isreal(Op const &op);\nisreal(op::Op)\nauto op = \"T\" * Op(\"Hop\", {0, 1});\nXDIAG_SHOW(op);\n\nop = 1.23 * Op(\"Hop\", {0, 1});\nXDIAG_SHOW(op);\n\narma::cx_mat m(arma::mat(\"0 0; 0 0\"), arma::mat(\"0 -1; 1 0\"));\n op = Op(\"Matrix\", 0, m);\nXDIAG_SHOW(op);\nXDIAG_SHOW(isreal(op));\nop = \"T\" * Op(\"Hop\", [1, 2])\n@show op\n\nop = 1.23 * Op(\"Hop\", [1, 2])\n@show op\n\nop = Op(\"Matrix\", 1, [0 -im; im 0])\n@show op\n@show isreal(op)\nGeneric operators in XDiag are represented as OpSum objects made up of Coupling and Op objects. Every Op is defined by a type. Here we list all the available types implemented in XDiag, their required number of sites, and the blocks for which they are available.
Hop A hopping term for \\(\\uparrow\\) and \\(\\downarrow\\) spins of the form $$ \\sum_{\\sigma=\\uparrow\\downarrow} (c^\\dagger_{i\\sigma}c_{j\\sigma} + \\textrm{h.c.})$$ 2 tJ, Electron, tJDistributed Hopup A hopping term for \\(\\uparrow\\) spins of the form $$ (c^\\dagger_{i\\uparrow}c_{j\\uparrow} + \\textrm{h.c.})$$ 2 tJ, Electron, tJDistributed Hopdn A hopping term for \\(\\downarrow\\) spins of the form $$ (c^\\dagger_{i\\downarrow}c_{j\\downarrow} + \\textrm{h.c.})$$ 2 tJ, Electron, tJDistributed HubbardU A uniform Hubbard interaction across the full lattice of the form $$ \\sum_i n_{i\\uparrow}n_{i\\downarrow}$$ 0 Electron Cdagup A fermionic creation operator for an \\(\\uparrow\\) spin \\(c^\\dagger_{i\\uparrow}\\) 1 tJ, Electron, tJDistributed Cdagdn A fermionic creation operator for an \\(\\downarrow\\) spin \\(c^\\dagger_{i\\downarrow}\\) 1 tJ, Electron, tJDistributed Cup A fermionic annihilation operator for an \\(\\uparrow\\) spin \\(c_{i\\uparrow}\\) 1 tJ, Electron, tJDistributed Cdn A fermionic annihilation operator for an \\(\\downarrow\\) spin \\(c_{i\\downarrow}\\) 1 tJ, Electron, tJDistributed Nup A number operator for an \\(\\uparrow\\) spin \\(n_{i\\uparrow}\\) 1 tJ, Electron, tJDistributed Ndn A number operator for an \\(\\downarrow\\) spin \\(n_{i\\downarrow}\\) 1 tJ, Electron, tJDistributed Ntot A number operator \\(n_i = n_{i\\uparrow} + n_{i\\downarrow}\\) 1 tJ, Electron, tJDistributed NtotNtot A density-density interaction \\(n_i n_j\\) 2 tJ, Electron, tJDistributed SdotS A Heisenberg interaction of the form $$ \\mathbf{S}_i \\cdot \\mathbf{S}_j = S^x_iS^x_j + S^y_iS^y_j + S^z_iS^z_j$$ 2 Spinhalf, tJ, Electron, SpinhalfDistributed, tJDistributed SzSz An Ising interaction of the form $ S^z_i S^z_j $ 2 Spinhalf, tJ, Electron, SpinhalfDistributed, tJDistributed Exchange A spin exchange interaction of the form $$ \\frac{1}{2}(S^+_i S^-_j + S^-_iS^+_j) = S^x_iS^x_j + S^y_iS^y_j j$$ 2 Spinhalf, tJ, Electron, SpinhalfDistributed, tJDistributed Sz A local magnetic moment in the \\(z\\)-direction $ S^z_i$ 1 Spinhalf, tJ, Electron, SpinhalfDistributed, tJDistributed S+ A local spin raising operator \\(S^+_i\\) 1 Spinhalf, SpinhalfDistributed S- A local spin lowering operator \\(S^-_i\\) 1 Spinhalf, SpinhalfDistributed ScalarChirality A scalar chirality interaction of the form $$ \\mathbf{S}_i \\cdot ( \\mathbf{S}_j \\times \\mathbf{S}_k)$$ 3 Spinhalf tJSzSz An Ising interaction as encountered in the \\(t-J\\) model of the form $$ S^z_i S^z_j - \\frac{n_i n_j}{4}$$ 2 tJ, tJDistributed tJSdotS An Heisenberg interaction as encountered in the \\(t-J\\) model of the form $$ \\mathbf{S}_i \\cdot \\mathbf{S}_j - \\frac{n_i n_j}{4}$$ 2 tJ, tJDistributed Matrix A generic spin interaction no an arbitrary number of sites defined via a coupling matrix arbitrary Spinhalf"},{"location":"documentation/operators/operator_types/#matrix-type","title":"Matrix type","text":"The Matrix interaction type is a special type with whom one can define generic interactions for the Spinhalf block. In addition to the type and sites argument, also a numerical matrix is provided when constructing the Op object. The matrix describes the operator acting on the \\(2^n\\) dimensional space spanned by the \\(n\\) sites of the operator. For example, we can represent a \\(S^x\\) spin operator as,
auto sx = arma::mat({{0, 1},{1, 0}});\nauto op = Op(\"Matrix\", 0, sx);\nMore generically, we can use this mechanism to construct arbitary spin interactions, e.g.
C++auto sx = arma::mat({{0, 1},{1, 0}});\nauto sz = arma::mat({{0.5, 1},{0, -0.5}});\n\narma::mat sxsz = arma::kron(sx, sz);\narma::mat sxszsxsz = arma::kron(sxsz, sxsz);\n\nauto op_sxsz = Op(\"Matrix\", {0, 1}, sxsz);\nauto op_sxszsxsz = Op(\"Matrix\", {0, 1, 2, 3}, sxsz);\nHere we have been using the Kronecker product function kron.
Object representing a generic many-body operator by a sum of operators of the form
\\[ \\mathcal{O} = \\sum_i c_i \\mathcal{O}_i. \\]Sources opsum.hpp, opsum.cpp
An OpSum is consists of a sum of pairs given by
A coupling constant \\(c_i\\) which is defined by a Coupling object. The coupling can either be a string name or a real/complex number.
An operator \\(\\mathcal{O}_i\\) defined by an Op object.
Generically, an OpSum can thus have coupling constants defined by either strings or numerical real/complex numbers. We call an OpSum plain if its couplings are only numerical numbers, and not strings. String couplings can be defined by using the access operator[]. If all string coupling constants are defined, the OpSum can be converted to a plain OpSum using the plain method shown below.
Thus, OpSums can be defined independently of the numerical values of the coupling constants, e.g. in an input file. Upon execution of the code, these constants can then be set. Most operations in XDiag require the OpSum to be convertible to a plain OpSum.
OpSums can be added and subtracted, as well as multiplied with and divided by a Scalar. Hence, OpSums carry the mathematical structure of a vector space.
"},{"location":"documentation/operators/opsum/#constructors","title":"Constructors","text":"The following constructors create an OpSum with a single pair of coupling and operator. Additional terms can be added using the + and += operators explained further below. If no Coupling is given, a numerical coefficient of 1.0 is assumed.
OpSum(Op const &op);\nOpSum(Coupling const &cpl, Op const &op);\nOpSum(op::Op)\nOpSum(cpl::Coupling, op::Op)\nAlternatively, an OpSum can also be constructed via the * operator, for example:
auto ops = OpSum();\nfor (int i = 0; i<N; ++i) {\n ops += J * Op(\"SzSz\", {i, (i + 1) % N});\n}\nXDiag allows all couplings to be complex. Depending on the operator type a complex coupling can have two meanings:
A complex prefactor \\(c\\) which upon hermitian conjugation with hc gets conjugated to \\(c^\\star\\). This is the case for the following interaction types: HubbardU, Cdagup, Cdagdn, Cup, Cdn, Nup, Ndn, Ntot, NtotNtot, SdotS, SzSz, Sz, S+, S-, ScalarChirality, tJSzSz, tJSdotS, Matrix Thus, a complex coupling can turn a Hermitian operator to a non-Hermitian operator.
The coupling is part of the definition of the operator. For, example a hopping operator of the form $$ ( t c^\\dagger_{i\\sigma}c_{j\\sigma} + \\textrm{h.c.}) = ( t c^\\dagger_{i\\sigma}c_{j\\sigma} + t^\\star c^\\dagger_{j\\sigma}c_{i\\sigma}) $$ A complex coupling \\(t\\) gives the hopping a phase, but the overall operator remains Hermitian and, thus, invariant under hc. This holds for the types Hop, Hopup, Hopdn, Exchange. In the latter case, complex spin exchange Exchange is defined as, $$ \\frac{1}{2}( J S^+_i S^-_j + J^\\star S^-_iS^+_j)$$
Converts an OpSum with possible string couplings to an OpSum with purely numerical real/complex couplings.
C++JuliaOpSum plain(OpSum const &ops) const;\nplain(ops::OpSum)\nCreates an OpSum with a single pair of Coupling and Op objects.
C++OpSum operator*(std::string cpl, Op const &op);\nOpSum operator*(double cpl, Op const &op);\nOpSum operator*(complex cpl, Op const &op);\nOpSum operator*(Scalar cpl, Op const &op);\nOpSum operator*(Coupling const& cpl, Op const &op);\nAdds two OpSum objects \\(\\mathcal{A} = \\sum_i a_i \\mathcal{A}_i\\) and \\(\\mathcal{B} = \\sum_i b_i \\mathcal{B}_i\\) to for the sum of the two operators, $$ \\mathcal{A} + \\mathcal{B} = \\sum_i a_i \\mathcal{A}_i + \\sum_i b_i \\mathcal{B}_i$$
C++JuliaOpSum &operator+=(OpSum const &ops);\nOpSum operator+(OpSum const &ops) const;\nBase.:+(ops::OpSum, ops2::OpSum)\nSubtracts to OpSum objects.
C++JuliaOpSum &operator-=(OpSum const &ops);\nOpSum operator-(OpSum const &ops) const;\nBase.:-(ops::OpSum, ops2::OpSum)\nMultiplies an OpSum \\(\\mathcal{A} = \\sum_i a_i \\mathcal{A}_i\\) with a scalar \\(b\\) to form
\\[\\mathcal{B} = b \\sum_i a_i \\mathcal{A}_i\\] C++OpSum &operator*=(Scalar const &cpl);\nOpSum &operator/=(Scalar const &cpl);\nOpSum operator*(Scalar const &cpl, OpSum const &op);\nOpSum operator*(OpSum const &op, Scalar const &cpl);\nOpSum operator/(OpSum const &op, Scalar const &cpl);\nSets a coupling constant defined as a string to a numerical value.
C++Scalar &operator[](std::string name);\nScalar const &operator[](std::string name) const;\nReturns a vector of strings with the coupling constants defined, i.e. the strings that define some of the Coupling objects.
C++std::vector<std::string> constants(OpSum const &ops) const;\n// Define the 1D transverse-field Ising chain\nint N = 12;\ndouble J = 1.0;\ndouble h = 0.5;\nauto Sx = arma::mat(\"0 1; 1 0\");\n\n// Option 1: coupling constants as numbers\nauto ops1 = OpSum();\nfor (int i = 0; i<N; ++i) {\n ops1 += J * Op(\"SzSz\", {i, (i+1)%N});\n ops1 += h * Op(\"Matrix\", i, Sx);\n}\n\n// Option 2: coupling constants as strings\nauto ops2 = OpSum();\nfor (int i = 0; i<N; ++i) {\n ops2 += \"J\" * Op(\"SzSz\", {i, (i+1)%N});\n ops2 += \"h\" * Op(\"Matrix\", i, Sx);\n}\nops2[\"J\"] = J;\nops2[\"h\"] = h;\n\nXDIAG_SHOW(isapprox(ops1, ops2));\nXDIAG_SHOW(isapprox(ops1 + ops2, 2.0 * ops1));\nA variable which can be either a real or complex double precision number. The algebraic operators +,-,*,\\ and comparison operators ==, != are defined.
Sources scalar.hpp, scalar.cpp
"},{"location":"documentation/operators/scalar/#constructors","title":"Constructors","text":"C++Scalar(double value);\nScalar(complex value);\nReturns true if Scalar is a real number, false otherwise.
bool isreal(Scalar const &s)\nReturns the real part of the Scalar.
C++double real(Scalar const &s)\nReturns the imaginary part of the Scalar.
C++double imag(Scalar const &s)\nReturns the complex conjugate of the Scalar.
C++Scalar conj(Scalar const &s)\nReturns the absolute value of the Scalar.
C++double abs(Scalar const &s)\nCompares two Scalars to whether they are approximately equal. Return true if the following condition holds.
\\[ | a - b | < atol + rtol*|b| \\] C++bool isapprox(Scalar const &a, Scalar const &b, double rtol = 1e-12, double atol = 1e-12)\nParameters
Name Description Default a first Scalar b second Scalar rtol relative tolerance 1e-12 atol absolute tolerance 0"},{"location":"documentation/operators/scalar/#to_string","title":"to_string","text":"Converts the scalar into a string.
C++std::string to_string(Scalar const &v);\nSymmetrizes an operator with respect to a PermutationGroup or a Representation.
Sources symmetrize.hpp, symmetrize.cpp
Symmetrization in this context means the following. In general, we are given an OpSum of the form,
\\[ O = \\sum_{A\\subseteq \\mathcal{L}} O_A,\\]where \\(O_A\\) denotes a local operator acting on sites \\(A=\\{a_1, \\ldots, a_{l_A}\\}\\) and \\(L\\) denotes the lattice. A PermutationGroup \\(\\mathcal{G}\\) is defined through its permutations \\(\\pi_1, \\ldots, \\pi_M\\). The symmetrized operator returned by this function is then
\\[ O^\\mathcal{G} = \\frac{1}{M}\\sum_{A\\subseteq \\mathcal{L}} \\sum_{\\pi \\in \\mathcal{G}} O_{\\pi(A)},\\]where \\(\\pi(A) = \\{\\pi(a_1), \\ldots,\\pi(a_{l_A})\\}\\) denotes the permutated set of sites of the local operator \\(O_A\\). If a Representation called \\(\\rho\\) is given in addition, the following operator is constructed,
\\[ O^\\mathcal{G, \\rho} = \\frac{1}{M}\\sum_{A\\subseteq \\mathcal{L}} \\sum_{\\pi \\in \\mathcal{G}} \\chi_\\rho(\\pi) O_{\\pi(A)},\\]where \\(\\chi_\\rho(\\pi)\\) denotes the characters of the representation \\(\\rho\\). This routine is useful to evaluate observables in symmetrized blocks.
"},{"location":"documentation/operators/symmetrize/#definition","title":"Definition","text":"C++JuliaOpSum symmetrize(Op const &op, PermutationGroup const &group);\nOpSum symmetrize(Op const &op, Representation const &irrep);\nOpSum symmetrize(OpSum const &ops, PermutationGroup const &group);\nOpSum symmetrize(OpSum const &ops, Representation const &irrep);\nsymmetrize(op::Op, group::PermutationGroup)\nsymmetrize(op::Op, irrep::Representation)\nsymmetrize(ops::OpSum, group::PermutationGroup)\nsymmetrize(ops::OpSum, irrep::Representation)\nint N = 4;\nint nup = 2;\nauto block = Spinhalf(N, nup);\nauto p1 = Permutation({0, 1, 2, 3});\nauto p2 = Permutation({1, 2, 3, 0});\nauto p3 = Permutation({2, 3, 0, 1});\nauto p4 = Permutation({3, 0, 1, 2});\nauto group = PermutationGroup({p1, p2, p3, p4});\nauto rep = Representation(group);\nauto block_sym = Spinhalf(N, rep);\n\nauto ops = OpSum();\nfor (int i=0; i<N; ++i) {\n ops += Op(\"SdotS\", {i, (i+1)%N});\n}\nauto [e0, psi] = eig0(ops, block);\nauto [e0s, psi_sym] = eig0(ops, block_sym);\n\nauto corr = Op(\"SdotS\", {0, 1});\nauto nn_corr = inner(corr, psi);\nauto corr_sym = symmetrize(corr, group);\nauto nn_corr_sym = innerC(corr_sym, psi_sym);\nXDIAG_SHOW(nn_corr);\nXDIAG_SHOW(nn_corr_sym);\nlet\n N = 4\n nup = 2\n block = Spinhalf(N, nup)\n p1 = Permutation([1, 2, 3, 4])\n p2 = Permutation([2, 3, 4, 1])\n p3 = Permutation([3, 4, 1, 2])\n p4 = Permutation([4, 1, 2, 3])\n group = PermutationGroup([p1, p2, p3, p4])\n rep = Representation([1, 1, 1, 1])\n block_sym = Spinhalf(N, group, rep)\n\n ops = OpSum()\n for i in 1:N\n ops += Op(\"HB\", 1.0, [i, mod1(i+1, N)])\n end\n\n e0, psi = eig0(ops, block);\n e0, psi_sym = eig0(ops, block_sym);\n\n corr = Op(\"HB\", 1.0, [1, 2])\n nn_corr = inner(corr, psi)\n corr_sym = symmetrize(corr, group)\n nn_corr_sym = inner(corr_sym, psi_sym)\n @show nn_corr, nn_corr_sym\nend\nSources create_state.hpp, create_state.cpp
"},{"location":"documentation/states/create_state/#product_state","title":"product_state","text":"Creates a filled product state.
C++JuliaState product_state(Block const &block, std::vector<std::string> const &local_state, bool real = true);\nproduct_state(block::Block, local_states::Vector{String}, real::Bool=true)\nCreate a filled random state with normal \\(\\mathcal{N}(0, 1)\\) distributed coefficients.
C++JuliaState random_state(Block const &block, bool real = true, int64_t seed = 42, bool normalized = true);\nrandom_state(block::Block, real::Bool=true, seed::Int64=42, normalized::Bool=true\nCreate a filled state with all zero entries.
C++JuliaState zero_state(Block const &block, bool real = true, int64_t n_cols = 1);\nzero_state(block::Block, real::Bool=true, n_col::Int64=1)\nSet all coefficients of a given state to zero.
C++Juliavoid zero(State &state);\nzero(state::State)\nauto block = Spinhalf(2);\nauto state = product_state(block, {\"Up\", \"Dn\"});\nXDIAG_SHOW(state.vector());\n\nzero(state);\nXDIAG_SHOW(state.vector());\n\nstate = random_state(block, false, 1234, true);\nXDIAG_SHOW(state.vectorC());\n\nstate = zero_state(block, true, 2);\nXDIAG_SHOW(state.vector());\nblock = Spinhalf(2)\nstate = product(block, [\"Up\", \"Dn\"])\ndisplay(vector(state))\n\nzero(state)\ndisplay(vector(state))\n\nstate = rand(block, false, 1234, true)\ndisplay(vector(state))\n\nstate = zeros(block, true, 2)\ndisplay(matrix(state))\nFills a State with a given model state, e.g. a ProductState or a RandomState.
Sources fill.hpp, fill.cpp
"},{"location":"documentation/states/fill/#definition","title":"Definition","text":"C++Juliavoid fill(State &state, ProductState const &pstate, int64_t ncol = 0);\nvoid fill(State &state, RandomState const &rstate, int64_t ncol = 0);\nfill(state::State, pstate::ProductState, ncol::Int64 = 1)\nfill(state::State, rstate::RandomState, ncol::Int64 = 1)\nauto block = Spinhalf(2);\nauto state = State(block); \nauto pstate = ProductState({\"Up\", \"Dn\"});\nfill(state, pstate);\nXDIAG_SHOW(state.vector());\n\nauto rstate = RandomState(1234);\nfill(state, rstate);\nXDIAG_SHOW(state.vector());\nblock = Spinhalf(2)\nstate = State(block)\npstate = ProductState([\"Up\", \"Dn\"])\nfill(state, pstate)\ndisplay(vector(state))\n\nrstate = RandomState(1234)\nfill(state, rstate)\ndisplay(vector(state))\nA product state of local configurations
Sources product_state.hpp, product_state.cpp
"},{"location":"documentation/states/product_state/#constructors","title":"Constructors","text":"C++JuliaProductState(int64_t nsites);\nProductState(std::vector<std::string> const &local_states);\nProductState(nsites::Int64)\nProductState(local_states::Vector{String})\nA ProductState can be iterated over, where at each iteration the string of the local configuration is retured. Here is an example:
C++Juliaauto pstate = ProductState({\"Up\", \"Dn\", \"Emp\", \"UpDn\"});\nfor (auto s : pstate) {\n Log(\"{}\", s);\n}\npstate = ProductState([\"Up\", \"Dn\", \"Emp\", \"UpDn\"])\nfor s in pstate\n @show s\nend\nReturns the number of sites of the product state
C++Juliaint64_t nsites(ProductState const &p);\nnsites(p::ProductState)\nReturns the number of sites of the product state. Same as \"nsites\".
C++Juliaint64_t size(ProductState const &p);\nsize(state::ProductState)\nSets the local configuration at the given site index to the given string.
C++Juliastd::string &operator[](int64_t i);\nsetindex!(state::ProductState, local_state::String, idx::Int64)\nReturns the string of the local configuration at the given site index.
C++Juliastd::string const &operator[](int64_t i) const;\ngetindex(state::ProductState, idx::Int64)\nAdds a local configuration add the end of the product state.
C++Juliavoid push_back(std::string l);\npush!(state::ProductState, local_state::String\nauto pstate = ProductState({\"Up\", \"Dn\", \"Emp\", \"UpDn\"});\nfor (auto s : pstate) {\n Log(\"{}\", s);\n}\nXDIAG_SHOW(to_string(pstate));\n\npstate = ProductState();\npstate.push_back(\"Dn\");\npstate.push_back(\"Up\");\npstate.push_back(\"Dn\");\nXDIAG_SHOW(pstate.nsites());\nfor (auto s : pstate) {\n Log(\"{}\", s);\n}\nXDIAG_SHOW(to_string(pstate));\npstate = ProductState([\"Up\", \"Dn\", \"Emp\", \"UpDn\"])\nfor s in pstate\n @show s\nend\n@show pstate\n\npstate = ProductState()\npush!(pstate, \"Dn\")\npush!(pstate, \"Up\")\npush!(pstate, \"Dn\")\n@show nsites(pstate)\nfor s in pstate\n @show s\nend\n@show pstate\nA random state with \\(\\mathcal{N}(0, 1)\\) normal distributed coefficients.
Sources random_state.hpp, random_state.cpp
"},{"location":"documentation/states/random_state/#constructors","title":"Constructors","text":"C++JuliaRandomState(int64_t seed = 42, bool normalized = true);\nRandomState(seed::Int64 = 42, normalized::Bool = true)\nblock = Spinhalf(2)\nstate = State(block, real=false) # complex State\nrstate1 = RandomState(1234)\nfill(state, rstate1)\ndisplay(vector(state))\n\nrstate2 = RandomState(4321)\nfill(state, rstate2)\ndisplay(vector(state))\n\nfill(state, rstate1)\ndisplay(vector(state))\nauto block = Spinhalf(2);\nauto state = State(block, false); // complex State\nauto rstate1 = RandomState(1234);\nfill(state, rstate1);\nXDIAG_SHOW(state.vectorC());\n\nauto rstate2 = RandomState(4321);\nfill(state, rstate2);\nXDIAG_SHOW(state.vectorC());\n\nfill(state, rstate1);\nXDIAG_SHOW(state.vectorC());\nA generic state describing a quantum wave function \\(|\\psi \\rangle\\).
Sources state.hpp state.cpp
"},{"location":"documentation/states/state/#constructors","title":"Constructors","text":"A state can be constructed in two ways:
By only specifying the block. In this case the state is initialized with all coefficients zero.
C++JuliaState(Block const &block, bool real = true, int64_t n_cols = 1);\nState(block::Block; real::Bool = true, n_cols::Int64 = 1)\nBy handing a vector or matrix of coefficients.
C++Juliatemplate <typename block_t, typename coeff_t>\nState(block_t const &block, arma::Col<coeff_t> const &vector);\n\ntemplate <typename block_t, typename coeff_t>\nState(block_t const &block, arma::Mat<coeff_t> const &matrix);\nState(block::Block, vec::Vector{Float64})\nState(block::Block, vec::Vector{ComplexF64})\nState(block::Block, mat::Matrix{Float64})\nState(block::Block, mat::Matrix{ComplexF64})\nReturns the number of sites of the block the state is defined on.
C++Juliaint64_t nsites(State const &s) const\nnsites(state::State)\nReturns whether two states are approximately equal.
C++Juliabool isapprox(State const &v, State const &w, double rtol = 1e-12,\n double atol = 1e-12);\nisapprox(v::State, w::State, rtol::Float64, atol::Float64)\nReturns whether the state is real.
C++Juliaint64_t isreal(State const &s) const;\nisreal(state::State)\nReturns whether the real part of the State.
C++JuliaState real(State const &s) const;\nreal(state::State)\nReturns whether the imaginary part of the State.
C++JuliaState imag(State const &s) const;\nimag(state::State)\nTurns a real State into a complex State. Does nothing if the state is already complex
C++Juliavoid make_complex(State &s);\nmake_complex!(state::State)\nReturns the dimension of the block the state is defined on.
C++Juliaint64_t dim(State const &s) const;\ndim(block::Spinhalf)\nReturns the size of the block the state is defined on. locally. Same as \"dim\" for non-distributed Blocks but different for distributed blocks.
C++Juliaint64_t size(State const &s);\nsize(s::State)\nReturns number of rows of the local storage. Same as \"size\"
C++Juliaint64_t nrows(State const &s);\nnrows(s::State)\nReturns number of columns of the local storage.
C++Juliaint64_t ncols(State const &s);\nncols(s::State)\nReturns a state created from the n-th column of the storage. Whether or not the storage is copied can be specified by setting the flag \"copy\".
C++JuliaState col(State const &s, int64_t n, bool copy = true);\ncol(s::State, n::Int64 = 1; copy::Bool = true)\nReturns a vector from the n-th column of the storage. In C++ use \"vector\"/\"vectorC\" to either get a real or complex vector.
C++Juliaarma::vec vector(State const &s, int64_t n = 0, bool copy = true);\narma::cx_vec vectorC(State const &s, int64_t n = 0, bool copy = true);\nvector(state::State; n::Int64 = 1, copy::Bool=true)\n# no vectorC method in julia\nReturns matrix representing the storage. In C++ use \"matrix\"/\"matrixC\" to either get a real or complex matrix.
C++Juliaarma::vec matrix(State const &s, bool copy = true);\narma::cx_vec matrixC(State const &s, bool copy = true);\nmatrix(state::State, copy::Bool=true)\n# no matrixC method in julia\nauto block = Spinhalf(2);\nauto psi1 = State(block, arma::vec(\"1.0 2.0 3.0 4.0\"));\nXDIAG_SHOW(psi1);\nXDIAG_SHOW(vector(psi1));\nmake_complex(psi1);\nXDIAG_SHOW(vectorC(psi1));\n\nauto psi2 = State(block, false, 3);\nXDIAG_SHOW(psi2);\nXDIAG_SHOW(matrixC(psi2));\n\nauto psi3 = State(block, arma::cx_vec(arma::vec(\"1.0 2.0 3.0 4.0\"),\n arma::vec(\"4.0 3.0 2.0 1.0\")));\nXDIAG_SHOW(vectorC(psi3));\nXDIAG_SHOW(vector(real(psi3)));\nXDIAG_SHOW(vector(imag(psi3)));\nblock = Spinhalf(2)\npsi1 = State(block, [1.0, 2.0, 3.0, 4.0])\n@show psi1\ndisplay(vector(psi1))\nmake_complex!(psi1)\ndisplay(vector(psi1))\n\npsi2 = State(block, real=false, n_cols=3)\n@show psi2\ndisplay(matrix(psi2))\n\npsi3 = State(block, [1.0+4.0im, 2.0+3.0im, 3.0+2.0im, 4.0+1.0im])\ndisplay(vector(psi3))\ndisplay(vector(real(psi3)))\ndisplay(vector(imag(psi3)))\nPermutations of indices or lattice sites. Basic building block of a PermutationGroup. Permutations can be multiplied, inverted and raised to a power.
Sources permutation.hpp, permutation.cpp
"},{"location":"documentation/symmetries/permutation/#constructors","title":"Constructors","text":""},{"location":"documentation/symmetries/permutation/#from-an-array","title":"From an array","text":"Creates an Permutation out of an array of integers, e.g. {0, 2, 1, 3}. If the input array is of size N then every number between 0 and N-1 must occur exactly once, otherwise the Permutation is invalid.
1-indexing in Julia / 0-indexing in C++
To enumerate the sites of a Permutation, we start counting at 1 in Julia and 0 in C++.
C++JuliaPermutation(std::initializer_list<int64_t> list);\nPermutation(std::vector<int32_t> const &array);\nPermutation(std::vector<int64_t> const &array);\nPermutation(arma::Col<int64_t> const &array);\nPermutation(int64_t *array, int64_t size);\nPermutation(array::Vector{Int64})\nConstructs an identity permutation of a given size, e.g. {0, 1, 2, 3}.
Permutation(int64_t size);\nPermutation(size::Integer)\nComputes the inverse permutation.
C++JuliaPermutation inverse(Permutation const &p);\ninverse(perm::Permutation)\nConcatenates two permutations by overloading the * operator.
Permutation operator*(Permutation const &p1, Permutation const &p2);\nBase.:*(p1::Permutation, p2::Permutation)\nRaises a permutation to an integer power.
C++JuliaPermutation pow(Permutation const &p, int64_t power);\nBase.:^(p::Permutation, power::Integer)\nReturns the size of a Permutation.
C++Juliaint64_t size(Permutation const &p);\nsize(p::Permutation)\nPermutation p1 = {0, 2, 1, 3};\nPermutation p2 = {2, 0, 1, 3};\n\nXDIAG_SHOW(inverse(p1));\nXDIAG_SHOW(p1*p2);\np1 = Permutation([1, 3, 2, 4])\np2 = Permutation([3, 1, 2, 4])\n\n@show inverse(p1)\n@show p1 * p2\nA group of permutations. Group axioms are verified during construction.
Sources permutation_group.hpp, permutation_group.cpp
"},{"location":"documentation/symmetries/permutation_group/#constructor","title":"Constructor","text":""},{"location":"documentation/symmetries/permutation_group/#from-permutations","title":"From Permutations","text":"Creates an PermutationGroup out of a vector of Permutation objects.
C++JuliaPermutationGroup(std::vector<Permutation> const &permutations);\nPermutationGroup(permutations::Vector{Permutation})\nCreates a PermutationGroup out of a matrix whose rows specify the individual permutations. If a raw pointer is handed, the matrix is assumed to be in column-major form.
C++PermutationGroup(arma::Mat<int64_t> const &matrix);\nPermutationGroup(int64_t *ptr, int64_t n_permutations, int64_t nsites);\nReturns the number of sites on which the permutations of the group acts.
C++Juliaint64_t nsites(PermutationGroup const &group);\nnsites(group::PermutationGroup)\nReturns the size of the permutation group, i.e. the number permutations.
C++Juliaint64_t size(PermutationGroup const &group);\nsize(group::PermutationGroup)\n// Define a cyclic group of order 3\nPermutation p1 = {0, 1, 2};\nPermutation p2 = {1, 2, 0};\nPermutation p3 = {2, 0, 1};\nauto C3 = PermutationGroup({p1, p2, p3});\n\nXDIAG_SHOW(C3.size());\nXDIAG_SHOW(C3.nsites());\nXDIAG_SHOW(C3.inverse(1)); // = 2\n# Define a cyclic group of order 3\np1 = Permutation([1, 2, 3])\np2 = Permutation([2, 3, 1])\np3 = Permutation([3, 1, 2])\nC3 = PermutationGroup([p1, p2, p3])\n\n@show size(C3)\n@show nsites(C3)\n@show inverse(C3, 1) # = 2\nA (1D) irreducible representation of a finite group. Upon creation, the group homomorphism properties are verified.
Source representation.hpp
"},{"location":"documentation/symmetries/representation/#constructors","title":"Constructors","text":""},{"location":"documentation/symmetries/representation/#trivial-representation","title":"Trivial representation","text":"Creates the trivial representation (all characters equal to 1) of a PermutationGroup
C++JuliaRepresentation(PermutationGroup const &group);\nRepresentation(group::PermutationGroup)\nCreates a 1D representation of a PermutationGroup with given real or complex characters.
C++template <typename T>\nRepresentation(PermutationGroup const &group, std::vector<T> const &characters);\ntemplate <typename T>\nRepresentation(PermutationGroup const &group, arma::Col<T> const &characters);\ntemplate <typename T>\nRepresentation(PermutationGroup const &group, T *characters, int64_t n_characters);\nThe template parameter T in C++ can either be double or complex.
Returns the size of the Representation, i.e. the number of group elements represented.
C++Juliaint64_t size(Representation const &irrep);\nsize(irrep::Representation)\nReturns the whether or not the Representation is real, I.E. the characters are real numbers and do not have an imaginary part.
C++Juliabool isreal(Representation const &irrep) const;\nisreal(irrep::Representation)\nMultiplies two Representations by overloading the * operator.
Representation operator*(Representation const &r1, Representation const &r2);\nBase.:*(r1::Representation, r2::Representation)\nPermutation p = {1, 2, 3, 0};\nauto C4 = PermutationGroup({pow(p, 0), pow(p, 1), pow(p, 2), pow(p, 3)});\nRepresentation r1(C4, arma::vec{1.0, -1.0, 1.0, -1.0});\nRepresentation r2(C4, arma::cx_vec{1.0, 1.0i, -1.0, -1.0i});\nXDIAG_SHOW(r1 * r2);\np = Permutation([2, 3, 4, 1])\nC4 = PermutationGroup(pow(p, 0), pow(p, 1), pow(p, 2), pow(p, 3))\nr1 = Representation(C4, [1, -1, 1, -1])\nr2 = Representation(C4, [1, 1im, -1, -1im])\n\n@show r1 * r2\nSources logger.hpp, logger.cpp
"},{"location":"documentation/utilities/logging/#setting-the-verbosity","title":"Setting the verbosity","text":"Algorithms implemented in XDiag do not output anything during their execution by default. However, it is typically useful to get some information on how the code is performing and even intermediary results at runtime. For this, the verbosity of the internal XDiag logging can be set using the function set_verbosity, which is defined as
void set_verbosity(int64_t level);\nset_verbosity(level::Integer);\nThere are several levels of verbosity, defining how much information is shown.
level outputed information 0 no information 1 some information 2 detailed informationFor example, when computing a ground state energy using the eigval0 function, we can set a higher verbosity level using
C++Juliaset_verbosity(2);\ndouble e0 = eigval0(bonds, block);\nset_verbosity(2);\ne0 = eigval0(bonds, block);\nThis will print detailed information, which can look like this
Lanczos iteration 1\nMVM: 0.00289 secs\nalpha: -0.2756971549998545\nbeta: 1.7639347562074059\neigs: -0.2756971549998545\nLanczos iteration 2\nMVM: 0.00244 secs\nalpha: -0.7116140394927443\nbeta: 2.3044797637130743\neigs: -2.2710052270892791 1.2836940325966804\nLanczos iteration 3\nMVM: 0.00210 secs\nalpha: -1.2772539678430306\nbeta: 2.6627870395174456\neigs: -3.7522788386927637 -0.6474957945455240 2.1352094709026579\nProducing nicely formatted output is unfortunately a bit cumbersome in standard C++. For this, the Log mechanism in XDiag can help. To simply write out a line of information you can call,
Log(\"hello from the logger\");\nBy default, a new line is added. It is also possible to set verbosity by handing the level as the first argument,
Log(2, \"hello from the logger only if global verbosity is set to >= 2\");\nThis message will only appear if the global verbosity level is set to a value \\(\\geq 2\\). Finally, XDiag also supports formatted output by using the fmtlib library. For example, numbers can be formated this way
Log(\"pi is around {:.4f} and the answer is {}\", 3.141592, 42);\nIn standard C++ measuring time is a bit awkward. To quickly monitor the CPU time spent by XDiag by simple functions.
"},{"location":"documentation/utilities/timing/#simple-timing-using-tic-toc","title":"Simple timing using tic() / toc()","text":"Similar as in Matlab one can use tic() and toc() to measure the time spent between two points in the code.
tic();\ndouble e0 = eigval0(bonds, block);\ntoc();\ntoc() will output the time spent since the last time tic() has been called.
To get the present time, simply call
auto time = rightnow();\nA timing (in second) between two time points can be written to output using
timing(begin, end);\nThis can even be accompanied by a message about what is being timed and a verbosity level (see Logging) can also be set. The full call signature is
timing(begin, end, message, level);\nrightnow() end end time computed using rightnow() message message string to be prepended to timing \"\" level verbosity level at which timing is printed 0"},{"location":"documentation/utilities/utils/","title":"Utility functions","text":""},{"location":"documentation/utilities/utils/#set_verbosity","title":"set_verbosity","text":"Set how much logging is generated my XDiag to monitor the progress and behaviour of the code. There are three verbosity levels that can be set:
This can be useful, e.g. to monitor the progress of an iterative algorithm
C++Juliavoid set_verbosity(int64_t level);\nset_verbosity(level::Int64)\nPrints a nice welcome message containing the version number and git commit used.
C++Juliavoid say_hello()\nsay_hello()\nIf say_hello is too much flower power for you, one can also just have a boring print-out of the version number using this function.
void print_version()\nprint_version()\nFor quick debugging in C++, XDiag features a simple macro which outputs the name and content of a variable calles XDIAG_SHOW(x). For example
Spinhalf block(16, 8);\nXDIAG_SHOW(block);\nwill write an output similar to
block:\n nsites : 16\n nup : 8\n dimension: 12,870\n ID : 0xa9127434d66b9878\nThe XDIAG_SHOW(x) macro can be used on any XDiag object and several other standard C++ objects as well.
cmake_minimum_required(VERSION 3.19)\n\nproject(\n tj_distributed_time_evolve\n)\n\nfind_package(xdiag_distributed REQUIRED HINTS ../../../install)\nadd_executable(main main.cpp)\ntarget_link_libraries(main PUBLIC xdiag::xdiag_distributed)\ncmake_minimum_required(VERSION 3.19)\n\nproject(\n hello_world\n)\n\nfind_package(xdiag REQUIRED HINTS \"/path/to/xdiag/install\")\nadd_executable(main main.cpp)\ntarget_link_libraries(main PRIVATE xdiag::xdiag)\n using XDiag\n say_hello()\n #include <xdiag/all.hpp>\n\n using namespace xdiag;\n\n int main() try {\n say_hello();\n } catch (Error e) {\n error_trace(e);\n }\nusing XDiag\n\nlet \n N = 16\n nup = N \u00f7 2\n block = Spinhalf(N, nup)\n\n # Define the nearest-neighbor Heisenberg model\n ops = OpSum()\n for i in 1:N\n ops += Op(\"HB\", \"J\", [i-1, i % N])\n end\n ops[\"J\"] = 1.0\n\n set_verbosity(2) # set verbosity for monitoring progress\n e0 = eigval0(ops, block) # compute ground state energy\n\n println(\"Ground state energy: $e0\")\nend\n#include <xdiag/all.hpp>\n\nusing namespace xdiag;\n\nint main() try {\n int N = 16;\n int nup = N / 2;\n Spinhalf block(N, nup);\n\n // Define the nearest-neighbor Heisenberg model\n OpSum ops;\n for (int i = 0; i < N; ++i) {\n ops += Op(\"HB\", \"J\", {i, (i + 1) % N});\n }\n ops[\"J\"] = 1.0;\n\n set_verbosity(2); // set verbosity for monitoring progress\n double e0 = eigval0(ops, block); // compute ground state energy\n\n Log(\"Ground state energy: {:.12f}\", e0);\n\n} catch (Error e) {\n error_trace(e);\n}\n#include <xdiag/all.hpp>\n\nusing namespace xdiag;\n\nvoid measure_density(int nsites, State const &v) {\n int rank;\n MPI_Comm_rank(MPI_COMM_WORLD, &rank);\n for (int i = 0; i < nsites; ++i) {\n complex sz = innerC(Bond(\"Ntot\", i), v);\n if (rank == 0) {\n printf(\"%.6f \", std::real(sz));\n }\n }\n if (rank == 0) {\n printf(\"\\n\");\n }\n}\n\nint main(int argc, char **argv) try {\n MPI_Init(&argc, &argv);\n\n int L = 6;\n int W = 4;\n double t = 1.0;\n double J = 0.1;\n double mu_0 = 10;\n\n int nsites = L * W;\n double precision = 1e-12;\n\n // Create square lattice t-J model\n OpSum ops;\n for (int x = 0; x < L-1; ++x) {\n for (int y = 0; y < W; ++y) {\n int nx = (x + 1) % L;\n int ny = (y + 1) % W;\n\n int site = x * W + y;\n int right = nx * W + y;\n int top = x * W + ny;\n ops += Op(\"Hop\", \"T\", {site, right});\n ops += Op(\"Exchange\", \"J\", {site, right});\n ops += Op(\"Hop\", \"T\", {site, top});\n ops += Op(\"Exchange\", \"J\", {site, top});\n\n\n\n if (x < L / 2) {\n Log(\"x {} y {} site {} t {} r {} +\", x, y, site, top, right);\n ops += Op(\"Ntot\", \"MUPLUS\", site);\n } else {\n Log(\"x {} y {} site {} t {} r {} -\", x, y, site, top, right);\n ops += Op(\"Ntot\", \"MUNEG\", site);\n }\n }\n }\n ops[\"T\"] = t;\n ops[\"J\"] = J;\n ops[\"MUPLUS\"] = mu_0;\n ops[\"MUNEG\"] = mu_0;\n\n auto block = tJDistributed(nsites, nsites / 2 - 1, nsites / 2 - 1);\n\n XDIAG_SHOW(block);\n\n Log.set_verbosity(2);\n tic();\n auto [e0, v] = eig0(ops, block);\n toc(\"gs\");\n\n ops[\"MUPLUS\"] = 0;\n ops[\"MUNEG\"] = 0;\n\n measure_density(nsites, v);\n\n // Do the time evolution with a step size tau\n double tau = 0.1;\n for (int i = 0; i < 40; ++i) {\n tic();\n v = time_evolve(ops, v, tau, precision);\n toc(\"time evolve\");\n tic();\n measure_density(nsites, v);\n toc(\"measure\");\n }\n\n MPI_Finalize();\n return EXIT_SUCCESS;\n} catch (std::exception const &e) {\n traceback(e);\n}\nLet us set up our first program using the xdiag library.
using XDiag\nsay_hello()\n#include <xdiag/all.hpp>\n\nusing namespace xdiag;\n\nint main() try {\n say_hello();\n} catch (Error e) {\n error_trace(e);\n}\nThe function say_hello() prints out a welcome message, which also contains information which exact XDiag version is used. In Julia this is all there is to it.
For the C++ code we need to create two files to compile the program. The first is the actual C++ code. What is maybe a bit unfamiliar is the try / catch block. XDiag implements a traceback mechanism for runtime errors, which is activated by this idiom. While not stricly necessary here, it is a good practice to make use of this.
Now that the application program is written, we next need to set up the compilation instructions using CMake. To do so we create a second file called CMakeLists.txt in the same directory.
cmake_minimum_required(VERSION 3.19)\n\nproject(\n hello_world\n)\n\nfind_package(xdiag REQUIRED HINTS \"/path/to/xdiag/install\")\nadd_executable(main main.cpp)\ntarget_link_libraries(main PRIVATE xdiag::xdiag)\nYou should replace \"/path/to/xdiag/install\" with the appropriate directory where your XDiag library is installed after compilation. This exact CMakeLists.txt file can be used to compile any XDiag application.
Info
For using the distributed XDiag library the last line of the above CMakeLists.txt should be changed to
target_link_libraries(main PUBLIC xdiag::xdiag_distributed)\nWe then compile the application code,
cmake -S . -B build\ncmake --build build\n./build/main\nWe are now ready to run our first actual calculation using XDiag. Our immediate goal will be to determine the ground state energy of a \\(S=1/2\\) Heisenberg model on a 1D chain lattice with periodic boundary conditions,
\\[ H = J \\sum_{\\langle i, j\\rangle} \\mathbf{S}_i \\cdot \\mathbf{S}_j,\\]where \\(\\mathbf{S}_i = (S^x_i, S^y_i, S^z_i)\\) denotes the vector of spin matrices at a given site \\(i\\). The notation \\(\\langle i, j\\rangle\\) refers to summatation over neighboring sites \\(i\\) and \\(j\\).
The first thing to define before any computation, is the Hilbert space our model will be defined on. For the present example, we use the Hilbert space class Spinhalf. We consider a chain lattice with \\(N=8\\) sites and create a \\(S=1/2\\) Hilbert space:
C++using namespace xdiag;\nint N = 8;\nauto hspace = Spinhalf(N);\nWe would like to know which spin configurations, the Hilbert space is made up of. To do so, we can iterate over the Hilbert space and print out the spin configurations.
C++for (auto spins : hspace) {\n Log(\"{}\", to_string(spins));\n}\nThis produces an output similar to the following:
\u2193\u2193\u2193\u2193\u2193\u2193\u2193\u2193\n\u2193\u2193\u2193\u2193\u2193\u2193\u2193\u2191\n\u2193\u2193\u2193\u2193\u2193\u2193\u2191\u2193\n\u2193\u2193\u2193\u2193\u2193\u2193\u2191\u2191\n\u2193\u2193\u2193\u2193\u2193\u2191\u2193\u2193\n\u2193\u2193\u2193\u2193\u2193\u2191\u2193\u2191\n\u2193\u2193\u2193\u2193\u2193\u2191\u2191\u2193\n...\nHere we already see several things at work. XDiag features a convenient way to write logs, with the Log class. The first argument to Log() is a format string. In C++ we use the fmt library, to be able to write structured format and format our output. The second argument turns our spins into a string. spins is of type ProductState, whose configuration on each site can be individually addressed.
Further, we notice that all \\(2^N\\) spin configurations are included in this Hilbert space. However, the Heisenberg model conserves the total \\(S^z = \\sum_i S^z_i\\), and thus we could limit ourselves to a block of the Hilbert space, which only contains configurations of a certain magnetization:
C++int nup = 4;\nauto block = Spinhalf(N, nup);\nfor (auto spins : block) {\n Log(\"{}\", to_string(spins));\n}\nThis produces an output similar to:
\u2193\u2193\u2193\u2193\u2191\u2191\u2191\u2191\n\u2193\u2193\u2193\u2191\u2193\u2191\u2191\u2191\n\u2193\u2193\u2193\u2191\u2191\u2193\u2191\u2191\n\u2193\u2193\u2193\u2191\u2191\u2191\u2193\u2191\n\u2193\u2193\u2193\u2191\u2191\u2191\u2191\u2193\n\u2193\u2193\u2191\u2193\u2193\u2191\u2191\u2191\n\u2193\u2193\u2191\u2193\u2191\u2193\u2191\u2191\n...\nWe see that the block now only contains configurations, where the total number of spins pointing up is 4. We can now quickly check the dimension of the Hilbert spaces, and confirm that the dimension of the block is reduced from \\(2^8=256\\) to \\(\\begin{pmatrix} 8 \\\\ 4 \\end{pmatrix} = 70\\),
C++XDIAG_SHOW(hspace.size());\nXDIAG_SHOW(block.size());\nwhich should print:
hspace.size():\n256\nblock.size():\n70\nHere, we introduced another functionality of XDiag, the XDIAG_SHOW macro which can be used for quick debug printing of XDiag objects.
"},{"location":"user_guide/input_output/","title":"Input/Output","text":"Julia features a variety of packages facilitating input and output of data. For C++, we provide convenient functionality for TOML and HDF5 files.
"},{"location":"user_guide/input_output/#toml","title":"TOML","text":"For simulations is can often be useful to read input paramters from a file. The
"},{"location":"user_guide/installation/","title":"Installation","text":""},{"location":"user_guide/installation/#julia-installation","title":"Julia installation","text":"Enter the package mode using ] in the Julia REPL and type:
add XDiag\nThat's it!
"},{"location":"user_guide/installation/#c-compilation","title":"C++ compilation","text":"Using XDiag with C++ is a two-step process. First the xdiag library needs to be compiled and installed. Therafter, application codes are compiled in a second step. The library can be compiled in two different versions:
g++, clang, or Intel's icpx)Download the source code using git
cd /path/to/where/xdiag/should/be\ngit clone https://github.com/awietek/xdiag.git\nCompile the default library
cd xdiag\ncmake -S . -B build\ncmake --build build\ncmake --install build\ninstall. Compile the distributed library
To use the distributed computing features of xdiag, the distributed library has to be built which requires MPI.
cd xdiag\ncmake -S . -B build -D XDIAG_DISTRIBUTED=On\ncmake --build build\ncmake --install build\nInfo
It might be necessary to explicitly define MPI compiler, e.g. mpicxx like this
cmake -S . -B build -D XDIAG_DISTRIBUTED=On -D CMAKE_CXX_COMPILER=mpicxx\nThe compilation process can be modified and also allows for further optimizations. We collect several common scenarios in the Advanced Compilation guide.
"}]} \ No newline at end of file +{"config":{"lang":["en"],"separator":"[\\s\\-]+","pipeline":["stopWordFilter"]},"docs":[{"location":"","title":"Home","text":"Quick Start Code on GitHub
"},{"location":"#overview","title":"Overview","text":"XDiag is a library for performing Exact Diagonalizations of quantum many-body systems. Key features include optimized combinatorical algorithms for navigating Hilbert spaces, iterative linear algebra algorithms, shared and distributed memory parallelization. It consist of two packages:
Please support our work by citing XDiag and the implemented algorithms if it is used in your published research.
@article{Wietek2018,\n title = {Sublattice coding algorithm and distributed memory parallelization for large-scale exact diagonalizations of quantum many-body systems},\n author = {Wietek, Alexander and L\\\"auchli, Andreas M.},\n journal = {Phys. Rev. E},\n volume = {98},\n issue = {3},\n pages = {033309},\n numpages = {10},\n year = {2018},\n month = {Sep},\n publisher = {American Physical Society},\n doi = {10.1103/PhysRevE.98.033309},\n url = {https://link.aps.org/doi/10.1103/PhysRevE.98.033309}\n}\nSupporting material for lecture held at Quant24 master's school at MPI PKS. Consists of a Jupyter notebook and a sample lattice file describing the \\(N=12\\) site triangular lattice Heisenberg model:
This notebook uses the Julia verision of XDiag and covers the basic functionality:
Hello World
Prints out a greeting containing information on the version of the code.
source
Groundstate energy
Computes the ground state energy of a simple Heisenberg spin \\(S=1/2\\) chain
source
\\(t\\)-\\(J\\) time evolution
Computes the time evolution of a state in the \\(t\\)-\\(J\\) model with distributed parallelization
source
Normal XDiag
Template CMakeLists.txt which can be used to compile applications with the normal XDiag library.
source
Distributed XDiag
Template CMakeLists.txt which can be used to compile applications with the distributed XDiag library.
source
Enter the package mode using ] in the Julia REPL and type:
add XDiag\nThat's it!
"},{"location":"installation/#c-compilation","title":"C++ Compilation","text":"Using XDiag with C++ is a two-step process. First the xdiag library needs to be compiled and installed. Therafter, application codes are compiled in a second step. Here we explain how to compile the library.
g++, clang, or Intel's icpx)Download the source code using git
cd /path/to/where/xdiag/should/be\ngit clone https://github.com/awietek/xdiag.git\nCompile the default library
cd xdiag\ncmake -S . -B build\ncmake --build build\ncmake --install build\ninstall. Compile the distributed library
To use the distributed computing features of xdiag, the distributed library has to be built which requires MPI.
cd xdiag\ncmake -S . -B build -D XDIAG_DISTRIBUTED=On\ncmake --build build\ncmake --install build\nInfo
It might be necessary to explicitly define MPI compiler, e.g. mpicxx like this
cmake -S . -B build -D XDIAG_DISTRIBUTED=On -D CMAKE_CXX_COMPILER=mpicxx\nParallel compilation To speed up the compilation process, the build step can be performed in parallel using the -j flag
cmake --build build -j\nListing compile options
The available compilation options can be displayed using
cmake -L .\nChoosing a certain compiler
The compiler (e.g. icpx) can be specified using
cmake -S . -B build -D CMAKE_CXX_COMPILER=icpx\nWarning
If the xdiag library is compiled with a certain compiler, it is advisable to also compile the application codes with the same compiler.
Setting the install path
In the installation step, the install directory can be set in the following way
cmake --install build --prefix /my/install/prefix\nDisabling HDF5/OpenMP
To disable support for HDF5 or OpenMP support, use
cmake -S . -B build -D XDIAG_DISABLE_OPENMP=On -D XDIAG_DISABLE_HDF5=On\nBuilding and running tests
To compile and run the testing programs, use
cmake -S . -B build -D BUILD_TESTING=On\ncmake --build build\nbuild/tests/tests\nBuilding the Julia wrapper locally
First, get the path to the CxxWrap package of julia. To do so, enter the Julia REPL,
julia\nusing CxxWrap\nCxxWrap.prefix_path()\n/path/to/libcxxwrap-julia-prefix. This is then used to configure the cmake compilation. cmake -S . -B build -D XDIAG_JULIA_WRAPPER=On -D CMAKE_PREFIX_PATH=/path/to/libcxxwrap-julia-prefix\ncmake --build build\ncmake --install build\nlibxdiagjl.so, (or the corresponding library format on non-Linux systems). The source files for the documentation can be found in the directory docs. The documentation is built using Material for MKDocs. To work on it locally, it can be served using
mkdocs serve\nfrom the xdiag root source directory. A local build of the documentation can then be accessed in a webbrowser at the adress
127.0.0.1:8000\nLet us set up our first program using the xdiag library.
using XDiag\nsay_hello()\n#include <xdiag/all.hpp>\n\nusing namespace xdiag;\n\nint main() try {\n say_hello();\n} catch (Error e) {\n error_trace(e);\n}\nThe function say_hello() prints out a welcome message, which also contains information which exact XDiag version is used. In Julia this is all there is to it.
For the C++ code we need to create two files to compile the program. The first is the actual C++ code. What is maybe a bit unfamiliar is the try / catch block. XDiag implements a traceback mechanism for runtime errors, which is activated by this idiom. While not stricly necessary here, it is a good practice to make use of this.
Now that the application program is written, we next need to set up the compilation instructions using CMake. To do so we create a second file called CMakeLists.txt in the same directory.
cmake_minimum_required(VERSION 3.19)\n\nproject(\n hello_world\n)\n\nfind_package(xdiag REQUIRED HINTS \"/path/to/xdiag/install\")\nadd_executable(main main.cpp)\ntarget_link_libraries(main PRIVATE xdiag::xdiag)\nYou should replace \"/path/to/xdiag/install\" with the appropriate directory where your XDiag library is installed after compilation. This exact CMakeLists.txt file can be used to compile any XDiag application.
Info
For using the distributed XDiag library the last line of the above CMakeLists.txt should be changed to
target_link_libraries(main PUBLIC xdiag::xdiag_distributed)\nWe then compile the application code,
cmake -S . -B build\ncmake --build build\nand finally run our first xdiag application.
./build/main\nWe compute the ground state energy of the \\(S=1/2\\) Heisenberg chain on a periodic chain lattice in one dimension. The Hamiltonian is given by
\\[ H = J\\sum_{\\langle i,j \\rangle} \\mathbf{S}_i \\cdot \\mathbf{S}_j\\]where \\(\\mathbf{S}_i = (S_i^x, S_i^y, S_i^z)\\) are the spin \\(S=1/2\\) operators and \\(\\langle i,j \\rangle\\) denotes summation over nearest-meighbor sites \\(i\\) and \\(j\\).
The following code, sets up the Hilbert space, defines the Hamiltonian and finally calls an iterative eigenvalue solver to compute the ground state energy.
JuliaC++using XDiag\n\nlet \n N = 16\n nup = N \u00f7 2\n block = Spinhalf(N, nup)\n\n # Define the nearest-neighbor Heisenberg model\n ops = OpSum()\n for i in 1:N\n ops += Op(\"HB\", \"J\", [i-1, i % N])\n end\n ops[\"J\"] = 1.0\n\n set_verbosity(2) # set verbosity for monitoring progress\n e0 = eigval0(ops, block) # compute ground state energy\n\n println(\"Ground state energy: $e0\")\nend\n#include <xdiag/all.hpp>\n\nusing namespace xdiag;\n\nint main() try {\n int N = 16;\n int nup = N / 2;\n Spinhalf block(N, nup);\n\n // Define the nearest-neighbor Heisenberg model\n OpSum ops;\n for (int i = 0; i < N; ++i) {\n ops += Op(\"HB\", \"J\", {i, (i + 1) % N});\n }\n ops[\"J\"] = 1.0;\n\n set_verbosity(2); // set verbosity for monitoring progress\n double e0 = eigval0(ops, block); // compute ground state energy\n\n Log(\"Ground state energy: {:.12f}\", e0);\n\n} catch (Error e) {\n error_trace(e);\n}\nSep. 9, 2024
Introduced 1-indexing everywhere in Julia version
Aug. 27, 2024
Lanczos routines and multicolumn States
Aug. 16, 2024
Small patch release providing small utility functions
Aug. 15, 2024
Basic functionality for three Hilbert space types, Spinhalf, tJ, and Electron, has been implemented. Features are:
A step-by-step guide to using XDiag
Installation Language Julia installation How to install the XDiag in Julia C++ compilation How to compile the XDiag C++ library First steps Language Writing code How to write an application using XDiag Compilation How to compile an application using XDiag in C++ Hilbert spaces How to define a Hilbert space Operators How to define an operator Ground state How compute a ground state (energy) Measurements How to measure an observable on a state Input/Output Language TOML How to read input from a TOML file"},{"location":"documentation/","title":"Documentation","text":"XDiag uses the C++ library Armadillo as a linear algebra backend. Documentation for linear algebra operations can, therefore, be found in the Armadillo Documentation.
"},{"location":"documentation/#algebra","title":"Algebra","text":"Name Description Language matrix Creates the full matrix representation of an operator on a block apply Applies an operator to a state \\(\\vert \\phi \\rangle = O \\vert \\psi\\rangle\\) dot Computes the dot product between two states inner Computes an expectation value \\(\\langle \\psi \\vert O \\vert \\psi \\rangle\\) norm Computes the 2-norm of a state norm1 Computes the 1-norm of a state norminf Computes the \\(\\infty\\)-norm of a state"},{"location":"documentation/#algorithms","title":"Algorithms","text":"Diagonalization
Name Description Language eigval0 Computes the lowest lying eigenvalue of an operator eig0 Computes the lowest lying eigenvalue and eigenvector of an operator eigvals_lanczos Performs an iterative eigenvalue calculation using the Lanczos algorithm eigs_lanczos Performs an iterative eigenvalue calculation building eigenvectors using the Lanczos algorithmTime evolution
Name Description Language time_evolve Performs a real-time evolution \\(e^{ -iHt} \\vert \\psi \\rangle\\) of a State with a given Hermitian operator \\(H\\) imaginary_time_evolve Performs a imaginary-time evolution \\(e^{ -\\tau H}\\vert\\psi\\rangle\\) of a State with a given Hermitian operator \\(H\\) evolve_lanczos Computes the exponential $e^{z H}\\vert\\psi\\rangle $ of a Hermitian operator times a State for a real or complex \\(z\\) using the Lanczos algorithm time_evolve_expokit Performs a real-time evolution \\(e^{ -iHt} \\vert \\psi \\rangle\\) using a highly accurate Lanczos algorithm"},{"location":"documentation/#blocks","title":"Blocks","text":"Shared memory
Name Description Language Spinhalf Block of a spin \\(S=1/2\\) type Hilbert space tJ Block of a \\(t-J\\) type Hilbert space Electron Block of a Electron type Hilbert spaceDistributed memory
Name Description Language SpinhalfDistributed Block of a spin \\(S=1/2\\) type Hilbert space (distributed computing) tJDistributed Block of a \\(t-J\\) type Hilbert space (distributed computing)"},{"location":"documentation/#building","title":"Building","text":"Name Description Language Compilation Advanced settings for compilation of the C++ library using CMake Documentation How to build and work on the documentation locally Julia Wrapper How to build and develop the Julia wrapper locally"},{"location":"documentation/#input-output","title":"Input / Output","text":"TOML
Name Description Language FileToml A file handler for TOML files read_opsum reads an OpSum from a TOML file read_permutation_group reads a PermutationGroup from a TOML file read_representation reads a Representation from a TOML fileHDF5
Name Description Language FileH5 A file handler for hdf5 files"},{"location":"documentation/#operators","title":"Operators","text":"Operator types A summary of all the operator types defined in XDiag Name Description Language Op A linear operator acting on the Hilbert space OpSum Sum of couplings times operators hc Returns the hermitian conjugate of an Op or OpSum. symmetrize Symmetrizes an operator with a PermutationGroup or Representation Scalar A scalar number which can be either real or complex Coupling Describes a coupling associated with an operator, either string or scalar"},{"location":"documentation/#states","title":"States","text":"Name Description Language State A generic state describing a quantum wave function ProductState A product state of local configurations RandomState A random state with normal distributed coefficients fill Fill a state with a given model state product_state Creates a filled product state random_state Create a filled random state with normal distributed coefficients zero_state Create a filled state with all zero entries zero Set all coefficients of a given state to zero"},{"location":"documentation/#symmetries","title":"Symmetries","text":"Name Description Language Permutation Permutations of indices or lattice sites PermutationGroup A group of permutations Representation A (1D) irreducible representation of a finite group"},{"location":"documentation/#utilities","title":"Utilities","text":"Name Description Language set_verbosity Sets how much information is printed during computations say_hello Prints a nice welcome message with version number print_version Prints the plain version number Logging Controling what is written to standard output Timing Measurng wall time straightforwardly XDIAG_SHOW Macro for printing debugging information"},{"location":"documentation/algebra/algebra/","title":"Algebra routines","text":"Several basic algebraic operations for states and operators.
Sources algebra.hpp, algebra.cpp
"},{"location":"documentation/algebra/algebra/#dot","title":"dot","text":"Computes the dot product \\(\\langle v \\vert w \\rangle\\) between two states \\(\\vert v \\rangle\\) and \\(\\vert w \\rangle\\). In C++, please use the dotC function if one of the two states is expected to be complex.
C++Juliadouble dot(State const &v, State const &w);\ncomplex dotC(State const &v, State const &w);\ndot(v::State, w::State)\nComputes the expectation value \\(\\langle v | O |v \\rangle\\) of an operator \\(O\\) and a state \\(|v\\rangle\\). The operator can either be an Op or an OpSum object. In C++, please use the innerC function if either the operator or the state are complex.
C++Juliadouble inner(Op const &op, State const &v);\ndouble inner(OpSum const &ops, State const &v);\ncomplex innerC(Op const &op, State const &v);\ncomplex innerC(OpSum const &ops, State const &v);\ninner(op::Op, v::State)\ninner(ops::OpSum, v::State)\nComputes the 2-norm \\(\\parallel |v \\rangle \\parallel_2\\) of a state \\(|v \\rangle\\) defined as
\\[ \\parallel |v \\rangle \\parallel_2 = \\sum_n |\\langle n | v \\rangle |^2, \\]where \\(\\{ |n\\rangle \\}\\) denotes an orthonormal basis of the block.
C++Juliadouble norm(State const &v);\nnorm(state::State)\nComputes the 1-norm \\(\\parallel |v \\rangle \\parallel_1\\) of a state \\(|v \\rangle\\) defined as
\\[ \\parallel |v \\rangle \\parallel_1 = \\sum_n |\\langle n | v \\rangle |, \\]where \\(\\{ |n\\rangle \\}\\) denotes an orthonormal basis of the block.
C++Juliadouble norm1(State const &v);\nnorm1(state::State)\nComputes the \\(\\infty\\)-norm \\(\\parallel |v \\rangle \\parallel_\\infty\\) of a state \\(|v \\rangle\\) defined as
\\[ \\parallel |v \\rangle \\parallel_\\infty = \\max_n |\\langle n | v \\rangle |, \\]where \\(\\{ |n\\rangle \\}\\) denotes an orthonormal basis of the block.
C++Juliadouble norminf(State const &v);\nnorminf(state::State)\nint N = 8;\nauto block = Spinhalf(N, N / 2);\nauto ops = OpSum();\nfor (int i=0; i<N; ++i) {\n ops += Op(\"SdotS\", {i, (i+1)%N});\n}\nauto [e0, psi] = eig0(ops, block);\n\nXDIAG_SHOW(norm(psi));\nXDIAG_SHOW(norm1(psi));\nXDIAG_SHOW(norminf(psi));\n\nXDIAG_SHOW(dot(psi, psi));\nXDIAG_SHOW(e0);\nXDIAG_SHOW(inner(ops, psi));\n\nauto phi = random_state(block);\nXDIAG_SHOW(phi.vector());\nXDIAG_SHOW(psi.vector());\nXDIAG_SHOW((psi + 2.0*phi).vector());\nXDIAG_SHOW((psi*complex(0,3.0) + phi/2.0).vectorC());\nlet \n N = 8\n block = Spinhalf(N, N \u00f7 2)\n ops = OpSum()\n for i in 1:N\n ops += Op(\"HB\", 1.0, [i, mod1(i+1, N)])\n end\n e0, psi = eig0(ops, block);\n\n @show norm(psi)\n @show norm1(psi)\n @show norminf(psi)\n\n @show dot(psi, psi)\n @show e0, inner(ops, psi)\n\n phi = rand(block)\n display(vector(phi))\n display(vector(psi))\n display(vector(psi + 2.0*phi))\n display(vector(psi*3.0im + phi/2.0))\nend\nApplies an operator given as an Op or OpSum to a State \\(\\vert w \\rangle = \\mathcal{O} \\vert v\\rangle\\).
Sources apply.hpp, apply.cpp
"},{"location":"documentation/algebra/apply/#definition","title":"Definition","text":"An operator \\(\\mathcal{O}\\) can be applied to a state \\(\\vert v\\rangle\\) in two ways:
Only the input state on which the operator acts is defined is given. The block of the output state is calculated and eventually created automatically.
C++JuliaState apply(Op const &op, State const &v);\nState apply(OpSum const &ops, State const &v);\napply(op::Op, v::State)\napply(ops::OpSum, v::State)\nThe output state is also handed as an argument which is overwritten. The compatibility of quantum numbers is checked. This way the output block is not created automatically and, thus, can be used to save computation time if the output block appears repeatedly in the computation.
C++Juliavoid apply(Op const &op, State const &v, State &w);\nvoid apply(OpSum const &ops, State const &v, State &w);\napply(op::Op, v::State, w::State)\napply(ops::OpSum, v::State, w::State)\nint N = 8;\nauto block = Spinhalf(N, N / 2);\nauto ops = OpSum();\nfor (int i=0; i<N; ++i){\n ops += Op(\"SdotS\", {i, (i+1)%N});\n}\nauto [e0, psi] = eig0(ops, block);\nauto phi = apply(Op(\"S+\", 2), psi);\nXDIAG_SHOW(inner(ops, psi));\nXDIAG_SHOW(inner(ops, phi));\nlet \n N = 8\n block = Spinhalf(N, N \u00f7 2)\n ops = OpSum()\n for i in 1:N\n ops += Op(\"HB\", 1.0, [i, mod1(i+1, N)])\n end\n e0, psi = eig0(ops, block);\n\n blockp = Spinhalf(N, N \u00f7 2 + 1)\n phi = zeros(blockp)\n apply(Op(\"S+\", 1.0, 2), psi, phi)\n @show inner(ops, psi)\n @show inner(ops, phi)\nend\nCreates a numerical matrix with real (matrix) or complex (matrixC) coefficients given an Op or OpSum on a certain block.
Sources matrix.hpp, matrix.cpp
"},{"location":"documentation/algebra/matrix/#definition","title":"Definition","text":"A matrix can be created in two ways:
Only the input block on which the operator is defined is given. The output block is calculated and eventually created automatically.
C++Juliatemplate <class block_t>\narma::mat matrix(Op const &op, block_t const &block);\n\ntemplate <class block_t>\narma::mat matrix(OpSum const &ops, block_t const &block);\n\ntemplate <class block_t>\narma::cx_mat matrixC(Op const &op, block_t const &block);\n\ntemplate <typename block_t>\narma::cx_mat matrixC(OpSum const &ops, block_t const &block);\nmatrix(op::Op, block::Block)\nmatrix(ops::OpSum, block::Block)\nThe output block is also handed as an argument. The compatibility of quantum numbers is checked. This way the output block is not created automatically and, thus, can be used to save computation time if the output block appears repeatedly in the computation.
C++Juliatemplate <class block_t>\narma::mat matrix(Op const &op, block_t const &block_in, \n block_t const &block_out);\n\ntemplate <class block_t>\narma::mat matrix(OpSum const &ops, block_t const &block_in, \n block_t const &block_out);\n\ntemplate <class block_t>\narma::cx_mat matrixC(Op const &op, block_t const &block_in, \n block_t const &block_out);\n\ntemplate <typename block_t>\narma::cx_mat matrixC(OpSum const &ops, block_t const &block_in, \n block_t const &block_out);\nmatrix(op::Op, block_in::Block, block_out::Block)\nmatrix(ops::OpSum, block_in::Block, block_out::Block)\nComment: In Julia, depending on whether a real/complex matrix is generated also a real/complex matrix is returned. The C++ version has to return a fixed type. If a real matrix is desired, use the function matrix. If a complex matrix is desired, use the function matrixC.
// Creates matrix H_{k=2} in Eq (18.23) of https://link.springer.com/content/pdf/10.1007/978-3-540-74686-7_18.pdf\nint N = 4;\nint nup = 3;\nint ndn = 2;\n\n// Define a Hubbard chain model\nauto ops = OpSum();\nfor (int i=0; i< N; ++i){\n ops += \"T\" * Op(\"Hop\", {i, (i+1) % N});\n}\nops+= \"U\" * Op(\"HubbardU\");\nops[\"T\"] = 1.0;\nops[\"U\"] = 5.0;\n\n// Create the a permutation group\nauto p1 = Permutation({0, 1, 2, 3});\nauto p2 = Permutation({1, 2, 3, 0});\nauto p3 = Permutation({2, 3, 0, 1});\nauto p4 = Permutation({3, 0, 1, 2});\nauto group = PermutationGroup({p1, p2, p3, p4});\nauto irrep = Representation(group, arma::vec{1.0, -1.0, 1.0, -1.0});\nauto block = Electron(N, nup, ndn, irrep);\nauto H = matrix(ops, block);\nH.print();\nlet\n # Creates matrix H_{k=2} in Eq (18.23) of https://link.springer.com/content/pdf/10.1007/978-3-540-74686-7_18.pdf\n N = 4\n nup = 3\n ndn = 2\n\n # Define a Hubbard chain model\n ops = OpSum()\n for i in 1:N\n ops += Op(\"Hop\", \"T\", [i, mod1(i+1, N)])\n end\n ops[\"T\"] = 1.0;\n ops[\"U\"] = 5.0;\n\n # Create the a permutation group\n p1 = Permutation([1, 2, 3, 4])\n p2 = Permutation([2, 3, 4, 1])\n p3 = Permutation([3, 4, 1, 2])\n p4 = Permutation([4, 1, 2, 3])\n group = PermutationGroup([p1, p2, p3, p4])\n irrep = Representation([1, -1, 1, -1])\n block = Electron(N, nup, ndn, group, irrep)\n\n H = matrix(ops, block)\n display(H)\nend\nComputes the groud state energy and the ground state of a Hermitian operator on a block by using an iterative Lanczos algorithm. This function is a shortcut for the eigs_lanczos function. We refer to eigs_lanczos for further details on the algorithm and the convergence criterion.
Sources sparse_diag.hpp, sparse_diag.cpp
"},{"location":"documentation/algorithms/eig0/#definition","title":"Definition","text":"C++Juliastd::tuple<double, State> eig0(OpSum const &ops, Block const &block,\n double precision = 1e-12, int64_t max_iterations = 1000,\n bool force_complex = false, int64_t random_seed = 42);\nfunction eig0(\n ops::OpSum,\n block::Block;\n precision::Real = 1e-12,\n maxiter::Int64 = 1000,\n force_complex::Bool = false,\n seed::Int64 = 42,\n)\nops State groundstate"},{"location":"documentation/algorithms/eig0/#usage-example","title":"Usage Example","text":"C++Julia int N = 8;\nint nup = N / 2;\nauto block = Spinhalf(N, nup);\n\n// Define the nearest-neighbor Heisenberg model\nauto ops = OpSum();\nfor (int i=0; i<N; ++i) {\n ops += \"J\" * Op(\"SdotS\", {i, (i+1) % N});\n}\nops[\"J\"] = 1.0;\nauto [e0, gs] = eig0(ops, block);\nlet \n N = 8\n nup = N \u00f7 2\n block = Spinhalf(N, nup)\n\n # Define the nearest-neighbor Heisenberg model\n ops = OpSum()\n for i in 1:N\n ops += Op(\"HB\", \"J\", [i, mod1(i+1, N)])\n end\n ops[\"J\"] = 1.0;\n\n e0, gs = eig0(ops, block);\nend\nPerforms an iterative eigenvalue calculation building eigenvectors using the Lanczos algorithm. Returns the tridiagonal matrix, eigenvalues, number of iterations and the stopping criterion. The Lanczos interations are performed twice, where at the second run the eigenvectors are built.
Sources eigs_lanczos.hpp, eigs_lanczos.cpp
"},{"location":"documentation/algorithms/eigs_lanczos/#definition","title":"Definition","text":"The Lanczos algorithm can be run in two distinct ways:
A random intial state \\(|\\psi_0\\rangle = |r\\rangle\\) with normal distributed entries is used.
C++Juliaeigs_lanczos_result_t\neigs_lanczos(OpSum const &ops, Block const &block, int64_t neigvals = 1,\n double precision = 1e-12, int64_t max_iterations = 1000,\n bool force_complex = false, double deflation_tol = 1e-7,\n int64_t random_seed = 42);\nfunction eigs_lanczos(\n ops::OpSum,\n block::Block;\n neigvals::Int64 = 1,\n precision::Float64 = 1e-12,\n max_iterations::Int64 = 1000,\n force_complex::Bool = false,\n deflation_tol::Float64 = 1e-7,\n random_seed::Int64 = 42,\n)\nThe initial state \\(|\\psi_0\\rangle\\) is explicitly specified
C++eigs_lanczos_result_t \neigs_lanczos(OpSum const &ops, State &state0, int64_t neigvals = 1,\n double precision = 1e-12, int64_t max_iterations = 1000,\n bool force_complex = false, double deflation_tol = 1e-7);\nA struct with the following entries
Entry Description alphas diagonal elements of the tridiagonal matrix betas off-diagonal elements of the tridiagonal matrix eigenvalues the computed Ritz eigenvalues of the tridiagonal matrix eigenvectors State of shape $D \\times $neigvals holding all low-lying eigenvalues up to neigvals niterations number of iterations performed criterion string denoting the reason why the algorithm stopped"},{"location":"documentation/algorithms/eigs_lanczos/#convergence-criterion","title":"Convergence criterion","text":"The algorithm terminates if the \\(k\\)-th (\\(k\\) is the argument neigvals) approximate eigenvalue changes only by a fraction smaller than \\(\\epsilon\\) (\\(k\\) is the argument precision), i.e.
Here, \\(\\tilde{e}_k^{(n)}\\) denotes the Lanczos approximation to the \\(k\\)-th eigenvalue after \\(n\\) iterations.
"},{"location":"documentation/algorithms/eigs_lanczos/#usage-example","title":"Usage Example","text":"C++int N = 8;\nint nup = N / 2;\nauto block = Spinhalf(N, nup);\n\n// Define the nearest-neighbor Heisenberg model\nauto ops = OpSum();\nfor (int i=0; i<N; ++i) {\n ops += \"J\" * Op(\"SdotS\", {i, (i+1) % N});\n}\nops[\"J\"] = 1.0;\n\n// With random intial state\nauto res = eigs_lanczos(ops, block);\nXDIAG_SHOW(res.alphas);\nXDIAG_SHOW(res.betas);\nXDIAG_SHOW(res.eigenvalues);\nXDIAG_SHOW(res.eigenvectors);\nComputes the groud state energy of a Hermitian operator on a block by using an iterative Lanczos algorithm. This function is a shortcut for the eigvals_lanczos function. We refer to eigvals_lanczos for further details on the algorithm and the convergence criterion.
Sources sparse_diag.hpp, sparse_diag.cpp
"},{"location":"documentation/algorithms/eigval0/#definition","title":"Definition","text":"C++Juliadouble eigval0(OpSum const &ops, Block const &block, double precision = 1e-12,\n int64_t max_iterations = 1000, bool force_complex = false,\n int64_t random_seed = 42);\nfunction eigval0(\n ops::OpSum,\n block::Block;\n precision::Real = 1e-12,\n maxiter::Integer = 1000,\n force_complex::Bool = false,\n seed::Integer = 42,\n)\nops"},{"location":"documentation/algorithms/eigval0/#usage-example","title":"Usage Example","text":"C++Julia int N = 8;\nint nup = N / 2;\nauto block = Spinhalf(N, nup);\n\n// Define the nearest-neighbor Heisenberg model\nauto ops = OpSum();\nfor (int i=0; i<N; ++i) {\n ops += \"J\" * Op(\"SdotS\", {i, (i+1) % N});\n}\nops[\"J\"] = 1.0;\ndouble e0 = eigval0(ops, block);\nlet \n N = 8\n nup = N \u00f7 2\n block = Spinhalf(N, nup)\n\n # Define the nearest-neighbor Heisenberg model\n ops = OpSum()\n for i in 1:N\n ops += Op(\"HB\", \"J\", [i, mod1(i+1, N)])\n end\n ops[\"J\"] = 1.0\n\n e0 = eigval0(ops, block);\nend\nPerforms an iterative eigenvalue calculation using the Lanczos algorithm. Returns the tridiagonal matrix, eigenvalues, number of iterations and the stopping criterion.
Sources eigvals_lanczos.hpp, eigvals_lanczos.cpp
"},{"location":"documentation/algorithms/eigvals_lanczos/#definition","title":"Definition","text":"The Lanczos algorithm can be run in thre distinct ways:
A random intial state \\(|\\psi_0\\rangle = |r\\rangle\\) with normal distributed entries is used.
C++Juliaeigvals_lanczos_result_t\neigvals_lanczos(OpSum const &ops, Block const &block, int64_t neigvals = 1,\n double precision = 1e-12, int64_t max_iterations = 1000,\n bool force_complex = false, double deflation_tol = 1e-7,\n int64_t random_seed = 42);\nfunction eigvals_lanczos(\n ops::OpSum,\n block::Block;\n neigvals::Int64 = 1,\n precision::Float64 = 1e-12,\n max_iterations::Int64 = 1000,\n force_complex::Bool = false,\n deflation_tol::Float64 = 1e-7,\n random_seed::Int64 = 42,\n)\nThe initial state \\(|\\psi_0\\rangle\\) is explicitly specified.
C++eigvals_lanczos_result_t \neigvals_lanczos(OpSum const &ops, State psi0, int64_t neigvals = 1,\n double precision = 1e-12, int64_t max_iterations = 1000,\n bool force_complex = false, double deflation_tol = 1e-7);\nNotice this version copies the initial state, which requires memory but keeps the orginal state intact.
The initial state \\(|\\psi_0\\rangle\\) is explicitly specified and overwritten in the process. This version can save memory, but the initial state \\(|\\psi_0\\rangle\\) cannot be used later.
C++eigvals_lanczos_result_t \neigvals_lanczos_inplace(OpSum const &ops, State &psi0, int64_t neigvals = 1,\n double precision = 1e-12, int64_t max_iterations = 1000,\n bool force_complex = false, double deflation_tol = 1e-7);\nA struct of type eigvals_lanczos_result_t with the following entries.
The algorithm terminates if the \\(k\\)-th (\\(k\\) is the argument neigvals) approximate eigenvalue changes only by a fraction smaller than \\(\\epsilon\\) (\\(k\\) is the argument precision), i.e.
Here, \\(\\tilde{e}_k^{(n)}\\) denotes the Lanczos approximation to the \\(k\\)-th eigenvalue after \\(n\\) iterations.
"},{"location":"documentation/algorithms/eigvals_lanczos/#usage-example","title":"Usage Example","text":"C++int N = 8;\nint nup = N / 2;\nauto block = Spinhalf(N, nup);\n\n// Define the nearest-neighbor Heisenberg model\nauto ops = OpSum();\nfor (int i=0; i<N; ++i) {\n ops += \"J\" * Op(\"SdotS\", {i, (i+1) % N});\n}\nops[\"J\"] = 1.0;\n\n// With random intial state\nauto res = eigvals_lanczos(ops, block);\nXDIAG_SHOW(res.alphas);\nXDIAG_SHOW(res.betas);\nXDIAG_SHOW(res.eigenvalues);\n\n// With specific initial state\nauto psi0 = product_state(block, {\"Up\", \"Dn\", \"Up\", \"Dn\", \"Up\", \"Dn\", \"Up\", \"Dn\"});\nauto res2 = eigvals_lanczos(ops, psi0);\nXDIAG_SHOW(res.alphas);\nXDIAG_SHOW(res.betas);\nXDIAG_SHOW(res.eigenvalues);\nComputes the exponential of a Hermitian operator \\(H\\) with an arbitrary real or complex prefactor \\(z\\) applied to a State \\(\\vert \\psi_0\\rangle\\),
\\[\\vert \\psi(z) \\rangle = e^{z(H - \\delta)} \\vert \\psi_0\\rangle.\\]Here, \\(\\delta\\) denotes a real number shifting the spectrum of \\(H\\). The algorithm implemented is described in the following publication.
On Krylov Subspace Approximations to the Matrix Exponential Operator Marlis Hochbruck and Christian Lubich SIAM Journal on Numerical Analysis, Vol. 34, Iss. 5 (1997) DOI: 10.1137/S0036142995280572
Sources evolve_lanczos.hpp, evolve_lanczos.cpp
"},{"location":"documentation/algorithms/evolve_lanczos/#definition","title":"Definition","text":"The method is provided in two variants:
Returning a new state while the input state remains untouched. This variant is safe to use and simple to code.
C++evolve_lanczos_result_t\nevolve_lanczos(OpSum const &H, State psi, double z, double precision = 1e-12,\n double shift = 0., bool normalize = false,\n int64_t max_iterations = 1000, double deflation_tol = 1e-7);\n\nevolve_lanczos_result_t\nevolve_lanczos(OpSum const &H, State psi, complex z, double precision = 1e-12,\n double shift = 0., bool normalize = false,\n int64_t max_iterations = 1000, double deflation_tol = 1e-7);\nAn inplace variant evolve_lanczos_inplace, where the input state is overwritten and contains the time evolved state upon exit. This version is more memory efficient than evolve_lanczos.
evolve_lanczos_inplace_result_t\nevolve_lanczos_inplace(OpSum const &H, State &psi, double z, \n double precision = 1e-12, double shift = 0.,\n bool normalize = false, int64_t max_iterations = 1000, \n double deflation_tol = 1e-7);\n\nevolve_lanczos_inplace_result_t\nevolve_lanczos_inplace(OpSum const &H, State &psi, complex z, \n double precision = 1e-12, double shift = 0.,\n bool normalize = false, int64_t max_iterations = 1000, \n double deflation_tol = 1e-7);\nThe parameter shift can be used to turn all eigenvalues of the matrix \\(H - \\delta \\;\\textrm{Id}\\) positive whenever \\(\\delta < E_0\\), where \\(E_0\\) denotes the ground state energy of \\(H\\).
A struct with the following entries
Entry Description alphas diagonal elements of the Lanczos tridiagonal matrix betas off-diagonal elements of the Lanczos tridiagonal matrix eigenvalues the computed Ritz eigenvalues of the tridiagonal matrix niterations number of iterations performed criterion string denoting the reason why the Lanczosalgorithm stopped state time-evolved State \\(\\vert \\psi(t)\\rangle\\) (not defined for inplace variant)"},{"location":"documentation/algorithms/evolve_lanczos/#convergence-criterion","title":"Convergence criterion","text":"The algorithm is estimating the following error,
\\[ \\varepsilon = \\parallel \\vert \\tilde{\\psi}(t)\\rangle - e^{z(H - \\delta)} \\vert \\psi_0\\rangle \\parallel_2, \\]where \\(\\vert \\tilde{\\psi}(t) \\rangle\\) denotes the approximation computed during the algorithm. As the exact solution is not available this error is estimated using the method described by Algorithm 2 in
Expokit: A Software Package for Computing Matrix Exponentials Roger B. Sidje ACM Trans. Math. Softw., 24(1):130-156, 1998. (1998) DOI: 10.1145/285861.285868
"},{"location":"documentation/algorithms/evolve_lanczos/#usage-example","title":"Usage Example","text":"C++int N = 8;\nint nup = N / 2;\nauto block = Spinhalf(N, nup);\n\n// Define the nearest-neighbor Heisenberg model\nauto ops = OpSum();\nfor (int i=0; i<N; ++i) {\n ops += Op(\"SdotS\", {i, (i+1) % N});\n}\n\n// Compute ground state energy\ndouble e0 = eigval0(ops, block);\n\nauto psi0 = product_state(block, {\"Up\", \"Dn\", \"Up\", \"Dn\", \"Up\", \"Dn\", \"Up\", \"Dn\"});\ndouble time = 1.0;\ndouble precision = 1e-12;\nauto res = evolve_lanczos(ops, psi0, time, precision, e0, true, 500);\nXDIAG_SHOW(res.alphas);\nXDIAG_SHOW(res.betas);\nComputes the imaginary-time evolution,
\\[\\vert \\psi(\\tau) \\rangle = e^{-(H - \\delta) \\tau} \\vert \\psi_0\\rangle,\\]of a State \\(\\vert \\psi_0 \\rangle\\) and a Hermitian operator \\(H\\) using an iterative algorithm. \\(\\delta\\) here denotes a real number which can be chosen as the ground state energy \\(\\delta=E_0\\) of \\(H\\).
Sources imaginary_time_evolve.hpp, imaginary_time_evolve.cpp
"},{"location":"documentation/algorithms/imaginary_time_evolve/#definition","title":"Definition","text":"The method is provided in two variants:
Returning a new state while the input state remains untouched. This variant is safe to use and simple to code.
C++State imaginary_time_evolve(OpSum const &H, State psi0, double time,\n double precision = 1e-12, double shift = 0.);\nAn inplace variant imaginary_time_evolve_inplace, where the input state is overwritten and contains the time evolved state upon exit. This version is more memory efficient than imaginary_time_evolve.
void imaginary_time_evolve_inplace(OpSum const &H, State &psi0, double time,\n double precision = 1e-12, shift = 0.);\nThe routine calls the subroutine evolve_lanczos implementing a Lanczos algorithm to perform the evolution. This routine can also be called explicitly if more control is desired. Please also confer to the page evolve_lanczos for further details on the specifics of the algorithm. The parameter shift can be used to turn all eigenvalues of the matrix \\(H - \\delta \\;\\textrm{Id}\\) positive whenever \\(\\delta < E_0\\), where \\(E_0\\) denotes the ground state energy of \\(H\\).
int N = 8;\nint nup = N / 2;\nauto block = Spinhalf(N, nup);\n\n// Define the nearest-neighbor Heisenberg model\nauto ops = OpSum();\nfor (int i=0; i<N; ++i) {\n ops += Op(\"SdotS\", {i, (i+1) % N});\n}\n\n// Compute ground state energy\ndouble e0 = eigval0(ops, block);\n\nauto psi0 = product_state(block, {\"Up\", \"Dn\", \"Up\", \"Dn\", \"Up\", \"Dn\", \"Up\", \"Dn\"});\ndouble time = 1.0;\ndouble precision = 1e-12;\nauto psi = imaginary_time_evolve(ops, psi0, time, precision, e0);\nimaginary_time_evolve_inplace(ops, psi0, time, precision, e0);\nXDIAG_SHOW(isapprox(psi0, psi));\nComputes the real-time evolution,
\\[\\vert \\psi(t) \\rangle = e^{-iHt} \\vert \\psi_0\\rangle,\\]of a State \\(\\vert \\psi_0 \\rangle\\) and a Hermitian operator \\(H\\) using an iterative algorithm.
Sources time_evolve.hpp, time_evolve.cpp
"},{"location":"documentation/algorithms/time_evolve/#definition","title":"Definition","text":"The method is provided in two variants:
Returning a new state while the input state remains untouched. This variant is safe to use and simple to code.
C++State time_evolve(OpSum const &H, State psi0, double time,\n double precision = 1e-12,\n std::string algorithm = \"lanczos\");\nAn inplace variant time_evolve_inplace, where the input state is overwritten and contains the time evolved state upon exit. This version is more memory efficient than time_evolve.
void time_evolve_inplace(OpSum const &H, State &psi0, double time,\n double precision = 1e-12,\n std::string algorithm = \"lanczos\");\nlanczos or expokit lanczos The algorithm parameter decised which backend is run. If lanczos is chosen, the evolve_lanczos routine is called with the standard arguments. Alternatively, expokit chooses the time_evolve_expokit routine. For a detailed documentation of the algorithms we refer to the evolve_lanczos and time_evolve_expokit pages. Broadly speaking, the expokit can yield higher precision states at arbitrarily long times at the cost of increased memory and computing time. In practice, we recommend analysing the effect of the precision parameters on the time evolution series obtained in both cases.
int N = 8;\nint nup = N / 2;\nauto block = Spinhalf(N, nup);\n\n// Define the nearest-neighbor Heisenberg model\nauto ops = OpSum();\nfor (int i=0; i<N; ++i) {\n ops += Op(\"SdotS\", {i, (i+1) % N});\n}\n\nauto psi0 = product_state(block, {\"Up\", \"Dn\", \"Up\", \"Dn\", \"Up\", \"Dn\", \"Up\", \"Dn\"});\ndouble time = 1.0;\nauto psi = time_evolve(ops, psi0, time);\ntime_evolve_inplace(ops, psi0, time);\nXDIAG_SHOW(isapprox(psi0, psi));\nComputes the real-time evolution,
\\[\\vert \\psi(t) \\rangle = e^{-iHt} \\vert \\psi_0\\rangle,\\]of a State \\(\\vert \\psi_0 \\rangle\\) and a Hermitian operator \\(H\\) using the iterative algorithm implemented by Expokit
Expokit: A Software Package for Computing Matrix Exponentials Roger B. Sidje ACM Trans. Math. Softw., 24(1):130-156, 1998. (1998) DOI: 10.1145/285861.285868
The algorithm features automatic stepsize control and computes approximate solutions with high precision according to our tests. Yet, the evolve_lanczos implementation is currently faster and more memory efficient.
Sources time_evolve_expokit.hpp, time_evolve_expokit.cpp
"},{"location":"documentation/algorithms/time_evolve_expokit/#definition","title":"Definition","text":"The method is provided in two variants:
Returning a new state while the input state remains untouched. This variant is safe to use and simple to code.
C++time_evolve_expokit_return_t time_evolve_expokit(\n OpSum const &ops, State state, double time, double precision = 1e-12,\n int64_t m = 30, double anorm = 0., int64_t nnorm = 2);\nAn inplace variant time_evolve_expokit_inplace, where the input state is overwritten and contains the time evolved state upon exit. This version is more memory efficient than time_evolve_expokit.
time_evolve_expokit_inplace_return_t time_evolve_expokit_inplace(\n OpSum const &ops, State &state, double time, double precision = 1e-12,\n int64_t m = 30, double anorm = 0., int64_t nnorm = 2);\nA struct with the following entries
Entry Description error the computed error estimate during evolution hump the \"hump\" as defined in Expokit 10.1145/285861.285868 state time-evolved State \\(\\vert \\psi(t)\\rangle\\) (not defined for inplace variant)"},{"location":"documentation/algorithms/time_evolve_expokit/#convergence-criterion","title":"Convergence criterion","text":"The algorithm is estimating the following error,
\\[ \\varepsilon = \\parallel \\vert \\tilde{\\psi}(t)\\rangle - e^{z(H - \\delta)} \\vert \\psi_0\\rangle \\parallel_2, \\]where \\(\\vert \\tilde{\\psi}(t) \\rangle\\) denotes the approximation computed during the algorithm. As the exact solution is not available this error is estimated using the method described by Algorithm 2 in
Expokit: A Software Package for Computing Matrix Exponentials Roger B. Sidje ACM Trans. Math. Softw., 24(1):130-156, 1998. (1998) DOI: 10.1145/285861.285868
"},{"location":"documentation/algorithms/time_evolve_expokit/#usage-example","title":"Usage Example","text":"C++int N = 8;\nint nup = N / 2;\nauto block = Spinhalf(N, nup);\n\n// Define the nearest-neighbor Heisenberg model\nauto ops = OpSum();\nfor (int i=0; i<N; ++i) {\n ops += Op(\"SdotS\", {i, (i+1) % N});\n}\n\nauto psi0 = product_state(block, {\"Up\", \"Dn\", \"Up\", \"Dn\", \"Up\", \"Dn\", \"Up\", \"Dn\"});\ndouble time = 1.0;\nauto res1 = time_evolve_expokit(ops, psi0, time, 20);\nauto res2 = time_evolve_expokit_inplace(ops, psi0, time, 20);\nXDIAG_SHOW(isapprox(psi0, res1.state));\nXDIAG_SHOW(res1.error);\nXDIAG_SHOW(res1.hump);\nA block in an Electron (fermions with \\(\\uparrow, \\downarrow\\) spin) Hilbert space.
Sources electron.hpp, electron.cpp
"},{"location":"documentation/blocks/electron/#constructors","title":"Constructors","text":"C++JuliaElectron(int64_t nsites, std::string backend = \"auto\");\nElectron(int64_t nsites, int64_t nup, int64_t ndn, std::string backend = \"auto\");\nElectron(int64_t nsites, Representation irrep, std::string backend = \"auto\");\nElectron(int64_t nsites, int64_t nup, int64_t ndn, Representation irrep, std::string backend = \"auto\");\nElectron(nsites::Integer, backend::String)\nElectron(nsites::Integer, nup::Integer, ndn::Integer, backend::String)\nElectron(nsites::Integer, irrep::Representation, backend::String)\nElectron(nsites::Integer, nup::Integer, ndn::Integer, irrep::Representation, backend::String)\nauto The parameter backend chooses how the block is coded internally. By using the default parameter auto the backend is chosen automatically. Alternatives are 32bit, 64bit
An Electron block can be iterated over, where at each iteration a ProductState representing the corresponding basis state is returned.
C++Juliaauto block = Electron(4, 2, 2);\nfor (auto pstate : block) {\n Log(\"{} {}\", to_string(pstate), index(block, pstate));\n}\nblock = Electron(4, 2, 2)\nfor pstate in block\n @show pstate, index(block, pstate) \nend\nReturns the index of a given ProductState in the basis of the tJ block.
C++Juliaint64_t index(tJ const &block, ProductState const &pstate);\nindex(block::tJ, pstate::ProductState)\n1-indexing
In the C++ version, the index count starts from \"0\" whereas in Julia the index count starts from \"1\".
"},{"location":"documentation/blocks/electron/#nsites","title":"nsites","text":"Returns the number of sites of the block.
C++Juliaint64_t nsites(tJ const &block);\nnsites(block::tJ)\nReturns the size of the block, i.e. its dimension.
C++Juliaint64_t size(tJ const &block) const;\nsize(block::tJ)\nReturns the dimension of the block, same as \"size\" for non-distributed blocks.
C++Juliaint64_t dim(tJ const &block) const;\ndim(block::tJ)\nReturns whether the block can be used with real arithmetic. Complex arithmetic is needed when a Representation is genuinely complex.
C++Juliabool isreal(tJ const &block);\nisreal(block::tJ)\nint N = 4;\nint nup = 2;\nint ndn = 1;\n\n// without number conservation\nauto block = Electron(N);\nXDIAG_SHOW(block);\n\n// with number conservation\nauto block_np = Electron(N, nup, ndn);\nXDIAG_SHOW(block_np);\n\n// with symmetries, without number conservation\nauto p1 = Permutation({0, 1, 2, 3});\nauto p2 = Permutation({1, 2, 3, 0});\nauto p3 = Permutation({2, 3, 0, 1});\nauto p4 = Permutation({3, 0, 1, 2});\nauto group = PermutationGroup({p1, p2, p3, p4});\nauto irrep = Representation(group, arma::vec{1, -1, 1, -1});\nauto block_sym = Electron(N, irrep);\nXDIAG_SHOW(block_sym);\n\n// with symmetries and number conservation\nauto block_sym_np = Electron(N, nup, ndn, irrep);\nXDIAG_SHOW(block_sym_np);\nXDIAG_SHOW(block_sym_np.nsites());\nXDIAG_SHOW(block_sym_np.size());\n\n// Iteration\nfor (auto pstate : block_sym_np) {\n Log(\"{} {}\", to_string(pstate), index(block_sym_np, pstate));\n}\nN = 4\nnup = 2\nndn = 1\n\n# without number conservation\nblock = Electron(N)\n@show block\n\n# with number conservation\nblock_np = Electron(N, nup, ndn)\n@show block_np\n\n# with symmetries, without number conservation\np1 = Permutation([1, 2, 3, 4])\np2 = Permutation([2, 3, 4, 1])\np3 = Permutation([3, 4, 1, 2])\np4 = Permutation([4, 1, 2, 3])\ngroup = PermutationGroup([p1, p2, p3, p4])\nrep = Representation([1, -1, 1, -1])\nblock_sym = Electron(N, group, rep)\n@show block_sym\n\n# with symmetries and number conservation\nblock_sym_np = Electron(N, nup, ndn, group, rep)\n@show block_sym_np\n\n@show nsites(block_sym_np)\n@show size(block_sym_np)\n\n# Iteration\nfor pstate in block_sym_np\n @show pstate, index(block_sym_np, pstate)\nend\n@show permutation_group(block_sym_np)\n@show irrep(block_sym_np)\nA block in a spin \\(S=1/2\\) Hilbert space.
Sources spinhalf.hpp, spinhalf.cpp
"},{"location":"documentation/blocks/spinhalf/#constructors","title":"Constructors","text":"C++JuliaSpinhalf(int64_t nsites, std::string backend = \"auto\");\nSpinhalf(int64_t nsites, int64_t nup, std::string backend = \"auto\");\nSpinhalf(int64_t nsites, Representation const &irrep, std::string backend = \"auto\");\nSpinhalf(int64_t nsites, int64_t nup, Representation const &irrep, std::string backend = \"auto\");\nSpinhalf(nsites::Integer, backend::String)\nSpinhalf(nsites::Integer, nup::Integer, backend::String)\nSpinhalf(nsites::Integer, irrep::Representation, backend::String)\nSpinhalf(nsites::Integer, nup::Integer, irrep::Representation, backend::String)\nauto The parameter backend chooses how the block is coded internally. By using the default parameter auto the backend is chosen automatically. Alternatives are 32bit, 64bit, 1sublattice, 2sublattice, 3sublattice, 4sublattice, and 5sublattice. The backends xsublattice implement the sublattice coding algorithm described in Wietek, L\u00e4uchli, Phys. Rev. E 98, 033309 (2018). The sublattice coding algorithms impose certain constraints on the symmetries used, as described in the reference.
An Spinhalf block can be iterated over, where at each iteration a ProductState representing the corresponding basis state is returned.
C++Juliaauto block = Spinhalf(4, 2);\nfor (auto pstate : block) {\n Log(\"{} {}\", to_string(pstate), block.index(pstate));\n}\nblock = Spinhalf(4, 2)\nfor pstate in block\n @show pstate, index(block, pstate) \nend\nReturns the index of a given ProductState in the basis of the Spinhalf block.
C++Juliaint64_t index(Spinhalf const &block, ProductState const &pstate);\nindex(block::Spinhalf, pstate::ProductState)\n1-indexing
In the C++ version, the index count starts from \"0\" whereas in Julia the index count starts from \"1\".
"},{"location":"documentation/blocks/spinhalf/#nsites","title":"nsites","text":"Returns the number of sites of the block.
C++Juliaint64_t nsites(Spinhalf const &block);\nnsites(block::Spinhalf)\nReturns the size of the block, i.e. its dimension.
C++Juliaint64_t size(Spinhalf const &block) const;\nsize(block::Spinhalf)\nReturns the dimension of the block, same as \"size\" for non-distributed blocks.
C++Juliaint64_t dim(Spinhalf const &block) const;\ndim(block::Spinhalf)\nReturns whether the block can be used with real arithmetic. Complex arithmetic is needed when a Representation is genuinely complex.
C++Juliabool isreal(Spinhalf const &block);\nisreal(block::Spinhalf)\nint N = 4;\nint nup = 2;\n\n// without Sz conservation\nauto block = Spinhalf(N);\nXDIAG_SHOW(block);\n\n// with Sz conservation\nauto block_sz = Spinhalf(N, nup);\nXDIAG_SHOW(block_sz);\n\n// with symmetries, without Sz\nPermutation p1 = {0, 1, 2, 3};\nPermutation p2 = {1, 2, 3, 0};\nPermutation p3 = {2, 3, 0, 1};\nPermutation p4 = {3, 0, 1, 2};\nauto group = PermutationGroup({p1, p2, p3, p4});\nauto irrep = Representation(group, arma::vec({1, -1, 1, -1}));\nauto block_sym = Spinhalf(N, irrep);\nXDIAG_SHOW(block_sym);\n\n// with symmetries and Sz\nauto block_sym_sz = Spinhalf(N, nup, irrep);\nXDIAG_SHOW(block_sym_sz);\n\nXDIAG_SHOW(block_sym_sz.nsites());\nXDIAG_SHOW(block_sym_sz.size());\n\n// Iteration\nfor (auto pstate : block_sym_sz) {\n Log(\"{} {}\", to_string(pstate), index(block_sym_sz, pstate));\n}\nN = 4\nnup = 2\n\n# without Sz conservation\nblock = Spinhalf(N)\n@show block\n\n\n# with Sz conservation\nblock_sz = Spinhalf(N, nup)\n@show block_sz\n\n# with symmetries, without Sz\np1 = Permutation([1, 2, 3, 4])\np2 = Permutation([2, 3, 4, 1])\np3 = Permutation([3, 4, 1, 2])\np4 = Permutation([4, 1, 2, 3])\ngroup = PermutationGroup([p1, p2, p3, p4])\nrep = Representation([1, -1, 1, -1])\nblock_sym = Spinhalf(N, group, rep)\n@show block_sym\n\n# with symmetries and Sz\nblock_sym_sz = Spinhalf(N, nup, group, rep)\n@show block_sym_sz\n\n@show nsites(block_sym_sz)\n@show size(block_sym_sz)\n\n# Iteration\nfor pstate in block_sym_sz\n @show pstate, index(block_sym_sz, pstate)\nend\n@show permutation_group(block_sym_sz)\n@show irrep(block_sym_sz)\nA block in a spin \\(S=1/2\\) Hilbert space with distributed computing capabilities.
Sources spinhalf_distributed.hpp, spinhalf_distributed.cpp
"},{"location":"documentation/blocks/spinhalf_distributed/#constructors","title":"Constructors","text":"C++SpinhalfDistributed(int64_t nsites, int64_t nup, std::string backend = \"auto\");\nauto The parameter backend chooses how the block is coded internally. By using the default parameter auto the backend is chosen automatically. Alternatives are 32bit, 64bit.
An SpinhalfDistributed block can be iterated over, where at each iteration a ProductState representing the corresponding basis state is returned.
C++auto block = SpinhalfDistributed(4, 2);\nfor (auto pstate : block) {\n Log(\"{} {}\", to_string(pstate), block.index(pstate));\n}\nReturns the index of a given ProductState in the basis of the SpinhalfDistributed block.
C++int64_t index(SpinhalfDistributed const &block, ProductState const &pstate);\nReturns the number of sites of the block.
C++int64_t nsites(SpinhalfDistributed const &block);\nReturns the size of the block on a local process.
C++int64_t size(SpinhalfDistributed const &block) const;\nReturns the dimension of the block, i.e. the sum of all sizes across all processes.
C++int64_t dim(SpinhalfDistributed const &block) const;\nReturns whether the block can be used with real arithmetic. Complex arithmetic is needed when a Representation is genuinely complex.
C++bool isreal(SpinhalfDistributed const &block);\nA block in a \\(t-J\\) type Hilbert space, i.e. fermions with \\(\\uparrow, \\downarrow\\) spin excluding doubly occupied sites.
Sources tj.hpp, tj.cpp
"},{"location":"documentation/blocks/tJ/#constructors","title":"Constructors","text":"C++JuliatJ(int64_t nsites, int64_t nup, int64_t ndn, std::string backend = \"auto\");\ntJ(int64_t nsites, int64_t nup, int64_t ndn, Representation const &irrep, std::string backend = \"auto\");\ntJ(nsites::Integer, nup::Integer, ndn::Integer, backend::String)\ntJ(nsites::Integer, nup::Integer, ndn::Integer, irrep::Representation, backend::String)\nauto The parameter backend chooses how the block is coded internally. By using the default parameter auto the backend is chosen automatically. Alternatives are 32bit, 64bit.
An tJ block can be iterated over, where at each iteration a ProductState representing the corresponding basis state is returned.
C++Juliaauto block = tJ(4, 2, 1);\nfor (auto pstate : block) {\n Log(\"{} {}\", to_string(pstate), index(block, pstate));\n}\nblock = tJ(4, 2, 1)\nfor pstate in block\n @show pstate, index(block, pstate) \nend\nReturns the index of a given ProductState in the basis of the tJ block.
C++Juliaint64_t index(tJ const &block, ProductState const &pstate);\nindex(block::tJ, pstate::ProductState)\n1-indexing
In the C++ version, the index count starts from \"0\" whereas in Julia the index count starts from \"1\".
"},{"location":"documentation/blocks/tJ/#nsites","title":"nsites","text":"Returns the number of sites of the block.
C++Juliaint64_t nsites(tJ const &block);\nnsites(block::tJ)\nReturns the size of the block, i.e. its dimension.
C++Juliaint64_t size(tJ const &block) const;\nsize(block::tJ)\nReturns the dimension of the block, same as \"size\" for non-distributed blocks.
C++Juliaint64_t dim(tJ const &block) const;\ndim(block::tJ)\nReturns whether the block can be used with real arithmetic. Complex arithmetic is needed when a Representation is genuinely complex.
C++Juliabool isreal(tJ const &block);\nisreal(block::tJ)\nint N = 4;\nint nup = 2;\nint ndn = 1;\n\n// without permutation symmetries\nauto block = tJ(N, nup, ndn);\nXDIAG_SHOW(block);\n\n// with permutation symmetries\nauto p1 = Permutation({0, 1, 2, 3});\nauto p2 = Permutation({1, 2, 3, 0});\nauto p3 = Permutation({2, 3, 0, 1});\nauto p4 = Permutation({3, 0, 1, 2});\nauto group = PermutationGroup({p1, p2, p3, p4});\nauto irrep = Representation(group, arma::vec{1, -1, 1, -1});\nauto block_sym = tJ(N, nup, ndn, irrep);\nXDIAG_SHOW(block_sym);\nXDIAG_SHOW(block_sym.nsites());\nXDIAG_SHOW(block_sym.size());\n\n// Iteration\nfor (auto pstate : block_sym) {\n Log(\"{} {}\", to_string(pstate), index(block_sym, pstate));\n}\nN = 4\nnup = 2\nndn = 1\n\n# without permutation symmetries\nblock = tJ(N, nup, ndn)\n@show block\n\n# with permutation symmetries\np1 = Permutation([1, 2, 3, 4])\np2 = Permutation([2, 3, 4, 1])\np3 = Permutation([3, 4, 1, 2])\np4 = Permutation([4, 1, 2, 3])\ngroup = PermutationGroup([p1, p2, p3, p4])\nrep = Representation([1, -1, 1, -1])\nblock_sym = tJ(N, nup, ndn, group, rep)\n@show block_sym\n\n@show nsites(block_sym)\n@show size(block_sym)\n\n# Iteration\nfor pstate in block_sym\n @show pstate, index(block_sym, pstate)\nend\n@show permutation_group(block_sym)\n@show irrep(block_sym)\nA block in a \\(t-J\\) type Hilbert space, i.e. fermions with \\(\\uparrow, \\downarrow\\) spin excluding doubly occupied sites with distributed computing capabilities.
Sources tj_distributed.hpp, tj_distributed.cpp
"},{"location":"documentation/blocks/tJ_distributed/#constructors","title":"Constructors","text":"C++tJDistributed(int64_t nsites, int64_t nup, int64_t ndn, std::string backend = \"auto\");\nauto The parameter backend chooses how the block is coded internally. By using the default parameter auto the backend is chosen automatically. Alternatives are 32bit, 64bit.
An tJDistributed block can be iterated over, where at each iteration a ProductState representing the corresponding basis state is returned.
C++auto block = tJDistributed(4, 2, 1);\nfor (auto pstate : block) {\n Log(\"{} {}\", to_string(pstate), index(block, pstate));\n}\nReturns the index of a given ProductState in the basis of the tJDistributed block.
C++int64_t index(tJDistributed const &block, ProductState const &pstate);\nReturns the number of sites of the block.
C++int64_t nsites(tJDistributed const &block);\nReturns the size of the block on a local process.
C++int64_t size(tJDistributed const &block) const;\nReturns the dimension of the block, i.e. the sum of all sizes across all processes.
C++int64_t dim(tJDistributed const &block) const;\nReturns whether the block can be used with real arithmetic. Complex arithmetic is needed when a Representation is genuinely complex.
C++bool isreal(tJDistributed const &block);\nDownload the source code using git
cd /path/to/where/xdiag/should/be\ngit clone https://github.com/awietek/xdiag.git\nCompile the default library
cd xdiag\ncmake -S . -B build\ncmake --build build\ncmake --install build\ninstall. Compile the distributed library
To use the distributed computing features of xdiag, the distributed library has to be built which requires MPI.
cd xdiag\ncmake -S . -B build -D XDIAG_DISTRIBUTED=On\ncmake --build build\ncmake --install build\nInfo
It might be necessary to explicitly define MPI compiler, e.g. mpicxx like this
cmake -S . -B build -D XDIAG_DISTRIBUTED=On -D CMAKE_CXX_COMPILER=mpicxx\nParallel compilation To speed up the compilation process, the build step can be performed in parallel using the -j flag
cmake --build build -j\nListing compile options
The available compilation options can be displayed using
cmake -L .\nShared / Static libraries
XDiag can be built either as a static or shared library. By default, a static library is built. To build a shared library, use the option XDIAG_SHARED_LIBS, e.g.
cmake -S . -B build -D XDIAG_SHARED_LIBS=On\nChoosing a certain compiler
The compiler (e.g. icpx) can be specified using
cmake -S . -B build -D CMAKE_CXX_COMPILER=icpx\nWarning
If the xdiag library is compiled with a certain compiler, it is advisable to also compile the application codes with the same compiler.
Setting the install path
In the installation step, the install directory can be set in the following way
cmake --install build --prefix /my/install/prefix\nDisabling HDF5/OpenMP
To disable support for HDF5 or OpenMP support, use
cmake -S . -B build -D XDIAG_DISABLE_OPENMP=On -D XDIAG_DISABLE_HDF5=On\nBuilding and running tests
To compile and run the testing programs, use
cmake -S . -B build -D BUILD_TESTING=On\ncmake --build build\nbuild/tests/tests\nmissing documentation
"},{"location":"documentation/compilation/building_documentation/","title":"Building Documentation","text":"The source files for the documentation can be found in the directory docs. The documentation is built using Material for MKDocs. To install it please run the following commands,
pip install mkdocs\npip install mkdocs-material\nTo work on the documentation locally, it can be served using
mkdocs serve\nfrom the xdiag root source directory. A local build of the documentation can then be accessed in a webbrowser at the adress
127.0.0.1:8000\nIn order to develop and extend the julia wrapper, one should work locally and build a local version of the xdiag Julia binaries. First, get the path to the CxxWrap package of julia. To do so, enter the Julia REPL,
julia\nusing CxxWrap\nCxxWrap.prefix_path()\n/path/to/libcxxwrap-julia-prefix. This is then used to configure the cmake compilation. cmake -S . -B build -D XDIAG_JULIA_WRAPPER=On -D CMAKE_PREFIX_PATH=/path/to/libcxxwrap-julia-prefix\ncmake --build build\ncmake --install build\nlibxdiagjl.so, (or the corresponding library format on non-Linux systems)."},{"location":"documentation/io/file_h5/","title":"FileH5","text":"A file handler for hdf5 files. The proper tool to write results of XDiag simulations to disk. Only provided for the C++ version as Julia already provides good tools handling hdf5 files.
Sources file_h5.hpp, file_h5.cpp
"},{"location":"documentation/io/file_h5/#constructors","title":"Constructors","text":"C++FileH5(std::string filename, std::string iomode = \"r\");\nr There are four possible values of iomode
r: read-only modew: secure write mode, new file is created if it does not existw!: forced write mode, existing file will be overwrittena: append mode, new datasets in an existing file can be createdReturns a handler to a value to be read or written from or to the hdf5 file.
C++hdf5::FileH5Handler operator[](std::string key);\nstd::string filename = XDIAG_DIRECTORY \"/misc/data/hdf5/write.h5\";\nauto fl = FileH5(filename, \"w!\");\n\n// Write output to the hdf5 file\nfl[\"val\"] = 12;\nfl[\"test/to\"] = 22;\nfl[\"test/to2/group\"] = 32;\nfl[\"test/to3/group2/asdf\"] = 42;\n\nauto mat = arma::cx_mat(3, 5, arma::fill::randn);\nfl[\"a/b/c/mat\"] = mat;\nA file handler for TOML files. Most funtionality is only provided for the C++ version as Julia already provides good tools handling TOML files.
Sources file_toml.hpp, file_toml.cpp
"},{"location":"documentation/io/file_toml/#constructors","title":"Constructors","text":"C++JuliaFileToml(const char *filename);\nFileToml(std::string filename);\nFileToml(filename::String)\nReturns whether or not the TOML file has a certain key defined.
C++bool defined(FileToml const &fl, std::string key);\nReturns a handler to a value to be read from the TOML file.
C++io::FileTomlHandler operator[](std::string key);\nauto fl = FileToml(XDIAG_DIRECTORY \"/misc/data/toml/input.toml\");\nXDIAG_SHOW(defined(fl, \"N\"));\n\nint N = fl[\"N\"].as<int>();\nint nup = fl[\"nup\"].as<int>();\ndouble J1 = fl[\"J1\"].as<double>();\ndouble J2 = fl[\"J2\"].as<double>();\n\nauto block = Spinhalf(N, nup);\nauto H = OpSum();\nfor (int i=0; i<N; ++i){\n H += J1 * Op(\"SdotS\", {i, (i+1)%N});\n H += J2 * Op(\"SdotS\", {i, (i+2)%N});\n}\ndouble e0 = eigval0(H, block);\nXDIAG_SHOW(e0);\nReads an OpSum object from a TOML file.
Sources read.hpp, read.cpp
"},{"location":"documentation/io/read_opsum/#definition","title":"Definition","text":"C++JuliaOpSum read_opsum(FileToml file, std::string tag);\nread_opsum(file::FileToml, tag::String);\nAn OpSum can be defined in a TOML file as a simple list. The entries of the list are themselves also lists, which contain two or more entries:
Here is typical example of a OpSum specification in a TOML file:
Interactions = [\n ['J1', 'SdotS', 0, 1],\n ['J1', 'SdotS', 1, 2],\n ['J1', 'SdotS', 2, 0],\n ['Jchi', 'ScalarChirality', 0, 1, 2],\n]\nstd::string file = XDIAG_DIRECTORY \"/misc/data/triangular.9.hop.sublattices.tsl.toml\";\nauto fl = FileToml(file);\nauto ops = read_opsum(fl, \"Interactions\");\nXDIAG_SHOW(ops);\n\nReads an PermutationGroup object from a TOML file.
Sources read.hpp, read.cpp
"},{"location":"documentation/io/read_permutation_group/#definition","title":"Definition","text":"C++JuliaOpSum read_permutation_group(FileToml file, std::string tag);\nread_permutation_group(file::FileToml, tag::String);\nA PermutationGroup can be defined in a TOML file as an integer matrix, where the rows are the integers of the permutation.
A typical specification of a \\(C_4\\) PermutationGroup is shown here:
Symmetries = [\n [0, 1, 2, 3],\n [1, 2, 3, 0],\n [2, 3, 0, 1],\n [3, 0, 1, 2]\n]\n1-indexing in Julia / 0-indexing in C++
To enumerate the sites of a Permutation, we start counting at 1 in Julia and 0 in C++.
"},{"location":"documentation/io/read_permutation_group/#usage-example","title":"Usage Example","text":"C++Juliastd::string file = XDIAG_DIRECTORY \"/misc/data/triangular.9.hop.sublattices.tsl.toml\";\nauto fl = FileToml(file);\nauto group = read_permutation_group(fl, \"Symmetries\");\nXDIAG_SHOW(group);\n\nReads an Representation object from a TOML file.
Sources read.hpp, read.cpp
"},{"location":"documentation/io/read_representation/#definition","title":"Definition","text":"C++JuliaOpSum read_representation(FileToml file, std::string irrep_tag, \n std::string group_tag = \"Symmetries\");\nread_representation(file::FileToml, irrep_tag::String, \n group_tag::String = \"Symmetries\");\nA Representation can be defined in a TOML file by specifying up two or three things:
characters of the representation allowed_symmetries of the representation, i.e. a list of the number of symmetries used in the irrep. By default all symmetries of the group are used in the representation.A typical specification of several Representations is shown here:
Symmetries = [\n [0, 1, 2, 3],\n [1, 2, 3, 0],\n [2, 3, 0, 1],\n [3, 0, 1, 2]\n ]\n\n # real irrep, momentum 0\n [k_0]\n characters=[1.0, 1.0, 1.0, 1.0]\n\n # complex irrep, momentum pi/2\n [k_pi2]\n characters=[\n [1.0, 0.0], \n [0.0, 1.0], \n [-1.0, 0.0],\n [0.0, -1.0]\n ]\n\n # real irrep, momentum pi\n [k_pi]\n characters=[1.0, -1.0, 1.0, -1.0]\n\n # real irrep on subgroup\n [k_pi2_half]\n characters=[1.0, -1.0]\n allowed_symmetries=[0, 2]\n1-indexing in Julia / 0-indexing in C++
To enumerate the sites of a Permutation, we start counting at 1 in Julia and 0 in C++.
"},{"location":"documentation/io/read_representation/#usage-example","title":"Usage Example","text":"The example reads the representation defined in the irreps.toml file, whose contents are shown up in the section Data format.
std::string file = XDIAG_DIRECTORY \"/misc/data/irreps.toml\";\nauto fl = FileToml(file);\n\nauto k_0 = read_representation(fl, \"k_0\");\nXDIAG_SHOW(k_0);\nXDIAG_SHOW(isreal(k_0));\n\nauto k_pi2 = read_representation(fl, \"k_pi2\");\nXDIAG_SHOW(k_pi2);\nXDIAG_SHOW(isreal(k_pi2));\n\nauto k_pi = read_representation(fl, \"k_pi\");\nXDIAG_SHOW(k_pi);\nXDIAG_SHOW(isreal(k_pi));\n\nauto k_pi2_half = read_representation(fl, \"k_pi2_half\");\nXDIAG_SHOW(k_pi2_half);\nXDIAG_SHOW(isreal(k_pi2_half));\n\nDescribes the coupling of an operator. A coupling can either be a string or a Scalar, which is either a real or complex double precision floating point number.
Sources coupling.hpp, coupling.cpp
"},{"location":"documentation/operators/coupling/#constructors","title":"Constructors","text":"C++Coupling(std::string value);\nCoupling(double value);\nCoupling(complex value);\nCoupling(Scalar value););\nReturns whether or not the coupling is a Scalar, i.e. a real or complex number.
C++bool isscalar(Coupling const& c)\nReturns whether or not the coupling is a string.
C++bool isstring(Coupling const& c)\nReturns the Scalar if the Coupling holds a scalar. Otherwise throws an xdiag::Error
Scalar scalar(Coupling const& c)\nReturns the string if the Coupling holds a string. Otherwise throws an xdiag::Error
std::string string(Coupling const& c)\nReturns the hermitian conjugate \\(\\mathcal{O}^\\dagger\\) of an operator \\(\\mathcal{O}\\) represented by an Op or OpSum object. Please note the details when conjugating complex couplings, outlined in OpSum # Complex couplings.
Sources hc.hpp, hc.cpp
"},{"location":"documentation/operators/hc/#definition","title":"Definition","text":"C++Op hc(Op const &op)\nOpSum hc(OpSum const &ops)\nauto cdagup = Op(\"Cdagup\", 0);\nauto sdots = Op(\"SdotS\", {0, 1});\nauto hop = (1.0 + 1.0i) * Op(\"Hop\", {0, 1});\nXDIAG_SHOW(cdagup == hc(cdagup));\nXDIAG_SHOW(sdots == hc(sdots));\nXDIAG_SHOW(hop == hc(hop));\nObject describing a single linear operator acting on a Hilbert space.
Sources op.hpp, op.cpp
Every operator is defined by up to three paramaters:
The type of the operator. This is a string argument which determines what kind of operator is represented. A detailed overview of the available types can be found at Operator types
The sites of the operator. This defines which physical sites (or orbitals) the operator acts upon. While most operator types require sites, there are also operator types (e.g. HubbardU) which do not need to define sites.
For special interactions, it can be necessary to additionally specify a numerical matrix, which can either be real or complex. An example is the operator type Matrix defining generic spin interactions.
Op(std::string type);\nOp(std::string type, int64_t site);\nOp(std::string type, std::vector<int64_t> const &sites);\nOp(std::string type, int64_t site, arma::mat const &matrix);\nOp(std::string type, std::vector<int64_t> const &sites, arma::mat const &matrix);\nOp(std::string type, int64_t site, arma::cx_mat const &matrix);\nOp(std::string type, std::vector<int64_t> const &sites, arma::cx_mat const &matrix);\nOp(type::String)\nOp(type::String, site::Int64)\nOp(type::String, sites::Vector{Int64})\nOp(type::String, site::Int64, matrix::Matrix{Float64})\nOp(type::String, sites::Vector{Int64}, matrix::Matrix{Float64})\nOp(type::String, site::Int64, matrix::Matrix{ComplexF64})\nOp(type::String, sites::Vector{Int64}, matrix::Matrix{ComplexF64})\n1-indexing in Julia / 0-indexing in C++
To enumerate the sites of an Op, we start counting at 1 in Julia and 0 in C++.
"},{"location":"documentation/operators/op/#methods","title":"Methods","text":""},{"location":"documentation/operators/op/#isreal","title":"isreal","text":"Returns whether an Op is a real operator.
C++Juliabool isreal(Op const &op);\nisreal(op::Op)\nauto op = \"T\" * Op(\"Hop\", {0, 1});\nXDIAG_SHOW(op);\n\nop = 1.23 * Op(\"Hop\", {0, 1});\nXDIAG_SHOW(op);\n\narma::cx_mat m(arma::mat(\"0 0; 0 0\"), arma::mat(\"0 -1; 1 0\"));\n op = Op(\"Matrix\", 0, m);\nXDIAG_SHOW(op);\nXDIAG_SHOW(isreal(op));\nop = \"T\" * Op(\"Hop\", [1, 2])\n@show op\n\nop = 1.23 * Op(\"Hop\", [1, 2])\n@show op\n\nop = Op(\"Matrix\", 1, [0 -im; im 0])\n@show op\n@show isreal(op)\nGeneric operators in XDiag are represented as OpSum objects made up of Coupling and Op objects. Every Op is defined by a type. Here we list all the available types implemented in XDiag, their required number of sites, and the blocks for which they are available.
Hop A hopping term for \\(\\uparrow\\) and \\(\\downarrow\\) spins of the form $$ \\textcolor{red}{-}\\sum_{\\sigma=\\uparrow\\downarrow} (c^\\dagger_{i\\sigma}c_{j\\sigma} + \\textrm{h.c.})$$ 2 tJ, Electron, tJDistributed Hopup A hopping term for \\(\\uparrow\\) spins of the form $$ \\textcolor{red}{-}(c^\\dagger_{i\\uparrow}c_{j\\uparrow} + \\textrm{h.c.})$$ 2 tJ, Electron, tJDistributed Hopdn A hopping term for \\(\\downarrow\\) spins of the form $$ \\textcolor{red}{-}(c^\\dagger_{i\\downarrow}c_{j\\downarrow} + \\textrm{h.c.})$$ 2 tJ, Electron, tJDistributed HubbardU A uniform Hubbard interaction across the full lattice of the form $$ \\sum_i n_{i\\uparrow}n_{i\\downarrow}$$ 0 Electron Cdagup A fermionic creation operator for an \\(\\uparrow\\) spin \\(c^\\dagger_{i\\uparrow}\\) 1 tJ, Electron, tJDistributed Cdagdn A fermionic creation operator for an \\(\\downarrow\\) spin \\(c^\\dagger_{i\\downarrow}\\) 1 tJ, Electron, tJDistributed Cup A fermionic annihilation operator for an \\(\\uparrow\\) spin \\(c_{i\\uparrow}\\) 1 tJ, Electron, tJDistributed Cdn A fermionic annihilation operator for an \\(\\downarrow\\) spin \\(c_{i\\downarrow}\\) 1 tJ, Electron, tJDistributed Nup A number operator for an \\(\\uparrow\\) spin \\(n_{i\\uparrow}\\) 1 tJ, Electron, tJDistributed Ndn A number operator for an \\(\\downarrow\\) spin \\(n_{i\\downarrow}\\) 1 tJ, Electron, tJDistributed Ntot A number operator \\(n_i = n_{i\\uparrow} + n_{i\\downarrow}\\) 1 tJ, Electron, tJDistributed NtotNtot A density-density interaction \\(n_i n_j\\) 2 tJ, Electron, tJDistributed SdotS A Heisenberg interaction of the form $$ \\mathbf{S}_i \\cdot \\mathbf{S}_j = S^x_iS^x_j + S^y_iS^y_j + S^z_iS^z_j$$ 2 Spinhalf, tJ, Electron, SpinhalfDistributed, tJDistributed SzSz An Ising interaction of the form $ S^z_i S^z_j $ 2 Spinhalf, tJ, Electron, SpinhalfDistributed, tJDistributed Exchange A spin exchange interaction of the form $$ \\frac{1}{2}(S^+_i S^-_j + S^-_iS^+_j) = S^x_iS^x_j + S^y_iS^y_j j$$ 2 Spinhalf, tJ, Electron, SpinhalfDistributed, tJDistributed Sz A local magnetic moment in the \\(z\\)-direction $ S^z_i$ 1 Spinhalf, tJ, Electron, SpinhalfDistributed, tJDistributed S+ A local spin raising operator \\(S^+_i\\) 1 Spinhalf, SpinhalfDistributed S- A local spin lowering operator \\(S^-_i\\) 1 Spinhalf, SpinhalfDistributed ScalarChirality A scalar chirality interaction of the form $$ \\mathbf{S}_i \\cdot ( \\mathbf{S}_j \\times \\mathbf{S}_k)$$ 3 Spinhalf tJSzSz An Ising interaction as encountered in the \\(t-J\\) model of the form $$ S^z_i S^z_j - \\frac{n_i n_j}{4}$$ 2 tJ, tJDistributed tJSdotS An Heisenberg interaction as encountered in the \\(t-J\\) model of the form $$ \\mathbf{S}_i \\cdot \\mathbf{S}_j - \\frac{n_i n_j}{4}$$ 2 tJ, tJDistributed Matrix A generic spin interaction no an arbitrary number of sites defined via a coupling matrix arbitrary Spinhalf"},{"location":"documentation/operators/operator_types/#matrix-type","title":"Matrix type","text":"The Matrix interaction type is a special type with whom one can define generic interactions for the Spinhalf block. In addition to the type and sites argument, also a numerical matrix is provided when constructing the Op object. The matrix describes the operator acting on the \\(2^n\\) dimensional space spanned by the \\(n\\) sites of the operator. For example, we can represent a \\(S^x\\) spin operator as,
auto sx = arma::mat({{0, 1},{1, 0}});\nauto op = Op(\"Matrix\", 0, sx);\nMore generically, we can use this mechanism to construct arbitary spin interactions, e.g.
C++auto sx = arma::mat({{0, 1},{1, 0}});\nauto sz = arma::mat({{0.5, 1},{0, -0.5}});\n\narma::mat sxsz = arma::kron(sx, sz);\narma::mat sxszsxsz = arma::kron(sxsz, sxsz);\n\nauto op_sxsz = Op(\"Matrix\", {0, 1}, sxsz);\nauto op_sxszsxsz = Op(\"Matrix\", {0, 1, 2, 3}, sxsz);\nHere we have been using the Kronecker product function kron.
Object representing a generic many-body operator by a sum of operators of the form
\\[ \\mathcal{O} = \\sum_i c_i \\mathcal{O}_i. \\]Sources opsum.hpp, opsum.cpp
An OpSum is consists of a sum of pairs given by
A coupling constant \\(c_i\\) which is defined by a Coupling object. The coupling can either be a string name or a real/complex number.
An operator \\(\\mathcal{O}_i\\) defined by an Op object.
Generically, an OpSum can thus have coupling constants defined by either strings or numerical real/complex numbers. We call an OpSum plain if its couplings are only numerical numbers, and not strings. String couplings can be defined by using the access operator[]. If all string coupling constants are defined, the OpSum can be converted to a plain OpSum using the plain method shown below.
Thus, OpSums can be defined independently of the numerical values of the coupling constants, e.g. in an input file. Upon execution of the code, these constants can then be set. Most operations in XDiag require the OpSum to be convertible to a plain OpSum.
OpSums can be added and subtracted, as well as multiplied with and divided by a Scalar. Hence, OpSums carry the mathematical structure of a vector space.
"},{"location":"documentation/operators/opsum/#constructors","title":"Constructors","text":"The following constructors create an OpSum with a single pair of coupling and operator. Additional terms can be added using the + and += operators explained further below. If no Coupling is given, a numerical coefficient of 1.0 is assumed.
OpSum(Op const &op);\nOpSum(Coupling const &cpl, Op const &op);\nOpSum(op::Op)\nOpSum(cpl::Coupling, op::Op)\nAlternatively, an OpSum can also be constructed via the * operator, for example:
auto ops = OpSum();\nfor (int i = 0; i<N; ++i) {\n ops += J * Op(\"SzSz\", {i, (i + 1) % N});\n}\nXDiag allows all couplings to be complex. Depending on the operator type a complex coupling can have two meanings:
A complex prefactor \\(c\\) which upon hermitian conjugation with hc gets conjugated to \\(c^\\star\\). This is the case for the following interaction types: HubbardU, Cdagup, Cdagdn, Cup, Cdn, Nup, Ndn, Ntot, NtotNtot, SdotS, SzSz, Sz, S+, S-, ScalarChirality, tJSzSz, tJSdotS, Matrix Thus, a complex coupling can turn a Hermitian operator to a non-Hermitian operator.
The coupling is part of the definition of the operator. For, example a hopping operator of the form $$ ( t c^\\dagger_{i\\sigma}c_{j\\sigma} + \\textrm{h.c.}) = ( t c^\\dagger_{i\\sigma}c_{j\\sigma} + t^\\star c^\\dagger_{j\\sigma}c_{i\\sigma}) $$ A complex coupling \\(t\\) gives the hopping a phase, but the overall operator remains Hermitian and, thus, invariant under hc. This holds for the types Hop, Hopup, Hopdn, Exchange. In the latter case, complex spin exchange Exchange is defined as, $$ \\frac{1}{2}( J S^+_i S^-_j + J^\\star S^-_iS^+_j)$$
Converts an OpSum with possible string couplings to an OpSum with purely numerical real/complex couplings.
C++JuliaOpSum plain(OpSum const &ops) const;\nplain(ops::OpSum)\nCreates an OpSum with a single pair of Coupling and Op objects.
C++OpSum operator*(std::string cpl, Op const &op);\nOpSum operator*(double cpl, Op const &op);\nOpSum operator*(complex cpl, Op const &op);\nOpSum operator*(Scalar cpl, Op const &op);\nOpSum operator*(Coupling const& cpl, Op const &op);\nAdds two OpSum objects \\(\\mathcal{A} = \\sum_i a_i \\mathcal{A}_i\\) and \\(\\mathcal{B} = \\sum_i b_i \\mathcal{B}_i\\) to for the sum of the two operators, $$ \\mathcal{A} + \\mathcal{B} = \\sum_i a_i \\mathcal{A}_i + \\sum_i b_i \\mathcal{B}_i$$
C++JuliaOpSum &operator+=(OpSum const &ops);\nOpSum operator+(OpSum const &ops) const;\nBase.:+(ops::OpSum, ops2::OpSum)\nSubtracts to OpSum objects.
C++JuliaOpSum &operator-=(OpSum const &ops);\nOpSum operator-(OpSum const &ops) const;\nBase.:-(ops::OpSum, ops2::OpSum)\nMultiplies an OpSum \\(\\mathcal{A} = \\sum_i a_i \\mathcal{A}_i\\) with a scalar \\(b\\) to form
\\[\\mathcal{B} = b \\sum_i a_i \\mathcal{A}_i\\] C++OpSum &operator*=(Scalar const &cpl);\nOpSum &operator/=(Scalar const &cpl);\nOpSum operator*(Scalar const &cpl, OpSum const &op);\nOpSum operator*(OpSum const &op, Scalar const &cpl);\nOpSum operator/(OpSum const &op, Scalar const &cpl);\nSets a coupling constant defined as a string to a numerical value.
C++Scalar &operator[](std::string name);\nScalar const &operator[](std::string name) const;\nReturns a vector of strings with the coupling constants defined, i.e. the strings that define some of the Coupling objects.
C++std::vector<std::string> constants(OpSum const &ops) const;\n// Define the 1D transverse-field Ising chain\nint N = 12;\ndouble J = 1.0;\ndouble h = 0.5;\nauto Sx = arma::mat(\"0 1; 1 0\");\n\n// Option 1: coupling constants as numbers\nauto ops1 = OpSum();\nfor (int i = 0; i<N; ++i) {\n ops1 += J * Op(\"SzSz\", {i, (i+1)%N});\n ops1 += h * Op(\"Matrix\", i, Sx);\n}\n\n// Option 2: coupling constants as strings\nauto ops2 = OpSum();\nfor (int i = 0; i<N; ++i) {\n ops2 += \"J\" * Op(\"SzSz\", {i, (i+1)%N});\n ops2 += \"h\" * Op(\"Matrix\", i, Sx);\n}\nops2[\"J\"] = J;\nops2[\"h\"] = h;\n\nXDIAG_SHOW(isapprox(ops1, ops2));\nXDIAG_SHOW(isapprox(ops1 + ops2, 2.0 * ops1));\nA variable which can be either a real or complex double precision number. The algebraic operators +,-,*,\\ and comparison operators ==, != are defined.
Sources scalar.hpp, scalar.cpp
"},{"location":"documentation/operators/scalar/#constructors","title":"Constructors","text":"C++Scalar(double value);\nScalar(complex value);\nReturns true if Scalar is a real number, false otherwise.
bool isreal(Scalar const &s)\nReturns the real part of the Scalar.
C++double real(Scalar const &s)\nReturns the imaginary part of the Scalar.
C++double imag(Scalar const &s)\nReturns the complex conjugate of the Scalar.
C++Scalar conj(Scalar const &s)\nReturns the absolute value of the Scalar.
C++double abs(Scalar const &s)\nCompares two Scalars to whether they are approximately equal. Return true if the following condition holds.
\\[ | a - b | < atol + rtol*|b| \\] C++bool isapprox(Scalar const &a, Scalar const &b, double rtol = 1e-12, double atol = 1e-12)\nParameters
Name Description Default a first Scalar b second Scalar rtol relative tolerance 1e-12 atol absolute tolerance 0"},{"location":"documentation/operators/scalar/#to_string","title":"to_string","text":"Converts the scalar into a string.
C++std::string to_string(Scalar const &v);\nSymmetrizes an operator with respect to a PermutationGroup or a Representation.
Sources symmetrize.hpp, symmetrize.cpp
Symmetrization in this context means the following. In general, we are given an OpSum of the form,
\\[ O = \\sum_{A\\subseteq \\mathcal{L}} O_A,\\]where \\(O_A\\) denotes a local operator acting on sites \\(A=\\{a_1, \\ldots, a_{l_A}\\}\\) and \\(L\\) denotes the lattice. A PermutationGroup \\(\\mathcal{G}\\) is defined through its permutations \\(\\pi_1, \\ldots, \\pi_M\\). The symmetrized operator returned by this function is then
\\[ O^\\mathcal{G} = \\frac{1}{M}\\sum_{A\\subseteq \\mathcal{L}} \\sum_{\\pi \\in \\mathcal{G}} O_{\\pi(A)},\\]where \\(\\pi(A) = \\{\\pi(a_1), \\ldots,\\pi(a_{l_A})\\}\\) denotes the permutated set of sites of the local operator \\(O_A\\). If a Representation called \\(\\rho\\) is given in addition, the following operator is constructed,
\\[ O^\\mathcal{G, \\rho} = \\frac{1}{M}\\sum_{A\\subseteq \\mathcal{L}} \\sum_{\\pi \\in \\mathcal{G}} \\chi_\\rho(\\pi) O_{\\pi(A)},\\]where \\(\\chi_\\rho(\\pi)\\) denotes the characters of the representation \\(\\rho\\). This routine is useful to evaluate observables in symmetrized blocks.
"},{"location":"documentation/operators/symmetrize/#definition","title":"Definition","text":"C++JuliaOpSum symmetrize(Op const &op, PermutationGroup const &group);\nOpSum symmetrize(Op const &op, Representation const &irrep);\nOpSum symmetrize(OpSum const &ops, PermutationGroup const &group);\nOpSum symmetrize(OpSum const &ops, Representation const &irrep);\nsymmetrize(op::Op, group::PermutationGroup)\nsymmetrize(op::Op, irrep::Representation)\nsymmetrize(ops::OpSum, group::PermutationGroup)\nsymmetrize(ops::OpSum, irrep::Representation)\nint N = 4;\nint nup = 2;\nauto block = Spinhalf(N, nup);\nauto p1 = Permutation({0, 1, 2, 3});\nauto p2 = Permutation({1, 2, 3, 0});\nauto p3 = Permutation({2, 3, 0, 1});\nauto p4 = Permutation({3, 0, 1, 2});\nauto group = PermutationGroup({p1, p2, p3, p4});\nauto rep = Representation(group);\nauto block_sym = Spinhalf(N, rep);\n\nauto ops = OpSum();\nfor (int i=0; i<N; ++i) {\n ops += Op(\"SdotS\", {i, (i+1)%N});\n}\nauto [e0, psi] = eig0(ops, block);\nauto [e0s, psi_sym] = eig0(ops, block_sym);\n\nauto corr = Op(\"SdotS\", {0, 1});\nauto nn_corr = inner(corr, psi);\nauto corr_sym = symmetrize(corr, group);\nauto nn_corr_sym = innerC(corr_sym, psi_sym);\nXDIAG_SHOW(nn_corr);\nXDIAG_SHOW(nn_corr_sym);\nlet\n N = 4\n nup = 2\n block = Spinhalf(N, nup)\n p1 = Permutation([1, 2, 3, 4])\n p2 = Permutation([2, 3, 4, 1])\n p3 = Permutation([3, 4, 1, 2])\n p4 = Permutation([4, 1, 2, 3])\n group = PermutationGroup([p1, p2, p3, p4])\n rep = Representation([1, 1, 1, 1])\n block_sym = Spinhalf(N, group, rep)\n\n ops = OpSum()\n for i in 1:N\n ops += Op(\"HB\", 1.0, [i, mod1(i+1, N)])\n end\n\n e0, psi = eig0(ops, block);\n e0, psi_sym = eig0(ops, block_sym);\n\n corr = Op(\"HB\", 1.0, [1, 2])\n nn_corr = inner(corr, psi)\n corr_sym = symmetrize(corr, group)\n nn_corr_sym = inner(corr_sym, psi_sym)\n @show nn_corr, nn_corr_sym\nend\nSources create_state.hpp, create_state.cpp
"},{"location":"documentation/states/create_state/#product_state","title":"product_state","text":"Creates a filled product state.
C++JuliaState product_state(Block const &block, std::vector<std::string> const &local_state, bool real = true);\nproduct_state(block::Block, local_states::Vector{String}, real::Bool=true)\nCreate a filled random state with normal \\(\\mathcal{N}(0, 1)\\) distributed coefficients.
C++JuliaState random_state(Block const &block, bool real = true, int64_t seed = 42, bool normalized = true);\nrandom_state(block::Block, real::Bool=true, seed::Int64=42, normalized::Bool=true\nCreate a filled state with all zero entries.
C++JuliaState zero_state(Block const &block, bool real = true, int64_t n_cols = 1);\nzero_state(block::Block, real::Bool=true, n_col::Int64=1)\nSet all coefficients of a given state to zero.
C++Juliavoid zero(State &state);\nzero(state::State)\nauto block = Spinhalf(2);\nauto state = product_state(block, {\"Up\", \"Dn\"});\nXDIAG_SHOW(state.vector());\n\nzero(state);\nXDIAG_SHOW(state.vector());\n\nstate = random_state(block, false, 1234, true);\nXDIAG_SHOW(state.vectorC());\n\nstate = zero_state(block, true, 2);\nXDIAG_SHOW(state.vector());\nblock = Spinhalf(2)\nstate = product(block, [\"Up\", \"Dn\"])\ndisplay(vector(state))\n\nzero(state)\ndisplay(vector(state))\n\nstate = rand(block, false, 1234, true)\ndisplay(vector(state))\n\nstate = zeros(block, true, 2)\ndisplay(matrix(state))\nFills a State with a given model state, e.g. a ProductState or a RandomState.
Sources fill.hpp, fill.cpp
"},{"location":"documentation/states/fill/#definition","title":"Definition","text":"C++Juliavoid fill(State &state, ProductState const &pstate, int64_t ncol = 0);\nvoid fill(State &state, RandomState const &rstate, int64_t ncol = 0);\nfill(state::State, pstate::ProductState, ncol::Int64 = 1)\nfill(state::State, rstate::RandomState, ncol::Int64 = 1)\nauto block = Spinhalf(2);\nauto state = State(block); \nauto pstate = ProductState({\"Up\", \"Dn\"});\nfill(state, pstate);\nXDIAG_SHOW(state.vector());\n\nauto rstate = RandomState(1234);\nfill(state, rstate);\nXDIAG_SHOW(state.vector());\nblock = Spinhalf(2)\nstate = State(block)\npstate = ProductState([\"Up\", \"Dn\"])\nfill(state, pstate)\ndisplay(vector(state))\n\nrstate = RandomState(1234)\nfill(state, rstate)\ndisplay(vector(state))\nA product state of local configurations
Sources product_state.hpp, product_state.cpp
"},{"location":"documentation/states/product_state/#constructors","title":"Constructors","text":"C++JuliaProductState(int64_t nsites);\nProductState(std::vector<std::string> const &local_states);\nProductState(nsites::Int64)\nProductState(local_states::Vector{String})\nA ProductState can be iterated over, where at each iteration the string of the local configuration is retured. Here is an example:
C++Juliaauto pstate = ProductState({\"Up\", \"Dn\", \"Emp\", \"UpDn\"});\nfor (auto s : pstate) {\n Log(\"{}\", s);\n}\npstate = ProductState([\"Up\", \"Dn\", \"Emp\", \"UpDn\"])\nfor s in pstate\n @show s\nend\nReturns the number of sites of the product state
C++Juliaint64_t nsites(ProductState const &p);\nnsites(p::ProductState)\nReturns the number of sites of the product state. Same as \"nsites\".
C++Juliaint64_t size(ProductState const &p);\nsize(state::ProductState)\nSets the local configuration at the given site index to the given string.
C++Juliastd::string &operator[](int64_t i);\nsetindex!(state::ProductState, local_state::String, idx::Int64)\nReturns the string of the local configuration at the given site index.
C++Juliastd::string const &operator[](int64_t i) const;\ngetindex(state::ProductState, idx::Int64)\nAdds a local configuration add the end of the product state.
C++Juliavoid push_back(std::string l);\npush!(state::ProductState, local_state::String\nauto pstate = ProductState({\"Up\", \"Dn\", \"Emp\", \"UpDn\"});\nfor (auto s : pstate) {\n Log(\"{}\", s);\n}\nXDIAG_SHOW(to_string(pstate));\n\npstate = ProductState();\npstate.push_back(\"Dn\");\npstate.push_back(\"Up\");\npstate.push_back(\"Dn\");\nXDIAG_SHOW(pstate.nsites());\nfor (auto s : pstate) {\n Log(\"{}\", s);\n}\nXDIAG_SHOW(to_string(pstate));\npstate = ProductState([\"Up\", \"Dn\", \"Emp\", \"UpDn\"])\nfor s in pstate\n @show s\nend\n@show pstate\n\npstate = ProductState()\npush!(pstate, \"Dn\")\npush!(pstate, \"Up\")\npush!(pstate, \"Dn\")\n@show nsites(pstate)\nfor s in pstate\n @show s\nend\n@show pstate\nA random state with \\(\\mathcal{N}(0, 1)\\) normal distributed coefficients.
Sources random_state.hpp, random_state.cpp
"},{"location":"documentation/states/random_state/#constructors","title":"Constructors","text":"C++JuliaRandomState(int64_t seed = 42, bool normalized = true);\nRandomState(seed::Int64 = 42, normalized::Bool = true)\nblock = Spinhalf(2)\nstate = State(block, real=false) # complex State\nrstate1 = RandomState(1234)\nfill(state, rstate1)\ndisplay(vector(state))\n\nrstate2 = RandomState(4321)\nfill(state, rstate2)\ndisplay(vector(state))\n\nfill(state, rstate1)\ndisplay(vector(state))\nauto block = Spinhalf(2);\nauto state = State(block, false); // complex State\nauto rstate1 = RandomState(1234);\nfill(state, rstate1);\nXDIAG_SHOW(state.vectorC());\n\nauto rstate2 = RandomState(4321);\nfill(state, rstate2);\nXDIAG_SHOW(state.vectorC());\n\nfill(state, rstate1);\nXDIAG_SHOW(state.vectorC());\nA generic state describing a quantum wave function \\(|\\psi \\rangle\\).
Sources state.hpp state.cpp
"},{"location":"documentation/states/state/#constructors","title":"Constructors","text":"A state can be constructed in two ways:
By only specifying the block. In this case the state is initialized with all coefficients zero.
C++JuliaState(Block const &block, bool real = true, int64_t n_cols = 1);\nState(block::Block; real::Bool = true, n_cols::Int64 = 1)\nBy handing a vector or matrix of coefficients.
C++Juliatemplate <typename block_t, typename coeff_t>\nState(block_t const &block, arma::Col<coeff_t> const &vector);\n\ntemplate <typename block_t, typename coeff_t>\nState(block_t const &block, arma::Mat<coeff_t> const &matrix);\nState(block::Block, vec::Vector{Float64})\nState(block::Block, vec::Vector{ComplexF64})\nState(block::Block, mat::Matrix{Float64})\nState(block::Block, mat::Matrix{ComplexF64})\nReturns the number of sites of the block the state is defined on.
C++Juliaint64_t nsites(State const &s) const\nnsites(state::State)\nReturns whether two states are approximately equal.
C++Juliabool isapprox(State const &v, State const &w, double rtol = 1e-12,\n double atol = 1e-12);\nisapprox(v::State, w::State, rtol::Float64, atol::Float64)\nReturns whether the state is real.
C++Juliaint64_t isreal(State const &s) const;\nisreal(state::State)\nReturns whether the real part of the State.
C++JuliaState real(State const &s) const;\nreal(state::State)\nReturns whether the imaginary part of the State.
C++JuliaState imag(State const &s) const;\nimag(state::State)\nTurns a real State into a complex State. Does nothing if the state is already complex
C++Juliavoid make_complex(State &s);\nmake_complex!(state::State)\nReturns the dimension of the block the state is defined on.
C++Juliaint64_t dim(State const &s) const;\ndim(block::Spinhalf)\nReturns the size of the block the state is defined on. locally. Same as \"dim\" for non-distributed Blocks but different for distributed blocks.
C++Juliaint64_t size(State const &s);\nsize(s::State)\nReturns number of rows of the local storage. Same as \"size\"
C++Juliaint64_t nrows(State const &s);\nnrows(s::State)\nReturns number of columns of the local storage.
C++Juliaint64_t ncols(State const &s);\nncols(s::State)\nReturns a state created from the n-th column of the storage. Whether or not the storage is copied can be specified by setting the flag \"copy\".
C++JuliaState col(State const &s, int64_t n, bool copy = true);\ncol(s::State, n::Int64 = 1; copy::Bool = true)\nReturns a vector from the n-th column of the storage. In C++ use \"vector\"/\"vectorC\" to either get a real or complex vector.
C++Juliaarma::vec vector(State const &s, int64_t n = 0, bool copy = true);\narma::cx_vec vectorC(State const &s, int64_t n = 0, bool copy = true);\nvector(state::State; n::Int64 = 1, copy::Bool=true)\n# no vectorC method in julia\nReturns matrix representing the storage. In C++ use \"matrix\"/\"matrixC\" to either get a real or complex matrix.
C++Juliaarma::vec matrix(State const &s, bool copy = true);\narma::cx_vec matrixC(State const &s, bool copy = true);\nmatrix(state::State, copy::Bool=true)\n# no matrixC method in julia\nauto block = Spinhalf(2);\nauto psi1 = State(block, arma::vec(\"1.0 2.0 3.0 4.0\"));\nXDIAG_SHOW(psi1);\nXDIAG_SHOW(vector(psi1));\nmake_complex(psi1);\nXDIAG_SHOW(vectorC(psi1));\n\nauto psi2 = State(block, false, 3);\nXDIAG_SHOW(psi2);\nXDIAG_SHOW(matrixC(psi2));\n\nauto psi3 = State(block, arma::cx_vec(arma::vec(\"1.0 2.0 3.0 4.0\"),\n arma::vec(\"4.0 3.0 2.0 1.0\")));\nXDIAG_SHOW(vectorC(psi3));\nXDIAG_SHOW(vector(real(psi3)));\nXDIAG_SHOW(vector(imag(psi3)));\nblock = Spinhalf(2)\npsi1 = State(block, [1.0, 2.0, 3.0, 4.0])\n@show psi1\ndisplay(vector(psi1))\nmake_complex!(psi1)\ndisplay(vector(psi1))\n\npsi2 = State(block, real=false, n_cols=3)\n@show psi2\ndisplay(matrix(psi2))\n\npsi3 = State(block, [1.0+4.0im, 2.0+3.0im, 3.0+2.0im, 4.0+1.0im])\ndisplay(vector(psi3))\ndisplay(vector(real(psi3)))\ndisplay(vector(imag(psi3)))\nPermutations of indices or lattice sites. Basic building block of a PermutationGroup. Permutations can be multiplied, inverted and raised to a power.
Sources permutation.hpp, permutation.cpp
"},{"location":"documentation/symmetries/permutation/#constructors","title":"Constructors","text":""},{"location":"documentation/symmetries/permutation/#from-an-array","title":"From an array","text":"Creates an Permutation out of an array of integers, e.g. {0, 2, 1, 3}. If the input array is of size N then every number between 0 and N-1 must occur exactly once, otherwise the Permutation is invalid.
1-indexing in Julia / 0-indexing in C++
To enumerate the sites of a Permutation, we start counting at 1 in Julia and 0 in C++.
C++JuliaPermutation(std::initializer_list<int64_t> list);\nPermutation(std::vector<int32_t> const &array);\nPermutation(std::vector<int64_t> const &array);\nPermutation(arma::Col<int64_t> const &array);\nPermutation(int64_t *array, int64_t size);\nPermutation(array::Vector{Int64})\nConstructs an identity permutation of a given size, e.g. {0, 1, 2, 3}.
Permutation(int64_t size);\nPermutation(size::Integer)\nComputes the inverse permutation.
C++JuliaPermutation inverse(Permutation const &p);\ninverse(perm::Permutation)\nConcatenates two permutations by overloading the * operator.
Permutation operator*(Permutation const &p1, Permutation const &p2);\nBase.:*(p1::Permutation, p2::Permutation)\nRaises a permutation to an integer power.
C++JuliaPermutation pow(Permutation const &p, int64_t power);\nBase.:^(p::Permutation, power::Integer)\nReturns the size of a Permutation.
C++Juliaint64_t size(Permutation const &p);\nsize(p::Permutation)\nPermutation p1 = {0, 2, 1, 3};\nPermutation p2 = {2, 0, 1, 3};\n\nXDIAG_SHOW(inverse(p1));\nXDIAG_SHOW(p1*p2);\np1 = Permutation([1, 3, 2, 4])\np2 = Permutation([3, 1, 2, 4])\n\n@show inverse(p1)\n@show p1 * p2\nA group of permutations. Group axioms are verified during construction.
Sources permutation_group.hpp, permutation_group.cpp
"},{"location":"documentation/symmetries/permutation_group/#constructor","title":"Constructor","text":""},{"location":"documentation/symmetries/permutation_group/#from-permutations","title":"From Permutations","text":"Creates an PermutationGroup out of a vector of Permutation objects.
C++JuliaPermutationGroup(std::vector<Permutation> const &permutations);\nPermutationGroup(permutations::Vector{Permutation})\nCreates a PermutationGroup out of a matrix whose rows specify the individual permutations. If a raw pointer is handed, the matrix is assumed to be in column-major form.
C++PermutationGroup(arma::Mat<int64_t> const &matrix);\nPermutationGroup(int64_t *ptr, int64_t n_permutations, int64_t nsites);\nReturns the number of sites on which the permutations of the group acts.
C++Juliaint64_t nsites(PermutationGroup const &group);\nnsites(group::PermutationGroup)\nReturns the size of the permutation group, i.e. the number permutations.
C++Juliaint64_t size(PermutationGroup const &group);\nsize(group::PermutationGroup)\n// Define a cyclic group of order 3\nPermutation p1 = {0, 1, 2};\nPermutation p2 = {1, 2, 0};\nPermutation p3 = {2, 0, 1};\nauto C3 = PermutationGroup({p1, p2, p3});\n\nXDIAG_SHOW(C3.size());\nXDIAG_SHOW(C3.nsites());\nXDIAG_SHOW(C3.inverse(1)); // = 2\n# Define a cyclic group of order 3\np1 = Permutation([1, 2, 3])\np2 = Permutation([2, 3, 1])\np3 = Permutation([3, 1, 2])\nC3 = PermutationGroup([p1, p2, p3])\n\n@show size(C3)\n@show nsites(C3)\n@show inverse(C3, 1) # = 2\nA (1D) irreducible representation of a finite group. Upon creation, the group homomorphism properties are verified.
Source representation.hpp
"},{"location":"documentation/symmetries/representation/#constructors","title":"Constructors","text":""},{"location":"documentation/symmetries/representation/#trivial-representation","title":"Trivial representation","text":"Creates the trivial representation (all characters equal to 1) of a PermutationGroup
C++JuliaRepresentation(PermutationGroup const &group);\nRepresentation(group::PermutationGroup)\nCreates a 1D representation of a PermutationGroup with given real or complex characters.
C++template <typename T>\nRepresentation(PermutationGroup const &group, std::vector<T> const &characters);\ntemplate <typename T>\nRepresentation(PermutationGroup const &group, arma::Col<T> const &characters);\ntemplate <typename T>\nRepresentation(PermutationGroup const &group, T *characters, int64_t n_characters);\nThe template parameter T in C++ can either be double or complex.
Returns the size of the Representation, i.e. the number of group elements represented.
C++Juliaint64_t size(Representation const &irrep);\nsize(irrep::Representation)\nReturns the whether or not the Representation is real, I.E. the characters are real numbers and do not have an imaginary part.
C++Juliabool isreal(Representation const &irrep) const;\nisreal(irrep::Representation)\nMultiplies two Representations by overloading the * operator.
Representation operator*(Representation const &r1, Representation const &r2);\nBase.:*(r1::Representation, r2::Representation)\nPermutation p = {1, 2, 3, 0};\nauto C4 = PermutationGroup({pow(p, 0), pow(p, 1), pow(p, 2), pow(p, 3)});\nRepresentation r1(C4, arma::vec{1.0, -1.0, 1.0, -1.0});\nRepresentation r2(C4, arma::cx_vec{1.0, 1.0i, -1.0, -1.0i});\nXDIAG_SHOW(r1 * r2);\np = Permutation([2, 3, 4, 1])\nC4 = PermutationGroup(pow(p, 0), pow(p, 1), pow(p, 2), pow(p, 3))\nr1 = Representation(C4, [1, -1, 1, -1])\nr2 = Representation(C4, [1, 1im, -1, -1im])\n\n@show r1 * r2\nSources logger.hpp, logger.cpp
"},{"location":"documentation/utilities/logging/#setting-the-verbosity","title":"Setting the verbosity","text":"Algorithms implemented in XDiag do not output anything during their execution by default. However, it is typically useful to get some information on how the code is performing and even intermediary results at runtime. For this, the verbosity of the internal XDiag logging can be set using the function set_verbosity, which is defined as
void set_verbosity(int64_t level);\nset_verbosity(level::Integer);\nThere are several levels of verbosity, defining how much information is shown.
level outputed information 0 no information 1 some information 2 detailed informationFor example, when computing a ground state energy using the eigval0 function, we can set a higher verbosity level using
C++Juliaset_verbosity(2);\ndouble e0 = eigval0(bonds, block);\nset_verbosity(2);\ne0 = eigval0(bonds, block);\nThis will print detailed information, which can look like this
Lanczos iteration 1\nMVM: 0.00289 secs\nalpha: -0.2756971549998545\nbeta: 1.7639347562074059\neigs: -0.2756971549998545\nLanczos iteration 2\nMVM: 0.00244 secs\nalpha: -0.7116140394927443\nbeta: 2.3044797637130743\neigs: -2.2710052270892791 1.2836940325966804\nLanczos iteration 3\nMVM: 0.00210 secs\nalpha: -1.2772539678430306\nbeta: 2.6627870395174456\neigs: -3.7522788386927637 -0.6474957945455240 2.1352094709026579\nProducing nicely formatted output is unfortunately a bit cumbersome in standard C++. For this, the Log mechanism in XDiag can help. To simply write out a line of information you can call,
Log(\"hello from the logger\");\nBy default, a new line is added. It is also possible to set verbosity by handing the level as the first argument,
Log(2, \"hello from the logger only if global verbosity is set to >= 2\");\nThis message will only appear if the global verbosity level is set to a value \\(\\geq 2\\). Finally, XDiag also supports formatted output by using the fmtlib library. For example, numbers can be formated this way
Log(\"pi is around {:.4f} and the answer is {}\", 3.141592, 42);\nIn standard C++ measuring time is a bit awkward. To quickly monitor the CPU time spent by XDiag by simple functions.
"},{"location":"documentation/utilities/timing/#simple-timing-using-tic-toc","title":"Simple timing using tic() / toc()","text":"Similar as in Matlab one can use tic() and toc() to measure the time spent between two points in the code.
tic();\ndouble e0 = eigval0(bonds, block);\ntoc();\ntoc() will output the time spent since the last time tic() has been called.
To get the present time, simply call
auto time = rightnow();\nA timing (in second) between two time points can be written to output using
timing(begin, end);\nThis can even be accompanied by a message about what is being timed and a verbosity level (see Logging) can also be set. The full call signature is
timing(begin, end, message, level);\nrightnow() end end time computed using rightnow() message message string to be prepended to timing \"\" level verbosity level at which timing is printed 0"},{"location":"documentation/utilities/utils/","title":"Utility functions","text":""},{"location":"documentation/utilities/utils/#set_verbosity","title":"set_verbosity","text":"Set how much logging is generated my XDiag to monitor the progress and behaviour of the code. There are three verbosity levels that can be set:
This can be useful, e.g. to monitor the progress of an iterative algorithm
C++Juliavoid set_verbosity(int64_t level);\nset_verbosity(level::Int64)\nPrints a nice welcome message containing the version number and git commit used.
C++Juliavoid say_hello()\nsay_hello()\nIf say_hello is too much flower power for you, one can also just have a boring print-out of the version number using this function.
void print_version()\nprint_version()\nFor quick debugging in C++, XDiag features a simple macro which outputs the name and content of a variable calles XDIAG_SHOW(x). For example
Spinhalf block(16, 8);\nXDIAG_SHOW(block);\nwill write an output similar to
block:\n nsites : 16\n nup : 8\n dimension: 12,870\n ID : 0xa9127434d66b9878\nThe XDIAG_SHOW(x) macro can be used on any XDiag object and several other standard C++ objects as well.
cmake_minimum_required(VERSION 3.19)\n\nproject(\n tj_distributed_time_evolve\n)\n\nfind_package(xdiag_distributed REQUIRED HINTS ../../../install)\nadd_executable(main main.cpp)\ntarget_link_libraries(main PUBLIC xdiag::xdiag_distributed)\ncmake_minimum_required(VERSION 3.19)\n\nproject(\n hello_world\n)\n\nfind_package(xdiag REQUIRED HINTS \"/path/to/xdiag/install\")\nadd_executable(main main.cpp)\ntarget_link_libraries(main PRIVATE xdiag::xdiag)\n using XDiag\n say_hello()\n #include <xdiag/all.hpp>\n\n using namespace xdiag;\n\n int main() try {\n say_hello();\n } catch (Error e) {\n error_trace(e);\n }\nusing XDiag\n\nlet \n N = 16\n nup = N \u00f7 2\n block = Spinhalf(N, nup)\n\n # Define the nearest-neighbor Heisenberg model\n ops = OpSum()\n for i in 1:N\n ops += Op(\"HB\", \"J\", [i-1, i % N])\n end\n ops[\"J\"] = 1.0\n\n set_verbosity(2) # set verbosity for monitoring progress\n e0 = eigval0(ops, block) # compute ground state energy\n\n println(\"Ground state energy: $e0\")\nend\n#include <xdiag/all.hpp>\n\nusing namespace xdiag;\n\nint main() try {\n int N = 16;\n int nup = N / 2;\n Spinhalf block(N, nup);\n\n // Define the nearest-neighbor Heisenberg model\n OpSum ops;\n for (int i = 0; i < N; ++i) {\n ops += Op(\"HB\", \"J\", {i, (i + 1) % N});\n }\n ops[\"J\"] = 1.0;\n\n set_verbosity(2); // set verbosity for monitoring progress\n double e0 = eigval0(ops, block); // compute ground state energy\n\n Log(\"Ground state energy: {:.12f}\", e0);\n\n} catch (Error e) {\n error_trace(e);\n}\n#include <xdiag/all.hpp>\n\nusing namespace xdiag;\n\nvoid measure_density(int nsites, State const &v) {\n int rank;\n MPI_Comm_rank(MPI_COMM_WORLD, &rank);\n for (int i = 0; i < nsites; ++i) {\n complex sz = innerC(Bond(\"Ntot\", i), v);\n if (rank == 0) {\n printf(\"%.6f \", std::real(sz));\n }\n }\n if (rank == 0) {\n printf(\"\\n\");\n }\n}\n\nint main(int argc, char **argv) try {\n MPI_Init(&argc, &argv);\n\n int L = 6;\n int W = 4;\n double t = 1.0;\n double J = 0.1;\n double mu_0 = 10;\n\n int nsites = L * W;\n double precision = 1e-12;\n\n // Create square lattice t-J model\n OpSum ops;\n for (int x = 0; x < L-1; ++x) {\n for (int y = 0; y < W; ++y) {\n int nx = (x + 1) % L;\n int ny = (y + 1) % W;\n\n int site = x * W + y;\n int right = nx * W + y;\n int top = x * W + ny;\n ops += Op(\"Hop\", \"T\", {site, right});\n ops += Op(\"Exchange\", \"J\", {site, right});\n ops += Op(\"Hop\", \"T\", {site, top});\n ops += Op(\"Exchange\", \"J\", {site, top});\n\n\n\n if (x < L / 2) {\n Log(\"x {} y {} site {} t {} r {} +\", x, y, site, top, right);\n ops += Op(\"Ntot\", \"MUPLUS\", site);\n } else {\n Log(\"x {} y {} site {} t {} r {} -\", x, y, site, top, right);\n ops += Op(\"Ntot\", \"MUNEG\", site);\n }\n }\n }\n ops[\"T\"] = t;\n ops[\"J\"] = J;\n ops[\"MUPLUS\"] = mu_0;\n ops[\"MUNEG\"] = mu_0;\n\n auto block = tJDistributed(nsites, nsites / 2 - 1, nsites / 2 - 1);\n\n XDIAG_SHOW(block);\n\n Log.set_verbosity(2);\n tic();\n auto [e0, v] = eig0(ops, block);\n toc(\"gs\");\n\n ops[\"MUPLUS\"] = 0;\n ops[\"MUNEG\"] = 0;\n\n measure_density(nsites, v);\n\n // Do the time evolution with a step size tau\n double tau = 0.1;\n for (int i = 0; i < 40; ++i) {\n tic();\n v = time_evolve(ops, v, tau, precision);\n toc(\"time evolve\");\n tic();\n measure_density(nsites, v);\n toc(\"measure\");\n }\n\n MPI_Finalize();\n return EXIT_SUCCESS;\n} catch (std::exception const &e) {\n traceback(e);\n}\nLet us set up our first program using the xdiag library.
using XDiag\nsay_hello()\n#include <xdiag/all.hpp>\n\nusing namespace xdiag;\n\nint main() try {\n say_hello();\n} catch (Error e) {\n error_trace(e);\n}\nThe function say_hello() prints out a welcome message, which also contains information which exact XDiag version is used. In Julia this is all there is to it.
For the C++ code we need to create two files to compile the program. The first is the actual C++ code. What is maybe a bit unfamiliar is the try / catch block. XDiag implements a traceback mechanism for runtime errors, which is activated by this idiom. While not stricly necessary here, it is a good practice to make use of this.
Now that the application program is written, we next need to set up the compilation instructions using CMake. To do so we create a second file called CMakeLists.txt in the same directory.
cmake_minimum_required(VERSION 3.19)\n\nproject(\n hello_world\n)\n\nfind_package(xdiag REQUIRED HINTS \"/path/to/xdiag/install\")\nadd_executable(main main.cpp)\ntarget_link_libraries(main PRIVATE xdiag::xdiag)\nYou should replace \"/path/to/xdiag/install\" with the appropriate directory where your XDiag library is installed after compilation. This exact CMakeLists.txt file can be used to compile any XDiag application.
Info
For using the distributed XDiag library the last line of the above CMakeLists.txt should be changed to
target_link_libraries(main PUBLIC xdiag::xdiag_distributed)\nWe then compile the application code,
cmake -S . -B build\ncmake --build build\n./build/main\nWe are now ready to run our first actual calculation using XDiag. Our immediate goal will be to determine the ground state energy of a \\(S=1/2\\) Heisenberg model on a 1D chain lattice with periodic boundary conditions,
\\[ H = J \\sum_{\\langle i, j\\rangle} \\mathbf{S}_i \\cdot \\mathbf{S}_j,\\]where \\(\\mathbf{S}_i = (S^x_i, S^y_i, S^z_i)\\) denotes the vector of spin matrices at a given site \\(i\\). The notation \\(\\langle i, j\\rangle\\) refers to summatation over neighboring sites \\(i\\) and \\(j\\).
The first thing to define before any computation, is the Hilbert space our model will be defined on. For the present example, we use the Hilbert space class Spinhalf. We consider a chain lattice with \\(N=8\\) sites and create a \\(S=1/2\\) Hilbert space:
C++using namespace xdiag;\nint N = 8;\nauto hspace = Spinhalf(N);\nWe would like to know which spin configurations, the Hilbert space is made up of. To do so, we can iterate over the Hilbert space and print out the spin configurations.
C++for (auto spins : hspace) {\n Log(\"{}\", to_string(spins));\n}\nThis produces an output similar to the following:
\u2193\u2193\u2193\u2193\u2193\u2193\u2193\u2193\n\u2193\u2193\u2193\u2193\u2193\u2193\u2193\u2191\n\u2193\u2193\u2193\u2193\u2193\u2193\u2191\u2193\n\u2193\u2193\u2193\u2193\u2193\u2193\u2191\u2191\n\u2193\u2193\u2193\u2193\u2193\u2191\u2193\u2193\n\u2193\u2193\u2193\u2193\u2193\u2191\u2193\u2191\n\u2193\u2193\u2193\u2193\u2193\u2191\u2191\u2193\n...\nHere we already see several things at work. XDiag features a convenient way to write logs, with the Log class. The first argument to Log() is a format string. In C++ we use the fmt library, to be able to write structured format and format our output. The second argument turns our spins into a string. spins is of type ProductState, whose configuration on each site can be individually addressed.
Further, we notice that all \\(2^N\\) spin configurations are included in this Hilbert space. However, the Heisenberg model conserves the total \\(S^z = \\sum_i S^z_i\\), and thus we could limit ourselves to a block of the Hilbert space, which only contains configurations of a certain magnetization:
C++int nup = 4;\nauto block = Spinhalf(N, nup);\nfor (auto spins : block) {\n Log(\"{}\", to_string(spins));\n}\nThis produces an output similar to:
\u2193\u2193\u2193\u2193\u2191\u2191\u2191\u2191\n\u2193\u2193\u2193\u2191\u2193\u2191\u2191\u2191\n\u2193\u2193\u2193\u2191\u2191\u2193\u2191\u2191\n\u2193\u2193\u2193\u2191\u2191\u2191\u2193\u2191\n\u2193\u2193\u2193\u2191\u2191\u2191\u2191\u2193\n\u2193\u2193\u2191\u2193\u2193\u2191\u2191\u2191\n\u2193\u2193\u2191\u2193\u2191\u2193\u2191\u2191\n...\nWe see that the block now only contains configurations, where the total number of spins pointing up is 4. We can now quickly check the dimension of the Hilbert spaces, and confirm that the dimension of the block is reduced from \\(2^8=256\\) to \\(\\begin{pmatrix} 8 \\\\ 4 \\end{pmatrix} = 70\\),
C++XDIAG_SHOW(hspace.size());\nXDIAG_SHOW(block.size());\nwhich should print:
hspace.size():\n256\nblock.size():\n70\nHere, we introduced another functionality of XDiag, the XDIAG_SHOW macro which can be used for quick debug printing of XDiag objects.
"},{"location":"user_guide/input_output/","title":"Input/Output","text":"Julia features a variety of packages facilitating input and output of data. For C++, we provide convenient functionality for TOML and HDF5 files.
"},{"location":"user_guide/input_output/#toml","title":"TOML","text":"For simulations is can often be useful to read input paramters from a file. The
"},{"location":"user_guide/installation/","title":"Installation","text":""},{"location":"user_guide/installation/#julia-installation","title":"Julia installation","text":"Enter the package mode using ] in the Julia REPL and type:
add XDiag\nThat's it!
"},{"location":"user_guide/installation/#c-compilation","title":"C++ compilation","text":"Using XDiag with C++ is a two-step process. First the xdiag library needs to be compiled and installed. Therafter, application codes are compiled in a second step. The library can be compiled in two different versions:
g++, clang, or Intel's icpx)Download the source code using git
cd /path/to/where/xdiag/should/be\ngit clone https://github.com/awietek/xdiag.git\nCompile the default library
cd xdiag\ncmake -S . -B build\ncmake --build build\ncmake --install build\ninstall. Compile the distributed library
To use the distributed computing features of xdiag, the distributed library has to be built which requires MPI.
cd xdiag\ncmake -S . -B build -D XDIAG_DISTRIBUTED=On\ncmake --build build\ncmake --install build\nInfo
It might be necessary to explicitly define MPI compiler, e.g. mpicxx like this
cmake -S . -B build -D XDIAG_DISTRIBUTED=On -D CMAKE_CXX_COMPILER=mpicxx\nThe compilation process can be modified and also allows for further optimizations. We collect several common scenarios in the Advanced Compilation guide.
"}]} \ No newline at end of file