N-fermion sector views and exact diagonalization with Eigen 3
In this example we show how to partially diagonalize a large model of lattice fermions. More specifically, we are going to use the single band Fermi-Hubbard model on a 2D square lattice with periodic boundary conditions.
The number of atoms in the considered lattice (cluster) is set to 16, which makes for an intractably large Hilbert space of dimension \(2^{32}\). Since the Hamiltonian of the Hubbard model preserves total numbers of spin-up and spin-down electrons, one can diagonalize the Hamiltonian within a single \(N\)-fermion sector or within a \((N_\uparrow, N_\downarrow)\)-multisector.
In libcommute’s terms, an \(N\)-fermion
sector is a subspace of a full
Hilbert space, which is spanned by all basis states with a fixed total
occupation of fermionic degrees of freedom (FDOF). Similarly, a
multisector is
spanned by all basis states with a fixed occupation \(N_1\) of a subset of
the FDOF \(\{S_1\}\), occupation \(N_2\) of another subset
\(\{S_2\}\) and so on. There can be any number of pairs
\((\{S_i\}, N_i)\) (sectors contributing to the multisector) as long as
all the subsets \(\{S_i\}\) are disjoint.
In our example we consider a moderately sized sector with \(N = 2\) (\(\dim = {32 \choose 2} = 496\)) and a multisector with \(N_\uparrow = 1, N_\downarrow = 1\) (\(\dim = {16 \choose 1}{16 \choose 1} = 256\)).
Note
In general, the N-fermion (multi)sector views and functions do not require a purely fermionic system. The definitions of a sector and a multisector stand in presence of bosons, spins and any other non-fermionic degrees of freedom.
We use Eigen’s real vector container to store
elements of state vectors within the (multi)sector, and the
SelfAdjointEigenSolver class to diagonalize Hamiltonian matrices.
1// 'state_vector_eigen3.hpp' must be included before 'libcommute.hpp'.
2// This way libcommute will know how to deal with Eigen's vector types.
3#include <libcommute/loperator/state_vector_eigen3.hpp>
4
5#include <libcommute/libcommute.hpp>
6
7// For Eigen::SelfAdjointEigenSolver
8#include <Eigen/Eigenvalues>
9
10using namespace libcommute;
11
12int main() {
13
14 //
15 // Parameters of the system
16 //
17
18 // Linear sizes of the lattice
19 int const Nx = 4;
20 int const Ny = 4;
21
22 // Hopping constant
23 double const t = 0.5;
24 // Chemical potential
25 double const mu = 1.0;
26 // Coulomb repulsion
27 double const U = 2.0;
28
29 // Fermi-Hubbard Hamiltonian to be constructed.
30 // Every creation and annihilation operator met in the Hamiltonian
31 // carry two integer indices (coordinates of a lattice site) and
32 // one string index (spin).
33 static_indices::expr_real<int, int, std::string> H;
34
35 // The following functions make quantum operators with
36 // statically typed indices.
37 using static_indices::c_dag; // Create an electron
38 using static_indices::c; // Destroy an electron
39 using static_indices::n; // Number of electrons
40
41 // Are two sites neighbors along the x-axis with periodicity?
42 auto neighbors_x = [](int ix, int jx) {
43 return std::abs(ix - jx) == 1 || std::abs(ix - jx) == Nx - 1;
44 };
45 // Are two sites neighbors along the y-axis with periodicity?
46 auto neighbors_y = [](int iy, int jy) {
47 return std::abs(iy - jy) == 1 || std::abs(iy - jy) == Ny - 1;
48 };
49
50 // Hopping terms of H
51 for(auto spin : {"up", "down"}) {
52 for(int ix = 0; ix < Nx; ++ix) {
53 for(int iy = 0; iy < Ny; ++iy) {
54 for(int jx = 0; jx < Nx; ++jx) {
55 for(int jy = 0; jy < Ny; ++jy) {
56 // Skip all pairs of lattice sites (ix,iy) and (jx,jy)
57 // that are not nearest neighbors.
58 if((neighbors_x(ix, jx) && iy == jy) ||
59 (ix == jx && neighbors_y(iy, jy))) {
60 // Add a hopping term
61 H += -t * c_dag(ix, iy, spin) * c(jx, jy, spin);
62 }
63 }
64 }
65 }
66 }
67 }
68
69 // Chemical potential terms
70 for(int ix = 0; ix < Nx; ++ix)
71 for(int iy = 0; iy < Ny; ++iy) {
72 H += -mu * (n(ix, iy, "up") + n(ix, iy, "down"));
73 }
74
75 // Coulomb repulsion terms
76 for(int ix = 0; ix < Nx; ++ix)
77 for(int iy = 0; iy < Ny; ++iy) {
78 H += U * n(ix, iy, "up") * n(ix, iy, "down");
79 }
80
81 // Automatically analyze structure of 'H' and construct a Hilbert space
82 auto hs = make_hilbert_space(H);
83 std::cout << "Full Hilbert space dimension: " << hs.dim() << '\n';
84
85 // Construct a 'loperator' object that represents action of expression 'H' on
86 // state vectors in the Hilbert space 'hs'.
87 auto Hop = make_loperator(H, hs);
88
89 //
90 // Diagonalize the N = 2 sector of the model using Eigen 3 and the
91 // 'n_fermion_sector_view' class.
92 //
93 {
94 int const N = 2;
95
96 auto const sector_size = n_fermion_sector_size(hs, N);
97 std::cout << "Size of the N = 2 sector: " << sector_size << '\n';
98
99 // Preallocate a Hamiltonian matrix to be filled and diagonalized.
100 // Note that we have to store matrix elements only within the small sector.
101 Eigen::MatrixXd Hmat(sector_size, sector_size);
102
103 // Preallocate an input state vector
104 Eigen::VectorXd psi(sector_size);
105
106 // Create a constant 2-fermion sector view of the source vector 'psi'
107 auto psi_view = make_const_nfs_view(psi, hs, N);
108
109 // Fill 'Hmat'
110 for(int i = 0; i < int(sector_size); ++i) {
111 // Reset vector |\psi> to zero.
112 psi = Eigen::VectorXd::Zero(int(sector_size));
113 psi[i] = 1; // 'psi' is the i-th basis vector now
114
115 // Create a 2-fermion sector view of the i-th column of 'Hmat'.
116 auto HMat_col_view = make_nfs_view(Hmat.col(i), hs, N);
117
118 // Act with Hop on the i-th sector basis state and write the result into
119 // the i-th column of Hmat.
120 Hop(psi_view, HMat_col_view);
121 }
122
123 // Diagonalize 'Hmat'
124 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> diag(Hmat,
125 Eigen::EigenvaluesOnly);
126
127 std::cout << "10 lowest eigenvalues of the N = 2 sector: "
128 << diag.eigenvalues().head(10).transpose() << '\n';
129 }
130
131 //
132 // Diagonalize the N_up = 1, N_down = 1 multisector of the model using
133 // Eigen 3 and the 'n_fermion_multisector_view' class.
134 //
135 {
136 int const N_up = 1, N_down = 1;
137
138 // Define sectors contributing to the multisector
139 sector_descriptor<decltype(hs)> sector_up{{}, N_up};
140 sector_descriptor<decltype(hs)> sector_down{{}, N_down};
141
142 // Fill index sets {S_up} and {S_down}
143 for(int ix = 0; ix < Nx; ++ix) {
144 for(int iy = 0; iy < Ny; ++iy) {
145 sector_up.indices.emplace(ix, iy, "up");
146 sector_down.indices.emplace(ix, iy, "down");
147 }
148 }
149
150 std::vector<sector_descriptor<decltype(hs)>> sectors = {sector_up,
151 sector_down};
152
153 auto const multisector_size = n_fermion_multisector_size(hs, sectors);
154 std::cout << "Size of the N_up = 1, N_down = 1 multisector: "
155 << multisector_size << '\n';
156
157 // Preallocate a Hamiltonian matrix to be filled and diagonalized.
158 // Note that we have to store matrix elements only within the multisector.
159 Eigen::MatrixXd Hmat(multisector_size, multisector_size);
160
161 // Preallocate an input state vector
162 Eigen::VectorXd psi(multisector_size);
163
164 // Create a constant multisector view of the source vector 'psi'
165 auto psi_view = make_const_nfms_view(psi, hs, sectors);
166
167 // Fill 'Hmat'
168 for(int i = 0; i < int(multisector_size); ++i) {
169 // Reset vector |\psi> to zero.
170 psi = Eigen::VectorXd::Zero(int(multisector_size));
171 psi[i] = 1; // 'psi' is the i-th basis vector now
172
173 // Create a multisector view of the i-th column of 'Hmat'.
174 auto HMat_col_view = make_nfms_view(Hmat.col(i), hs, sectors);
175
176 // Act with Hop on the i-th multisector basis state and write the result
177 // into the i-th column of Hmat.
178 Hop(psi_view, HMat_col_view);
179 }
180
181 // Diagonalize 'Hmat'
182 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> diag(Hmat,
183 Eigen::EigenvaluesOnly);
184
185 std::cout << "10 lowest eigenvalues of the multisector: "
186 << diag.eigenvalues().head(10).transpose() << '\n';
187 }
188
189 return 0;
190}