# 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.

$\hat H = -t\sum_{\langle i,j \rangle, \sigma} c^\dagger_{i,\sigma} c_{j,\sigma} -\mu \sum_{i, \sigma} n_{i,\sigma} + U \sum_i n_{i,\uparrow} n_{i,\downarrow}.$

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
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'
125                                                        Eigen::EigenvaluesOnly);
126
127    std::cout << "10 lowest eigenvalues of the N = 2 sector: "
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'