关于Fermi Hubbard模型严格对角化的算法和C++实现

目标：使用严格对角化求解一维Hubbard模型基态

难点：理解算法过程以及程序实现

算法在前面已经介绍过，包括三个步骤

1. 多体波函数State的构建，需要用到位运算。我们直接构建固定粒子数的态。理解二进制表示State。 

2. 哈密顿量矩阵的构建，动能项在实空间，如何对应到多体希尔伯特空间，以及Hubbard项双占构型的构建。这里直接用稀疏矩阵来表示。

3. 基态的求解，这里不需要手动实现算法，直接调用第三方库求解即可。

python实现，numpy + scipy足够。理解python中位运算。然后后两步都很简单，外加Numba优化，很容易获得很强的性能，难点只是Python中的位运算。由于Python是动态语言，一些位运算操作比较有技巧。优化的不好，性能很差。

Fortran实现，理论上这是性能最高的实现方式之一，稀疏矩阵有古老的ARPACK，但是要熟悉Fortran及第三方库的使用。

C++实现，位运算非常好实现，稀疏矩阵使用Eigen3，但是Eigen3暂时不支持对角化求基态操作，使用第三方库，还有一种方案是调用APRACK。那么性能将会非常优秀。

还可以加openMP优化，充分利用多核CPU并行。

#include <vector>
#include <iostream>
#include <unordered_map>
#include <algorithm>
#include <tuple>
#include <chrono>

#include "Eigen/Core"
#include "Eigen/SparseCore"
#include "Spectra/SymEigsSolver.h"
#include "Spectra/MatOp/SparseSymMatProd.h"

using namespace Spectra;

typedef Eigen::SparseMatrix<double> SpMat; 
typedef Eigen::Triplet<double> T;
typedef std::tuple<u_int32_t, u_int32_t, double> Hopping;
typedef std::vector<Hopping> Hoppings;
typedef std::unordered_map<u_int32_t, u_int32_t> StateIdxMap;
typedef std::vector<u_int32_t> Uint32Array;

typedef std::chrono::system_clock::time_point time_point;

const int L = 12;
const int Nup = 6;
const int Ndown = 6;
const double U = 10.0;
const double t = -1;

int numberOf1(u_int32_t n) {
    int cnt = 0;
    while (n != 0) {
        cnt++;
        n &= (n - 1);
    }
    return cnt;
}

uint32_t factorial(uint32_t n) {
    uint32_t res = 1;
    while (n != 0) {
        res *= (n--);
    }
    return res;
}

Hoppings getNNHoppings1D(int n) {
    Hoppings hoppings;
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            if (abs(j - i) == 1) {
                hoppings.push_back({i, j, t});
            }
        }
    }
    return hoppings;
}

Uint32Array getStateWithFixSpinNumber(u_int32_t n) {
    Uint32Array states;
    states.reserve(1 << L);
    for (unsigned int i = 0; i < (1 << L); i++) {
        if (numberOf1(i) == n) {
            states.emplace_back(i);
        }
    }
    return states;
}

Uint32Array getStateWithFixNumber(Uint32Array& spinUpStates, Uint32Array& spinDownStates) {
    Uint32Array states;
    states.reserve(1 << L);
    for (auto spinUpState : spinUpStates) {
        for (auto spinDownState : spinDownStates) {
            states.emplace_back((spinUpState << L) + spinDownState);
        }
    }
    return states;
}

StateIdxMap getStateIdxMap(Uint32Array states) {
    StateIdxMap stateIdxMap;
    stateIdxMap.reserve(1 << L);
    for (u_int32_t i = 0; i < states.size(); i++) {
        stateIdxMap[states[i]] = i;
    }
    return stateIdxMap;
}

u_int32_t doubleOccNum(u_int32_t state) {
    u_int32_t spinUpState = state >> L;
    u_int32_t spinDownState = state & ((1 << L) - 1);
    return numberOf1(spinUpState & spinDownState);
}

int perm(u_int32_t state, u_int32_t siteA, u_int32_t siteB) {
    if (siteA < siteB) std::swap(siteA, siteB);
    return (numberOf1(state & ((1 << siteA) - 1) & ~((1 << (siteB + 1)) - 1)) % 2 ? -1 : 1);
}

SpMat Hamiltonian(Uint32Array& states, StateIdxMap& stateIdxMap, std::vector<Hopping>& hoppings) {
    u_int32_t n = states.size();
    std::vector<T> tripletList;
    tripletList.reserve(L * n);
    for (u_int32_t i = 0; i < n; i++) {
        tripletList.emplace_back(T(i,i,U * doubleOccNum(states[i])));
        for (auto [a, b, t] : hoppings) {
            for (u_int32_t s = 0; s <= 1; s++) {
                if ((states[i] & (1 << (s * L) << b)) && !(states[i] & (1 << (s * L) << a))) {
                    u_int32_t state = states[i] ^ (1 << (s * L) << b) ^ (1 << (s * L) << a);
                    u_int32_t j = stateIdxMap[state];
                    tripletList.emplace_back(T(j, i, t * perm(states[i], a, b)));
                }
            }
        }
    }
    SpMat H(n, n);
    H.setFromTriplets(tripletList.begin(), tripletList.end());
    return H;
}


void printTimeCost(std::string taskName, time_point t0, time_point t1) {
    std::chrono::microseconds duration = std::chrono::duration_cast<std::chrono::microseconds>(t1 - t0);
    double second = double(duration.count()) * std::chrono::microseconds::period::num / std::chrono::microseconds::period::den;
    std::cout << taskName << " Time Cost: " << second << " seconds." << std::endl;
}


int main()
{
    // Constructing State of Hilbert Space with Fix N and Sz
    std::cout <<"Start Constructing State" << std::endl;
    time_point t0 = std::chrono::system_clock::now();

    Uint32Array spinUPStates = getStateWithFixSpinNumber(Nup);
    Uint32Array spinDownStates = getStateWithFixSpinNumber(Ndown);
    Uint32Array states = getStateWithFixNumber(spinUPStates, spinDownStates);
    StateIdxMap statesIdxMap = getStateIdxMap(states);

    time_point t1 = std::chrono::system_clock::now();
    printTimeCost("Constructing State", t0, t1);
    std::cout << "All  States Num: " << (1 << (2 * L)) << std :: endl;
    std::cout << "FixN States Num: " << states.size() << std::endl; 
    std::cout << "-----------------------------" << std::endl;
    
    // Add hopping in real space
    Hoppings hoppings = getNNHoppings1D(L);

    // Constructing Hamiltonian
    std::cout << "Start Constructing Hamiltonian" << std::endl;
    t0 = std::chrono::system_clock::now();
    SpMat H = Hamiltonian(states, statesIdxMap, hoppings);
    t1 = std::chrono::system_clock::now();
    printTimeCost("Constructing Hamiltonian", t0, t1);
    std::cout << "-----------------------------" << std::endl;

    // Solve Ground State of Hamiltonian
    std::cout << "Start Solve Eigenvalue" << std::endl;
    t0 = std::chrono::system_clock::now();
    SparseSymMatProd<double> op(H);
    
    SymEigsSolver<SparseSymMatProd<double>> eigs(op, 1, 4);

    eigs.init();
    int nconv = eigs.compute(SortRule::SmallestAlge);

    Eigen::VectorXd evalues;
    if (eigs.info() == CompInfo::Successful) {
        evalues = eigs.eigenvalues();
    }
    t1 = std::chrono::system_clock::now();
    printTimeCost("Eig", t0, t1);
    std::cout << "Eigenvalues found:\n" << evalues << std::endl;
}