matlib科学计算平均场能带技巧

最新推荐文章于 2021-11-22 00:18:07 发布

wwxy261

最新推荐文章于 2021-11-22 00:18:07 发布

阅读量274

点赞数

分类专栏：算法

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

版权

算法专栏收录该内容

3633 篇文章 117 订阅

订阅专栏

1. 注意使用向量化编程技巧，使得代码中尽可能少用循环，而用向量化计算实现。

2. 注意同时在主计算程序中，少用逻辑控制，这样能提高计算效率，后期也能改成GPU计算

3. 还可以优化的点是把spin也做成多维数组，这样代码写起来简单

nk = 200;                        
k = linspace(0,2*pi,nk);
t = 1;          % 最近邻hooping
U = 0.5;        % Hubbard系数
mu = -sqrt(3);  % 化学势1/4填充取值
N = 40;         % N/2是对应的y方向原子个数
m = 4*N;        % 总的原子的个数
T = 0.001;      % 温度
N1avg = 0.25*eye(m);
N2avg = 0.25*eye(m);
N3avg = 0.25*eye(m);
N4avg = 0.25*eye(m);

N1avg(1,1) = 0;
N2avg(1,1) = 0;
N3avg(1,1) = 0;
N4avg(1,1) = 0;
N1avg(2,2) = 0.7;
N2avg(2,2) = 0.1;
N3avg(2,2) = 0.1;
N4avg(2,2) = 0.1;

Ncc = 6;        % 自洽计算迭代次数
band1=zeros(m,nk);
band2=zeros(m,nk);
band3=zeros(m,nk);
band4=zeros(m,nk);

for ncc = 1:Ncc
    N1avg0=zeros(m,m);                   %初始化N1avg的临时变量用于k空间的积分
    N2avg0=zeros(m,m);                   %初始化N2avg的临时变量用于k空间的积分
    N3avg0=zeros(m,m); 
    N4avg0=zeros(m,m);
    m1avg = (N1avg-N2avg)/2;
    m2avg = (N1avg+N2avg-2*N3avg)/(2*sqrt(3));
    m3avg = (N1avg+N2avg+N3avg-3*N4avg)/(2*sqrt(6));
    C = 8/5*U*(m1avg.^2+m2avg.^2+m3avg.^2);
    U1 = m1avg/2+m2avg/(2*sqrt(3)) + m3avg/(2*sqrt(6));
    U2 = -m1avg/2+m2avg/(2*sqrt(3))+ m3avg/(2*sqrt(6));
    U3 =  -2*m2avg/(2*sqrt(3)) + m3avg/(2*sqrt(6));
    U4 =  -3*m3avg/(2*sqrt(6));
    for i = 1:nk
        Hk =  Hamiltonian_H0(k(i),t,N,mu);
        [A1,E1] = eig(Hk-16/5*U*U1+C);
        [A2,E2] = eig(Hk-16/5*U*U2+C);
        [A3,E3] = eig(Hk-16/5*U*U3+C);
        [A4,E4] = eig(Hk-16/5*U*U4+C);
        f1 = Fermi_funtion(diag(E1),0,T)';
        f2 = Fermi_funtion(diag(E2),0,T)';
        f3 = Fermi_funtion(diag(E3),0,T)';
        f4 = Fermi_funtion(diag(E4),0,T)';
        fermi_1 = f1(ones(m,1),:);
        fermi_2 = f2(ones(m,1),:);
        fermi_3 = f3(ones(m,1),:);
        fermi_4 = f4(ones(m,1),:);
        N1 = sum(A1.*conj(A1).*fermi_1,2);
        N2 = sum(A2.*conj(A2).*fermi_2,2);
        N3 = sum(A3.*conj(A3).*fermi_3,2);
        N4 = sum(A4.*conj(A4).*fermi_4,2);
        N1avg0 = N1avg0+diag(N1);
        N2avg0 = N2avg0+diag(N2);
        N3avg0 = N3avg0+diag(N3);
        N4avg0 = N4avg0+diag(N4);
    end
    N1avg = N1avg0/nk;
    N2avg = N2avg0/nk;
    N3avg = N3avg0/nk;
    N4avg = N4avg0/nk;
    fprintf('已完成第%d次迭代\n',ncc)
end
fprintf("自洽计算完成\n")
fprintf("能带计算开始\n")

%%{
for i = 1:nk
    Hk =  Hamiltonian_H0(k(i),t,N,mu);
    [~,E1] = eig(Hk-16/5*U*U1+C);
    [~,E2] = eig(Hk-16/5*U*U2+C);
    [~,E3] = eig(Hk-16/5*U*U3+C);
    [~,E4] = eig(Hk-16/5*U*U4+C);
    band1(:,i) = diag(E1);
    band2(:,i) = diag(E2);
    band3(:,i) = diag(E3);
    band4(:,i) = diag(E4);
end
plot(k,band1);
hold on
plot(k,band2)
set(gca,'YLim',[-0.5 0.5]);
set(gca,'XLim',[0 2*pi]);
%%}