mirror of https://github.com/pfloos/quack
rename UGW routines
This commit is contained in:
parent
6b56acf66b
commit
eee5f6bac8
|
|
@ -118,7 +118,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons
|
|||
! Excitation densities !
|
||||
!----------------------!
|
||||
|
||||
call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
|
||||
call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
|
||||
|
||||
!------------------------------------------------!
|
||||
! Compute self-energy and renormalization factor !
|
||||
|
|
@ -131,8 +131,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons
|
|||
|
||||
else
|
||||
|
||||
call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC,EcGM)
|
||||
call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,Z)
|
||||
call UGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC,Z,EcGM)
|
||||
|
||||
end if
|
||||
|
||||
|
|
@ -154,8 +153,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons
|
|||
! Find graphical solution of the QP equation
|
||||
|
||||
do is=1,nspin
|
||||
call unrestricted_QP_graph(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eta,eHF(:,is), &
|
||||
OmRPA,rho_RPA(:,:,:,is),eGWlin(:,is),eGW(:,is))
|
||||
call UGW_QP_graph(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eta,eHF(:,is), &
|
||||
OmRPA,rho_RPA(:,:,:,is),eGWlin(:,is),eGW(:,is))
|
||||
end do
|
||||
|
||||
end if
|
||||
|
|
@ -177,8 +176,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons
|
|||
|
||||
if(BSE) then
|
||||
|
||||
call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,S, &
|
||||
ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF,eGW,EcBSE)
|
||||
call UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,S, &
|
||||
ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF,eGW,EcBSE)
|
||||
|
||||
if(exchange_kernel) then
|
||||
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
subroutine unrestricted_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW)
|
||||
subroutine UGW_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW)
|
||||
|
||||
! Compute the graphical solution of the QP equation
|
||||
|
||||
|
|
@ -80,4 +80,4 @@ subroutine unrestricted_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eG
|
|||
|
||||
end do
|
||||
|
||||
end subroutine unrestricted_QP_graph
|
||||
end subroutine
|
||||
|
|
@ -1,4 +1,4 @@
|
|||
subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,rho)
|
||||
subroutine UGW_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,rho)
|
||||
|
||||
! Compute excitation densities for unrestricted reference
|
||||
|
||||
|
|
@ -103,4 +103,4 @@ subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aaaa,ER
|
|||
enddo
|
||||
enddo
|
||||
|
||||
end subroutine unrestricted_excitation_density
|
||||
end subroutine
|
||||
|
|
@ -95,7 +95,7 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s
|
|||
|
||||
call phULR(isp_W,.true.,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,eW, &
|
||||
ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
|
||||
call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
|
||||
call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
|
||||
|
||||
! Spin-conserved manifold
|
||||
|
||||
|
|
@ -122,7 +122,7 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s
|
|||
|
||||
call phULR(isp_W,.true.,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, &
|
||||
ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
|
||||
call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
|
||||
call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
|
||||
|
||||
end if
|
||||
|
||||
|
|
@ -176,7 +176,7 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s
|
|||
|
||||
call phULR(isp_W,.true.,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, &
|
||||
ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
|
||||
call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
|
||||
call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
|
||||
|
||||
end if
|
||||
|
||||
|
|
|
|||
|
|
@ -1,6 +1,6 @@
|
|||
subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta, &
|
||||
nBas,nC,nO,nV,nR,nS,S,ERI_aaaa,ERI_aabb,ERI_bbbb, &
|
||||
dipole_int_aa,dipole_int_bb,cW,eW,eGW,EcBSE)
|
||||
subroutine UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta, &
|
||||
nBas,nC,nO,nV,nR,nS,S,ERI_aaaa,ERI_aabb,ERI_bbbb, &
|
||||
dipole_int_aa,dipole_int_bb,cW,eW,eGW,EcBSE)
|
||||
|
||||
! Compute the Bethe-Salpeter excitation energies
|
||||
|
||||
|
|
@ -82,7 +82,7 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_f
|
|||
call phULR(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, &
|
||||
eW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
|
||||
|
||||
call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
|
||||
call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
|
||||
|
||||
|
||||
!----------------------------!
|
||||
|
|
@ -1,5 +1,5 @@
|
|||
subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,eGW, &
|
||||
ERI_aaaa,ERI_aabb,ERI_bbbb,Omega,rho,A_lr)
|
||||
subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,eGW, &
|
||||
ERI_aaaa,ERI_aabb,ERI_bbbb,Omega,rho,A_lr)
|
||||
|
||||
! Compute the extra term for Bethe-Salpeter equation for linear response in the unrestricted formalism
|
||||
|
||||
|
|
@ -149,4 +149,4 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n
|
|||
|
||||
end if
|
||||
|
||||
end subroutine unrestricted_Bethe_Salpeter_A_matrix
|
||||
end subroutine
|
||||
|
|
@ -1,5 +1,5 @@
|
|||
subroutine unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, &
|
||||
ERI_aaaa,ERI_aabb,ERI_bbbb,Omega,rho,B_lr)
|
||||
subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, &
|
||||
ERI_aaaa,ERI_aabb,ERI_bbbb,Omega,rho,B_lr)
|
||||
|
||||
! Compute the extra term for Bethe-Salpeter equation for linear response
|
||||
|
||||
|
|
@ -150,4 +150,4 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n
|
|||
|
||||
end if
|
||||
|
||||
end subroutine unrestricted_Bethe_Salpeter_B_matrix
|
||||
end subroutine
|
||||
|
|
@ -0,0 +1,134 @@
|
|||
subroutine UGW_self_energy(eta,nBas,nC,nO,nV,nR,nSt,e,Om,rho,Sig,Z,EcGM)
|
||||
|
||||
! Compute diagonal of the correlation part of the self-energy
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
double precision,intent(in) :: eta
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC(nspin)
|
||||
integer,intent(in) :: nO(nspin)
|
||||
integer,intent(in) :: nV(nspin)
|
||||
integer,intent(in) :: nR(nspin)
|
||||
integer,intent(in) :: nSt
|
||||
double precision,intent(in) :: e(nBas,nspin)
|
||||
double precision,intent(in) :: Om(nSt)
|
||||
double precision,intent(in) :: rho(nBas,nBas,nSt,nspin)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: i,a,p,q,m
|
||||
double precision :: num,eps
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: Sig(nBas,nBas,nspin)
|
||||
double precision,intent(out) :: Z(nBas,nspin)
|
||||
double precision :: EcGM(nspin)
|
||||
|
||||
! Initialize
|
||||
|
||||
Sig(:,:,:) = 0d0
|
||||
Z(:,:) = 0d0
|
||||
EcGM(:) = 0d0
|
||||
|
||||
!--------------!
|
||||
! Spin-up part !
|
||||
!--------------!
|
||||
|
||||
! Occupied part of the correlation self-energy
|
||||
|
||||
do p=nC(1)+1,nBas-nR(1)
|
||||
do q=nC(1)+1,nBas-nR(1)
|
||||
do i=nC(1)+1,nO(1)
|
||||
do m=1,nSt
|
||||
eps = e(p,1) - e(i,1) + Om(m)
|
||||
num = rho(p,i,m,1)*rho(q,i,m,1)
|
||||
Sig(p,q,1) = Sig(p,q,1) + num*eps/(eps**2 + eta**2)
|
||||
if(p == q) Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! Virtual part of the correlation self-energy
|
||||
|
||||
do p=nC(1)+1,nBas-nR(1)
|
||||
do q=nC(1)+1,nBas-nR(1)
|
||||
do a=nO(1)+1,nBas-nR(1)
|
||||
do m=1,nSt
|
||||
eps = e(p,1) - e(a,1) - Om(m)
|
||||
num = rho(p,a,m,1)*rho(q,a,m,1)
|
||||
Sig(p,q,1) = Sig(p,q,1) + num*eps/(eps**2 + eta**2)
|
||||
if(p == q) Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! GM correlation energy
|
||||
|
||||
do i=nC(1)+1,nO(1)
|
||||
do a=nO(1)+1,nBas-nR(1)
|
||||
do m=1,nSt
|
||||
eps = e(a,1) - e(i,1) + Om(m)
|
||||
num = rho(a,i,m,1)**2
|
||||
EcGM(1) = EcGM(1) - num*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
!----------------!
|
||||
! Spin-down part !
|
||||
!----------------!
|
||||
|
||||
! Occupied part of the correlation self-energy
|
||||
|
||||
do p=nC(2)+1,nBas-nR(2)
|
||||
do q=nC(2)+1,nBas-nR(2)
|
||||
do i=nC(2)+1,nO(2)
|
||||
do m=1,nSt
|
||||
eps = e(p,2) - e(i,2) + Om(m)
|
||||
num = rho(p,i,m,2)*rho(q,i,m,2)
|
||||
Sig(p,q,2) = Sig(p,q,2) + num*eps/(eps**2 + eta**2)
|
||||
if(p == q) Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! Virtual part of the correlation self-energy
|
||||
|
||||
do p=nC(2)+1,nBas-nR(2)
|
||||
do q=nC(2)+1,nBas-nR(2)
|
||||
do a=nO(2)+1,nBas-nR(2)
|
||||
do m=1,nSt
|
||||
eps = e(p,2) - e(a,2) - Om(m)
|
||||
num = rho(p,a,m,2)*rho(q,a,m,2)
|
||||
Sig(p,q,2) = Sig(p,q,2) + num*eps/(eps**2 + eta**2)
|
||||
if(p == q) Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! GM correlation energy
|
||||
|
||||
do i=nC(2)+1,nO(2)
|
||||
do a=nO(2)+1,nBas-nR(2)
|
||||
do m=1,nSt
|
||||
eps = e(a,2) - e(i,2) + Om(m)
|
||||
num = rho(a,i,m,2)**2
|
||||
EcGM(2) = EcGM(2) - num*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! Compute renormalization factor from derivative
|
||||
|
||||
Z(:,:) = 1d0/(1d0 - Z(:,:))
|
||||
|
||||
end subroutine
|
||||
|
|
@ -0,0 +1,126 @@
|
|||
subroutine UGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nSt,e,Om,rho,Sig,Z,EcGM)
|
||||
|
||||
! Compute diagonal of the correlation part of the self-energy
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
double precision,intent(in) :: eta
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC(nspin)
|
||||
integer,intent(in) :: nO(nspin)
|
||||
integer,intent(in) :: nV(nspin)
|
||||
integer,intent(in) :: nR(nspin)
|
||||
integer,intent(in) :: nSt
|
||||
double precision,intent(in) :: e(nBas,nspin)
|
||||
double precision,intent(in) :: Om(nSt)
|
||||
double precision,intent(in) :: rho(nBas,nBas,nSt,nspin)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: i,a,p,m
|
||||
double precision :: num,eps
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: Sig(nBas,nspin)
|
||||
double precision,intent(out) :: Z(nBas,nspin)
|
||||
double precision :: EcGM(nspin)
|
||||
|
||||
! Initialize
|
||||
|
||||
Sig(:,:) = 0d0
|
||||
Z(:,:) = 0d0
|
||||
EcGM(:) = 0d0
|
||||
|
||||
!--------------!
|
||||
! Spin-up part !
|
||||
!--------------!
|
||||
|
||||
! Occupied part of the correlation self-energy
|
||||
|
||||
do p=nC(1)+1,nBas-nR(1)
|
||||
do i=nC(1)+1,nO(1)
|
||||
do m=1,nSt
|
||||
eps = e(p,1) - e(i,1) + Om(m)
|
||||
num = rho(p,i,m,1)**2
|
||||
Sig(p,1) = Sig(p,1) + num*eps/(eps**2 + eta**2)
|
||||
Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! Virtual part of the correlation self-energy
|
||||
|
||||
do p=nC(1)+1,nBas-nR(1)
|
||||
do a=nO(1)+1,nBas-nR(1)
|
||||
do m=1,nSt
|
||||
eps = e(p,1) - e(a,1) - Om(m)
|
||||
num = rho(p,a,m,1)**2
|
||||
Sig(p,1) = Sig(p,1) + num*eps/(eps**2 + eta**2)
|
||||
Z(p,1) = Z(p,1) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! GM correlation energy
|
||||
|
||||
do i=nC(1)+1,nO(1)
|
||||
do a=nO(1)+1,nBas-nR(1)
|
||||
do m=1,nSt
|
||||
eps = e(a,1) - e(i,1) + Om(m)
|
||||
num = rho(a,i,m,1)**2
|
||||
EcGM(1) = EcGM(1) - num*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
!----------------!
|
||||
! Spin-down part !
|
||||
!----------------!
|
||||
|
||||
! Occupied part of the correlation self-energy
|
||||
|
||||
do p=nC(2)+1,nBas-nR(2)
|
||||
do i=nC(2)+1,nO(2)
|
||||
do m=1,nSt
|
||||
eps = e(p,2) - e(i,2) + Om(m)
|
||||
num = rho(p,i,m,2)**2
|
||||
Sig(p,2) = Sig(p,2) + num*eps/(eps**2 + eta**2)
|
||||
Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! Virtual part of the correlation self-energy
|
||||
|
||||
do p=nC(2)+1,nBas-nR(2)
|
||||
do a=nO(2)+1,nBas-nR(2)
|
||||
do m=1,nSt
|
||||
eps = e(p,2) - e(a,2) - Om(m)
|
||||
num = rho(p,a,m,2)**2
|
||||
Sig(p,2) = Sig(p,2) + num*eps/(eps**2 + eta**2)
|
||||
Z(p,2) = Z(p,2) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! GM correlation energy
|
||||
|
||||
do i=nC(2)+1,nO(2)
|
||||
do a=nO(2)+1,nBas-nR(2)
|
||||
do m=1,nSt
|
||||
eps = e(a,2) - e(i,2) + Om(m)
|
||||
num = rho(a,i,m,2)**2
|
||||
EcGM(2) = EcGM(2) - num*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! Compute renormalization factor from derivative
|
||||
|
||||
Z(:,:) = 1d0/(1d0 - Z(:,:))
|
||||
|
||||
end subroutine
|
||||
|
|
@ -135,7 +135,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
|
|||
! Excitation densities !
|
||||
!----------------------!
|
||||
|
||||
call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
|
||||
call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
|
||||
|
||||
!------------------------------------------------!
|
||||
! Compute self-energy and renormalization factor !
|
||||
|
|
@ -148,8 +148,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
|
|||
|
||||
else
|
||||
|
||||
call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC,EcGM)
|
||||
call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,Z)
|
||||
call UGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC,Z,EcGM)
|
||||
|
||||
end if
|
||||
|
||||
|
|
@ -222,8 +221,8 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
|
|||
|
||||
if(BSE) then
|
||||
|
||||
call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, &
|
||||
S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eGW,eGW,EcBSE)
|
||||
call UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, &
|
||||
S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eGW,eGW,EcBSE)
|
||||
|
||||
if(exchange_kernel) then
|
||||
|
||||
|
|
|
|||
|
|
@ -204,7 +204,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
|
|||
! Excitation densities !
|
||||
!----------------------!
|
||||
|
||||
call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
|
||||
call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
|
||||
|
||||
!------------------------------------------------!
|
||||
! Compute self-energy and renormalization factor !
|
||||
|
|
@ -217,8 +217,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
|
|||
|
||||
else
|
||||
|
||||
call unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC,EcGM)
|
||||
call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,Z)
|
||||
call UGW_self_energy(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC,Z,EcGM)
|
||||
|
||||
end if
|
||||
|
||||
|
|
@ -362,8 +361,8 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
|
|||
|
||||
if(BSE) then
|
||||
|
||||
call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, &
|
||||
S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,c,eGW,eGW,EcBSE)
|
||||
call UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, &
|
||||
S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,c,eGW,eGW,EcBSE)
|
||||
|
||||
if(exchange_kernel) then
|
||||
|
||||
|
|
|
|||
|
|
@ -1,90 +0,0 @@
|
|||
subroutine unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,Z)
|
||||
|
||||
! Compute the renormalization factor in the unrestricted formalism
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
double precision,intent(in) :: eta
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC(nspin)
|
||||
integer,intent(in) :: nO(nspin)
|
||||
integer,intent(in) :: nV(nspin)
|
||||
integer,intent(in) :: nR(nspin)
|
||||
integer,intent(in) :: nSt
|
||||
double precision,intent(in) :: e(nBas,nspin)
|
||||
double precision,intent(in) :: Omega(nSt)
|
||||
double precision,intent(in) :: rho(nBas,nBas,nSt,nspin)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: i,a,p,jb
|
||||
double precision :: eps
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: Z(nBas,nspin)
|
||||
|
||||
! Initialize
|
||||
|
||||
Z(:,:) = 0d0
|
||||
|
||||
!--------------!
|
||||
! Spin-up part !
|
||||
!--------------!
|
||||
|
||||
! Occupied part of the renormalization factor
|
||||
|
||||
do p=nC(1)+1,nBas-nR(1)
|
||||
do i=nC(1)+1,nO(1)
|
||||
do jb=1,nSt
|
||||
eps = e(p,1) - e(i,1) + Omega(jb)
|
||||
Z(p,1) = Z(p,1) + rho(p,i,jb,1)**2*(eps/(eps**2 + eta**2))**2
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! Virtual part of the correlation self-energy
|
||||
|
||||
do p=nC(1)+1,nBas-nR(1)
|
||||
do a=nO(1)+1,nBas-nR(1)
|
||||
do jb=1,nSt
|
||||
eps = e(p,1) - e(a,1) - Omega(jb)
|
||||
Z(p,1) = Z(p,1) + rho(p,a,jb,1)**2*(eps/(eps**2 + eta**2))**2
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
!----------------!
|
||||
! Spin-down part !
|
||||
!----------------!
|
||||
|
||||
! Occupied part of the correlation self-energy
|
||||
|
||||
do p=nC(2)+1,nBas-nR(2)
|
||||
do i=nC(2)+1,nO(2)
|
||||
do jb=1,nSt
|
||||
eps = e(p,2) - e(i,2) + Omega(jb)
|
||||
Z(p,2) = Z(p,2) + rho(p,i,jb,2)**2*(eps/(eps**2 + eta**2))**2
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! Virtual part of the correlation self-energy
|
||||
|
||||
do p=nC(2)+1,nBas-nR(2)
|
||||
do a=nO(2)+1,nBas-nR(2)
|
||||
do jb=1,nSt
|
||||
eps = e(p,2) - e(a,2) - Omega(jb)
|
||||
Z(p,2) = Z(p,2) + rho(p,a,jb,2)**2*(eps/(eps**2 + eta**2))**2
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! Final rescaling
|
||||
|
||||
Z(:,:) = 1d0/(1d0 + Z(:,:))
|
||||
|
||||
end subroutine unrestricted_renormalization_factor
|
||||
|
|
@ -1,118 +0,0 @@
|
|||
subroutine unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,SigC,EcGM)
|
||||
|
||||
! Compute diagonal of the correlation part of the self-energy
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
double precision,intent(in) :: eta
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC(nspin)
|
||||
integer,intent(in) :: nO(nspin)
|
||||
integer,intent(in) :: nV(nspin)
|
||||
integer,intent(in) :: nR(nspin)
|
||||
integer,intent(in) :: nSt
|
||||
double precision,intent(in) :: e(nBas,nspin)
|
||||
double precision,intent(in) :: Omega(nSt)
|
||||
double precision,intent(in) :: rho(nBas,nBas,nSt,nspin)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: i,a,p,q,jb
|
||||
double precision :: eps
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: SigC(nBas,nBas,nspin)
|
||||
double precision :: EcGM(nspin)
|
||||
|
||||
! Initialize
|
||||
|
||||
SigC(:,:,:) = 0d0
|
||||
EcGM(:) = 0d0
|
||||
|
||||
!--------------!
|
||||
! Spin-up part !
|
||||
!--------------!
|
||||
|
||||
! Occupied part of the correlation self-energy
|
||||
|
||||
do p=nC(1)+1,nBas-nR(1)
|
||||
do q=nC(1)+1,nBas-nR(1)
|
||||
do i=nC(1)+1,nO(1)
|
||||
do jb=1,nSt
|
||||
eps = e(p,1) - e(i,1) + Omega(jb)
|
||||
SigC(p,q,1) = SigC(p,q,1) + rho(p,i,jb,1)*rho(q,i,jb,1)*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! Virtual part of the correlation self-energy
|
||||
|
||||
do p=nC(1)+1,nBas-nR(1)
|
||||
do q=nC(1)+1,nBas-nR(1)
|
||||
do a=nO(1)+1,nBas-nR(1)
|
||||
do jb=1,nSt
|
||||
eps = e(p,1) - e(a,1) - Omega(jb)
|
||||
SigC(p,q,1) = SigC(p,q,1) + rho(p,a,jb,1)*rho(q,a,jb,1)*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! GM correlation energy
|
||||
|
||||
do i=nC(1)+1,nO(1)
|
||||
do a=nO(1)+1,nBas-nR(1)
|
||||
do jb=1,nSt
|
||||
eps = e(a,1) - e(i,1) + Omega(jb)
|
||||
EcGM(1) = EcGM(1) - rho(a,i,jb,1)**2*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
!----------------!
|
||||
! Spin-down part !
|
||||
!----------------!
|
||||
|
||||
! Occupied part of the correlation self-energy
|
||||
|
||||
do p=nC(2)+1,nBas-nR(2)
|
||||
do q=nC(2)+1,nBas-nR(2)
|
||||
do i=nC(2)+1,nO(2)
|
||||
do jb=1,nSt
|
||||
eps = e(p,2) - e(i,2) + Omega(jb)
|
||||
SigC(p,q,2) = SigC(p,q,2) + rho(p,i,jb,2)*rho(q,i,jb,2)*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! Virtual part of the correlation self-energy
|
||||
|
||||
do p=nC(2)+1,nBas-nR(2)
|
||||
do q=nC(2)+1,nBas-nR(2)
|
||||
do a=nO(2)+1,nBas-nR(2)
|
||||
do jb=1,nSt
|
||||
eps = e(p,2) - e(a,2) - Omega(jb)
|
||||
SigC(p,q,2) = SigC(p,q,2) + rho(p,a,jb,2)*rho(q,a,jb,2)*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! GM correlation energy
|
||||
|
||||
do i=nC(2)+1,nO(2)
|
||||
do a=nO(2)+1,nBas-nR(2)
|
||||
do jb=1,nSt
|
||||
eps = e(a,2) - e(i,2) + Omega(jb)
|
||||
EcGM(2) = EcGM(2) - rho(a,i,jb,2)**2*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
end subroutine unrestricted_self_energy_correlation
|
||||
|
|
@ -1,111 +0,0 @@
|
|||
subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,SigC,EcGM)
|
||||
|
||||
! Compute diagonal of the correlation part of the self-energy
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
double precision,intent(in) :: eta
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC(nspin)
|
||||
integer,intent(in) :: nO(nspin)
|
||||
integer,intent(in) :: nV(nspin)
|
||||
integer,intent(in) :: nR(nspin)
|
||||
integer,intent(in) :: nSt
|
||||
double precision,intent(in) :: e(nBas,nspin)
|
||||
double precision,intent(in) :: Omega(nSt)
|
||||
double precision,intent(in) :: rho(nBas,nBas,nSt,nspin)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: i,a,p,q,jb
|
||||
double precision :: eps
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: SigC(nBas,nspin)
|
||||
double precision :: EcGM(nspin)
|
||||
|
||||
! Initialize
|
||||
|
||||
SigC(:,:) = 0d0
|
||||
EcGM(:) = 0d0
|
||||
|
||||
!--------------!
|
||||
! Spin-up part !
|
||||
!--------------!
|
||||
|
||||
! Occupied part of the correlation self-energy
|
||||
|
||||
do p=nC(1)+1,nBas-nR(1)
|
||||
do i=nC(1)+1,nO(1)
|
||||
do jb=1,nSt
|
||||
eps = e(p,1) - e(i,1) + Omega(jb)
|
||||
SigC(p,1) = SigC(p,1) + rho(p,i,jb,1)**2*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! Virtual part of the correlation self-energy
|
||||
|
||||
do p=nC(1)+1,nBas-nR(1)
|
||||
do a=nO(1)+1,nBas-nR(1)
|
||||
do jb=1,nSt
|
||||
eps = e(p,1) - e(a,1) - Omega(jb)
|
||||
SigC(p,1) = SigC(p,1) + rho(p,a,jb,1)**2*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! GM correlation energy
|
||||
|
||||
do i=nC(1)+1,nO(1)
|
||||
do a=nO(1)+1,nBas-nR(1)
|
||||
do jb=1,nSt
|
||||
eps = e(a,1) - e(i,1) + Omega(jb)
|
||||
EcGM(1) = EcGM(1) - rho(a,i,jb,1)**2*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
!----------------!
|
||||
! Spin-down part !
|
||||
!----------------!
|
||||
|
||||
! Occupied part of the correlation self-energy
|
||||
|
||||
do p=nC(2)+1,nBas-nR(2)
|
||||
do i=nC(2)+1,nO(2)
|
||||
do jb=1,nSt
|
||||
eps = e(p,2) - e(i,2) + Omega(jb)
|
||||
SigC(p,2) = SigC(p,2) + rho(p,i,jb,2)**2*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! Virtual part of the correlation self-energy
|
||||
|
||||
do p=nC(2)+1,nBas-nR(2)
|
||||
do a=nO(2)+1,nBas-nR(2)
|
||||
do jb=1,nSt
|
||||
eps = e(p,2) - e(a,2) - Omega(jb)
|
||||
SigC(p,2) = SigC(p,2) + rho(p,a,jb,2)**2*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! GM correlation energy
|
||||
|
||||
do i=nC(2)+1,nO(2)
|
||||
do a=nO(2)+1,nBas-nR(2)
|
||||
do jb=1,nSt
|
||||
eps = e(a,2) - e(i,2) + Omega(jb)
|
||||
EcGM(2) = EcGM(2) - rho(a,i,jb,2)**2*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
|
||||
end subroutine unrestricted_self_energy_correlation_diag
|
||||
|
|
@ -57,8 +57,8 @@ subroutine phULR(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,e,
|
|||
call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
|
||||
|
||||
if(BSE) &
|
||||
call unrestricted_Bethe_Salpeter_A_matrix(ispin,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,e, &
|
||||
ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,Aph)
|
||||
call UGW_phBSE_static_kernel_A(ispin,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,e, &
|
||||
ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,Aph)
|
||||
|
||||
! Tamm-Dancoff approximation
|
||||
|
||||
|
|
@ -75,8 +75,8 @@ subroutine phULR(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,e,
|
|||
call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
|
||||
|
||||
if(BSE) &
|
||||
call unrestricted_Bethe_Salpeter_B_matrix(ispin,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, &
|
||||
ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,Bph)
|
||||
call UGW_phBSE_static_kernel_B(ispin,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, &
|
||||
ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,Bph)
|
||||
|
||||
! Build A + B and A - B matrices
|
||||
|
||||
|
|
|
|||
|
|
@ -94,7 +94,7 @@ subroutine phUACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spin
|
|||
|
||||
call phULR(isp_W,.true.,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,eW, &
|
||||
ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
|
||||
call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
|
||||
call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
|
||||
|
||||
! Spin-conserved manifold
|
||||
|
||||
|
|
@ -121,7 +121,7 @@ subroutine phUACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spin
|
|||
|
||||
call phULR(isp_W,.true.,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, &
|
||||
ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
|
||||
call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
|
||||
call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
|
||||
|
||||
end if
|
||||
|
||||
|
|
@ -175,7 +175,7 @@ subroutine phUACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spin
|
|||
|
||||
call phULR(isp_W,.true.,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, &
|
||||
ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
|
||||
call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
|
||||
call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
|
||||
|
||||
end if
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue