第一性原理计算二维/三维电子磁化率方法

1. 用第一性原理完成标准的结构优化，静态自洽计算

2.  可以先计算出沿着高对称线的能带。

3. 在布里渊区均匀打点

这里有一个技巧是直接使用到格式分数坐标，然后沿着倒格矢等分就行。这样便于后面的磁化率计算，Python脚本如下：

import numpy as np
import matplotlib.pyplot as plt

G1 = np.array([0  ,1])        # Reciprocal Lattice 1
G2 = np.array([1, 0])        # Reciprocal Lattice 2


n1 = 100                            # Reciprocal Lattice 1 上的格点数
n2 = 100                            # Reciprocal Lattice 2 上的格点数
n = n1*n2                           # 总的格点数

e1 = G1/n1
e2 = G2/n2

KPOS = np.zeros((n1,n2,2))

for i in range(n1):
    for j in range(n2):
        KPOS[i][j][:] = i*e1 + j*e2

KPOS = KPOS.reshape(n,2)

output = open('KPOINTS.txt','w')
for i in range(n):
    #plt.scatter(KPOS[i][0],KPOS[i][1])
    output.write("%f %f %f %f\n" %(KPOS[i][0],KPOS[i][1],0.0,1.0))
#plt.show()
output.closed()

4. 需要用一个脚本将数据取出来，画二维能带，这里用vaspkit 三维能带功能。

5. 最后是一个C++程序，计算磁化率，二维体系O(n^4)，三维体系O(n^6)的时间复杂度，所以用C++

#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;
    double mu = 0;
    double T = 0.001;
    double delta = 0.01;
    ifstream infile;
    infile.open("Energy.txt");
    infile>>nx>>ny;
    //vector<double> kx = linspace(0,2*PI,nx);
    //vector<double> ky = linspace(0,2*PI,ny);
    vector<vector<double>> E(nx,vector<double>(ny,0));
    vector<vector<complex<double>>> chi(nx,vector<complex<double>>(ny));
    vector<vector<double>> real_chi(nx,vector<double>(ny,0));
    
    for(int i=0;i<nx;i++){
        for(int j=0;j<ny;j++){
            infile>>E[i][j];
        }
    }
    infile.close();
    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<nx;i++){
        for(int j=0;j<ny;j++){
            
            for(int k=0;k<nx;k++){
                for(int l=0;l<ny;l++){
                    int index_kq_x = (i+k)%nx;
                    int index_kq_y = (j+l)%ny;
                    double Ek = E[k][l];
                    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][j] += (f1-f2)/(Ek-Ekq+1i*delta);
                }
            }
        }
    }

    ofstream out("Datachi.txt");
    for(int i=0;i<nx;i++){
        for(int j=0;j<ny;j++){
            out<<fixed<<setprecision(4)<<chi[i][j].real()<<" ";
        }
        out<<endl;
    }
    endTime = clock();//计时结束
    cout << "The run time is: " <<(double)(endTime - startTime) / CLOCKS_PER_SEC << "s" << endl;
}