#include <cstdlib>
#include <iostream>

#include <libcommute/libcommute.hpp>

using namespace libcommute;

int main() {

  //
  // Parameters of the system
  //

  // Linear sizes of the lattice
  int const Nx = 10;
  int const Ny = 10;

  // Hopping constant
  double const t = 0.5;
  // Coulomb repulsion
  double const U = 2.0;
  // Electron-phonon coupling constant
  double const g = 0.1;

  // Expression with real coefficients 'H' will represent the Hamiltonian.
  // It is initially set to zero by its default-constructor.
  // Every creation and annihilation operator met in the expression must
  // carry two integer (coordinates of a lattice site) and one string
  // index.
  expression<double,               // type of expression's coefficients
             int, int, std::string // types of operator indices
             > H;

  // The following 'factory' functions make quantum operators with
  // statically typed indices and real coefficients.
  using static_indices::c_dag; // Create an electron
  using static_indices::c;     // Destroy an electron
  using static_indices::n;     // Number of electrons
  using static_indices::a_dag; // Create a phonon
  using static_indices::a;     // Destroy a phonon

  // Are two sites neighbors along the x-axis with periodicity?
  auto neighbors_x = [Nx](int ix, int jx) {
    return std::abs(ix - jx) == 1 || std::abs(ix - jx) == Nx - 1;
  };
  // Are two sites neighbors along the y-axis with periodicity?
  auto neighbors_y = [Ny](int iy, int jy) {
    return std::abs(iy - jy) == 1 || std::abs(iy - jy) == Ny - 1;
  };

  // Hopping terms of H
  for(auto spin : {"up", "down"}) {
    for(int ix = 0; ix < Nx;   ix) {
      for(int iy = 0; iy < Ny;   iy) {
        for(int jx = 0; jx < Nx;   jx) {
          for(int jy = 0; jy < Ny;   jy) {
            // Skip all pairs of lattice sites (ix,iy) and (jx,jy) that are
            // not nearest-neighbors.
            if((neighbors_x(ix, jx) && iy == jy) ||
               (ix == jx && neighbors_y(iy, jy))
            ) {
              // Add a hopping term
              H  = -t * c_dag(ix, iy, spin) * c(jx, jy, spin);
            }
          }
        }
      }
    }
  }

  // Coulomb repulsion terms
  for(int ix = 0; ix < Nx;   ix)
    for(int iy = 0; iy < Ny;   iy) {
      H  = U * n(ix, iy, "up") * n(ix, iy, "down");
    }

  // Energy of phonons
  for(int ix = 0; ix < Nx;   ix)
    for(int iy = 0; iy < Ny;   iy) {
      // The spin index is left blank for bosonic operators
      H  = a_dag(ix, iy, "") * a(ix, iy, "");
    }

  // Electron-phonon coupling
  for(auto spin : {"up", "down"}) {
    for(int ix = 0; ix < Nx;   ix) {
      for(int iy = 0; iy < Ny;   iy) {
        H  = g * n(ix, iy, spin) * (a_dag(ix, iy, "")   a(ix, iy, ""));
      }
    }
  }

  // Total number of terms (monomials) in 'H'.
  std::cout << "Total number of terms in H: " << H.size() << '\n';
  // Is H Hermitian?
  std::cout << "H^\\dagger - H = " << (conj(H) - H) << '\n';

  // Does H commute with N and S_z?
  decltype(H) N, S_z;
  for(int ix = 0; ix < Nx;   ix) {
    for(int iy = 0; iy < Ny;   iy) {
      N  = n(ix, iy, "up")   n(ix, iy, "down");
      S_z  = 0.5 * (n(ix, iy, "up") - n(ix, iy, "down"));
    }
  }
  std::cout << "[H, N] = " << (H * N - N * H) << '\n';
  std::cout << "[H, S_z] = " << (H * S_z - S_z * H) << '\n';

  // Iterate over all terms in 'H' and print those of degree 3.
  //
  // Monomials of degree 3 come from the electron-phonon coupling and
  // are products of two fermionic and one bosonic operators.
  for(auto const& term: H) {
    if(term.monomial.size() == 3) {
      // term.coeff is coefficient in front of the monomial
      std::cout << term.monomial << " => " << term.coeff << '\n';
    }
  }

  return 0;
}

#include <algorithm>
#include <cstdlib>
#include <iostream>
#include <vector>

#include <libcommute/libcommute.hpp>

using namespace libcommute;

int main() {

  // The following 'factory' functions make spin operators with statically typed
  // indices and real coefficients.
  using static_indices::S_p;  // Spin-1/2 raising operator S_ 
  using static_indices::S_m;  // Spin-1/2 lowering operator S_-
  using static_indices::S_z;  // Spin-1/2 operator S_z

  // Expression 'H' will represent the exchange interaction term.
  // Our spin operators will carry one integer index (site 1 or 2).
  auto H = 0.5 * (S_p(1) * S_m(2)   S_m(1) * S_p(2))   S_z(1) * S_z(2);

  // Print 'H'
  std::cout << "H = " << H << '\n';

  // Automatically analyze structure of 'H' and construct a 4-dimensional
  // Hilbert space (direct product of two spin-1/2 spaces).
  auto hs = make_hilbert_space(H);
  std::cout << "dim(hs) = " << hs.dim() << '\n';

  // Construct a 'loperator' object that represents action of expression 'H' on
  // state vectors in the Hilbert space 'hs'.
  auto Hop = make_loperator(H, hs);

  // Here, we will act with 'Hop' on each of the 4 basis states |\psi> in 'hs',
  // |\phi> = Hop |\psi>, and print components of |\phi>. In other words,
  // we are going to construct the matrix representation <\phi|Hop|\psi>.

  // Preallocate state vectors.
  // Other containers, such as Eigen::VectorXd could be used instead.
  std::vector<double> phi(4), psi(4);
  // Iterate over basis states
  for(int i = 0; i < 4;   i) {
    std::fill(phi.begin(), phi.end(), 0);
    std::fill(psi.begin(), psi.end(), 0);
    psi[i] = 1; // 'psi' is i-th basis vector now

    phi = Hop(psi);
    // NB.: It is generally recommended to use the in-place syntax
    //  Hop(psi, phi);
    // as it eliminates a memory allocation needed to store the result.

    std::cout << "H|" << i << "> = ";
    for(int j = 0; j < 4;   j) {
      std::cout << " (" << phi[j] << ")|" << j << ">";
    }
    std::cout << '\n';
  }

  return 0;
}