石墨烯能带结构（考虑三近邻），两种方法对比 （1）构造哈密顿 （2）matlab 后处理wannier90hr.dat （晶体规范）

已于 2024-01-14 23:30:48 修改
阅读量675 收藏 11
点赞数 4
文章标签： matlab 数据库 人工智能
于 2024-01-11 02:57:04 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/w126126/article/details/135516698
版权

                                      单层石墨烯能带图。方法一用+表示 ，方法2红色实线表示

(1)根据原胞构建哈密顿 

%%Nie-Wei Wang,2024-1-8 (分解版，有基础)
clear all
%%%%%%%%晶格矢量 正格矢量
a1=[2.4682263781246880    0.0000000000000000    0.0000000000000000];           
a2=[-1.2341131890623440    2.1375467457960933    0.0000000000000000];            
a3=[0.0000000000000000    0.0000000000000000   10.0000000000000000];
% C1=[0.0000018561517373 -0.0000018561604471  0.5000000000000000];
% C2=[0.6666682727452335  0.3333317562640553  0.5000000000000000];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
V=dot(a1,cross(a2,a3));%  %%%晶胞体积

b1=2*pi*cross(a2,a3)/V; %%%%倒格矢量
b2=2*pi*cross(a3,a1)/V;
b3=2*pi*cross(a1,a2)/V;

a1=a1(1:2);                     %晶格矢量a1 二维
a2=a2(1:2);                     %晶格矢量a2
 
b1=b1(1:2);                     %%%倒空间晶格矢量b1，二维
b2=b2(1:2);                     %%%倒格矢b2
 
gamma=[0 0];                    %%%高对称点gamma 
M=1/2*b1;                       %%%高对称点M    
K=1/3*b1+1/3*b2;                %%%高对称点K   


kn=20;     %两个高对称点之间分成有kn个点
% Gamma-M段  每个高对称点都有两个 150和151都是M点
kx(1:kn)=linspace(gamma(1),M(1),kn);
ky(1:kn)=linspace(gamma(2),M(2),kn);
% M-K段
kx(kn+1:2*kn)=linspace(M(1),K(1),kn);
ky(kn+1:2*kn)=linspace(M(2),K(2),kn);
% K-Gamma段
kx(2*kn+1:3*kn)=linspace(K(1),0,kn);
ky(2*kn+1:3*kn)=linspace(K(2),0,kn);
k=zeros(3*kn,2);
k(:,1)=kx;                            %%存储kx
k(:,2)=ky; 
Nw=2; %%%%能带条数需要修改
k_new=zeros(1,3*kn);    %%%为了沿着高对称点画出能带需要将k做调整,按照路径的模长走
energy=zeros(Nw,kn);    %%存放本征值
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 1:3*kn   %%%%高对称路径
    if i<=kn
        kv=[kx(i)-gamma(1),ky(i)-gamma(2)];
        k_new(i)=norm(kv);                        %gamma->M      
    elseif i>kn && i<=2*kn
        kv=[kx(i)-M(1),ky(i)-M(2)];
        k_new(i)=norm(M-gamma)+norm(kv);           %M->K 
    else 
        kv=[kx(i)-K(1),ky(i)-K(2)];
        k_new(i)=norm(kv)+norm(M-gamma)+norm(K-M); %%%K->gamma 
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%路径不要动
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
efermi=0;
e1=-0.91;
e2=e1;
t1=2.877065;  %跃迁积分
t2=0.235328;
t3=0.227017;

Nw=2; %两条能带

h0=[e1,t1;t1,e2];
h1=[0,t1;0 0];
h2=[t3,0;0,t3];
%%%%%%%%%%%%%%%%最近邻
R1=a1;
R2=a1+a2;
R3=a2;
R4=-a1;
R5=-(a1+a2);
R6=-a2;
R7=-a2-2*a1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%k=[k(1,:);k(kn,:);k(2*kn,:)];
for i = 1:3*kn
%%%%%%%%%%%%%%%normal   
%最近邻hR4,hR5  12
%次近邻 hR1-6 11 
%三近邻 hR3 hR7 hR6 12
hNN=(exp(1i*dot(k(i,:),R4))...
   +exp(1i*dot(k(i,:),R5)))*t1;   %%12
h2NN=(exp(1i*dot(k(i,:),R1))...
     +exp(1i*dot(k(i,:),R2))...
     +exp(1i*dot(k(i,:),R3))...
     +exp(1i*dot(k(i,:),R4))...
     +exp(1i*dot(k(i,:),R5))...
     +exp(1i*dot(k(i,:),R6)))*t2; %%11
h3NN=(exp(1i*dot(k(i,:),R3))...
     +exp(1i*dot(k(i,:),R6))...
     +exp(1i*dot(k(i,:),R7)))*t3; %%11
h=hNN+h3NN;
hsum=h0+[0,h;h',0]+[h2NN,0;0,h2NN];
[V,D] = eig(hsum);             %%%计算 A 的特征值D和右特征向量。%https://ww2.mathworks.cn/help/matlab/ref/eig.html
[d,ind]=sort(diag(real(D)));%使用 diag(D) 从 D 的对角线上提取特征值，然后按升序对得到的向量进行排序。
                            % sort 的第二个输出返回索引的置换向量。
Vs = V(:,ind);              %%%对本征矢量进行排序,本征值是一列一列的。
Ds = D(ind,ind);            %使用 ind 对 D 的对角线元素进行重新排序。
                            % 由于 D 中的特征值对应于 V 的各列中的特征向量，因此您还必须使用相同的索引对 V 的列进行重新排序。
energy(:,i)=d;

end

%%%%%%%plot 画图
hold on
for m=1:Nw                                 %%%  原胞中有几个轨道就有几条能带

    plot(k_new,energy(m,:)-efermi,'blue','LineWidth',1)  %%
    hold on
end

 
%%%%设置xticks布里渊区高对称点
xticks([k_new(1) k_new(kn+1)  k_new(2*kn+1)  k_new(3*kn)]);
xticklabels({'\Gamma','M','K','\Gamma'})
 
%%%辅助线
xline(k_new(kn+1),'--');
xline(k_new(2*kn+1),'--');
yline(0,'--');
%grid on
 
%%%
xlim([k_new(1) k_new(3*kn)])
ylim([-10,15])
ylabel("Energy")
box on

（2） 第二种：matlab 后处理 wannier_hr.dat 。

%%%2023/4/18_wnw_hebust
 
% $$\begin{aligned}
% h(\vec{k})=\sum_mH_{nm}\text{e}^{\text{i}\vec{k}\cdot(\vec{d}_m-\vec{d}_n)}
% \end{aligned}$$
 
clear all
 
efermi=0;  %费米能 在求解能带时grep fermi OUTCAR 
 
%在POSCAR 读取晶格矢量a1,a2,a3    
a1=[2.4682263781246880    0.0000000000000000    0.0000000000000000];           
a2=[-1.2341131890623440    2.1375467457960933    0.0000000000000000];            
a3=[0.0000000000000000    0.0000000000000000   10.0000000000000000];
 
b1=2*pi*cross(a2,a3)/dot(a1,cross(a2,a3));  %%%倒空间晶格矢量b1，三维
b2=2*pi*cross(a3,a1)/dot(a2,cross(a3,a1));  %%%倒格矢b2
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 正格矢和倒格矢由三维变成二维
a1=a1(1:2);                     %晶格矢量a1 二维
a2=a2(1:2);                     %晶格矢量a2
 
b1=b1(1:2);                     %%%倒空间晶格矢量b1，二维
b2=b2(1:2);                     %%%倒格矢b2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
gamma=[0,0];                    %%%高对称点gamma 
M=1/2*b1;                       %%%高对称点M    
K=1/3*b1+1/3*b2;                %%%高对称点K    
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
kn=20;     %两个高对称点之间分成有kn个点
% Gamma-M段  每个高对称点都有两个 150和151都是M点
kx(1:kn)=linspace(gamma(1),M(1),kn);
ky(1:kn)=linspace(gamma(2),M(2),kn);
% M-K段
kx(kn+1:2*kn)=linspace(M(1),K(1),kn);
ky(kn+1:2*kn)=linspace(M(2),K(2),kn);
% K-Gamma段
kx(2*kn+1:3*kn)=linspace(K(1),0,kn);
ky(2*kn+1:3*kn)=linspace(K(2),0,kn);
k=zeros(3*kn,2);
k(:,1)=kx;                            %%存储kx
k(:,2)=ky;                            %%存储ky 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%读取'wannier90_hr.dat'文件
S = importdata('wannier90_hr.dat');  %将wannier90_hr.dat放入此代码的同一文件夹内
Data = S.data;                 %按照7个为一组进行排列，竖排。
Nw = Data(1);                  %num_wannier=22
Data(1) = [];                  %这个命令表示从数组Data中删除第一个元素。
Nrpts = Data(1);               %共取了331个R点，存在简并
Data(1) = [];                  %这个命令表示从数组Data中删除第一个元素。
Ndegen = Data(1:Nrpts)';       %提取每个R点的简并
Data(1:Nrpts)=[];              %将每个R点的简并值删掉
H = reshape(Data,7,[])';       % 则 reshape(A,2,2,[]) 将 Data 重构为一个 [[],7] 的数组
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k_new=zeros(1,3*kn);    %%%为了沿着高对称点画出能带需要将k做调整,按照路径的模长走
energy=zeros(Nw,kn);    %%存放本征值
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 1:3*kn 
    if i<=kn
        kv=[kx(i)-gamma(1),ky(i)-gamma(2)];
        k_new(i)=norm(kv);                        %gamma->M      
    elseif i>kn && i<=2*kn
        kv=[kx(i)-M(1),ky(i)-M(2)];
        k_new(i)=norm(M-gamma)+norm(kv);           %M->K 
    else 
        kv=[kx(i)-K(1),ky(i)-K(2)];
        k_new(i)=norm(kv)+norm(M-gamma)+norm(K-M); %%%K->gamma 
    end
   h=zeros(Nw);
   %%%%%%%%%%%%%%a1,a2,a3的系数范围R=l*a1+m*a2+n*a3，n3=0;
for ii=1:Nw^2:length(H(:,1))  %%H的行数
    l=H(ii,1);           %第一列l
    m=H(ii,2);           %第2列m
    temp= find(H(:,1)==l & H(:,2)==m);   %%%查找符合的行号l,m,n==0；
    Htemp= H(min(temp):max(temp),:);     %%%将符合的行写入临时变量Htemp
    Htemp(:,6:7)=Htemp(:,6:7)/Ndegen(max(temp)/Nw^2);%%%将hoping项除以R的简并数
    Hmat = reshape(Htemp(:,6),[Nw,Nw]);              %将每个R点的hopping 实步改成num_wannier^2的矩阵
    %%只考虑最近邻
    for i1=1:Nw
        for i2=1:Nw
             if abs(Hmat(i1,i2))<0.22
                 Hmat(i1,i2)=0;
             end
        end
    end
    h=h+exp(1i*dot(k(i,:),(l*a1+m*a2)))*Hmat;        %倒空间的哈密顿量表示
    
end 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[V,D] = eig(h);             %%%计算 A 的特征值D和右特征向量。%https://ww2.mathworks.cn/help/matlab/ref/eig.html
[d,ind]=sort(diag(real(D)));%使用 diag(D) 从 D 的对角线上提取特征值，然后按升序对得到的向量进行排序。
                            % sort 的第二个输出返回索引的置换向量。
Vs = V(:,ind);              %%%对本征矢量进行排序,本征值是一列一列的。
Ds = D(ind,ind);            %使用 ind 对 D 的对角线元素进行重新排序。
                            % 由于 D 中的特征值对应于 V 的各列中的特征向量，因此您还必须使用相同的索引对 V 的列进行重新排序。
energy(:,i)=d;
end
 
%%%%%%%plot 画图
for m=1:Nw                                 %%%  原胞中有几个轨道就有几条能带
    plot(k_new,energy(m,:)-efermi,'color','r','LineWidth',1)  %%
    hold on
end
 
 
%%%%设置xticks布里渊区高对称点
xticks([k_new(1) k_new(kn+1)  k_new(2*kn+1)  k_new(3*kn)]);
xticklabels({'\Gamma','M','K','\Gamma'})
 
%%%辅助线
xline(k_new(kn+1),'--');
xline(k_new(2*kn+1),'--');
yline(0,'--');
%grid on
 
%%%
xlim([k_new(1) k_new(3*kn)])
ylim([-10,15])
ylabel("Energy")
box on

补充：DFT 计算结果和wannier90的对比    方法：一文搞定VASP+wannier90构造紧束缚模型一文搞定VASP+wannier90构造紧束缚模型 - 知乎

补充   

#检查wannier90.wout的坐标时可能和POSCAR中的不一样，wannier90自动发现对称性了，选了另一种原子。按照右边那个原胞看wannier90hr.dat。（笔记）我用的晶体规范，没考虑原子坐标，上面结果没影响。

Wang Niewei
关注
  • 4
    点赞
  • 踩
  • 11
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
MATLAB计算石墨烯能带，判断最近邻，对角化哈密顿量matlab程序,matlab程序实例源码.zip
10-15
MATLAB计算石墨烯能带，判断最近邻，对角化哈密顿量matlab程序,matlab程序实例源码
Vasp第一性原理计算从零单排——石墨烯能带的计算
wwxy1995的博客
10-19 5529
本次计算借助二维材料2dmatpedia数据库给出输入文件，并结合自己的工具进行能带计算。 第一步，下载所需要计算体系的输入文件 理论上需要首先进行结构优化，比如这里的relax1，relax2。 其中结构优化的作用是重新得到POSCAR，结构优化形成的文件写在CONTCAR中 执行命令 cp CONTCAR POSCAR 然后开始静态计算，这里要特别注意，这里二维材料的静态计算...
MATLAB-wannier90_hr.dat后处理以及紧束缚能带拟合（晶体规范）
w126126的博客
04-19 1341
初次接触vasp,在学习的过程中发现没有完整的关于wannier画能带的Matlab程序，经过三周的学习编写了后处理wannier90_hr.dat的Matlab代码。参考了微信公众号“拾梦的星星”《Mathematica处理wannier90_hr.dat文件并画图》一文，增加了对代码框架理解，同时改进了倒空间高对称K点路径的设置，在提取wannier90_hr.dat数据过程中参考了...。Wannier90是一款强大的开源代码，可以生成局域wannier函数，并使用它们来高效精确地计算材料的电子性质。
锯齿石墨烯能带matlab计算程序
wwxy1995的博客
08-22 4690
这个程序是初步优化后的matlab版本，主要思路是先生成体系的格点坐标，再运用坐标生成体系的哈密顿量，然后进行对角化计算能带，能带的计算使用一维体系超原胞的处理方法。可以进一步优化 主程序 nx=3; % ny=100; % 体系宽度（y方向的长度） [x,y]=zigzag_graphene(nx,ny); %plot(x,y,'.','MarkerSize',20) t1=-2.7;...
【ASE学习】-测量石墨烯结构的碳碳键平均键长
袁家瑜的博客
10-30 806
本文使用ASE测量下面团簇结构的碳碳键键长，用于验证结构优化方法是否合理。即该团簇结构中的碳碳键平均键长为：1.427652873317778埃。实验上测得的碳碳键键长为1.42埃。因此计算误差为：0.5%左右。
Vasp 石墨烯能带计算
wwxy1995的博客
10-16 6797
这里暂时不做结构优化，使用数据库结构 首先做自洽计算产生正确的基态电子密度 INCAR SYSTEM = C PREC = accutat ISTART = 0 ICHARG = 2 EDIFF = 0.1E-03 NSW = 0 IBRION = 2 ISIF = 2 POSCAR C 2.45944 1.000 0.000 0.000 0.500 0.866 0.000 0.00...
文件的读写基本操作
热门推荐
范高伦的博客
09-07 1万+
一、文件是计算机中数据持久化存储的表现形式 读写文件标准操作格式1： 1、打开文件：file1 = open('文件名','读写模式') 2、操作文件 3、关闭文件：file1.close() 文件操作完毕后必须关闭，否则长期保持对文件的连接状态，造成内存溢出的现象发生 读写文件操作格式2： # 1、打开文件 file1 = open('demo.txt','w') # 2、操作文件 file1.write('hello world') # 3、关闭文件 file1.close() 1、打开文件：wi
石墨烯能带matlab代码-ttlg:用动量空间连续模型计算具有两个独立扭曲角的扭曲三层石墨烯的电子带结构和密度。
05-22
石墨烯能带matlab代码ttlg-多元三层石墨烯的电子结构计算 使用动量空间连续谱模型计算具有两个独立扭曲角的扭曲三层石墨烯的电子能带结构和密度。 所有脚本都在MATLAB中。 引文 有关ttlg模型的参考，请参见并引用...
MATLAB1.rar_matlab次短路_哈密顿回路_回路程序_最小顶点覆盖_覆盖路径
09-19
MATLAB经典算法程序 经典程序。顶点覆盖近似算法，哈密尔顿回路，画等温线，模拟退火应用，生成全排列矩阵，随机数的产生，最大流和最小截，最短路和次短路，最短路径，最小生成树Prim算法
基于MATLAB的哈密顿回路算法-TSP模拟退火 程序源代码.rar
04-23
基于MATLAB的哈密顿回路算法-TSP模拟退火 程序源代码.rar
基于模糊神经网络的时间序列预测（以hopkinsirandeath数据集为例，MATLAB）
weixin_39402231的博客
06-29 556
采用Takagi Sugeno Kang模糊神经网络对时间序列进行预测，以hopkinsirandeat数据为例进行说明，程序运行环境为MATLAB R2021B
利用MATLAB批量读取图像时出现名称排序错乱问题解决方法sort-nat函数
Zephyr Hou 的博客
06-26 423
利用MATLAB批量读取图像时出现名称排序错乱问题解决方法sort-nat函数
如何使用MATLAB写测试（4）重复劳动？TestMethodSetup！
爱编程的鱼的博客
06-24 284
类似的还有在test method运行后做收尾工作的TestMethodTeardown和TestClassTeardown。哟，不错，熟练掌握了TestParameter以及function handle，确实是名合格的实习生。我们来自俄罗斯的实习生是学Computer Vision的，最近他为我们展示了一个全新的程序。我看了后说，你老板要发飙的，上次学的TestParameter呢?，原来可以用TestMethodSetup来做一些重复的初始化工作。实习生翻阅了文档后，发现了？
基于MATLAB仿真设计无线充电系统
zzhao2580的博客
06-27 1156
通过学习无线充电相关课程知识，通过课程设计无线充电系统，将所学习的WPT，DC-DC，APFC进行整合得到整个无线充电系统，通过进行仿真研究其系统特性，完成我们预期系统功能和指标
matlab编辑稀疏单位方阵
qq_65800683的博客
06-24 190
创建 1000×1000 稀疏单位方阵，并查看稀疏模式。
数学建模系列（4/4）：Matlab建模实战
weidl001的博客
06-21 1364
在本篇文章中，我们详细介绍了如何在Matlab中进行数学模型的实现和优化。通过具体的案例和代码示例，读者可以更好地理解如何在Matlab中应用这些建模方法。高级主题部分的内容展示了Matlab在建模中的更多可能性，包括与Simulink的结合、并行计算和GUI设计。
数字图像处理实验报告小论文（Matlab语言）
m0_65963495的博客
06-24 1255
本文复现与测试图像处理技术的灰度图像处理、小波变换和图像分割处理的算法。首先，通过直方图均衡化技术增强了图像的对比度。接着，利用小波变换对图像进行了多尺度分析，并通过重构展示了不同频率成分对图像的影响。在图像分割处理部分，本文复现了固定阈值分割和大津算法全局阈值分割两种方法，并测试了它们在不同阈值下的分割效果。
4、matlab双目相机标定实验
最新发布
XU157303764的博客
06-30 800
matlab双目相机标定实验
matlab 石墨烯 armchair
02-06
MATLAB是一个非常强大的科学计算软件，可以用来进行各种数学模拟和数据分析。石墨烯是一种二维的碳材料，具有许多独特的性质，比如超高的导热性和电导率。armchair则是石墨烯晶格结构的一种特定排列方式。 在MATLAB中，可以使用各种数学模型和算法来模拟石墨烯的性质和行为。比如，可以使用MATLAB来进行石墨烯的电子结构计算，研究其电导率和能带结构。也可以使用MATLAB来模拟石墨烯的热传导行为，探究其在高温环境下的稳定性。 此外，MATLAB还可以用来进行石墨烯的可视化分析，将研究结果以图表或动画的形式展示出来，更直观地理解石墨烯的特性和行为。通过MATLAB，可以更深入地研究石墨烯在纳米电子器件、材料科学和纳米技术领域的应用潜力。 总之，MATLAB在石墨烯研究中具有重要的作用，可以帮助科学家们更深入地理解石墨烯的特性和行为，推动石墨烯材料科学的发展。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值