C++ Eigen3实现Hatree

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

Hatree-Fork计算过程

1void Hatree_Fork(std::vector<double> &ks,Eigen::MatrixXd N_up_avg, Eigen::MatrixXd N_down_avg, int ncc){
2    auto I = Eigen::MatrixXd::Identity(2*N,2*N);
3    int nk = ks.size();
4    std::ofstream out("N_avg.txt");
5    while(ncc--){
6        auto N_up_avg_tmp = Eigen::MatrixXd(N_up_avg);
7        auto N_down_avg_tmp = Eigen::MatrixXd(N_down_avg);
8        
9        for(auto k:ks){
10            auto H0 = Hamiltonian0(k);
11            auto Hk_up = H0 + U*N_down_avg -0.5*U*I;
12            auto Hk_down = H0 + U*N_up_avg -0.5*U*I;
13            Eigen::SelfAdjointEigenSolver<Eigen::MatrixXcd> eigensolver_up(Hk_up);
14            Eigen::SelfAdjointEigenSolver<Eigen::MatrixXcd> eigensolver_down(Hk_down);
15            auto Ek_up = eigensolver_up.eigenvalues();
16            auto Ak_up = eigensolver_up.eigenvectors();
17            auto Ek_down = eigensolver_down.eigenvalues();
18            auto Ak_down = eigensolver_down.eigenvectors();
19
20            
21            auto Fermi_up = Fermi(Ek_up).transpose().replicate<2*N,1>();
22            auto Fermi_down = Fermi(Ek_down).transpose().replicate<2*N,1>();
23            auto Nk_up_tmp = Ak_up.array()*Ak_up.conjugate().array()*Fermi_up.array();
24            auto Nk_down_tmp = Ak_down.array()*Ak_down.conjugate().array()*Fermi_down.array();
25            auto Nk_up_vector = Eigen::MatrixXd(Nk_up_tmp.rowwise().sum().real());
26            auto Nk_down_vector = Eigen::MatrixXd(Nk_down_tmp.rowwise().sum().real());
27
28            auto Nk_up = Eigen::MatrixXd(Nk_up_vector.asDiagonal());
29            auto Nk_down = Eigen::MatrixXd(Nk_down_vector.asDiagonal());
30            N_up_avg_tmp += Nk_up;
31            N_down_avg_tmp += Nk_down;
32        }
33        N_up_avg = Eigen::MatrixXd(N_up_avg_tmp/nk);
34        N_down_avg = Eigen::MatrixXd(N_down_avg_tmp/nk);
35    }
36    //out<<2*N<<std::endl;
37    out<<N_up_avg.diagonal()<<std::endl;
38    out<<N_down_avg.diagonal()<<std::endl;
39}

向量化的费米函数实现：

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

能带计算函数

1void calculated_band(std::vector<double>& ks,Eigen::MatrixXd N_up_avg, Eigen::MatrixXd N_down_avg){
2    std::ofstream out("Ek.txt");
3    auto I = Eigen::MatrixXd::Identity(2*N,2*N);
4    int kn = ks.size();
5    out<<kn<<std::endl;
6    for(auto k:ks){
7        auto H0 = Hamiltonian0(k);
8        auto Hk_up = H0 + U*N_down_avg-0.5*U*I;
9        auto Hk_down = H0 + U*N_up_avg -0.5*U*I;
10        Eigen::SelfAdjointEigenSolver<Eigen::MatrixXcd> eigensolver_up(Hk_up);
11        Eigen::SelfAdjointEigenSolver<Eigen::MatrixXcd> eigensolver_down(Hk_down);
12        out << k << std::endl;
13        out << eigensolver_up.eigenvalues().transpose() << std::endl;
14        out << eigensolver_down.eigenvalues().transpose() << std::endl;
15    }
16}

单粒子格林函数实现：

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