1% test5.m
2% reads and plots the output of Siesta/Src/SiestaXC/Testers/test5.f90,
3% of the fit of the Vydrov-vanVoorhis vdW functional kernel.
4% J.M.Soler, Jul.2012
5
6% Read data
7data = load('n1n2phi.table');
8r  = data(:,1);
9kf = data(:,2);
10kg = data(:,3);
11n1n2phi_exact = data(:,4);
12n1n2phi_interp = data(:,5);
13r = unique(r);
14kf = unique(kf);
15kg = unique(kg);
16nr = numel(r);
17nkf = numel(kf);
18nkg = numel(kg);
19n1n2phi_exact = reshape(n1n2phi_exact,nkg,nkf,nr);
20n1n2phi_interp = reshape(n1n2phi_interp,nkg,nkf,nr);
21
22% Plot data
23[kf,kg] = meshgrid(kf,kg);
24for ir = 1:nr
25    figure(ir)
26    surf(kf,kg,n1n2phi_exact(:,:,ir)*27.2e6), hold on
27    surf(kf,kg,n1n2phi_interp(:,:,ir)*27.2e6)
28%    surf( kf, kg, (n1n2phi_interp(:,:,ir)-n1n2phi_exact(:,:,ir))*27.2e6 )
29%    surf( kf, kg, n1n2phi_interp(:,:,ir)./n1n2phi_exact(:,:,ir) )
30    xlabel('k_F'), ylabel('\nabla(n) / n')
31    zlabel('n_1 n_2 \Phi (k_1,k_2,r_{12} (\mueV / bohr^6)')
32    title(['r_{12} = ',num2str(r(ir))])
33    hold off
34end
35