关于无相互作用极化率的计算

想让APP和微信一样，可以流畅运行小程序？ | 体验即送大疆、华为、樱桃键盘！>>>

首先，我们需要在布里渊区均匀采样，生成一系列均匀的点，由于均匀的封闭性，给定k和q后，k+q就可以由已经存在的点表示出来。需要注意的是，我们需要对超出范围的k+q的点做平移(用编号取模就行）

之后，就可以完成极化率的计算。相当于这里全部用了离散坐标。

对于二维体系时间复杂度是O(N^4), C++ 1秒可以算10^7次，200*200得格点，C++大概需要160秒的时间，实际用时303.249s

square lattice matlib代码

nx = 100;
ny = 100;
kx = linspace(0,2*pi,nx);
ky = linspace(0,2*pi,ny);
[KX,KY] = meshgrid(kx,ky);
E = -(cos(KX) + cos(KY));
%mesh(KX,KY,E);
mu = 0;
T = 0.001;
delta = 0.01*0;     % 小虚部
chi = zeros(nx,ny);
Ef=0;
%[C,h] = contour(KX,KY,E,[Ef,0,0]);
for i = 1:nx
    for j = 1:ny
        for k=1:nx
            for l=1:ny
                index_kq_x = mod(i+k,nx)+1;
                index_kq_y = mod(j+l,ny)+1;
                Ek = E(k,l);
                Ekq = E(index_kq_x,index_kq_y);
                chi(i,j) = chi(i,j) + (Fermi_funtion(Ek,mu,T)-Fermi_funtion(Ekq,mu,T))/(Ek-Ekq+1i*delta);
            end
        end
    end
end
mesh(KX,KY,-real(chi))

算得非常慢

 

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 = 200, ny = 200;
    double mu = 0;
    double T = 0.026;
    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<vector<complex<double>>> chi(nx,vector<complex<double>>(ny));
    vector<vector<double>> real_chi(nx,vector<double>(ny,0));

    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[i]));
        }
    }

    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<nx;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);
                    if(k!=index_kq_x&&l!=index_kq_y){
                        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;
}

 

Python画图脚本

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 = 10
ny = 10
chi = np.array(list_text)
kx = np.linspace(0,2*pi,nx)
ky = np.linspace(0,2*pi,ny)
KX, KY = np.meshgrid(kx, ky)
#plt.contourf(KX,KY,chi)
#plt.show()
fig = plt.figure(figsize=plt.figaspect(0.5))
ax = fig.add_subplot(1, 2, 1, projection='3d')
surf = ax.plot_surface(KX, KY, chi, rstride=1, cstride=1, cmap=cm.coolwarm,
                       linewidth=0, antialiased=False)
ax.set_zlim(-1.01, 1.01)
fig.colorbar(surf, shrink=0.5, aspect=10)
plt.show()