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}\nEnter 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 \"../../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:
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
Source algebra.hpp
"},{"location":"documentation/algebra/algebra/#norm","title":"norm","text":"Computes the 2-norm of a state
JuliaC++norm(state::State)\ndouble norm(State const &v);\nComputes the 1-norm of a state
JuliaC++norm1(state::State)\ndouble norm1(State const &v);\nComputes the \\(\\infty\\)-norm of a state
JuliaC++norminf(state::State)\ndouble norminf(State const &v);\nComputes the dot product between two states. In C++, please use the dotC function if one of the two states is expected to be complex.
JuliaC++dot(v::State, w::State)\ndouble dot(State const &v, State const &w);\ncomplex dotC(State const &v, State const &w);\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.
JuliaC++inner(ops::OpSum, v::State)\ninner(op::Op, v::State)\ndouble inner(OpSum const &ops, State const &v);\ndouble inner(Op const &op, State const &v);\ncomplex innerC(OpSum const &ops, State const &v);\ncomplex innerC(Op const &op, State const &v);\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\nint N = 8;\nauto block = Spinhalf(N, N / 2);\nauto ops = OpSum();\nfor (int i=0; i<N; ++i) {\n ops += Op(\"HB\", 1.0, {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 = rand(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());\nApplies an operator to a state \\(\\vert w \\rangle = O \\vert v\\rangle\\) or block of states \\(\\left( \\vert w_1 \\rangle \\dots \\vert w_M \\rangle \\right) = O \\left( \\vert v_1 \\rangle \\dots \\vert v_M \\rangle \\right)\\).
Source apply.hpp
JuliaC++apply(op::Op, v::State, w::State, precision::Float64 = 1e-12)\napply(ops::OpSum, v::State, w::State, precision::Float64 = 1e-12)\nvoid apply(Op const &op, State const &v, State &w, double precision = 1e-12);\nvoid apply(OpSum const &ops, State const &v, State &w, double precision = 1e-12);\nThe resulting state is handed as the third argument and is overwritten upon exit.
"},{"location":"documentation/algebra/apply/#parameters","title":"Parameters","text":"Name Description ops / op OpSum or Op defining the operator v input state $\\vert v\\rangle $ w output state \\(\\vert w \\rangle = O \\vert v\\rangle\\) precision precision with which checks for zero are performed (default \\(10^{-12}\\))"},{"location":"documentation/algebra/apply/#usage-example","title":"Usage Example","text":"JuliaC++let \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\nint N = 8;\nauto block = Spinhalf(N, N / 2);\nauto ops = OpSum();\nfor (int i=0; i<N; ++i){\n ops += Op(\"HB\", 1.0, {i, (i+1)%N});\n}\nauto [e0, psi] = eig0(ops, block);\n\nauto blockp = Spinhalf(N, N / 2 + 1);\nauto phi = zeros(blockp);\napply(Op(\"S+\", 1.0, 2), psi, phi);\nXDIAG_SHOW(inner(ops, psi));\nXDIAG_SHOW(inner(ops, phi));\nmatrix(ops, block; force_complex=false)\nmatrixC(ops, block) # c++ only\nCreates the full matrix representation of a given OpSum on a block.
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.
Source matrix.hpp
"},{"location":"documentation/algebra/matrix/#parameters","title":"Parameters","text":"Name Description ops OpSum defining the operator block block on which the operator is defined force_complex flag to determine if returned matrix is forced to be complex (Julia only)"},{"location":"documentation/algebra/matrix/#returns","title":"Returns","text":"Type Description matrix matrix representation of opsum on block"},{"location":"documentation/algebra/matrix/#definition","title":"Definition","text":"JuliaC++matrix(ops::Opsum, block::Block; force_complex::Bool=false)\ntemplate <typename block_t>\narma::mat matrix(OpSum const &ops, block_t const &block);\n\ntemplate <typename block_t>\narma::cx_mat matrixC(OpSum const &ops, block_t const &block);\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\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\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 += Op(\"HOP\", \"T\", {i, (i+1) % N});\n}\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({1, -1, 1, -1});\nauto block = Electron(N, nup, ndn, group, irrep);\n\nauto H = matrix(ops, block);\nH.print();\nComputes the groud state energy and the ground state of an operator on a block.
Source sparse_diag.hpp
JuliaC++function 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)\nstd::tuple<double, State> eig0(OpSum const &ops, Block const &block,\n double precision = 1e-12,\n int64_t max_iterations = 1000,\n bool force_complex = false,\n int64_t random_seed = 42);\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\nint 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(\"HB\", \"J\", {i, (i+1) % N});\n}\nops[\"J\"] = 1.0;\nauto [e0, gs] = eig0(ops, block);\nPerforms an iterative eigenvalue calculation building eigenvectors using the Lanczos algorithm
Source eigs_lanczos.hpp
JuliaC++C++ (explicit state initialization)function 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)\neigs_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);\neigs_lanczos_result_t \neigs_lanczos(OpSum const &ops, Block const &block,\n State const &state, int64_t neigenvalues = 1,\n double precision = 1e-12, int64_t max_iterations = 1000,\n bool force_complex = false);\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/eigval0/","title":"eigval0","text":"Computes the groud state energy of an operator on a block.
Source sparse_diag.hpp
JuliaC++function 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)\ndouble 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);\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\nint 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(\"HB\", \"J\", {i, (i+1) % N});\n}\nops[\"J\"] = 1.0;\ndouble e0 = eigval0(ops, block);\nPerforms an iterative eigenvalue calculation using the Lanczos algorithm.
Source eigvals_lanczos.hpp
JuliaC++C++ (explicit state initialization)function 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)\neigvals_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);\neigvals_lanczos_result_t\neigvals_lanczos(OpSum const &ops, Block const &block, \n State const &state, int64_t neigvals = 1,\n double precision = 1e-12, int64_t max_iterations = 1000,\n bool force_complex = false);\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 niterations number of iterations performed criterion string denoting the reason why the algorithm stopped"},{"location":"documentation/blocks/electron/","title":"Electron","text":"Representation of a block in a Electron (spinful fermion) Hilbert space.
Source electron.hpp
"},{"location":"documentation/blocks/electron/#constructors","title":"Constructors","text":"JuliaC++Electron(n_sites::Integer)\nElectron(n_sites::Integer, n_up::Integer, n_dn::Integer)\nElectron(n_sites::Integer, group::PermutationGroup, irrep::Representation)\nElectron(n_sites::Integer, n_up::Integer, n_dn::Integer, \n group::PermutationGroup, irrep::Representation)\nElectron(int64_t n_sites);\nElectron(int64_t n_sites, int64_t n_up, int64_t n_dn);\nElectron(int64_t n_sites, PermutationGroup permutation_group,\n Representation irrep);\nElectron(int64_t n_sites, int64_t n_up, int64_t n_dn, \n PermutationGroup group, Representation irrep);\nAn Electron block can be iterated over, where at each iteration a ProductState representing the corresponding basis state is returned.
JuliaC++block = Electron(4, 2, 2)\nfor pstate in block\n @show pstate, index(block, pstate) \nend\nauto block = Electron(4, 2, 2);\nfor (auto pstate : block) {\n Log(\"{} {}\", to_string(pstate), block.index(pstate));\n}\nindex
Returns the index of a given ProductState in the basis of the Electron block.
JuliaC++index(block::Electron, pstate::ProductState)\nint64_t index(ProductState const &pstate) const;\n1-indexing
In the C++ version, the index count starts from \"0\" whereas in Julia the index count starts from \"1\".
n_sites
Returns the number of sites of the block.
JuliaC++n_sites(block::Electron)\nint64_t n_sites() const;\nn_up
Returns the number of \"up\" electrons.
JuliaC++n_up(block::Electron)\nint64_t n_up() const;\nn_dn
Returns the number of \"down\" electrons.
JuliaC++n_dn(block::Electron)\nint64_t n_dn() const;\npermutation_group
Returns the PermutationGroup of the block, if defined.
JuliaC++permutation_group(block::Electron)\nPermutationGroup permutation_group() const;\nirrep
Returns the Representation of the block, if defined.
JuliaC++irrep(block::Electron)\nRepresentation irrep() const;\nsize
Returns the size of the block, i.e. its dimension.
JuliaC++size(block::Electron)\nint64_t size() const;\ndim
Returns the dimension of the block, same as \"size\" for non-distributed blocks.
JuliaC++dim(block::Electron)\nint64_tdim() const;\nisreal
Returns whether the block can be used with real arithmetic. Complex arithmetic is needed when a Representation is genuinely complex.
JuliaC++isreal(block::Electron; precision::Real=1e-12)\nint64_t isreal(double precision = 1e-12) const;\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 n_sites(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)\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({1, -1, 1, -1});\nauto block_sym = Electron(N, group, irrep);\nXDIAG_SHOW(block_sym);\n\n// with symmetries and number conservation\nauto block_sym_np = Electron(N, nup, ndn, group, irrep);\nXDIAG_SHOW(block_sym_np);\n\nXDIAG_SHOW(block_sym_np.n_sites());\nXDIAG_SHOW(block_sym_np.size());\n\n// Iteration\nfor (auto pstate : block_sym_np) {\n Log(\"{} {}\", to_string(pstate), block_sym_np.index(pstate));\n}\nXDIAG_SHOW(block_sym_np.permutation_group());\nXDIAG_SHOW(block_sym_np.irrep());\nRepresentation of a block in a spin \\(S=1/2\\) Hilbert space.
Source spinhalf.hpp
"},{"location":"documentation/blocks/spinhalf/#constructors","title":"Constructors","text":"JuliaC++Spinhalf(n_sites::Integer)\nSpinhalf(n_sites::Integer, n_up::Integer)\nSpinhalf(n_sites::Integer, group::PermutationGroup, irrep::Representation)\nSpinhalf(n_sites::Integer, n_up::Integer, group::PermutationGroup, \n irrep::Representation)\nSpinhalf(int64_t n_sites);\nSpinhalf(int64_t n_sites, int64_t n_up);\nSpinhalf(int64_t n_sites, PermutationGroup permutation_group,\n Representation irrep);\nSpinhalf(int64_t n_sites, int64_t n_up, PermutationGroup group,\n Representation irrep);\nAn Spinhalf block can be iterated over, where at each iteration a ProductState representing the corresponding basis state is returned.
JuliaC++block = Spinhalf(4, 2)\nfor pstate in block\n @show pstate, index(block, pstate) \nend\nauto block = Spinhalf(4, 2);\nfor (auto pstate : block) {\n Log(\"{} {}\", to_string(pstate), block.index(pstate));\n}\nindex
Returns the index of a given ProductState in the basis of the Spinhalf block.
JuliaC++index(block::Spinhalf, pstate::ProductState)\nint64_t index(ProductState const &pstate) const;\n1-indexing
In the C++ version, the index count starts from \"0\" whereas in Julia the index count starts from \"1\".
n_sites
Returns the number of sites of the block.
JuliaC++n_sites(block::Spinhalf)\nint64_t n_sites() const;\nn_up
Returns the number of \"up\" spins.
JuliaC++n_up(block::Spinhalf)\nint64_t n_up() const;\npermutation_group
Returns the PermutationGroup of the block, if defined.
JuliaC++permutation_group(block::Spinhalf)\nPermutationGroup permutation_group() const;\nirrep
Returns the Representation of the block, if defined.
JuliaC++irrep(block::Spinhalf)\nRepresentation irrep() const;\nsize
Returns the size of the block, i.e. its dimension.
JuliaC++size(block::Spinhalf)\nint64_t size() const;\ndim
Returns the dimension of the block, same as \"size\" for non-distributed blocks.
JuliaC++dim(block::Spinhalf)\nint64_tdim() const;\nisreal
Returns whether the block can be used with real arithmetic. Complex arithmetic is needed when a Representation is genuinely complex.
JuliaC++isreal(block::Spinhalf; precision::Real=1e-12)\nint64_t isreal(double precision = 1e-12) const;\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 n_sites(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)\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({1, -1, 1, -1});\nauto block_sym = Spinhalf(N, group, irrep);\nXDIAG_SHOW(block_sym);\n\n// with symmetries and Sz\nauto block_sym_sz = Spinhalf(N, nup, group, irrep);\nXDIAG_SHOW(block_sym_sz);\n\nXDIAG_SHOW(block_sym_sz.n_sites());\nXDIAG_SHOW(block_sym_sz.size());\n\n// Iteration\nfor (auto pstate : block_sym_sz) {\n Log(\"{} {}\", to_string(pstate), block_sym_sz.index(pstate));\n}\nXDIAG_SHOW(block_sym_sz.permutation_group());\nXDIAG_SHOW(block_sym_sz.irrep());\nRepresentation of a block in a \\(t-J\\) type Hilbert space.
Source tj.hpp
"},{"location":"documentation/blocks/tJ/#constructors","title":"Constructors","text":"JuliaC++tJ(n_sites::Integer, n_up::Integer, n_dn::Integer)\ntJ(n_sites::Integer, n_up::Integer, n_dn::Integer, \n group::PermutationGroup, irrep::Representation)\ntJ(int64_t n_sites, int64_t n_up, int64_t n_dn);\ntJ(int64_t n_sites, int64_t n_up, int64_t n_dn, \n PermutationGroup group, Representation irrep);\nAn tJ block can be iterated over, where at each iteration a ProductState representing the corresponding basis state is returned.
JuliaC++block = tJ(4, 2, 1)\nfor pstate in block\n @show pstate, index(block, pstate) \nend\nauto block = tJ(4, 2, 1);\nfor (auto pstate : block) {\n Log(\"{} {}\", to_string(pstate), block.index(pstate));\n}\nindex
Returns the index of a given ProductState in the basis of the tJ block.
JuliaC++index(block::tJ, pstate::ProductState)\nint64_t index(ProductState const &pstate) const;\n1-indexing
In the C++ version, the index count starts from \"0\" whereas in Julia the index count starts from \"1\".
n_sites
Returns the number of sites of the block.
JuliaC++n_sites(block::tJ)\nint64_t n_sites() const;\nn_up
Returns the number of \"up\" electrons.
JuliaC++n_up(block::tJ)\nint64_t n_up() const;\nn_dn
Returns the number of \"down\" electrons.
JuliaC++n_dn(block::tJ)\nint64_t n_dn() const;\npermutation_group
Returns the PermutationGroup of the block, if defined.
JuliaC++permutation_group(block::tJ)\nPermutationGroup permutation_group() const;\nirrep
Returns the Representation of the block, if defined.
JuliaC++irrep(block::tJ)\nRepresentation irrep() const;\nsize
Returns the size of the block, i.e. its dimension.
JuliaC++size(block::tJ)\nint64_t size() const;\ndim
Returns the dimension of the block, same as \"size\" for non-distributed blocks.
JuliaC++dim(block::tJ)\nint64_tdim() const;\nisreal
Returns whether the block can be used with real arithmetic. Complex arithmetic is needed when a Representation is genuinely complex.
JuliaC++isreal(block::tJ; precision::Real=1e-12)\nint64_t isreal(double precision = 1e-12) const;\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 n_sites(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)\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({1, -1, 1, -1});\nauto block_sym = tJ(N, nup, ndn, group, irrep);\nXDIAG_SHOW(block_sym);\n\nXDIAG_SHOW(block_sym.n_sites());\nXDIAG_SHOW(block_sym.size());\n\n// Iteration\nfor (auto pstate : block_sym) {\n Log(\"{} {}\", to_string(pstate), block_sym.index(pstate));\n}\nXDIAG_SHOW(block_sym.permutation_group());\nXDIAG_SHOW(block_sym.irrep());\nDescribes the coupling of a local operator. A coupling can either be a string, a real/complex number or even a real/complex matrix. It allows for converting to real/complex numbers or matrices as well as strings, whenever this conversion is sensible.
Source coupling.hpp
"},{"location":"documentation/operators/coupling/#constructors","title":"Constructors","text":"JuliaC++Coupling(name::String)\nCoupling(val::Float64)\nCoupling(val::ComplexF64)\nCoupling(mat::Matrix{Float64})\nCoupling(mat::Matrix{ComplexF64})\nCoupling(std::string value);\nCoupling(double value);\nCoupling(complex value);\nCoupling(arma::mat const &value);\nCoupling(arma::cx_mat const &value);\ntype
Returns the type of the Coupling, i.e. a string which either reads \"string\", \"double\", \"complex\", \"mat\", or \"cx_mat\"
JuliaC++type(cpl::Coupling)\nstd::string type() const;\nisreal
Returns whether or not the coupling is real. Throws an error if the coupling is given as a string, since then it cannot be determined whether the operator is real.
JuliaC++isreal(cpl::Coupling)\nbool isreal() const;\nismatrix
Returns whether or not the coupling is defined as a matrix. Throws an error if the coupling is given as a string, since then it cannot be determined whether the operator is real.
JuliaC++ismatrix(cpl::Coupling)\nbool ismatrix() const;\nisexplicit
Returns false if the coupling is defined as a string, otherwise true
JuliaC++isexplicit(cpl::Coupling)\nbool isexplicit() const;\nA Coupling can be converted to the values it represents, so a string, real/complex number or a real/complex matrix. Initially real values can be cast to complex.
JuliaC++convert(::Type{String}, cpl::Coupling)\nconvert(::Type{Float64}, cpl::Coupling)\nconvert(::Type{ComplexF64}, cpl::Coupling)\nconvert(::Type{Matrix{Float64}}, cpl::Coupling)\nconvert(::Type{Matrix{ComplexF64}}, cpl::Coupling)\ntemplate <typename coeff_t> coeff_t as() const;\ncpl = Coupling(\"J\")\n@show type(cpl)\n@show isexplicit(cpl)\n\ncpl = Coupling(1.23)\n@show ismatrix(cpl)\n@show convert(Float64, cpl)\n@show convert(ComplexF64, cpl)\n\ncpl = Coupling([1 2; -2 1])\n@show ismatrix(cpl)\n@show isreal(cpl)\n@show convert(Matrix{Float64}, cpl)\n@show convert(Matrix{ComplexF64}, cpl)\nauto cpl = Coupling(\"J\");\nXDIAG_SHOW(cpl.type());\nXDIAG_SHOW(cpl.isexplicit());\n\ncpl = Coupling(1.23);\nXDIAG_SHOW(cpl.ismatrix());\nXDIAG_SHOW(cpl.as<double>());\nXDIAG_SHOW(cpl.as<complex>());\n\ncpl = Coupling(arma::mat(\"1 2; -2 1\"));\nXDIAG_SHOW(cpl.ismatrix());\nXDIAG_SHOW(cpl.isreal());\nXDIAG_SHOW(cpl.as<arma::mat>());\nXDIAG_SHOW(cpl.as<arma::cx_mat>());\nA local operator acting on several lattice sites.
Source op.hpp
"},{"location":"documentation/operators/op/#constructors","title":"Constructors","text":"JuliaC++Op(type::String, coupling::String, sites::Vector{Int64})\nOp(type::String, coupling::String, site::Int64)\n\nOp(type::String, coupling::Float64, sites::Vector{Int64})\nOp(type::String, coupling::Float64, site::Int64)\n\nOp(type::String, coupling::ComplexF64, sites::Vector{Int64})\nOp(type::String, coupling::ComplexF64, site::Int64)\n\nOp(type::String, coupling::Matrix{Float64}, sites::Vector{Int64})\nOp(type::String, coupling::Matrix{Float64}, site::Int64)\n\nOp(type::String, coupling::Matrix{ComplexF64}, sites::Vector{Int64})\nOp(type::String, coupling::Matrix{ComplexF64}, site::Int64)\nOp(std::string type, std::string coupling, std::vector<int64_t> const &sites)\nOp(std::string type, std::string coupling, int64_t site)\n\nOp(std::string type, double coupling, std::vector<int64_t> const &sites)\nOp(std::string type, double coupling, int64_t site)\n\nOp(std::string type, complex coupling, std::vector<int64_t> const &sites)\nOp(std::string type, complex coupling, int64_t site)\n\nOp(std::string type, arma::mat const &coupling, std::vector<int64_t> const &sites)\nOp(std::string type, arma::mat const &coupling, int64_t site)\n\nOp(std::string type, arma::cx_mat const &coupling, std::vector<int64_t> const &sites)\nOp(std::string type, arma::cx_mat const &coupling, int64_t site)\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++.
The coupling can take on several types and allow some flexibility in defining operators.
type Description string the coupling is represented as a string, e.g. \"\\(t\\)\" or \"\\(J\\)\" in a \\(t-J\\) model real/cplx number the actual numerical value of the coupling real/cplx matrix more generic interactions can be specified as matrices"},{"location":"documentation/operators/op/#methods","title":"Methods","text":"type
Returns the type of the operator
JuliaC++type(op::Op)\nstd::string type() const;\ncoupling
Returns the coupling of the operator
JuliaC++coupling(op::Op)\nCoupling const &coupling() const;\nThis returns an object of type Coupling, which can then be converted to an appropriate type.
size
Returns how many sites the operator is defined on
JuliaC++size(op::Op)\nint64_t size() const;\ngetindex / operator[]
Returns the site with the given index.
JuliaC++getindex(op::Op, idx::Int64)\nint64_t operator[](int64_t idx) const;\nsites
Returns all the sites the operator is defined on
JuliaC++sites(op::Op)\nstd::vector<int64_t> const &sites() const;\nisreal
Returns whether or not the coupling is real. Throws an error if the coupling is given as a string, since then it cannot be determined whether the operator is real.
JuliaC++isreal(op::Op)\nbool isreal() const;\nismatrix
Returns whether or not the coupling is defined as a matrix. Throws an error if the coupling is given as a string, since then it cannot be determined whether the operator is real.
JuliaC++ismatrix(op::Op)\nbool ismatrix() const;\nisexplicit
Returns false if the coupling is defined as a string, otherwise true
JuliaC++isexplicit(op::Op)\nbool isexplicit() const;\nop = Op(\"HOP\", \"T\", [1, 2])\n@show op\n@show type(op)\n@show convert(String, coupling(op))\n@show size(op), op[1], op[2]\n@show sites(op) == [1, 2]\n@show isexplicit(op)\n\nop = Op(\"HOP\", 1.23, [1, 2])\n@show op\n@show isreal(op)\n@show ismatrix(op)\n@show isexplicit(op)\n\nop = Op(\"SY\", [0 -im; im 0], 1)\n@show op\n@show isreal(op)\n@show ismatrix(op)\n@show isexplicit(op)\nauto op = Op(\"HOP\", \"T\", {0, 1});\nXDIAG_SHOW(op);\nXDIAG_SHOW(op.type());\nXDIAG_SHOW(op.coupling().as<std::string>());\nXDIAG_SHOW(op.size());\nXDIAG_SHOW(op[0]);\nXDIAG_SHOW(op[1]);\nXDIAG_SHOW(op.isexplicit());\n\n op = Op(\"HOP\", 1.23, {0, 1});\nXDIAG_SHOW(op);\nXDIAG_SHOW(op.isreal());\nXDIAG_SHOW(op.ismatrix());\nXDIAG_SHOW(op.isexplicit());\n\narma::cx_mat m(arma::mat(\"0 0; 0 0\"), arma::mat(\"0 -1; 1 0\"));\nop = Op(\"SY\", m, 0);\nXDIAG_SHOW(op);\nXDIAG_SHOW(op.isreal());\nXDIAG_SHOW(op.ismatrix());\nXDIAG_SHOW(op.isexplicit());\nA sum of local operators acting on several lattice sites.
Source opsum.hpp
"},{"location":"documentation/operators/opsum/#constructor","title":"Constructor","text":"JuliaC++OpSum()\nOpSum(ops::Vector{Op})\nOpSum() = default;\nOpSum(std::vector<Op> const &ops);\nsize
Returns the number of local Op operators.
JuliaC++size(ops::OpSum)\nint64_t size() const;\ndefined
Returns bool whether a coupling (of type string) is defined
JuliaC++defined(ops::OpSum, name::String)\nbool defined(std::string name) const;\nsetindex! / operator[]
Sets a coupling given as a string to a certain numerical value or matrix
JuliaC++Base.setindex!(ops::OpSum, cpl, name::String)\nCoupling &operator[](std::string name);\ngetindex / operator[]
Returns the value of a Coupling defined as a string.
JuliaC++getindex(ops::OpSum, name::String)\nCoupling const &operator[](std::string name) const;\ncouplings
Returns all the possible names of Coupling as a vector of strings
JuliaC++couplings(ops::OpSum)\nstd::vector<std::string> couplings() const;\nisreal
Returns whether or not the OpSum is real. This will throw an error if some Coupling are only defined as a string.
JuliaC++isreal(ops::OpSum)\nbool isreal() const;\nisexplicit
Returns false if there exist a Coupling which is defined as a string, otherwise true.
JuliaC++isexplicit(ops::OpSum)\nbool isexplicit() const;\noperator +
Adds a single Op or a full OpSum
JuliaC+++(ops::OpSum, op2::Op) = OpSum(ops.cxx_opsum + op2.cxx_op)\n+(ops::OpSum, ops2::OpSum) = OpSum(ops.cxx_opsum + ops2.cxx_opsum)\nvoid operator+=(Op const &op);\nvoid operator+=(OpSum const &ops);\nOpSum operator+(Op const &op) const;\nOpSum operator+(OpSum const &ops) const;\n# Define the 1D transverse-field Ising chain\nlet \n N = 12\n J = 1.0\n h = 0.5\n Sx = [0 1; 1 0]\n\n ops = OpSum()\n for i in 1:N\n ops += Op(\"ISING\", \"J\", [i, mod1(i+1, N)])\n ops += Op(\"SX\", h * Sx, i)\n end\n\n ops[\"J\"] = 1.0;\n @show ops\n @show defined(ops, \"J\")\n @show isreal(ops)\n @show isexplicit(ops)\nend\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\nauto ops = OpSum();\nfor (int i = 0; i<N; ++i) {\n ops += Op(\"ISING\", \"J\", {i, (i+1)%N});\n ops += Op(\"SX\", arma::mat(h*Sx), i);\n}\nops[\"J\"] = 1.0;\nXDIAG_SHOW(ops);\nXDIAG_SHOW(ops.defined(\"J\"));\nXDIAG_SHOW(ops.isreal());\nXDIAG_SHOW(ops.isexplicit());\nSymmetrizes an operator with respect to a permutation symmetry group.
Source symmetrize.hpp
JuliaC++symmetrize(op::Op, group::PermutationGroup)\nsymmetrize(op::Op, group::PermutationGroup, irrep::Representation)\nsymmetrize(ops::OpSum, group::PermutationGroup)\nsymmetrize(ops::OpSum, group::PermutationGroup, irrep::Representation)\nOpSum symmetrize(Op const &op, PermutationGroup const &group);\nOpSum symmetrize(Op const &op, PermutationGroup const &group, Representation const &irrep);\nOpSum symmetrize(OpSum const &ops, PermutationGroup const &group);\nOpSum symmetrize(OpSum const &ops, PermutationGroup const &group, Representation const &irrep);\nSymmetrization 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/#usage-example","title":"Usage Example","text":"JuliaC++let\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\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({1, 1, 1, 1});\nauto block_sym = Spinhalf(N, group, rep);\n\nauto ops = OpSum();\nfor (int i=0; i<N; ++i) {\n ops += Op(\"HB\", 1.0, {i, (i+1)%N});\n}\nauto [e0, psi] = eig0(ops, block);\nauto [e0s, psi_sym] = eig0(ops, block_sym);\n\nauto corr = Op(\"HB\", 1.0, {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);\nSource create_state.hpp
"},{"location":"documentation/states/create_state/#product","title":"product","text":"Creates a filled product state.
JuliaC++product(block::Block, local_states::Vector{String}, real::Bool=true)\nState product(Block const &block, std::vector<std::string> const &local_state, bool real = true);\nCreate a filled random state with normal distributed coefficients.
JuliaC++rand(block::Block, real::Bool=true, seed::Int64=42, normalized::Bool=true\nState rand(Block const &block, bool real = true, int64_t seed = 42, bool normalized = true);\nCreate a filled state with all zero entries.
JuliaC++zeros(block::Block, real::Bool=true, n_col::Int64=1)\nState zeros(Block const &block, bool real = true, int64_t n_cols = 1);\nSet all coefficients of a given state to zero.
JuliaC++zero(state::State)\nvoid zero(State &state);\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))\nauto block = Spinhalf(2);\nauto state = product(block, {\"Up\", \"Dn\"});\nXDIAG_SHOW(state.vector());\n\nzero(state);\nXDIAG_SHOW(state.vector());\n\nstate = rand(block, false, 1234, true);\nXDIAG_SHOW(state.vectorC());\n\nstate = zeros(block, true, 2);\nXDIAG_SHOW(state.vector());\nFills a State with a given model state, e.g. a ProductState or a RandomState.
Source fill.hpp
"},{"location":"documentation/states/fill/#definition","title":"Definition","text":"JuliaC++fill(state::State, pstate::ProductState, ncol::Int64 = 1)\nfill(state::State, rstate::RandomState, ncol::Int64 = 1)\nvoid fill(State &state, ProductState const &pstate, int64_t ncol = 0);\nvoid fill(State &state, RandomState const &rstate, int64_t ncol = 0);\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))\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());\nA product state of local configurations
Source product_state.hpp
"},{"location":"documentation/states/product_state/#constructors","title":"Constructors","text":"JuliaC++ProductState(n_sites::Int64)\nProductState(local_states::Vector{String})\nProductState(int64_t n_sites);\nProductState(std::vector<std::string> const &local_states);\nA Product state can be iterated over, where at each iteration the string of the local configuration is retured. Here is an example:
JuliaC++pstate = ProductState([\"Up\", \"Dn\", \"Emp\", \"UpDn\"])\nfor s in pstate\n @show s\nend\nauto pstate = ProductState({\"Up\", \"Dn\", \"Emp\", \"UpDn\"});\nfor (auto s : pstate) {\n Log(\"{}\", s);\n}\nn_sites
Returns the number of sites of the product state
JuliaC++n_sites(state::ProductState)\nint64_t n_sites() const\nsize
Returns the number of sites of the product state. Same as \"n_sites\".
JuliaC++size(state::ProductState)\nint64_t size() const\nsetindex! / operator[]
Sets the local configuration at the given site index to the given string.
JuliaC++setindex!(state::ProductState, local_state::String, idx::Int64)\nstd::string &operator[](int64_t i);\ngetindex / operator[]
Returns the string of the local configuration at the given site index.
JuliaC++getindex(state::ProductState, idx::Int64)\nstd::string const &operator[](int64_t i) const;\npush! / push_back
Adds a local configuration add the end of the product state.
JuliaC++push!(state::ProductState, local_state::String\nvoid push_back(std::string l);\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 n_sites(pstate)\nfor s in pstate\n @show s\nend\n@show pstate\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.n_sites());\nfor (auto s : pstate) {\n Log(\"{}\", s);\n}\nXDIAG_SHOW(to_string(pstate));\nA random state with normal distributed coefficients
Source random_state.hpp
"},{"location":"documentation/states/random_state/#constructors","title":"Constructors","text":"JuliaC++RandomState(seed::Int64 = 42, normalized::Bool = true)\nRandomState(int64_t seed = 42, bool normalized = true);\nseed
Returns the seed of the random state
JuliaC++seed(state::RandomState)\nint64_t seed() const;\nsize
Returns whether the state is normalized.
JuliaC++normalized(state::RandomState)\nbool normalized() const;\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
Source state.hpp
"},{"location":"documentation/states/state/#constructors","title":"Constructors","text":"JuliaC++State(block::Block; real::Bool = true, n_cols::Int64 = 1)\nState(block::Block, vec::Vector{Float64})\nState(block::Block, vec::Vector{ComplexF64})\nState(block::Block, mat::Matrix{Float64})\nState(block::Block, mat::Matrix{ComplexF64})\nState(Block const &block, bool real = true, int64_t n_cols = 1);\n\ntemplate <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);\nn_sites
Returns the number of sites of the block the state is defined on.
JuliaC++n_sites(state::State)\nint64_t n_sites() const\nisreal
Returns whether the state is real.
JuliaC++isreal(state::State)\nint64_t isreal() const;\nreal
Returns whether the real part of the State.
JuliaC++real(state::State)\nState real() const;\nimag
Returns whether the imaginary part of the State.
JuliaC++imag(state::State)\nState imag() const;\nmake_complex! / make_complex
Turns a real State into a complex State. Does nothing if the state is already complex
JuliaC++make_complex!(state::State)\nvoid make_complex();\ndim
Returns the dimension of the block the state is defined on.
JuliaC++dim(block::Spinhalf)\nint64_t dim() const;\nsize
Returns the size of the block the state is defined on. locally. Same as \"dim\" for non-distributed Blocks but different for distributed blocks.
JuliaC++size(block::Spinhalf)\nint64_t size() const;\nn_rows
Returns number of rows of the local storage. Same as \"size\"
JuliaC++n_rows(block::Spinhalf)\nint64_t n_rows() const;\nn_cols
Returns number of columns of the local storage.
JuliaC++n_cols(block::Spinhalf)\nint64_t n_cols() const;\ncol
Returns 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\".
JuliaC++col(state::State, n::Int64 = 1; copy::Bool = true)\nState col(int64_t n, bool copy = true) const;\nvector/vectorC
Returns a vector from the n-th column of the storage. In C++ use \"vector\"/\"vectorC\" to either get a real or complex vector.
JuliaC++vector(state::State; n::Int64 = 1)\n# no vectorC method in julia\narma::vec vector(int64_t n = 0, bool copy = true) const;\narma::cx_vec vectorC(int64_t n = 0, bool copy = true) const;\nmatrix/matrixC
Returns matrix representing the storage. In C++ use \"matrix\"/\"matrixC\" to either get a real or complex matrix.
JuliaC++matrix(state::State)\n# no matrixC method in julia\narma::vec matrix(bool copy = true) const;\narma::cx_vec matrixC(bool copy = true) const;\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)))\nauto block = Spinhalf(2);\nauto psi1 = State(block, arma::vec(\"1.0 2.0 3.0 4.0\"));\nXDIAG_SHOW(psi1);\nXDIAG_SHOW(psi1.vector());\npsi1.make_complex();\nXDIAG_SHOW(psi1.vectorC());\n\nauto psi2 = State(block, false, 3);\nXDIAG_SHOW(psi2);\nXDIAG_SHOW(psi2.matrixC());\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(psi3.vectorC());\nXDIAG_SHOW(psi3.real().vector());\nXDIAG_SHOW(psi3.imag().vector());\nPermutations of indices or lattice sites
Source permutation.hpp
"},{"location":"documentation/symmetries/permutation/#constructors","title":"Constructors","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.
Permutation(array::Vector{Int64})\nPermutation(std::vector<int64_t> const &array);\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/symmetries/permutation/#methods","title":"Methods","text":"inverse
Computes the inverse permutation.
JuliaC++inverse(perm::Permutation)\n// As a member function\nPermutation inverse() const;\n\n// As a non-member function\nPermutation inverse(Permutation const &p);\n\"*\" operator
Concatenates two permutations by overloading the * operator.
Base.:*(p1::Permutation, p2::Permutation)\nPermutation operator*(Permutation const &p1, Permutation const &p2);\nsize
Returns the size of the permutation, i.e. the number of indices being permuted.
JuliaC++size(perm::Permutation)\n// As a member function\nint64_t size() const;\n\n// As a non-member function\nint64_t size(Permutation const &p);\np1 = Permutation([1, 3, 2, 4])\np2 = Permutation([3, 1, 2, 4])\n\n@show inverse(p1)\n@show p1 * p2\nPermutation p1 = {0, 2, 1, 3};\nPermutation p2 = {2, 0, 1, 3};\n\nXDIAG_SHOW(inverse(p1));\nXDIAG_SHOW(p1*p2);\nA group of permutations
Source permutation_group.hpp
"},{"location":"documentation/symmetries/permutation_group/#constructor","title":"Constructor","text":"Creates an PermutationGroup out of a vector of Permutation objects.
JuliaC++PermutationGroup(permutations::Vector{Permutation})\nPermutationGroup(std::vector<Permutation> const &permutations);\nn_sites
Returns the number of sites on which the permutations of the group acts.
JuliaC++n_sites(group::PermutationGroup)\nint64_t n_sites() const\nsize
Returns the size of the permutation group, i.e. the number permutations
JuliaC++size(group::PermutationGroup)\nint64_t size() const;\ninverse
Given an index of a permutation, it returns the index of the inverse permutation.
JuliaC++inverse(group::PermutationGroup, idx::Integer)\n// As a member function\nint64_t inverse(int64_t sym) const;\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 n_sites(C3)\n@show inverse(C3, 1) # = 2\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.n_sites());\nXDIAG_SHOW(C3.inverse(1)); // = 2\nA (1D) irreducible representation of a finite group.
Source representation.hpp
"},{"location":"documentation/symmetries/representation/#constructors","title":"Constructors","text":"Creates a Representation from a vector of complex numbers
JuliaC++Representation(characters::Vector{<:Number})\nRepresentation(std::vector<complex> const &characters);\nsize
Returns the size of the Representation, i.e. the number of characters.
JuliaC++size(irrep::Representation)\nint64_t size() const;\nisreal
Returns the whether or not the Representation is real, I.E. the characters are real numbers and do not have an imaginary part.
JuliaC++isreal(irrep::Representation; precision::Real=1e-12)\nbool isreal(double precision = 1e-12) const;\n\"*\" operator
Multiplies two Representations by overloading the * operator.
Base.:*(p1::Representation, p2::Representation)\nRepresentation operator*(Representation const &p1, Representation const &p2);\nr1 = Representation([1, -1, 1, -1])\nr2 = Representation([1, 1im, -1, -1im])\n\n@show r1 * r2\nRepresentation r1 = {1, -1, 1, -1};\nRepresentation r2 = {1, 1i, -1, -1i};\n\nXDIAG_SHOW(r1 * r2);\nAlgorithms 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
set_verbosity(level::Integer);\nvoid set_verbosity(int64_t level);\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
JuliaC++set_verbosity(2);\ne0 = eigval0(bonds, block);\nset_verbosity(2);\ndouble e0 = 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. For example, numbers can be formated this way
Log(\"pi is around {:.4f} and the answer is {}\", 3.141592, 42);\nlogger.hpp
"},{"location":"documentation/utilities/timing/","title":"Timing","text":"In 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":"Miscellaneous 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
JuliaC++set_verbosity(level::Int64)\nvoid set_verbosity(int64_t level);\nPrints a nice welcome message containing the version number and git commit used.
JuliaC++say_hello()\nvoid say_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.
print_version()\nvoid print_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 n_sites : 16\n n_up : 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 \"../../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 n_sites, State const &v) {\n int rank;\n MPI_Comm_rank(MPI_COMM_WORLD, &rank);\n for (int i = 0; i < n_sites; ++i) {\n complex sz = innerC(Bond(\"NUMBER\", 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 n_sites = 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(\"NUMBER\", \"MUPLUS\", site);\n } else {\n Log(\"x {} y {} site {} t {} r {} -\", x, y, site, top, right);\n ops += Op(\"NUMBER\", \"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(n_sites, n_sites / 2 - 1, n_sites / 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(n_sites, 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(n_sites, v);\n toc(\"measure\");\n }\n\n MPI_Finalize();\n return EXIT_SUCCESS;\n} catch (std::exception const &e) {\n traceback(e);\n}\nQuick 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}\nEnter 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 \"../../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:
Supporting 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
Source algebra.hpp
"},{"location":"documentation/algebra/algebra/#norm","title":"norm","text":"Computes the 2-norm of a state
JuliaC++norm(state::State)\ndouble norm(State const &v);\nComputes the 1-norm of a state
JuliaC++norm1(state::State)\ndouble norm1(State const &v);\nComputes the \\(\\infty\\)-norm of a state
JuliaC++norminf(state::State)\ndouble norminf(State const &v);\nComputes the dot product between two states. In C++, please use the dotC function if one of the two states is expected to be complex.
JuliaC++dot(v::State, w::State)\ndouble dot(State const &v, State const &w);\ncomplex dotC(State const &v, State const &w);\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.
JuliaC++inner(ops::OpSum, v::State)\ninner(op::Op, v::State)\ndouble inner(OpSum const &ops, State const &v);\ndouble inner(Op const &op, State const &v);\ncomplex innerC(OpSum const &ops, State const &v);\ncomplex innerC(Op const &op, State const &v);\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\nint N = 8;\nauto block = Spinhalf(N, N / 2);\nauto ops = OpSum();\nfor (int i=0; i<N; ++i) {\n ops += Op(\"HB\", 1.0, {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 = rand(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());\nApplies an operator to a state \\(\\vert w \\rangle = O \\vert v\\rangle\\) or block of states \\(\\left( \\vert w_1 \\rangle \\dots \\vert w_M \\rangle \\right) = O \\left( \\vert v_1 \\rangle \\dots \\vert v_M \\rangle \\right)\\).
Source apply.hpp
JuliaC++apply(op::Op, v::State, w::State, precision::Float64 = 1e-12)\napply(ops::OpSum, v::State, w::State, precision::Float64 = 1e-12)\nvoid apply(Op const &op, State const &v, State &w, double precision = 1e-12);\nvoid apply(OpSum const &ops, State const &v, State &w, double precision = 1e-12);\nThe resulting state is handed as the third argument and is overwritten upon exit.
"},{"location":"documentation/algebra/apply/#parameters","title":"Parameters","text":"Name Description ops / op OpSum or Op defining the operator v input state $\\vert v\\rangle $ w output state \\(\\vert w \\rangle = O \\vert v\\rangle\\) precision precision with which checks for zero are performed (default \\(10^{-12}\\))"},{"location":"documentation/algebra/apply/#usage-example","title":"Usage Example","text":"JuliaC++let \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\nint N = 8;\nauto block = Spinhalf(N, N / 2);\nauto ops = OpSum();\nfor (int i=0; i<N; ++i){\n ops += Op(\"HB\", 1.0, {i, (i+1)%N});\n}\nauto [e0, psi] = eig0(ops, block);\n\nauto blockp = Spinhalf(N, N / 2 + 1);\nauto phi = zeros(blockp);\napply(Op(\"S+\", 1.0, 2), psi, phi);\nXDIAG_SHOW(inner(ops, psi));\nXDIAG_SHOW(inner(ops, phi));\nmatrix(ops, block; force_complex=false)\nmatrixC(ops, block) # c++ only\nCreates the full matrix representation of a given OpSum on a block.
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.
Source matrix.hpp
"},{"location":"documentation/algebra/matrix/#parameters","title":"Parameters","text":"Name Description ops OpSum defining the operator block block on which the operator is defined force_complex flag to determine if returned matrix is forced to be complex (Julia only)"},{"location":"documentation/algebra/matrix/#returns","title":"Returns","text":"Type Description matrix matrix representation of opsum on block"},{"location":"documentation/algebra/matrix/#definition","title":"Definition","text":"JuliaC++matrix(ops::Opsum, block::Block; force_complex::Bool=false)\ntemplate <typename block_t>\narma::mat matrix(OpSum const &ops, block_t const &block);\n\ntemplate <typename block_t>\narma::cx_mat matrixC(OpSum const &ops, block_t const &block);\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\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\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 += Op(\"HOP\", \"T\", {i, (i+1) % N});\n}\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({1, -1, 1, -1});\nauto block = Electron(N, nup, ndn, group, irrep);\n\nauto H = matrix(ops, block);\nH.print();\nComputes the groud state energy and the ground state of an operator on a block.
Source sparse_diag.hpp
JuliaC++function 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)\nstd::tuple<double, State> eig0(OpSum const &ops, Block const &block,\n double precision = 1e-12,\n int64_t max_iterations = 1000,\n bool force_complex = false,\n int64_t random_seed = 42);\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\nint 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(\"HB\", \"J\", {i, (i+1) % N});\n}\nops[\"J\"] = 1.0;\nauto [e0, gs] = eig0(ops, block);\nPerforms an iterative eigenvalue calculation building eigenvectors using the Lanczos algorithm
Source eigs_lanczos.hpp
JuliaC++C++ (explicit state initialization)function 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)\neigs_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);\neigs_lanczos_result_t \neigs_lanczos(OpSum const &ops, Block const &block,\n State const &state, int64_t neigenvalues = 1,\n double precision = 1e-12, int64_t max_iterations = 1000,\n bool force_complex = false);\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/eigval0/","title":"eigval0","text":"Computes the groud state energy of an operator on a block.
Source sparse_diag.hpp
JuliaC++function 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)\ndouble 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);\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\nint 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(\"HB\", \"J\", {i, (i+1) % N});\n}\nops[\"J\"] = 1.0;\ndouble e0 = eigval0(ops, block);\nPerforms an iterative eigenvalue calculation using the Lanczos algorithm.
Source eigvals_lanczos.hpp
JuliaC++C++ (explicit state initialization)function 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)\neigvals_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);\neigvals_lanczos_result_t\neigvals_lanczos(OpSum const &ops, Block const &block, \n State const &state, int64_t neigvals = 1,\n double precision = 1e-12, int64_t max_iterations = 1000,\n bool force_complex = false);\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 niterations number of iterations performed criterion string denoting the reason why the algorithm stopped"},{"location":"documentation/blocks/electron/","title":"Electron","text":"Representation of a block in a Electron (spinful fermion) Hilbert space.
Source electron.hpp
"},{"location":"documentation/blocks/electron/#constructors","title":"Constructors","text":"JuliaC++Electron(n_sites::Integer)\nElectron(n_sites::Integer, n_up::Integer, n_dn::Integer)\nElectron(n_sites::Integer, group::PermutationGroup, irrep::Representation)\nElectron(n_sites::Integer, n_up::Integer, n_dn::Integer, \n group::PermutationGroup, irrep::Representation)\nElectron(int64_t n_sites);\nElectron(int64_t n_sites, int64_t n_up, int64_t n_dn);\nElectron(int64_t n_sites, PermutationGroup permutation_group,\n Representation irrep);\nElectron(int64_t n_sites, int64_t n_up, int64_t n_dn, \n PermutationGroup group, Representation irrep);\nAn Electron block can be iterated over, where at each iteration a ProductState representing the corresponding basis state is returned.
JuliaC++block = Electron(4, 2, 2)\nfor pstate in block\n @show pstate, index(block, pstate) \nend\nauto block = Electron(4, 2, 2);\nfor (auto pstate : block) {\n Log(\"{} {}\", to_string(pstate), block.index(pstate));\n}\nindex
Returns the index of a given ProductState in the basis of the Electron block.
JuliaC++index(block::Electron, pstate::ProductState)\nint64_t index(ProductState const &pstate) const;\n1-indexing
In the C++ version, the index count starts from \"0\" whereas in Julia the index count starts from \"1\".
n_sites
Returns the number of sites of the block.
JuliaC++n_sites(block::Electron)\nint64_t n_sites() const;\nn_up
Returns the number of \"up\" electrons.
JuliaC++n_up(block::Electron)\nint64_t n_up() const;\nn_dn
Returns the number of \"down\" electrons.
JuliaC++n_dn(block::Electron)\nint64_t n_dn() const;\npermutation_group
Returns the PermutationGroup of the block, if defined.
JuliaC++permutation_group(block::Electron)\nPermutationGroup permutation_group() const;\nirrep
Returns the Representation of the block, if defined.
JuliaC++irrep(block::Electron)\nRepresentation irrep() const;\nsize
Returns the size of the block, i.e. its dimension.
JuliaC++size(block::Electron)\nint64_t size() const;\ndim
Returns the dimension of the block, same as \"size\" for non-distributed blocks.
JuliaC++dim(block::Electron)\nint64_tdim() const;\nisreal
Returns whether the block can be used with real arithmetic. Complex arithmetic is needed when a Representation is genuinely complex.
JuliaC++isreal(block::Electron; precision::Real=1e-12)\nint64_t isreal(double precision = 1e-12) const;\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 n_sites(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)\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({1, -1, 1, -1});\nauto block_sym = Electron(N, group, irrep);\nXDIAG_SHOW(block_sym);\n\n// with symmetries and number conservation\nauto block_sym_np = Electron(N, nup, ndn, group, irrep);\nXDIAG_SHOW(block_sym_np);\n\nXDIAG_SHOW(block_sym_np.n_sites());\nXDIAG_SHOW(block_sym_np.size());\n\n// Iteration\nfor (auto pstate : block_sym_np) {\n Log(\"{} {}\", to_string(pstate), block_sym_np.index(pstate));\n}\nXDIAG_SHOW(block_sym_np.permutation_group());\nXDIAG_SHOW(block_sym_np.irrep());\nRepresentation of a block in a spin \\(S=1/2\\) Hilbert space.
Source spinhalf.hpp
"},{"location":"documentation/blocks/spinhalf/#constructors","title":"Constructors","text":"JuliaC++Spinhalf(n_sites::Integer)\nSpinhalf(n_sites::Integer, n_up::Integer)\nSpinhalf(n_sites::Integer, group::PermutationGroup, irrep::Representation)\nSpinhalf(n_sites::Integer, n_up::Integer, group::PermutationGroup, \n irrep::Representation)\nSpinhalf(int64_t n_sites);\nSpinhalf(int64_t n_sites, int64_t n_up);\nSpinhalf(int64_t n_sites, PermutationGroup permutation_group,\n Representation irrep);\nSpinhalf(int64_t n_sites, int64_t n_up, PermutationGroup group,\n Representation irrep);\nAn Spinhalf block can be iterated over, where at each iteration a ProductState representing the corresponding basis state is returned.
JuliaC++block = Spinhalf(4, 2)\nfor pstate in block\n @show pstate, index(block, pstate) \nend\nauto block = Spinhalf(4, 2);\nfor (auto pstate : block) {\n Log(\"{} {}\", to_string(pstate), block.index(pstate));\n}\nindex
Returns the index of a given ProductState in the basis of the Spinhalf block.
JuliaC++index(block::Spinhalf, pstate::ProductState)\nint64_t index(ProductState const &pstate) const;\n1-indexing
In the C++ version, the index count starts from \"0\" whereas in Julia the index count starts from \"1\".
n_sites
Returns the number of sites of the block.
JuliaC++n_sites(block::Spinhalf)\nint64_t n_sites() const;\nn_up
Returns the number of \"up\" spins.
JuliaC++n_up(block::Spinhalf)\nint64_t n_up() const;\npermutation_group
Returns the PermutationGroup of the block, if defined.
JuliaC++permutation_group(block::Spinhalf)\nPermutationGroup permutation_group() const;\nirrep
Returns the Representation of the block, if defined.
JuliaC++irrep(block::Spinhalf)\nRepresentation irrep() const;\nsize
Returns the size of the block, i.e. its dimension.
JuliaC++size(block::Spinhalf)\nint64_t size() const;\ndim
Returns the dimension of the block, same as \"size\" for non-distributed blocks.
JuliaC++dim(block::Spinhalf)\nint64_tdim() const;\nisreal
Returns whether the block can be used with real arithmetic. Complex arithmetic is needed when a Representation is genuinely complex.
JuliaC++isreal(block::Spinhalf; precision::Real=1e-12)\nint64_t isreal(double precision = 1e-12) const;\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 n_sites(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)\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({1, -1, 1, -1});\nauto block_sym = Spinhalf(N, group, irrep);\nXDIAG_SHOW(block_sym);\n\n// with symmetries and Sz\nauto block_sym_sz = Spinhalf(N, nup, group, irrep);\nXDIAG_SHOW(block_sym_sz);\n\nXDIAG_SHOW(block_sym_sz.n_sites());\nXDIAG_SHOW(block_sym_sz.size());\n\n// Iteration\nfor (auto pstate : block_sym_sz) {\n Log(\"{} {}\", to_string(pstate), block_sym_sz.index(pstate));\n}\nXDIAG_SHOW(block_sym_sz.permutation_group());\nXDIAG_SHOW(block_sym_sz.irrep());\nRepresentation of a block in a \\(t-J\\) type Hilbert space.
Source tj.hpp
"},{"location":"documentation/blocks/tJ/#constructors","title":"Constructors","text":"JuliaC++tJ(n_sites::Integer, n_up::Integer, n_dn::Integer)\ntJ(n_sites::Integer, n_up::Integer, n_dn::Integer, \n group::PermutationGroup, irrep::Representation)\ntJ(int64_t n_sites, int64_t n_up, int64_t n_dn);\ntJ(int64_t n_sites, int64_t n_up, int64_t n_dn, \n PermutationGroup group, Representation irrep);\nAn tJ block can be iterated over, where at each iteration a ProductState representing the corresponding basis state is returned.
JuliaC++block = tJ(4, 2, 1)\nfor pstate in block\n @show pstate, index(block, pstate) \nend\nauto block = tJ(4, 2, 1);\nfor (auto pstate : block) {\n Log(\"{} {}\", to_string(pstate), block.index(pstate));\n}\nindex
Returns the index of a given ProductState in the basis of the tJ block.
JuliaC++index(block::tJ, pstate::ProductState)\nint64_t index(ProductState const &pstate) const;\n1-indexing
In the C++ version, the index count starts from \"0\" whereas in Julia the index count starts from \"1\".
n_sites
Returns the number of sites of the block.
JuliaC++n_sites(block::tJ)\nint64_t n_sites() const;\nn_up
Returns the number of \"up\" electrons.
JuliaC++n_up(block::tJ)\nint64_t n_up() const;\nn_dn
Returns the number of \"down\" electrons.
JuliaC++n_dn(block::tJ)\nint64_t n_dn() const;\npermutation_group
Returns the PermutationGroup of the block, if defined.
JuliaC++permutation_group(block::tJ)\nPermutationGroup permutation_group() const;\nirrep
Returns the Representation of the block, if defined.
JuliaC++irrep(block::tJ)\nRepresentation irrep() const;\nsize
Returns the size of the block, i.e. its dimension.
JuliaC++size(block::tJ)\nint64_t size() const;\ndim
Returns the dimension of the block, same as \"size\" for non-distributed blocks.
JuliaC++dim(block::tJ)\nint64_tdim() const;\nisreal
Returns whether the block can be used with real arithmetic. Complex arithmetic is needed when a Representation is genuinely complex.
JuliaC++isreal(block::tJ; precision::Real=1e-12)\nint64_t isreal(double precision = 1e-12) const;\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 n_sites(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)\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({1, -1, 1, -1});\nauto block_sym = tJ(N, nup, ndn, group, irrep);\nXDIAG_SHOW(block_sym);\n\nXDIAG_SHOW(block_sym.n_sites());\nXDIAG_SHOW(block_sym.size());\n\n// Iteration\nfor (auto pstate : block_sym) {\n Log(\"{} {}\", to_string(pstate), block_sym.index(pstate));\n}\nXDIAG_SHOW(block_sym.permutation_group());\nXDIAG_SHOW(block_sym.irrep());\nDescribes the coupling of a local operator. A coupling can either be a string, a real/complex number or even a real/complex matrix. It allows for converting to real/complex numbers or matrices as well as strings, whenever this conversion is sensible.
Source coupling.hpp
"},{"location":"documentation/operators/coupling/#constructors","title":"Constructors","text":"JuliaC++Coupling(name::String)\nCoupling(val::Float64)\nCoupling(val::ComplexF64)\nCoupling(mat::Matrix{Float64})\nCoupling(mat::Matrix{ComplexF64})\nCoupling(std::string value);\nCoupling(double value);\nCoupling(complex value);\nCoupling(arma::mat const &value);\nCoupling(arma::cx_mat const &value);\ntype
Returns the type of the Coupling, i.e. a string which either reads \"string\", \"double\", \"complex\", \"mat\", or \"cx_mat\"
JuliaC++type(cpl::Coupling)\nstd::string type() const;\nisreal
Returns whether or not the coupling is real. Throws an error if the coupling is given as a string, since then it cannot be determined whether the operator is real.
JuliaC++isreal(cpl::Coupling)\nbool isreal() const;\nismatrix
Returns whether or not the coupling is defined as a matrix. Throws an error if the coupling is given as a string, since then it cannot be determined whether the operator is real.
JuliaC++ismatrix(cpl::Coupling)\nbool ismatrix() const;\nisexplicit
Returns false if the coupling is defined as a string, otherwise true
JuliaC++isexplicit(cpl::Coupling)\nbool isexplicit() const;\nA Coupling can be converted to the values it represents, so a string, real/complex number or a real/complex matrix. Initially real values can be cast to complex.
JuliaC++convert(::Type{String}, cpl::Coupling)\nconvert(::Type{Float64}, cpl::Coupling)\nconvert(::Type{ComplexF64}, cpl::Coupling)\nconvert(::Type{Matrix{Float64}}, cpl::Coupling)\nconvert(::Type{Matrix{ComplexF64}}, cpl::Coupling)\ntemplate <typename coeff_t> coeff_t as() const;\ncpl = Coupling(\"J\")\n@show type(cpl)\n@show isexplicit(cpl)\n\ncpl = Coupling(1.23)\n@show ismatrix(cpl)\n@show convert(Float64, cpl)\n@show convert(ComplexF64, cpl)\n\ncpl = Coupling([1 2; -2 1])\n@show ismatrix(cpl)\n@show isreal(cpl)\n@show convert(Matrix{Float64}, cpl)\n@show convert(Matrix{ComplexF64}, cpl)\nauto cpl = Coupling(\"J\");\nXDIAG_SHOW(cpl.type());\nXDIAG_SHOW(cpl.isexplicit());\n\ncpl = Coupling(1.23);\nXDIAG_SHOW(cpl.ismatrix());\nXDIAG_SHOW(cpl.as<double>());\nXDIAG_SHOW(cpl.as<complex>());\n\ncpl = Coupling(arma::mat(\"1 2; -2 1\"));\nXDIAG_SHOW(cpl.ismatrix());\nXDIAG_SHOW(cpl.isreal());\nXDIAG_SHOW(cpl.as<arma::mat>());\nXDIAG_SHOW(cpl.as<arma::cx_mat>());\nA local operator acting on several lattice sites.
Source op.hpp
"},{"location":"documentation/operators/op/#constructors","title":"Constructors","text":"JuliaC++Op(type::String, coupling::String, sites::Vector{Int64})\nOp(type::String, coupling::String, site::Int64)\n\nOp(type::String, coupling::Float64, sites::Vector{Int64})\nOp(type::String, coupling::Float64, site::Int64)\n\nOp(type::String, coupling::ComplexF64, sites::Vector{Int64})\nOp(type::String, coupling::ComplexF64, site::Int64)\n\nOp(type::String, coupling::Matrix{Float64}, sites::Vector{Int64})\nOp(type::String, coupling::Matrix{Float64}, site::Int64)\n\nOp(type::String, coupling::Matrix{ComplexF64}, sites::Vector{Int64})\nOp(type::String, coupling::Matrix{ComplexF64}, site::Int64)\nOp(std::string type, std::string coupling, std::vector<int64_t> const &sites)\nOp(std::string type, std::string coupling, int64_t site)\n\nOp(std::string type, double coupling, std::vector<int64_t> const &sites)\nOp(std::string type, double coupling, int64_t site)\n\nOp(std::string type, complex coupling, std::vector<int64_t> const &sites)\nOp(std::string type, complex coupling, int64_t site)\n\nOp(std::string type, arma::mat const &coupling, std::vector<int64_t> const &sites)\nOp(std::string type, arma::mat const &coupling, int64_t site)\n\nOp(std::string type, arma::cx_mat const &coupling, std::vector<int64_t> const &sites)\nOp(std::string type, arma::cx_mat const &coupling, int64_t site)\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++.
The coupling can take on several types and allow some flexibility in defining operators.
type Description string the coupling is represented as a string, e.g. \"\\(t\\)\" or \"\\(J\\)\" in a \\(t-J\\) model real/cplx number the actual numerical value of the coupling real/cplx matrix more generic interactions can be specified as matrices"},{"location":"documentation/operators/op/#methods","title":"Methods","text":"type
Returns the type of the operator
JuliaC++type(op::Op)\nstd::string type() const;\ncoupling
Returns the coupling of the operator
JuliaC++coupling(op::Op)\nCoupling const &coupling() const;\nThis returns an object of type Coupling, which can then be converted to an appropriate type.
size
Returns how many sites the operator is defined on
JuliaC++size(op::Op)\nint64_t size() const;\ngetindex / operator[]
Returns the site with the given index.
JuliaC++getindex(op::Op, idx::Int64)\nint64_t operator[](int64_t idx) const;\nsites
Returns all the sites the operator is defined on
JuliaC++sites(op::Op)\nstd::vector<int64_t> const &sites() const;\nisreal
Returns whether or not the coupling is real. Throws an error if the coupling is given as a string, since then it cannot be determined whether the operator is real.
JuliaC++isreal(op::Op)\nbool isreal() const;\nismatrix
Returns whether or not the coupling is defined as a matrix. Throws an error if the coupling is given as a string, since then it cannot be determined whether the operator is real.
JuliaC++ismatrix(op::Op)\nbool ismatrix() const;\nisexplicit
Returns false if the coupling is defined as a string, otherwise true
JuliaC++isexplicit(op::Op)\nbool isexplicit() const;\nop = Op(\"HOP\", \"T\", [1, 2])\n@show op\n@show type(op)\n@show convert(String, coupling(op))\n@show size(op), op[1], op[2]\n@show sites(op) == [1, 2]\n@show isexplicit(op)\n\nop = Op(\"HOP\", 1.23, [1, 2])\n@show op\n@show isreal(op)\n@show ismatrix(op)\n@show isexplicit(op)\n\nop = Op(\"SY\", [0 -im; im 0], 1)\n@show op\n@show isreal(op)\n@show ismatrix(op)\n@show isexplicit(op)\nauto op = Op(\"HOP\", \"T\", {0, 1});\nXDIAG_SHOW(op);\nXDIAG_SHOW(op.type());\nXDIAG_SHOW(op.coupling().as<std::string>());\nXDIAG_SHOW(op.size());\nXDIAG_SHOW(op[0]);\nXDIAG_SHOW(op[1]);\nXDIAG_SHOW(op.isexplicit());\n\n op = Op(\"HOP\", 1.23, {0, 1});\nXDIAG_SHOW(op);\nXDIAG_SHOW(op.isreal());\nXDIAG_SHOW(op.ismatrix());\nXDIAG_SHOW(op.isexplicit());\n\narma::cx_mat m(arma::mat(\"0 0; 0 0\"), arma::mat(\"0 -1; 1 0\"));\nop = Op(\"SY\", m, 0);\nXDIAG_SHOW(op);\nXDIAG_SHOW(op.isreal());\nXDIAG_SHOW(op.ismatrix());\nXDIAG_SHOW(op.isexplicit());\nA sum of local operators acting on several lattice sites.
Source opsum.hpp
"},{"location":"documentation/operators/opsum/#constructor","title":"Constructor","text":"JuliaC++OpSum()\nOpSum(ops::Vector{Op})\nOpSum() = default;\nOpSum(std::vector<Op> const &ops);\nsize
Returns the number of local Op operators.
JuliaC++size(ops::OpSum)\nint64_t size() const;\ndefined
Returns bool whether a coupling (of type string) is defined
JuliaC++defined(ops::OpSum, name::String)\nbool defined(std::string name) const;\nsetindex! / operator[]
Sets a coupling given as a string to a certain numerical value or matrix
JuliaC++Base.setindex!(ops::OpSum, cpl, name::String)\nCoupling &operator[](std::string name);\ngetindex / operator[]
Returns the value of a Coupling defined as a string.
JuliaC++getindex(ops::OpSum, name::String)\nCoupling const &operator[](std::string name) const;\ncouplings
Returns all the possible names of Coupling as a vector of strings
JuliaC++couplings(ops::OpSum)\nstd::vector<std::string> couplings() const;\nisreal
Returns whether or not the OpSum is real. This will throw an error if some Coupling are only defined as a string.
JuliaC++isreal(ops::OpSum)\nbool isreal() const;\nisexplicit
Returns false if there exist a Coupling which is defined as a string, otherwise true.
JuliaC++isexplicit(ops::OpSum)\nbool isexplicit() const;\noperator +
Adds a single Op or a full OpSum
JuliaC+++(ops::OpSum, op2::Op) = OpSum(ops.cxx_opsum + op2.cxx_op)\n+(ops::OpSum, ops2::OpSum) = OpSum(ops.cxx_opsum + ops2.cxx_opsum)\nvoid operator+=(Op const &op);\nvoid operator+=(OpSum const &ops);\nOpSum operator+(Op const &op) const;\nOpSum operator+(OpSum const &ops) const;\n# Define the 1D transverse-field Ising chain\nlet \n N = 12\n J = 1.0\n h = 0.5\n Sx = [0 1; 1 0]\n\n ops = OpSum()\n for i in 1:N\n ops += Op(\"ISING\", \"J\", [i, mod1(i+1, N)])\n ops += Op(\"SX\", h * Sx, i)\n end\n\n ops[\"J\"] = 1.0;\n @show ops\n @show defined(ops, \"J\")\n @show isreal(ops)\n @show isexplicit(ops)\nend\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\nauto ops = OpSum();\nfor (int i = 0; i<N; ++i) {\n ops += Op(\"ISING\", \"J\", {i, (i+1)%N});\n ops += Op(\"SX\", arma::mat(h*Sx), i);\n}\nops[\"J\"] = 1.0;\nXDIAG_SHOW(ops);\nXDIAG_SHOW(ops.defined(\"J\"));\nXDIAG_SHOW(ops.isreal());\nXDIAG_SHOW(ops.isexplicit());\nSymmetrizes an operator with respect to a permutation symmetry group.
Source symmetrize.hpp
JuliaC++symmetrize(op::Op, group::PermutationGroup)\nsymmetrize(op::Op, group::PermutationGroup, irrep::Representation)\nsymmetrize(ops::OpSum, group::PermutationGroup)\nsymmetrize(ops::OpSum, group::PermutationGroup, irrep::Representation)\nOpSum symmetrize(Op const &op, PermutationGroup const &group);\nOpSum symmetrize(Op const &op, PermutationGroup const &group, Representation const &irrep);\nOpSum symmetrize(OpSum const &ops, PermutationGroup const &group);\nOpSum symmetrize(OpSum const &ops, PermutationGroup const &group, Representation const &irrep);\nSymmetrization 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/#usage-example","title":"Usage Example","text":"JuliaC++let\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\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({1, 1, 1, 1});\nauto block_sym = Spinhalf(N, group, rep);\n\nauto ops = OpSum();\nfor (int i=0; i<N; ++i) {\n ops += Op(\"HB\", 1.0, {i, (i+1)%N});\n}\nauto [e0, psi] = eig0(ops, block);\nauto [e0s, psi_sym] = eig0(ops, block_sym);\n\nauto corr = Op(\"HB\", 1.0, {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);\nSource create_state.hpp
"},{"location":"documentation/states/create_state/#product","title":"product","text":"Creates a filled product state.
JuliaC++product(block::Block, local_states::Vector{String}, real::Bool=true)\nState product(Block const &block, std::vector<std::string> const &local_state, bool real = true);\nCreate a filled random state with normal distributed coefficients.
JuliaC++rand(block::Block, real::Bool=true, seed::Int64=42, normalized::Bool=true\nState rand(Block const &block, bool real = true, int64_t seed = 42, bool normalized = true);\nCreate a filled state with all zero entries.
JuliaC++zeros(block::Block, real::Bool=true, n_col::Int64=1)\nState zeros(Block const &block, bool real = true, int64_t n_cols = 1);\nSet all coefficients of a given state to zero.
JuliaC++zero(state::State)\nvoid zero(State &state);\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))\nauto block = Spinhalf(2);\nauto state = product(block, {\"Up\", \"Dn\"});\nXDIAG_SHOW(state.vector());\n\nzero(state);\nXDIAG_SHOW(state.vector());\n\nstate = rand(block, false, 1234, true);\nXDIAG_SHOW(state.vectorC());\n\nstate = zeros(block, true, 2);\nXDIAG_SHOW(state.vector());\nFills a State with a given model state, e.g. a ProductState or a RandomState.
Source fill.hpp
"},{"location":"documentation/states/fill/#definition","title":"Definition","text":"JuliaC++fill(state::State, pstate::ProductState, ncol::Int64 = 1)\nfill(state::State, rstate::RandomState, ncol::Int64 = 1)\nvoid fill(State &state, ProductState const &pstate, int64_t ncol = 0);\nvoid fill(State &state, RandomState const &rstate, int64_t ncol = 0);\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))\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());\nA product state of local configurations
Source product_state.hpp
"},{"location":"documentation/states/product_state/#constructors","title":"Constructors","text":"JuliaC++ProductState(n_sites::Int64)\nProductState(local_states::Vector{String})\nProductState(int64_t n_sites);\nProductState(std::vector<std::string> const &local_states);\nA Product state can be iterated over, where at each iteration the string of the local configuration is retured. Here is an example:
JuliaC++pstate = ProductState([\"Up\", \"Dn\", \"Emp\", \"UpDn\"])\nfor s in pstate\n @show s\nend\nauto pstate = ProductState({\"Up\", \"Dn\", \"Emp\", \"UpDn\"});\nfor (auto s : pstate) {\n Log(\"{}\", s);\n}\nn_sites
Returns the number of sites of the product state
JuliaC++n_sites(state::ProductState)\nint64_t n_sites() const\nsize
Returns the number of sites of the product state. Same as \"n_sites\".
JuliaC++size(state::ProductState)\nint64_t size() const\nsetindex! / operator[]
Sets the local configuration at the given site index to the given string.
JuliaC++setindex!(state::ProductState, local_state::String, idx::Int64)\nstd::string &operator[](int64_t i);\ngetindex / operator[]
Returns the string of the local configuration at the given site index.
JuliaC++getindex(state::ProductState, idx::Int64)\nstd::string const &operator[](int64_t i) const;\npush! / push_back
Adds a local configuration add the end of the product state.
JuliaC++push!(state::ProductState, local_state::String\nvoid push_back(std::string l);\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 n_sites(pstate)\nfor s in pstate\n @show s\nend\n@show pstate\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.n_sites());\nfor (auto s : pstate) {\n Log(\"{}\", s);\n}\nXDIAG_SHOW(to_string(pstate));\nA random state with normal distributed coefficients
Source random_state.hpp
"},{"location":"documentation/states/random_state/#constructors","title":"Constructors","text":"JuliaC++RandomState(seed::Int64 = 42, normalized::Bool = true)\nRandomState(int64_t seed = 42, bool normalized = true);\nseed
Returns the seed of the random state
JuliaC++seed(state::RandomState)\nint64_t seed() const;\nsize
Returns whether the state is normalized.
JuliaC++normalized(state::RandomState)\nbool normalized() const;\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
Source state.hpp
"},{"location":"documentation/states/state/#constructors","title":"Constructors","text":"JuliaC++State(block::Block; real::Bool = true, n_cols::Int64 = 1)\nState(block::Block, vec::Vector{Float64})\nState(block::Block, vec::Vector{ComplexF64})\nState(block::Block, mat::Matrix{Float64})\nState(block::Block, mat::Matrix{ComplexF64})\nState(Block const &block, bool real = true, int64_t n_cols = 1);\n\ntemplate <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);\nn_sites
Returns the number of sites of the block the state is defined on.
JuliaC++n_sites(state::State)\nint64_t n_sites() const\nisreal
Returns whether the state is real.
JuliaC++isreal(state::State)\nint64_t isreal() const;\nreal
Returns whether the real part of the State.
JuliaC++real(state::State)\nState real() const;\nimag
Returns whether the imaginary part of the State.
JuliaC++imag(state::State)\nState imag() const;\nmake_complex! / make_complex
Turns a real State into a complex State. Does nothing if the state is already complex
JuliaC++make_complex!(state::State)\nvoid make_complex();\ndim
Returns the dimension of the block the state is defined on.
JuliaC++dim(block::Spinhalf)\nint64_t dim() const;\nsize
Returns the size of the block the state is defined on. locally. Same as \"dim\" for non-distributed Blocks but different for distributed blocks.
JuliaC++size(block::Spinhalf)\nint64_t size() const;\nn_rows
Returns number of rows of the local storage. Same as \"size\"
JuliaC++n_rows(block::Spinhalf)\nint64_t n_rows() const;\nn_cols
Returns number of columns of the local storage.
JuliaC++n_cols(block::Spinhalf)\nint64_t n_cols() const;\ncol
Returns 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\".
JuliaC++col(state::State, n::Int64 = 1; copy::Bool = true)\nState col(int64_t n, bool copy = true) const;\nvector/vectorC
Returns a vector from the n-th column of the storage. In C++ use \"vector\"/\"vectorC\" to either get a real or complex vector.
JuliaC++vector(state::State; n::Int64 = 1)\n# no vectorC method in julia\narma::vec vector(int64_t n = 0, bool copy = true) const;\narma::cx_vec vectorC(int64_t n = 0, bool copy = true) const;\nmatrix/matrixC
Returns matrix representing the storage. In C++ use \"matrix\"/\"matrixC\" to either get a real or complex matrix.
JuliaC++matrix(state::State)\n# no matrixC method in julia\narma::vec matrix(bool copy = true) const;\narma::cx_vec matrixC(bool copy = true) const;\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)))\nauto block = Spinhalf(2);\nauto psi1 = State(block, arma::vec(\"1.0 2.0 3.0 4.0\"));\nXDIAG_SHOW(psi1);\nXDIAG_SHOW(psi1.vector());\npsi1.make_complex();\nXDIAG_SHOW(psi1.vectorC());\n\nauto psi2 = State(block, false, 3);\nXDIAG_SHOW(psi2);\nXDIAG_SHOW(psi2.matrixC());\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(psi3.vectorC());\nXDIAG_SHOW(psi3.real().vector());\nXDIAG_SHOW(psi3.imag().vector());\nPermutations of indices or lattice sites
Source permutation.hpp
"},{"location":"documentation/symmetries/permutation/#constructors","title":"Constructors","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.
Permutation(array::Vector{Int64})\nPermutation(std::vector<int64_t> const &array);\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/symmetries/permutation/#methods","title":"Methods","text":"inverse
Computes the inverse permutation.
JuliaC++inverse(perm::Permutation)\n// As a member function\nPermutation inverse() const;\n\n// As a non-member function\nPermutation inverse(Permutation const &p);\n\"*\" operator
Concatenates two permutations by overloading the * operator.
Base.:*(p1::Permutation, p2::Permutation)\nPermutation operator*(Permutation const &p1, Permutation const &p2);\nsize
Returns the size of the permutation, i.e. the number of indices being permuted.
JuliaC++size(perm::Permutation)\n// As a member function\nint64_t size() const;\n\n// As a non-member function\nint64_t size(Permutation const &p);\np1 = Permutation([1, 3, 2, 4])\np2 = Permutation([3, 1, 2, 4])\n\n@show inverse(p1)\n@show p1 * p2\nPermutation p1 = {0, 2, 1, 3};\nPermutation p2 = {2, 0, 1, 3};\n\nXDIAG_SHOW(inverse(p1));\nXDIAG_SHOW(p1*p2);\nA group of permutations
Source permutation_group.hpp
"},{"location":"documentation/symmetries/permutation_group/#constructor","title":"Constructor","text":"Creates an PermutationGroup out of a vector of Permutation objects.
JuliaC++PermutationGroup(permutations::Vector{Permutation})\nPermutationGroup(std::vector<Permutation> const &permutations);\nn_sites
Returns the number of sites on which the permutations of the group acts.
JuliaC++n_sites(group::PermutationGroup)\nint64_t n_sites() const\nsize
Returns the size of the permutation group, i.e. the number permutations
JuliaC++size(group::PermutationGroup)\nint64_t size() const;\ninverse
Given an index of a permutation, it returns the index of the inverse permutation.
JuliaC++inverse(group::PermutationGroup, idx::Integer)\n// As a member function\nint64_t inverse(int64_t sym) const;\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 n_sites(C3)\n@show inverse(C3, 1) # = 2\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.n_sites());\nXDIAG_SHOW(C3.inverse(1)); // = 2\nA (1D) irreducible representation of a finite group.
Source representation.hpp
"},{"location":"documentation/symmetries/representation/#constructors","title":"Constructors","text":"Creates a Representation from a vector of complex numbers
JuliaC++Representation(characters::Vector{<:Number})\nRepresentation(std::vector<complex> const &characters);\nsize
Returns the size of the Representation, i.e. the number of characters.
JuliaC++size(irrep::Representation)\nint64_t size() const;\nisreal
Returns the whether or not the Representation is real, I.E. the characters are real numbers and do not have an imaginary part.
JuliaC++isreal(irrep::Representation; precision::Real=1e-12)\nbool isreal(double precision = 1e-12) const;\n\"*\" operator
Multiplies two Representations by overloading the * operator.
Base.:*(p1::Representation, p2::Representation)\nRepresentation operator*(Representation const &p1, Representation const &p2);\nr1 = Representation([1, -1, 1, -1])\nr2 = Representation([1, 1im, -1, -1im])\n\n@show r1 * r2\nRepresentation r1 = {1, -1, 1, -1};\nRepresentation r2 = {1, 1i, -1, -1i};\n\nXDIAG_SHOW(r1 * r2);\nAlgorithms 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
set_verbosity(level::Integer);\nvoid set_verbosity(int64_t level);\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
JuliaC++set_verbosity(2);\ne0 = eigval0(bonds, block);\nset_verbosity(2);\ndouble e0 = 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. For example, numbers can be formated this way
Log(\"pi is around {:.4f} and the answer is {}\", 3.141592, 42);\nlogger.hpp
"},{"location":"documentation/utilities/timing/","title":"Timing","text":"In 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":"Miscellaneous 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
JuliaC++set_verbosity(level::Int64)\nvoid set_verbosity(int64_t level);\nPrints a nice welcome message containing the version number and git commit used.
JuliaC++say_hello()\nvoid say_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.
print_version()\nvoid print_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 n_sites : 16\n n_up : 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 \"../../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 n_sites, State const &v) {\n int rank;\n MPI_Comm_rank(MPI_COMM_WORLD, &rank);\n for (int i = 0; i < n_sites; ++i) {\n complex sz = innerC(Bond(\"NUMBER\", 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 n_sites = 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(\"NUMBER\", \"MUPLUS\", site);\n } else {\n Log(\"x {} y {} site {} t {} r {} -\", x, y, site, top, right);\n ops += Op(\"NUMBER\", \"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(n_sites, n_sites / 2 - 1, n_sites / 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(n_sites, 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(n_sites, v);\n toc(\"measure\");\n }\n\n MPI_Finalize();\n return EXIT_SUCCESS;\n} catch (std::exception const &e) {\n traceback(e);\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:
+