边缘态能带计算C++代码demo

#define EIGEN_USE_MKL_ALL
#define EIGEN_VECTORIZE_SSE4_2

const double PI = std::acos(-1);
const double rt3 = std::sqrt(3);
const std::complex<double> i1(0.0, 1.0);

Eigen::MatrixXcd Hamiltonian_H0_SU4(double k, int N, double t = 1){
    Eigen::Matrix4cd A;
    A << rt3,     0,       1,  std::exp(i1*k),
	      0,     rt3,      1,    1,
	      1,     1,       rt3,   0,
	     std::exp(i1*k), 1, 0, rt3;
    Eigen::Matrix4cd B = Eigen::Matrix4cd::Zero();
    B(0, 2) = 1;
    B(1, 3) = -1;
    Eigen::MatrixXcd H = Eigen::MatrixXcd::Zero(4*N,4*N);
    Eigen::Matrix4cd Bdiag = B.transpose();
    H.block<4,4>(0,0) = A;
    for(int i = 4; i < 4 * N; i += 4) {
        H.block<4,4>(i,i) = A;
        H.block<4,4>(i - 4,i) = Bdiag;
        H.block<4,4>(i, i - 4) = B;
    }
    return H.block(1, 1, 4 * N - 1, 4 * N - 1);
}


Eigen::MatrixXcd Hamiltonian_Graphene(double k, int N, double t = -1.0) {
    Eigen::Matrix2cd T = Eigen::Matrix2cd::Zero(2,2);
    Eigen::Matrix2cd Dup = Eigen::Matrix2cd::Zero(2,2);
    T(0,1) = 1.0 + std::exp(-i1*k);
    T(1,0) = 1.0 + std::exp(-i1*k);
    Dup(1,0) = 1.0 + 0.0 * i1;
    Eigen::Matrix2cd Ddown = Dup.transpose();
    Eigen::MatrixXcd H = Eigen::MatrixXcd::Zero(2*N,2*N);
    H.block<2,2>(0,0) = T;
    for(int i=2;i<2*N;i+=2){
        H.block<2,2>(i,i) = T;
        H.block<2,2>(i-2,i) = Dup;
        H.block<2,2>(i,i-2) = Ddown;
    }
    return t*H;
}

void solve_H0_band(std::vector<double>& ks, int N) {
    std::ofstream out("H0_band.txt");
    for(auto k : ks) {
        Eigen::MatrixXcd H0 = Hamiltonian_H0_SU4(k, N);
        Eigen::SelfAdjointEigenSolver<Eigen::MatrixXcd> eigensolver(H0);
        out << eigensolver.eigenvalues().transpose() << std::endl;
    }
}