关于计算电子磁化率以及RPA的总结

核心就是这个公式 

求和对整个布里渊区进行。要求只能计算一次能谱（因为能谱可能由第一性原理软件读出来）

目前想到的方法是均匀划分网格，这样给定k和q以后，k+q就一定是算过的。对于正方格子，很简单。对于六角格子，平移成平行四边形。

二维对整个布里渊区计算，时间复杂度O(N^4), 个人电脑C++能接受的大小是200*200。（优点对称性体现的很清楚，缺点是计算速度太慢，关键点看起来不容易）

沿着高对称线计算q，时间复杂度O(N^3), 个人电脑C++ 400*400，大约需要20s，可以接受。（推荐）

 

关于从第一性原理软件结果读取能量本征值，需要处理成画能带的数据格式，然后指定能带条数处理。这里需要k点和能量严格对应好。否则，程序就会错误。

核心计算的C++程序，从Energy,txt计算，输出高对称线上的chi

#include <iostream>
#include <cmath>
#include <vector>
#include <fstream>
#include <iomanip>
#include <complex>
#include <ctime>
using namespace std;

#define PI 3.1415926

vector<double> linspace(double min, double max, int n){
    vector<double> result;
    // vector iterator
    int iterator = 0;
    for (int i = 0; i <= n-2; i++){
        double temp = min + i*(max-min)/(floor((double)n) - 1);
        result.insert(result.begin() + iterator, temp);
        iterator += 1;
    }
    //iterator += 1;
    result.insert(result.begin() + iterator, max);
    return result;
}

double fermi_function(double E,double mu,double T){
    return 1/(exp((E-mu)/T)+1);
}

int main()
{
    int nx, ny, nq;
    double mu = 0;
    double T = 0.001;
    double delta = 0.01;
    //vector<double> kx = linspace(0,2*PI,nx);
    //vector<double> ky = linspace(0,2*PI,ny);
    ifstream infile;
    infile.open("Energy.txt");
    infile>>nx>>ny;
    nq = nx;

    vector<vector<double>> E(nx,vector<double>(ny,0));
    vector<complex<double>> chi(nq);
    vector<double> real_chi(nq);

    for(int i=0;i<nx;i++){
        for(int j=0;j<ny;j++){
            infile>>E[i][j];
        }
    }

    clock_t startTime,endTime;
    startTime = clock();              //计时开始
    cout<<"start run"<<endl;

    /*for(int i=0;i<nx;i++){
        for(int j=0;j<ny;j++){
            E[i][j] = -(cos(kx[i])+cos(ky[j]));
        }
    }*/

    for(int i=0;i<nq;i++){

        for(int j=0;j<nx;j++){
            for(int k=0;k<ny;k++){
                int index_kq_x = (i+j)%nx;
                int index_kq_y = (i*0+k)%ny;
                double Ek = E[j][k];
                double Ekq = E[index_kq_x][index_kq_y];
                double f1 = fermi_function(Ek,0,T);
                double f2 = fermi_function(Ekq,0,T);
                chi[i] += (f1-f2)/(Ek-Ekq+1i*delta);
            }
        }
    }


    ofstream out("Datachi.txt");
    for(int i=0;i<nq;i++){
        out<<fixed<<setprecision(4)<<chi[i].real()<<endl;
    }

    endTime = clock();
    cout << "The run time is: " <<(double)(endTime - startTime) / CLOCKS_PER_SEC << "s" << endl;
}

Python 读取画图并作RPA计算程序，做了RPA以后，峰值会放大

import numpy as np
import matplotlib.pyplot as plt
from math import pi

file = open("Datachi.txt", "r")
row = file.readlines()

list_text = []
for line in row:
    line = list(line.strip().split(' '))
    s = []
    for i in line:
        s.append(float(i))
    list_text.append(s)
#print(list_text)

nx = 400
chi = -np.array(list_text)/(nx*nx)
U = 1
chi = chi/(1-U*chi)
plt.plot(chi)
plt.show()