% Numerische und grafische Auswertung der Energie-Wellenzahl-Dispersion
% eines homogenen Elektronengases in der Hartree-Fock-Naeherung und in der
% RPA-Korrelations-Naeherung.
%
% Impulse in Einheiten 2*pi/a, a=Gitterkonstante
% Energien in Rydberg
%
% H. Sormann   17-10-2001
% Last update  20-12-2004
%

  clc

  clear all

  metal=input('sodium: 1  aluminum: 2  ');
  
  if(metal==1)
    fprintf(' SODIUM  bcc  Gitterkonstante (Bohr) = 8.003\n');
    a=8.003;
    fac=2.0*pi/a;
    gamma=1.0;
    a0kf=fac*(3*gamma/4/pi)^(1.0/3.0);
    rs=(9*pi/4)^(1/3) /a0kf;

    lapw=load('Na_energy.sor');
    ef_lapw = lapw(1,1);
    lapw(1,:)=[];

% Daten aus Lyo and Plummer, 1988:
    k_exp=[0 0.5 0.95 1.43 1.9 2.2];
    e_exp=[-4.55 -4.20 -3.70 -2.60 -1.10 0];
    k_exp=0.2818*k_exp;                        % in 2pi/a
    e_exp=0.5882*e_exp/13.6;                   % in Ry
    
% Rechnung Sormann: Energie-Dispersion inkl. Exchange und 
%                   RPA-Korrelation (12-11-03):
    rpa=load('corr_rpa.dat');
    ef_rpa=rpa(7,2);
    
  end

  if(metal==2)
    fprintf(' ALUMINUM  fcc  Gitterkonstante (Bohr) = 7.6395\n');
    a=7.6395;
    fac=2.0*pi/a;
    gamma=3.0;
    a0kf=fac*(3*gamma/2/pi)^(1.0/3.0);
    rs=(9*pi/4)^(1/3) /a0kf
    fprintf('ENTER\n');
    pause

    lapw=load('Al_energy.sor');
    ef_lapw = lapw(1,1);
    lapw(1,:)=[];
  end


  kfeld=0.001:0.01:1.00;   % k in units  2pi/a
  kfeld1=fac/a0kf * kfeld;   % k in units kf

  bandwidth_free=a0kf^2;
  e_free=a0kf^2 * kfeld1.^2 - bandwidth_free;

  kk=0.00001;
  e_bottom=a0kf^2*kk^2 - a0kf/pi * (2.0 + (1.0-kk^2)/kk * log(abs((1.0+kk)/(1.0-kk))));
  kk=0.99999;
  e_fermi=a0kf^2*kk^2 - a0kf/pi * (2.0 + (1.0-kk^2)/kk * log(abs((1.0+kk)/(1.0-kk))));

  bandwidth_x = e_fermi-e_bottom;

  e_x=a0kf^2*kfeld1.^2 - a0kf/pi * (2.0 + (1.0-kfeld1.^2)./kfeld1 .* log(abs((1.0+kfeld1)./(1.0-kfeld1)))) - e_fermi;

 % fprintf('ENTER\n');
 % pause

  if(metal==1)
    plot(kfeld',e_free','--')
    axis([0 0.7 -0.6 0.2])
    title('ENERGIE-K-DISPERSION  SODIUM   RS=3.9405 BOHR');
    xlabel('k (2*pi/a)');
    ylabel('ENERGIE [RYD]');
    legend('freies Elektron: 0.2372 Ry (0.0 Proz.)');
    grid on
    hold on
    pause
    plot(lapw(:,2),lapw(:,5)-ef_lapw,'ro')
    legend('freies Elektron: 0.2372 Ry (0.0 Proz.)','LAPW-Bandstruktur: 0.2401 Ry (+1.2 Proz.)');
    fprintf('ENTER\n');
    pause
    plot(kfeld',e_x')
    legend('freies Elektron: 0.2372 Ry (0.0 Proz.)','LAPW-Bandstruktur: 0.2401 Ry (+1.2 Proz.)','HF-Elektron: 0.5472 Ry (+130.7 Proz.)');
    fprintf('ENTER\n');
    pause
    plot(k_exp,e_exp,'k*')
    legend('freies Elektron: 0.2372 Ry (0.0 Proz.)','LAPW-Bandstruktur: 0.2401 Ry (+1.2 Proz.)','HF-Elektron: 0.5472 Ry (+130.7 Proz.)','exp: Lyo and Plummer 1988: 0.1968 Ry (-17.0 Proz.)');
    fprintf('ENTER\n');
    pause
    plot(0.62*rpa(:,1),rpa(:,2)-ef_rpa,'diamond-')
    legend('freies Elektron: 0.2372 Ry (0.0 Proz.)','LAPW-Bandstruktur: 0.2401 Ry (+1.2 Proz.)','HF-Elektron: 0.5472 Ry (+130.7 Proz.)','exp: Lyo and Plummer 1988: 0.1968 Ry (-17.0 Proz.)','RPA-Korrelation (Sormann, 2004): 0.1901 Ry (-19.9 Proz.)');
  end

  if(metal==2)
    plot(kfeld',e_free','--',kfeld',e_x')
%    axis([0 0.7 -0.6 0.2])
    title('ENERGIE-K-DISPERSION  ALUMINUM');
    xlabel('k (2*pi/a)');
    ylabel('ENERGIE [RYD]');
    legend('freies Elektron','HF-Elektron');
    grid on
    hold on
    fprintf('ENTER\n');
    pause
    plot(lapw(:,2),lapw(:,5)-ef_lapw,'ro')
    legend('freies Elektron','HF-Elektron','LAPW-Bandstruktur');
  end
  hold off
