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;
}