无相互作用极化率计算——沿着高对称线计算

最新推荐文章于 2024-07-20 22:41:11 发布

wwxy261

最新推荐文章于 2024-07-20 22:41:11 发布

阅读量237

点赞数

分类专栏：算法

本文链接：https://blog.csdn.net/wwxy1995/article/details/109562135

版权

算法专栏收录该内容

3633 篇文章 120 订阅

订阅专栏

求和需要对整个布里渊区求和（是否能通过对称性简化，这里暂时不考虑这个问题）

不需要将整个布里渊区的k点都计算出来，只需要计算沿着高对称线的k点，对于二维体系，此时的时间复杂度降低为O(n^3) ，现在对于200*200*200的体系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 = 400, ny = 400;
    int nq = 400;                                       // 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);
    vector<vector<double>> E(nx,vector<double>(ny,0));
    vector<complex<double>> chi(nq);
    vector<double> real_chi(nq);

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

画图代码

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 = 200
ny = 200
chi = -np.array(list_text)
plt.plot(chi)
plt.show()