Program Listing for File fermionBasis.hpp¶
↰ Return to documentation for file (models/fermionBasis.hpp)
#pragma once
#include "nrgcore/qOperator.hpp"
#include "utils/qmatrix.hpp"
#include <cmath>
#include <cstddef>
#include <iostream>
#include <map>
#include <nrgcore/qsymmetry.hpp>
#include <ostream>
#include <set>
class fermionBasis {
public:
enum modelSymmetry {
chargeOnly, // Only charge of the system is conserved.
spinOnly,
chargeAndSpin
}; // Add a key for the exact diagonisation
std::vector<qOperator> f_dag_operator; // f_dag operator
std::vector<std::vector<double>> fnParticle;
fermionBasis(size_t dof, modelSymmetry l) {
// create the basis in nstates x nstates
createFermionBasis(dof); // create the basis in nstates x nstates
// create the switch statement for l
switch (l) {
case chargeOnly:
std::cout << "chargeOnly" << std::endl;
create_QuantumNChargeonly(); // create the quantum numbers
break;
case spinOnly:
std::cout << "spinOnly" << std::endl;
{
std::vector<size_t> tm1;
std::vector<size_t> tm2;
for (size_t j = 0; j < dof; 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
} //
break;
case chargeAndSpin:
std::cout << "chargeAndSpin" << std::endl;
create_QuantumNspinCharge();
break;
}
create_Block_structure(); // create the block structure
set_f_dag_operators(); // set the fermion operators
DebugView = false; //
}
[[nodiscard]] auto get_raw_f_dag_operator() const { return fermionOprMat; }
[[nodiscard]] auto get_f_dag_operator() const { return f_dag_operator; }
void create_QuantumNChargeonly() {
// Here We assume that
// charge of the individual spin is conserved
nQuantum.clear();
for (size_t i = 0; i < fnParticle[0].size(); i++) {
int tm1{0};
for (auto &j : fnParticle) {
tm1 += static_cast<int>(j[i]);
// std::cout << fnpr[j][i] << " ";
}
nQuantum.push_back({tm1});
// std::cout << ":Q:" << std::vector({tm1, tm2});
// std::cout << std::endl;
}
}
void create_QuantumNspinCharge() {
// Here We assume that
// charge of the individual spin is conserved
nQuantum.clear();
for (size_t i = 0; i < fnParticle[0].size(); i++) {
int tm1{0};
int tm2{0};
for (size_t j = 0; j < fnParticle.size(); j++) {
if (j % 2 == 0) { // Spin up
tm1 += static_cast<int>(fnParticle[j][i]);
} else { // Spin down
tm2 += static_cast<int>(fnParticle[j][i]);
}
// std::cout << fnpr[j][i] << " ";
}
nQuantum.push_back({tm1, tm2});
// std::cout << ":Q:" << std::vector({tm1, tm2});
// std::cout << std::endl;
}
}
void create_QuantumSpinOnly(std::vector<std::vector<size_t>> &qsymmetry) {
// qsymmetry defines the symmetries of propblem
// qsymmetry.size = 2 . If there is four fermion flavours the qsymmetry is
// (Charge of individual spin is conserved)
// i.e qsymmetry = {{0,2},{1,3}}
// check for empty
if (qsymmetry.empty()) {
std::cout << "Warning: qsymmetry is empty" << std::endl;
// create odd even pair for spin up and down.
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);
}
}
qsymmetry = std::vector<std::vector<size_t>>({tm1, tm2});
}
if (qsymmetry.size() != 2) {
throw std::runtime_error("qsymmetry.size() != 2");
}
nQuantum.clear();
for (size_t i = 0; i < fnParticle[0].size(); i++) {
std::vector<int> tmVec;
{
size_t iq = 0;
int tm1{0}; // up spin
int tm2{0}; // down spin
for (auto jq : qsymmetry[iq]) {
tm1 += static_cast<int>(fnParticle[jq][i]);
}
iq = 1;
for (auto jq : qsymmetry[iq]) {
tm2 += static_cast<int>(fnParticle[jq][i]);
}
tmVec.push_back(tm1 - tm2);
}
nQuantum.push_back(tmVec);
// std::cout << ":Q:" << std::vector({tm1, tm2});
// std::cout << std::endl;
}
}
void createQNumbers(const std::vector<std::vector<size_t>> &qsymmetry) {
// qsymmetry defines the symmetries of propblem
// qsymmetry.size = 2 for a system with charge and conserved
// (Charge of individual spin is conserved)
// i.e qsymmetry = {{0,2},{1,3}}
nQuantum.clear();
for (size_t i = 0; i < fnParticle[0].size(); i++) {
std::vector<int> tmVec;
for (const auto &iq : qsymmetry) {
int tm1{0};
for (auto jq : iq) {
tm1 += static_cast<int>(fnParticle[jq][i]);
}
tmVec.push_back(tm1);
}
nQuantum.push_back(tmVec);
// std::cout << ":Q:" << std::vector({tm1, tm2});
// std::cout << std::endl;
}
}
void create_Block_structure() {
nQBlocks.clear();
unique_Qnumbers.clear();
// unique_Indices of the quantum numbers
std::vector<size_t> unique_Indices;
{
std::set<std::vector<int>> tm;
for (size_t i = 0; i < fnParticle[0].size(); i++) {
if (tm.insert(nQuantum[i]).second) {
unique_Indices.push_back(i);
}
}
}
// sub blocks indices where Quantum numbers are same.
for (auto i : unique_Indices) {
unique_Qnumbers.push_back(nQuantum[i]);
std::vector<size_t> tm1{i};
for (size_t j = i; j < fnParticle[0].size(); j++) {
if (nQuantum[i] == nQuantum[j]) {
if (i != j) {
tm1.push_back(j);
}
}
}
if (DebugView) {
std::cout << "Coupled : " << nQuantum[i] << " : ";
for (auto j : tm1) {
for (auto &jx : fnParticle) {
std::cout << jx[j] << " ";
}
std::cout << std::endl;
}
}
// << tm1 << std::endl;
nQBlocks.push_back(tm1);
}
// std::cout << "nQBlocks : " << nQBlocks.size() << std::endl;
//
}
auto get_unique_Qnumbers() {
if (unique_Qnumbers.empty()) {
std::cout << "Warning: unique_Qnumbers is empty !" << std::endl;
}
return unique_Qnumbers;
}
std::vector<qOperator>
get_block_operators(const std::vector<qmatrix<>> &sys_operators) {
// This function can be used to calculate any other
// operators which are a combinations of the f operators
std::vector<qOperator> block_operators(sys_operators.size(), qOperator());
for (size_t i = 0; i < nQBlocks.size(); i++) {
for (size_t j = 0; j < nQBlocks.size(); j++) {
size_t tdim = nQBlocks[i].size();
size_t tdim_p = nQBlocks[j].size();
for (size_t ipr = 0; ipr < sys_operators.size(); ipr++) {
qmatrix<> fmatrix(tdim, tdim_p, 0);
double tvalue{0};
for (size_t ik = 0; ik < nQBlocks[i].size(); ik++) {
for (size_t ik_p = 0; ik_p < nQBlocks[j].size(); ik_p++) {
fmatrix(ik, ik_p) =
sys_operators[ipr](nQBlocks[i][ik], nQBlocks[j][ik_p]);
tvalue += std::fabs(
sys_operators[ipr](nQBlocks[i][ik], nQBlocks[j][ik_p]));
}
}
if (tvalue > 0) {
block_operators[ipr].set(fmatrix, i, j);
}
}
}
} // End of operator
return block_operators;
}
qOperator get_block_Hamiltonian(const qmatrix<double> &sys_operators) {
// This function can be used to calculate any other
// check the Hamiltonian is Harmitian i.e H == H.T
{
double tvalue = 0;
for (size_t i = 0; i < sys_operators.getcolumn(); i++) {
for (size_t j = i + 1; j < sys_operators.getrow(); j++) {
tvalue += std::fabs(sys_operators(i, j) - sys_operators(j, i));
}
}
if (tvalue > 1e-10) {
throw std::runtime_error("Hamiltonian is not Hermitian\n ");
}
}
// operators which are a combinations of the f operators
// std::cout << "&sys_operators" << sys_operators;
qOperator block_operators;
for (size_t i = 0; i < nQBlocks.size(); i++) {
size_t tdim = nQBlocks[i].size();
qmatrix<> fmatrix(tdim, tdim, 0);
for (size_t ik = 0; ik < nQBlocks[i].size(); ik++) {
for (size_t ik_p = 0; ik_p < nQBlocks[i].size(); ik_p++) {
fmatrix(ik, ik_p) = sys_operators(nQBlocks[i][ik], nQBlocks[i][ik_p]);
}
}
block_operators.set(fmatrix, i, i);
} //
// test all the offdiagonal elements are zero
double tvalue = 0;
for (size_t i = 0; i < nQBlocks.size(); i++) {
for (size_t j = i + 1; j < nQBlocks.size(); j++) {
for (size_t ik = 0; ik < nQBlocks[i].size(); ik++) {
for (size_t ik_p = 0; ik_p < nQBlocks[j].size(); ik_p++) {
tvalue += sys_operators(nQBlocks[i][ik], nQBlocks[j][ik_p]);
}
}
}
}
if (std::abs(tvalue) > 1e-10) {
throw std::runtime_error("Warning: Hamiltonian is not block diagonal \n");
}
//
// / ///////////////////////////////////////
return block_operators;
}
// End of operator
[[nodiscard]] auto get_basis() const { return fnParticle; }
void setDebugMode(bool debug = true) { DebugView = debug; }
void createFermionBasis(size_t ldof) {
qmatrix<double> fdag({0, 0, 1, 0});
qmatrix<double> sigz({1, 0, 0, -1});
qmatrix<double> id2({1, 0, 0, 1});
fermionOprMat.clear();
for (size_t i = 0; i < ldof; i++) {
qmatrix<> prev(1, 1, 1);
for (size_t j = 0; j < ldof - i - 1; j++) {
prev = id2.krDot(prev);
}
qmatrix<> nxt(1, 1, 1);
for (size_t j = 0; j < i; j++) {
nxt = sigz.krDot(nxt);
}
fermionOprMat.push_back(prev.krDot(fdag.krDot(nxt)));
}
fnParticle.clear();
// std::cout << fup << fdw;
for (size_t i = 0; i < ldof; i++) {
fnParticle.push_back(
(fermionOprMat[i].dot(fermionOprMat[i].cTranspose())).getdiagonal());
}
//
}
std::vector<qmatrix<>> fermionOprMat;
void set_f_dag_operators() {
f_dag_operator = get_block_operators(fermionOprMat);
}
private:
bool DebugView{false};
// foperator in the full basis
// no of particles operator
// basis states -- and quantum numbers
// n_up n_down as a quntum number
std::vector<std::vector<int>> nQuantum;
std::vector<std::vector<int>> unique_Qnumbers;
// Blocked no of particles i.e.,
std::vector<std::vector<size_t>> nQBlocks;
};