结合之前Matlab设计出来的向量化算法，实现了Hatree-Fork算法

Hatree-Fork计算过程

void Hatree_Fork(std::vector<double> &ks,Eigen::MatrixXd N_up_avg, Eigen::MatrixXd N_down_avg, int ncc){
    auto I = Eigen::MatrixXd::Identity(2*N,2*N);
    int nk = ks.size();
    std::ofstream out("N_avg.txt");
    while(ncc--){
        auto N_up_avg_tmp = Eigen::MatrixXd(N_up_avg);
        auto N_down_avg_tmp = Eigen::MatrixXd(N_down_avg);
        
        for(auto k:ks){
            auto H0 = Hamiltonian0(k);
            auto Hk_up = H0 + U*N_down_avg -0.5*U*I;
            auto Hk_down = H0 + U*N_up_avg -0.5*U*I;
            Eigen::SelfAdjointEigenSolver<Eigen::MatrixXcd> eigensolver_up(Hk_up);
            Eigen::SelfAdjointEigenSolver<Eigen::MatrixXcd> eigensolver_down(Hk_down);
            auto Ek_up = eigensolver_up.eigenvalues();
            auto Ak_up = eigensolver_up.eigenvectors();
            auto Ek_down = eigensolver_down.eigenvalues();
            auto Ak_down = eigensolver_down.eigenvectors();

            
            auto Fermi_up = Fermi(Ek_up).transpose().replicate<2*N,1>();
            auto Fermi_down = Fermi(Ek_down).transpose().replicate<2*N,1>();
            auto Nk_up_tmp = Ak_up.array()*Ak_up.conjugate().array()*Fermi_up.array();
            auto Nk_down_tmp = Ak_down.array()*Ak_down.conjugate().array()*Fermi_down.array();
            auto Nk_up_vector = Eigen::MatrixXd(Nk_up_tmp.rowwise().sum().real());
            auto Nk_down_vector = Eigen::MatrixXd(Nk_down_tmp.rowwise().sum().real());

            auto Nk_up = Eigen::MatrixXd(Nk_up_vector.asDiagonal());
            auto Nk_down = Eigen::MatrixXd(Nk_down_vector.asDiagonal());
            N_up_avg_tmp += Nk_up;
            N_down_avg_tmp += Nk_down;
        }
        N_up_avg = Eigen::MatrixXd(N_up_avg_tmp/nk);
        N_down_avg = Eigen::MatrixXd(N_down_avg_tmp/nk);
    }
    //out<<2*N<<std::endl;
    out<<N_up_avg.diagonal()<<std::endl;
    out<<N_down_avg.diagonal()<<std::endl;
}

向量化的费米函数实现：

Eigen::VectorXd Fermi(Eigen::VectorXd E){
    auto x = beta*(E.array()-mu);
    auto f = (x.exp()+1).pow(-1);
    return Eigen::VectorXd(f);
}

能带计算函数

void calculated_band(std::vector<double>& ks,Eigen::MatrixXd N_up_avg, Eigen::MatrixXd N_down_avg){
    std::ofstream out("Ek.txt");
    auto I = Eigen::MatrixXd::Identity(2*N,2*N);
    int kn = ks.size();
    out<<kn<<std::endl;
    for(auto k:ks){
        auto H0 = Hamiltonian0(k);
        auto Hk_up = H0 + U*N_down_avg-0.5*U*I;
        auto Hk_down = H0 + U*N_up_avg -0.5*U*I;
        Eigen::SelfAdjointEigenSolver<Eigen::MatrixXcd> eigensolver_up(Hk_up);
        Eigen::SelfAdjointEigenSolver<Eigen::MatrixXcd> eigensolver_down(Hk_down);
        out << k << std::endl;
        out << eigensolver_up.eigenvalues().transpose() << std::endl;
        out << eigensolver_down.eigenvalues().transpose() << std::endl;
    }
}

单粒子格林函数实现：

Eigen::MatrixXcd Green_Function(double k, double omega, Eigen::MatrixXd N_avg){
    // 1/(z-H(k))
    auto I = Eigen::MatrixXd::Identity(2*N,2*N);
    auto H0 = Hamiltonian0(k);
    auto Hk = H0 + U*N_avg -0.5*U*I;
    auto z = omega + std::complex<double>(0,delta);
    auto Gk = Eigen::MatrixXcd((z*I-Hk).inverse());
    return Gk;
}