Program Listing for File kondoSC.hpp

Return to documentation for file (models/kondoSC.hpp)

#pragma once
#include "models/fermionBasis.hpp"
#include "nrgcore/qOperator.hpp"
#include "nrgcore/qsymmetry.hpp"
#include "utils/qmatrix.hpp"
#include <algorithm>
#include <cmath>
#include <complex>
#include <cstddef>
#include <iostream>
#include <map>
#include <optional>
#include <string>
#include <vector>
class kondoSC {
  fermionBasis localSCbabsis;
public:
  kondoSC(double JKondo, double spinS, double localDelta)
      : localSCbabsis(2, fermionBasis::spinOnly) {
    createBasis(JKondo, spinS,
                localDelta); // create the basis in nstates x nstates
  }
  [[nodiscard]] std::vector<std::vector<int>> get_basis() const { return n_Q; }
  [[nodiscard]] std::vector<std::vector<double>> get_eigenvaluesQ() const {
    return eigenvalues_Q;
  }
  [[nodiscard]] std::vector<double> get_chi_Q() const { return chi_Q; }
  //
  std::vector<std::vector<double>> eigenvalues_Q;
  std::vector<double>              chi_Q;
  std::vector<std::vector<int>>    n_Q;
  //    ########################################
private:
  void createBasis(double JKondo, double spinS, double localDelta) {
    // TODO(saurabh): check spinS is multiple of 2.
    auto spinDim = static_cast<size_t>(2. * spinS + 1.);
    std::cout << "spinDim: " << spinDim << std::endl;
    qmatrix<double> spinSz(spinDim, spinDim, 0);
    qmatrix<double> spinSplus(spinDim, spinDim, 0);
    // Spin Z
    for (size_t i = 0; i < spinDim; i++) {
      double m     = -spinS + i * 1.0;
      spinSz(i, i) = m;
      if (i != (spinDim - 1)) {
        std::cout << "spinDim: " << i << " " << spinDim << std::endl;
        spinSplus(i, i + 1) = std::sqrt(spinS * (spinS + 1.) - m * (m + 1));
      }
    }
    std::cout << "spinSz: \n" << spinSz << std::endl;
    std::cout << "spinSplus: \n" << spinSplus << std::endl;
    auto spinSminus = spinSplus.cTranspose();
    auto spinSx     = (spinSplus + spinSminus) * 0.5;
    auto spinSy     = (spinSplus - spinSminus) * 0.5; // -i is omitted here.
    // ########################################
    std::cout << "spinSx: " << spinSx << std::endl;
    std::cout << "spinSy: " << spinSy << std::endl;
    std::cout << "spinSz: " << spinSz << std::endl;
    // ########################################
    // create basis
    qmatrix<> fdag = {0, 0, 1, 0};
    qmatrix<> sigz = {1, 0, 0, -1};
    qmatrix<> id2  = {1, 0, 0, 1};
    // Create c_up and c_Down operator
    // Addting and additional operators should be done in the same way.
    auto c_up_dag   = fdag.krDot(id2);
    auto c_Down_dag = sigz.krDot(fdag);
    // Create Wilson Site spin operators
    auto wSpinx = (c_up_dag.dot(c_Down_dag.cTranspose()) +
                   c_Down_dag.dot(c_up_dag.cTranspose())) *
                  0.5;
    auto wSpiny = (c_Down_dag.dot(c_up_dag.cTranspose()) -
                   c_up_dag.dot(c_Down_dag.cTranspose())) *
                  0.5; // -i is omitted here.
    auto wSpinz = (c_up_dag.dot(c_up_dag.cTranspose()) -
                   c_Down_dag.dot(c_Down_dag.cTranspose())) *
                  0.5;
    auto wNtotal =
        (c_up_dag.dot(c_up_dag.cTranspose()) +
         c_Down_dag.dot(
             c_Down_dag.cTranspose())); // Create Hamiltonian spin operator
    // ########################################################
    std::cout << "wSpinx: " << wSpinx << std::endl;
    std::cout << "wSpiny: " << wSpiny << std::endl;
    std::cout << "wSpinz: " << wSpinz << std::endl;
    // ############################################
    auto Hamiltonian =
        (spinSz.krDot(wSpinz) + spinSx.krDot(wSpinx) //
         - spinSy.krDot(wSpiny)) * // Imaginary part is taken care here
            JKondo
        // SC term
        + spinSz.id().krDot(
              (c_up_dag.dot(c_Down_dag) +
               c_Down_dag.cTranspose().dot(c_Down_dag.cTranspose()))) *
              localDelta;
    // End
    // std::cout << "spinSz: " << spinSz << std::endl;
    spinSz = qmatrix(spinSz.krDot(wSpinx.id()));
    // std::cout << "spinSz: " << spinSz << std::endl;
    wSpinz  = spinSx.id().krDot(wSpinz);
    wNtotal = spinSx.id().krDot(wNtotal);
    std::cout << "wSpinz: " << wSpinz * 2 << std::endl;
    c_up_dag   = spinSx.id().krDot(c_up_dag);
    c_Down_dag = spinSx.id().krDot(c_Down_dag);
    fnParticle.clear();
    // Number of particles
    fnParticle.push_back(((spinSz.id() + spinSz * 2.) * 0.5).getdiagonal());
    fnParticle.push_back(((spinSz.id() - spinSz * 2.) * 0.5).getdiagonal());
    fnParticle.push_back(((wNtotal + wSpinz * 2.) * 0.5).getdiagonal());
    fnParticle.push_back(((wNtotal - wSpinz * 2.) * 0.5).getdiagonal());
    //
    std::cout << "fnParticle: " << fnParticle << std::endl;
    // set the symmetries of the system
    // create_QuantumNspinCharge();
    // spin is only quantum number that is conserved
    std::vector<size_t> tm1;
    std::vector<size_t> tm2;
    for (size_t j = 0; j < fnParticle.size(); j++) {
      if (j % 2 == 0) { // Spin up
        tm1.push_back(j);
      } else { // Spin down
        tm2.push_back(j);
      }
    }
    std::vector<std::vector<size_t>> SpinSz = {tm1, tm2};
    create_QuantumSpinOnly(SpinSz);
    // create the quantum numbers
    // create the quantum numbers
    create_Block_structure();
    // ####################################################################
    n_Q = get_unique_Qnumbers();
    // set chi_Q
    chi_Q.clear();
    for (auto ai : n_Q) {
      double t_charge = std::accumulate(ai.begin(), ai.end(), 0);
      chi_Q.push_back(std::pow(-1., t_charge));
    }
    //
    // set foperator
    auto h_blocked = get_block_Hamiltonian(Hamiltonian);
    std::cout << "h_blocked: " << h_blocked << std::endl;
    std::cout << "Hamiltonian: " << Hamiltonian << std::endl;
    // Diagonalize the hamilton
    eigenvalues_Q.clear();
    eigenvalues_Q.resize(n_Q.size(), {});
    for (size_t i = 0; i < n_Q.size(); i++) {
      eigenvalues_Q[i] = (h_blocked.get(i, i)).value()->diag();
    }
    // TODO: rotate the f operator
    // ####################################################################
    f_dag_operator = get_block_operators({c_up_dag, c_Down_dag});
    std::cout << "f_dag_operators: " << f_dag_operator.size() << std::endl;
    std::vector<qOperator> topr(f_dag_operator.size(), qOperator());
    for (size_t ip = 0; ip < f_dag_operator.size(); ip++) {
      for (size_t i = 0; i < n_Q.size(); i++) {
        for (size_t j = 0; j < n_Q.size(); j++) {
          auto tfopr = f_dag_operator[ip].get(i, j);
          if (tfopr) {
            topr[ip].set((h_blocked.get(i, i))
                             .value()
                             ->cTranspose()
                             .dot(*tfopr.value())
                             .dot(*(h_blocked.get(j, j)).value()),
                         i, j);
          }
        }
      }
    }
    f_dag_operator = topr;
  }
  //    ######################################
};