当前位置:网站首页>关于无相互作用极化率的计算
关于无相互作用极化率的计算
2020-11-09 11:37:00 【osc_22rhv8iu】
首先,我们需要在布里渊区均匀采样,生成一系列均匀的点,由于均匀的封闭性,给定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()
版权声明
本文为[osc_22rhv8iu]所创,转载请带上原文链接,感谢
https://my.oschina.net/u/4272135/blog/4709113
边栏推荐
猜你喜欢
range_sensor_layer
El table dynamic header
Visual Studio (MAC) installation process notes
Recommendation system, in-depth paper analysis gbdt + LR
android studio创建平板模拟器方法
Android Development - service application, timer implementation (thread + service)
nodejs学习笔记(慕课网nodejs从零开发web Server博客项目)
050_ object-oriented
配置交换机Trunk接口流量本地优先转发(集群/堆叠)
Learning notes of nodejs
随机推荐
Graph node classification and message passing
彩虹排序 | 荷兰旗问题
商品管理系统——商品新增本地保存实现部分
做用户,绕不开画像!
典型分布式系统分析:Dynamo
El table dynamic header
使用CopyMemory API出现 尝试读取或写入受保护的内存。这通常指示其他内存已损坏。
Complete set of linked list operations of data structure and algorithm series (3) (go)
Start learning discrete mathematics again
AI fresh student's annual salary has increased to 400000, you can still make a career change now!
From the practice, this paper discusses the problems caused by the inconsistent design of ruby syntax.
基于synchronized锁的深度解析
如何保证消息不被重复消费?(如何保证消息消费的幂等性)
用一种简单的方式实现终端文字粘贴板
Application of cloud gateway equipment on easynts in Xueliang project
Sql分组查询后取每组的前N条记录
搭建全分布式集群全过程
Capture bubbles? Is browser a fish?
1450. 在既定时间做作业的学生人数
nodejs学习笔记(慕课网nodejs从零开发web Server博客项目)