Se for utilizar o material, favor citar como:

1)J.P. Braga e F.S. Carvalho - Métodos Numéricos em Química Quântica, 2021.
2)J.P. Braga, Fundamentos de Química Quântica, Editora UFV, 2007.

PROGRAMA PRINCIPAL SEM EFEITO DE CORRELAÇÃO
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% HeKS-troca
clc,clear all,close all

Z=2;c=[10 7];r=0:0.01:20;n=0;
[p1,dp1,ddp1]=phi1(r);
[p2,dp2,ddp2]=phi2(r);

for w=1:10;    
n=n+1;

S11=matrizS(r,p1,p1);
S12=matrizS(r,p1,p2);
S21=S12;
S22=matrizS(r,p2,p2);

T11=matrizT(r,1,1);
T12=matrizT(r,1,2);
T21=T12;
T22=matrizT(r,2,2);

V11=matrizV(r,p1,p1,Z);
V12=matrizV(r,p1,p2,Z);
V21=V12;
V22=matrizV(r,p2,p2,Z);

r1=r;
r2=r;

J11=matrizJ(r,r1,r2,c,1,1);
J12=matrizJ(r,r1,r2,c,1,2);
J21=J12;
J22=matrizJ(r,r1,r2,c,2,2);

dens=rho(p1,p2,c);

E11=matrizE(r,c,p1,p1,dens);
E12=matrizE(r,c,p1,p2,dens);
E21=E12;
E22=matrizE(r,c,p2,p2,dens);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Matriz H
H(1,1)=T11+V11+E11+J11;
H(1,2)=T12+V12+E12+J12;
H(2,1)=T21+V21+E21+J21;
H(2,2)=T22+V22+E22+J22;

%Matriz S
S(1,1)=S11;
S(1,2)=S12;
S(2,1)=S12;
S(2,2)=S22;

[C,E]=eig(H,S);

c1=c(1);c2=c(2);
[n c1 c2]

c(1)=C(1,1);
c(2)=C(2,1);

Coef1(n)=C(1,1);
Coef2(n)=C(2,1);

end

%% calculando a densidade
dens=rho(p1,p2,c);

figure(1)
plot(r,dens,'-')
set(gca,'FontSize',18)
xlabel('r / u.a.','FontSize',18)
ylabel('\rho(r)','FontSize',18)
axis([0 5 0 2])


%% calculando a energia
eps_troca=2*E(1,1)-Jint(r,r1,r2,c)+kint(r,p1,p2,c)

%% calculando os potenciais

Vxc= -(3*dens/pi).^(1/3);
Ven= - Z./r;

r2=r';

[x,y]=meshgrid(r,r2);

[phi1r1,phi1r2]=phi1m(r,x,y);
[phi2r1,phi2r2]=phi2m(r,x,y);

rho2=2*(c(1)*phi1r2+c(2)*phi2r2).*conj(c(1)*phi1r2+c(2)*phi2r2);
kernel=(2*pi)*(y./x).*(x+y-abs(x-y)).*rho2;

Vee=trapz(r2,kernel,1);

figure(2)
plot(r,Vee,'-',r,Vxc,'--',r,Ven,'-.')
set(gca,'FontSize',18)
xlabel('r / u.a.','FontSize',18)
ylabel('V(r) / u.a.','FontSize',18)
axis([0 5 -4 3])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PROGRAMA PRINCIPAL COM EFEITO DE CORRELAÇÃO
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% HeKS-troca-corr
clc,clear all,close all

Z=2;c=[10 7];r=0:0.01:20;n=0;
[p1,dp1,ddp1]=phi1(r);
[p2,dp2,ddp2]=phi2(r);

for w=1:10;
    
n=n+1;

S11=matrizS(r,p1,p1);
S12=matrizS(r,p1,p2);
S21=S12;
S22=matrizS(r,p2,p2);

T11=matrizT(r,1,1);
T12=matrizT(r,1,2);
T21=T12;
T22=matrizT(r,2,2);

V11=matrizV(r,p1,p1,Z);
V12=matrizV(r,p1,p2,Z);
V21=V12;
V22=matrizV(r,p2,p2,Z);

r1=r;
r2=r;

J11=matrizJ(r,r1,r2,c,1,1);
J12=matrizJ(r,r1,r2,c,1,2);
J21=J12;
J22=matrizJ(r,r1,r2,c,2,2);

dens=rho(p1,p2,c);

E11=matrizE(r,c,p1,p1,dens);
E12=matrizE(r,c,p1,p2,dens);
E21=E12;
E22=matrizE(r,c,p2,p2,dens);

%% efeito de correlacao
C11=corr(r,p1,p1,dens);
C12=corr(r,p1,p2,dens);
C21=C12;
C22=corr(r,p2,p2,dens);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Matriz H
H(1,1)=T11+V11+E11+J11+C11;
H(1,2)=T12+V12+E12+J12+C12;
H(2,1)=T21+V21+E21+J21+C21;
H(2,2)=T22+V22+E22+J22+C22;

%Matriz S
S(1,1)=S11;
S(1,2)=S12;
S(2,1)=S12;
S(2,2)=S22;

[C,E]=eig(H,S);
    
[n c(1) c(2)]

c(1)=C(1,1);
c(2)=C(2,1);

Coef1(n)=C(1,1);
Coef2(n)=C(2,1);

end

%% calculando a densidade
dens=rho(p1,p2,c);

figure(1)
plot(r,dens,'-')
set(gca,'FontSize',18)
xlabel('r / u.a.','FontSize',18)
ylabel('\rho(r)','FontSize',18)
axis([0 5 0 2])
%print -deps -r600 rhoHeKScor

%% calculando a energia
c0=0.0252;r0=30;Y=(3*pi/4)^(1/3);

Bcorr=-c0*((1+(Y^3)./((r0^3)*dens)).*((1+(r0*dens.^(1/3))/Y).^(-1)).*(r0./(3*Y*dens.^(2/3)))-...
    ((Y^3)./((r0^3)*dens.^2)).*log(1+((r0*dens.^(1/3))/Y))-Y./(6*r0*dens.^(4/3))+(2*Y^2)./(3*(r0^2)*dens.^(5/3)));

fat=trapz(r,4*pi*r.^2.*Bcorr.*dens.^2);

eps_trocaEcorr=2*E(1,1)-Jint(r,r1,r2,c)+kint(r,p1,p2,c)+fat

%% calculando os potenciais

Vxc= -(3*dens/pi).^(1/3);
Ven= - Z./r;

r2=r';

[x,y]=meshgrid(r,r2);

[phi1r1,phi1r2]=phi1m(r,x,y);
[phi2r1,phi2r2]=phi2m(r,x,y);

rho2=2*(c(1)*phi1r2+c(2)*phi2r2).*conj(c(1)*phi1r2+c(2)*phi2r2);
kernel=(2*pi)*(y./x).*(x+y-abs(x-y)).*rho2;

Vee=trapz(r2,kernel,1);

epscorr=-c0*((1+(Y^3)./(dens*r0^3)).*log(1+(r0*dens.^(1/3))/Y)+...
    (Y./(2*r0*dens.^(1/3)))-((Y^2)./((r0^2)*dens.^(2/3)))-1/3);

Bcorr=-c0*((1+(Y^3)./((r0^3)*dens)).*((1+(r0*dens.^(1/3))/Y).^(-1)).*(r0./(3*Y*dens.^(2/3)))-...
    ((Y^3)./((r0^3)*dens.^2)).*log(1+((r0*dens.^(1/3))/Y))-Y./(6*r0*dens.^(4/3))+(2*Y^2)./(3*(r0^2)*dens.^(5/3)));

Vcor=epscorr+Bcorr.*dens;

figure(2)
plot(r,Vee,'-',r,Vxc,'--',r,Ven,'-.'),hold on,plot(r,Vcor,':','LineWidth',1.5)
set(gca,'FontSize',18)
xlabel('r / u.a.','FontSize',18)
ylabel('V(r) / u.a.','FontSize',18)
axis([0 5 -4 3])
%print -deps -r600 VsHeKScor

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

FUNÇÕES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [p,dp,ddp]=phi1(r)
psi=exp(-1.4536*r);
N1=1./(4*pi*trapz(r,r.^2.*psi.*psi))^0.5;
p=N1*exp(-1.4536*r);
dp=-1.4536*N1*exp(-1.4536*r);
ddp=N1*(1.4536^2)*exp(-1.4536*r);


function [p,dp,ddp]=phi2(r)
psi=exp(-2.9109*r);
N2=1./(4*pi*trapz(r,r.^2.*psi.*psi))^0.5;
p=N2*exp(-2.9109*r);
dp=-2.9109*N2*exp(-2.9109*r);
ddp=N2*(2.9109^2)*exp(-2.9109*r);


function f=rho(p,q,c)
f=2*(c(1)*p+c(2)*q).*conj((c(1)*p+c(2)*q));


function f=matrizT(r,p,q)
if p == 1
    [p,dp,ddp]=phi1(r);
else
    [p,dp,ddp]=phi2(r);
end
if q == 1
    [q,dq,ddq]=phi1(r);
else
    [q,dq,ddq]=phi2(r);
end    
kernel=-0.5*4*pi*(p.*(2*r.*dq + r.*r.*ddq));
f=trapz(r,kernel);


function f=matrizV(r,p,q,Z)
kernel=-Z*(4*pi*r).*p.*q;
f=trapz(r,kernel);


function f=matrizS(r,p,q)
kernel=4*pi*r.*r.*p.*q;
f=trapz(r,kernel);


function [f,g]=phi1m(r,x,y)
psi=exp(-1.4536*r);
N1=1./(4*pi*trapz(r,r.^2.*psi.*psi))^0.5;
f=N1*exp(-1.4536*x);
g=N1*exp(-1.4536*y);


function [f,g]=phi2m(r,x,y)
psi=exp(-2.9109*r);
N2=1./(4*pi*trapz(r,r.^2.*psi.*psi))^0.5;
f=N2*exp(-2.9109*x);
g=N2*exp(-2.9109*y);


function f=matrizJ(r,r1,r2,c,i,j)
[x,y]=meshgrid(r1,r2);
[phi1r1,phi1r2]=phi1m(r,x,y);
[phi2r1,phi2r2]=phi2m(r,x,y);
rho=2*(c(1)*phi1r2+c(2)*phi2r2).*conj(c(1)*phi1r2+c(2)*phi2r2);
if i==1
    phii=phi1r1;
else
    phii=phi1r2;
end
if j==1
    phij=phi1r1;
else
    phij=phi1r2;
end
kernel=(8*pi^2)*x.*y.*(x+y-abs(x-y)).*rho.*phii.*phij;
f=trapz(r2,trapz(r1,kernel,2));


function f=matrizE(r,c,p,q,dens)
kernel=-((3/pi)^(1/3))*4*pi*r.*r.*(dens.^(1/3)).*p.*q;
f=trapz(r,kernel);


function f=corr(r,p,q,dens)
c0=0.0252;r0=30;Y=(3*pi/4)^(1/3);
epscorr=-c0*((1+(Y^3)./(dens*r0^3)).*log(1+(r0*dens.^(1/3))/Y)+...
    (Y./(2*r0*dens.^(1/3)))-((Y^2)./((r0^2)*dens.^(2/3)))-1/3);
Bcorr=-c0*((1+(Y^3)./((r0^3)*dens)).*((1+(r0*dens.^(1/3))/Y).^(-1)).*(r0./(3*Y*dens.^(2/3)))-...
    ((Y^3)./((r0^3)*dens.^2)).*log(1+((r0*dens.^(1/3))/Y))-Y./(6*r0*dens.^(4/3))+(2*Y^2)./(3*(r0^2)*dens.^(5/3)));
kernel=4*pi*r.*r.*(epscorr+dens.*Bcorr).*p.*q;
f=trapz(r,kernel);


function f=Jint(r,r1,r2,c)
[x,y]=meshgrid(r1,r2);
[phi1r1,phi1r2]=phi1m(r,x,y);
[phi2r1,phi2r2]=phi2m(r,x,y);
rho1=2*(c(1)*phi1r1+c(2)*phi2r1).*conj(c(1)*phi1r1+c(2)*phi2r1);
rho2=2*(c(1)*phi1r2+c(2)*phi2r2).*conj(c(1)*phi1r2+c(2)*phi2r2);
kernel=0.5*(8*pi^2)*x.*y.*(x+y-abs(x-y)).*rho1.*rho2;
f=trapz(r2,trapz(r1,kernel,2));


function f=kint(r,p,q,c)
rho=2*(c(1)*p+c(2)*q).*conj((c(1)*p+c(2)*q));
kernel=4*pi*(1/4)*(3/pi)^(1/3)*r.^2.*rho.^(4/3);
f=trapz(r,kernel);