PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 12:23:42 +01:00
Code Releases Activity

more cleanup

Browse Source
This commit is contained in:
Pierre-Francois Loos 2024-09-10 09:40:22 +02:00
parent 65ba14706c
commit 826d80d594
6 changed files with 135 additions and 130 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

41
src/GW/RG0W0.f90
View File

@ -1,5 +1,5 @@
subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
! Perform G0W0 calculation
@ -28,6 +28,7 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
logical,intent(in) :: regularize
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
@ -35,9 +36,9 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart)
double precision,intent(in) :: eHF(nOrb)
! Local variables
@ -88,15 +89,15 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
! Memory allocation
allocate(Aph(nS,nS),Bph(nS,nS),SigC(nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS), &
eGW(nBas),eGWlin(nBas))
allocate(Aph(nS,nS),Bph(nS,nS),SigC(nOrb),Z(nOrb),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nOrb,nOrb,nS), &
eGW(nOrb),eGWlin(nOrb))
!-------------------!
! Compute screening !
!-------------------!
call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,Aph)
if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,lambda,eHF,ERI,Aph)
if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
@ -106,15 +107,15 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
! Compute spectral weights !
!--------------------------!
call RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho)
call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho)
!------------------------!
! Compute GW self-energy !
!------------------------!
if(regularize) call GW_regularization(nBas,nC,nO,nV,nR,nS,eHF,Om,rho)
if(regularize) call GW_regularization(nOrb,nC,nO,nV,nR,nS,eHF,Om,rho)
call RGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z)
call RGW_self_energy_diag(eta,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z)
!-----------------------------------!
! Solve the quasi-particle equation !
@ -136,24 +137,24 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
write(*,*) ' *** Quasiparticle energies obtained by root search *** '
write(*,*)
call RGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eHF,eGW,Z)
call RGW_QP_graph(eta,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eHF,eGW,Z)
end if
! Plot self-energy, renormalization factor, and spectral function
if(plot_self) call RGW_plot_self_energy(nBas,eta,nC,nO,nV,nR,nS,eHF,eHF,Om,rho)
if(plot_self) call RGW_plot_self_energy(nOrb,eta,nC,nO,nV,nR,nS,eHF,eHF,Om,rho)
!--------------------!
! Cumulant expansion !
!--------------------!
! call RGWC(dotest,eta,nBas,nC,nO,nV,nR,nS,Om,rho,eHF,eHF,eGW,Z)
! call RGWC(dotest,eta,nOrb,nC,nO,nV,nR,nS,Om,rho,eHF,eHF,eGW,Z)
! Compute the RPA correlation energy
call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,eGW,ERI,Aph)
if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,lambda,eGW,ERI,Aph)
if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,lambda,ERI,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
@ -161,14 +162,14 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
! Dump results !
!--------------!
call print_RG0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM)
call print_RG0W0(nOrb,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM)
! Perform BSE calculation
if(dophBSE) then
call RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta, &
nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE)
nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
@ -188,7 +189,7 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
write(*,*) '-------------------------------------------------------------'
write(*,*)
call RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,EcBSE)
call RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,eHF,eGW,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
@ -205,7 +206,7 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
if(doppBSE) then
call RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE)
call RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'

8
src/GW/RGW.f90
View File

@ -77,7 +77,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
call wall_time(start_GW)
call RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
@ -94,7 +94,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
call wall_time(start_GW)
call evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
singlet,triplet,linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
@ -150,7 +150,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
call wall_time(start_GW)
! TODO
call ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call ufG0W0(dotest,TDA_W,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
@ -167,7 +167,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
call wall_time(start_GW)
! TODO
call ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call ufRGW(dotest,TDA_W,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW

39
src/GW/evRGW.f90
View File

@ -1,5 +1,5 @@
subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
singlet,triplet,linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
! Perform self-consistent eigenvalue-only GW calculation
@ -32,14 +32,15 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
logical,intent(in) :: regularize
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart)
! Local variables
@ -89,8 +90,8 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
! Memory allocation
allocate(Aph(nS,nS),Bph(nS,nS),eGW(nBas),eOld(nBas),Z(nBas),SigC(nBas), &
Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS),error_diis(nBas,max_diis),e_diis(nBas,max_diis))
allocate(Aph(nS,nS),Bph(nS,nS),eGW(nOrb),eOld(nOrb),Z(nOrb),SigC(nOrb), &
Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nOrb,nOrb,nS),error_diis(nOrb,max_diis),e_diis(nOrb,max_diis))
! Initialization
@ -113,20 +114,20 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
! Compute screening
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
if(.not.TDA_W) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
! Compute spectral weights
call RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho)
call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho)
! Compute correlation part of the self-energy
if(regularize) call GW_regularization(nBas,nC,nO,nV,nR,nS,eGW,Om,rho)
if(regularize) call GW_regularization(nOrb,nC,nO,nV,nR,nS,eGW,Om,rho)
call RGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
call RGW_self_energy_diag(eta,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
! Solve the quasi-particle equation
@ -142,7 +143,7 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
write(*,*) ' *** Quasiparticle energies obtained by root search *** '
write(*,*)
call RGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eOld,eOld,eGW,Z)
call RGW_QP_graph(eta,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,eOld,eOld,eGW,Z)
end if
@ -152,7 +153,7 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
! Print results
call print_evRGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM)
call print_evRGW(nOrb,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM)
! Linear mixing or DIIS extrapolation
@ -164,7 +165,7 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
n_diis = min(n_diis+1,max_diis)
if(abs(rcond) > 1d-7) then
call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGW-eOld,eGW)
call DIIS_extrapolation(rcond,nOrb,nOrb,n_diis,error_diis,e_diis,eGW-eOld,eGW)
else
n_diis = 0
end if
@ -203,8 +204,8 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
!--------------------!
! TODO
!call RGWC(dotest, eta, nBas, nC, nO, nV, nR, nS, Om, rho, eHF, eGW, eGW, Z)
call RGWC(dotest, eta, nBas, nC, nO, nV, nR, nS, Om, rho, eHF, eHF, eGW, Z)
!call RGWC(dotest, eta, nOrb, nC, nO, nV, nR, nS, Om, rho, eHF, eGW, eGW, Z)
call RGWC(dotest, eta, nOrb, nC, nO, nV, nR, nS, Om, rho, eHF, eHF, eGW, Z)
! Deallocate memory
@ -215,7 +216,7 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
if(dophBSE) then
call RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta, &
nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGW,eGW,EcBSE)
nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eGW,eGW,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
@ -235,7 +236,7 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
write(*,*) '-----------------------------------------------------------'
write(*,*)
call RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcBSE)
call RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
@ -252,7 +253,7 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
if(doppBSE) then
call RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE)
call RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'

27
src/GW/ufBSE.f90
View File

@ -1,4 +1,4 @@
subroutine ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
subroutine ufBSE(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
! Unfold BSE@GW equations
@ -8,6 +8,7 @@ subroutine ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
@ -15,9 +16,9 @@ subroutine ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: eGW(nBas)
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: eGW(nOrb)
! Local variables
@ -84,12 +85,12 @@ subroutine ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
ia = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
do a=nO+1,nOrb-nR
ia = ia + 1
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
do b=nO+1,nOrb-nR
jb = jb + 1
H(ia,jb) = (eGW(a) - eGW(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
@ -107,14 +108,14 @@ subroutine ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
iajb=0
do i=nC+1,nO
do a=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do j=nC+1,nO
do b=nO+1,nBas-nR
do b=nO+1,nOrb-nR
iajb = iajb + 1
kc = 0
do k=nC+1,nO
do c=nO+1,nBas-nR
do c=nO+1,nOrb-nR
kc = kc + 1
tmp = sqrt(2d0)*Kronecker_delta(k,j)*ERI(b,a,c,i)
@ -139,16 +140,16 @@ subroutine ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
iajb = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do j=nC+1,nO
do b=nO+1,nBas-nR
do b=nO+1,nOrb-nR
iajb = iajb + 1
kcld = 0
do k=nC+1,nO
do c=nO+1,nBas-nR
do c=nO+1,nOrb-nR
do l=nC+1,nO
do d=nO+1,nBas-nR
do d=nO+1,nOrb-nR
kcld = kcld + 1
tmp = ((eHF(a) + eGW(b) - eHF(i) - eGW(j))*Kronecker_delta(i,k)*Kronecker_delta(a,c) &

45
src/GW/ufG0W0.f90
View File

@ -1,4 +1,4 @@
subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
subroutine ufG0W0(dotest,TDA_W,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
! Unfold G0W0 equations
@ -11,6 +11,7 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
logical,intent(in) :: TDA_W
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
@ -18,8 +19,8 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eHF(nOrb)
! Local variables
@ -93,10 +94,10 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
! Memory allocation
allocate(Om(nS),Aph(nS,nS),Bph(nS,nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS))
allocate(Om(nS),Aph(nS,nS),Bph(nS,nS),XpY(nS,nS),XmY(nS,nS),rho(nOrb,nOrb,nS))
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
@ -106,7 +107,7 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
! Compute spectral weights !
!--------------------------!
call RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho)
call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho)
deallocate(Aph,Bph,XpY,XmY)
@ -154,7 +155,7 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
ija = 0
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do a=nO+1,nOrb-nR
ija = ija + 1
H(1 ,1+ija) = sqrt(2d0)*ERI(p,a,i,j)
@ -170,8 +171,8 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
iab = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
iab = iab + 1
H(1 ,1+n2h1p+iab) = sqrt(2d0)*ERI(p,i,b,a)
@ -188,13 +189,13 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
ija = 0
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do a=nO+1,nOrb-nR
ija = ija + 1
klc = 0
do k=nC+1,nO
do l=nC+1,nO
do c=nO+1,nBas-nR
do c=nO+1,nOrb-nR
klc = klc + 1
H(1+ija,1+klc) &
@ -215,14 +216,14 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
iab = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
iab = iab + 1
kcd = 0
do k=nC+1,nO
do c=nO+1,nBas-nR
do d=nO+1,nBas-nR
do c=nO+1,nOrb-nR
do d=nO+1,nOrb-nR
kcd = kcd + 1
H(1+n2h1p+iab,1+n2h1p+kcd) &
@ -304,7 +305,7 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
iab = 0
do ia=1,nS
do b=nO+1,nBas-nR
do b=nO+1,nOrb-nR
iab = iab + 1
H(1+n2h1p+iab,1+n2h1p+iab) = eHF(b) + Om(ia)
@ -318,7 +319,7 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
iab = 0
do ia=1,nS
do b=nO+1,nBas-nR
do b=nO+1,nOrb-nR
iab = iab + 1
H(1 ,1+n2h1p+iab) = sqrt(2d0)*rho(p,b,ia)
@ -409,7 +410,7 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
ija = 0
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do a=nO+1,nOrb-nR
ija = ija + 1
if(abs(H(1+ija,s)) > cutoff2) &
@ -422,8 +423,8 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
iab = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
iab = iab + 1
if(abs(H(1+n2h1p+iab,s)) > cutoff2) &
@ -478,7 +479,7 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
iab = 0
do ia=1,nS
do b=nO+1,nBas-nR
do b=nO+1,nOrb-nR
iab = iab + 1
if(abs(H(1+n2h1p+iab,s)) > cutoff2) &

105
src/GW/ufGW.f90
View File

@ -1,4 +1,4 @@
subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
subroutine ufRGW(dotest,TDA_W,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
! Unfold GW equations
@ -11,6 +11,7 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
logical,intent(in) :: TDA_W
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
@ -18,8 +19,8 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eHF(nOrb)
! Local variables
@ -68,7 +69,7 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
n2h1p = nO*nO*nV
n2p1h = nV*nV*nO
nH = nBas + n2h1p + n2p1h
nH = nOrb + n2h1p + n2p1h
! Memory allocation
@ -89,14 +90,14 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
! Memory allocation
allocate(Om(nS),Aph(nS,nS),Bph(nS,nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS))
allocate(Om(nS),Aph(nS,nS),Bph(nS,nS),XpY(nS,nS),XmY(nS,nS),rho(nOrb,nOrb,nS))
! Spin manifold
ispin = 1
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
@ -106,7 +107,7 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
! Compute spectral weights !
!--------------------------!
call RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho)
call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho)
deallocate(Aph,Bph,XpY,XmY)
@ -141,7 +142,7 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
! Block F !
!---------!
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
H(p,p) = eHF(p)
end do
@ -149,16 +150,16 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
! Block V2h1p !
!-------------!
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
ija = 0
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do a=nO+1,nOrb-nR
ija = ija + 1
H(p ,nBas+ija) = sqrt(2d0)*ERI(p,a,i,j)
H(nBas+ija,p ) = sqrt(2d0)*ERI(p,a,i,j)
H(p ,nOrb+ija) = sqrt(2d0)*ERI(p,a,i,j)
H(nOrb+ija,p ) = sqrt(2d0)*ERI(p,a,i,j)
end do
end do
@ -170,16 +171,16 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
! Block V2p1h !
!-------------!
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
iab = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
iab = iab + 1
H(p ,nBas+n2h1p+iab) = sqrt(2d0)*ERI(p,i,b,a)
H(nBas+n2h1p+iab,p ) = sqrt(2d0)*ERI(p,i,b,a)
H(p ,nOrb+n2h1p+iab) = sqrt(2d0)*ERI(p,i,b,a)
H(nOrb+n2h1p+iab,p ) = sqrt(2d0)*ERI(p,i,b,a)
end do
end do
@ -194,16 +195,16 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
ija = 0
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do a=nO+1,nOrb-nR
ija = ija + 1
klc = 0
do k=nC+1,nO
do l=nC+1,nO
do c=nO+1,nBas-nR
do c=nO+1,nOrb-nR
klc = klc + 1
H(nBas+ija,nBas+klc) &
H(nOrb+ija,nOrb+klc) &
= ((eHF(i) + eHF(j) - eHF(a))*Kronecker_delta(j,l)*Kronecker_delta(a,c) &
- 2d0*ERI(j,c,a,l))*Kronecker_delta(i,k)
@ -221,17 +222,17 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
iab = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
iab = iab + 1
kcd = 0
do k=nC+1,nO
do c=nO+1,nBas-nR
do d=nO+1,nBas-nR
do c=nO+1,nOrb-nR
do d=nO+1,nOrb-nR
kcd = kcd + 1
H(nBas+n2h1p+iab,nBas+n2h1p+kcd) &
H(nOrb+n2h1p+iab,nOrb+n2h1p+kcd) &
= ((eHF(a) + eHF(b) - eHF(i))*Kronecker_delta(i,k)*Kronecker_delta(a,c) &
+ 2d0*ERI(a,k,i,c))*Kronecker_delta(b,d)
@ -275,7 +276,7 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
! Block F !
!---------!
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
H(p,p) = eHF(p)
end do
@ -283,15 +284,15 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
! Block W2h1p !
!-------------!
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
ija = 0
do i=nC+1,nO
do ja=1,nS
ija = ija + 1
H(p ,nBas+ija) = sqrt(2d0)*rho(p,i,ja)
H(nBas+ija,p ) = sqrt(2d0)*rho(p,i,ja)
H(p ,nOrb+ija) = sqrt(2d0)*rho(p,i,ja)
H(nOrb+ija,p ) = sqrt(2d0)*rho(p,i,ja)
end do
end do
@ -302,15 +303,15 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
! Block W2p1h !
!-------------!
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
iab = 0
do ia=1,nS
do b=nO+1,nBas-nR
do b=nO+1,nOrb-nR
iab = iab + 1
H(p ,nBas+n2h1p+iab) = sqrt(2d0)*rho(p,b,ia)
H(nBas+n2h1p+iab,p ) = sqrt(2d0)*rho(p,b,ia)
H(p ,nOrb+n2h1p+iab) = sqrt(2d0)*rho(p,b,ia)
H(nOrb+n2h1p+iab,p ) = sqrt(2d0)*rho(p,b,ia)
end do
end do
@ -326,7 +327,7 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
do ja=1,nS
ija = ija + 1
H(nBas+ija,nBas+ija) = eHF(i) - Om(ja)
H(nOrb+ija,nOrb+ija) = eHF(i) - Om(ja)
end do
end do
@ -337,10 +338,10 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
iab = 0
do ia=1,nS
do b=nO+1,nBas-nR
do b=nO+1,nOrb-nR
iab = iab + 1
H(nBas+n2h1p+iab,nBas+n2h1p+iab) = eHF(b) + Om(ia)
H(nOrb+n2h1p+iab,nOrb+n2h1p+iab) = eHF(b) + Om(ia)
end do
end do
@ -375,7 +376,7 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
Z(:) = 0d0
do s=1,nH
do p=nC+1,nBas-nR
do p=nC+1,nOrb-nR
Z(s) = Z(s) + H(p,s)**2
end do
end do
@ -425,7 +426,7 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
write(*,'(1X,A7,I3,A16,1X,F15.6,1X,F15.6)') &
' (',p,') ',H(p,s),H(p,s)**2
end do
do p=nO+1,nBas-nR
do p=nO+1,nOrb-nR
if(abs(H(p,s)) > cutoff2) &
write(*,'(1X,A16,I3,A7,1X,F15.6,1X,F15.6)') &
' (',p,') ',H(p,s),H(p,s)**2
@ -434,12 +435,12 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
ija = 0
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do a=nO+1,nOrb-nR
ija = ija + 1
if(abs(H(nBas+ija,s)) > cutoff2) &
if(abs(H(nOrb+ija,s)) > cutoff2) &
write(*,'(1X,A3,I3,A1,I3,A6,I3,A7,1X,F15.6,1X,F15.6)') &
' (',i,',',j,') -> (',a,') ',H(nBas+ija,s),H(nBas+ija,s)**2
' (',i,',',j,') -> (',a,') ',H(nOrb+ija,s),H(nOrb+ija,s)**2
end do
end do
@ -447,13 +448,13 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
iab = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
iab = iab + 1
if(abs(H(nBas+n2h1p+iab,s)) > cutoff2) &
if(abs(H(nOrb+n2h1p+iab,s)) > cutoff2) &
write(*,'(1X,A7,I3,A6,I3,A1,I3,A3,1X,F15.6,1X,F15.6)') &
' (',i,') -> (',a,',',b,') ',H(nBas+n2h1p+iab,s),H(nBas+n2h1p+iab,s)**2
' (',i,') -> (',a,',',b,') ',H(nOrb+n2h1p+iab,s),H(nOrb+n2h1p+iab,s)**2
end do
end do
@ -487,7 +488,7 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
write(*,'(1X,A7,I3,A16,1X,F15.6,1X,F15.6)') &
' (',p,' ) ',H(p,s),H(p,s)**2
end do
do p=nO+1,nBas-nR
do p=nO+1,nOrb-nR
if(abs(H(p,s)) > cutoff2) &
write(*,'(1X,A7,I3,A16,1X,F15.6,1X,F15.6)') &
' (',p,' ) ',H(p,s),H(p,s)**2
@ -498,21 +499,21 @@ subroutine ufRGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
do ja=1,nS
ija = ija + 1
if(abs(H(nBas+ija,s)) > cutoff2) &
if(abs(H(nOrb+ija,s)) > cutoff2) &
write(*,'(1X,A7,I3,A1,I3,A12,1X,F15.6,1X,F15.6)') &
' (',i,',',ja,') ',H(nBas+ija,s),H(nBas+ija,s)**2
' (',i,',',ja,') ',H(nOrb+ija,s),H(nOrb+ija,s)**2
end do
end do
iab = 0
do ia=1,nS
do b=nO+1,nBas-nR
do b=nO+1,nOrb-nR
iab = iab + 1
if(abs(H(nBas+n2h1p+iab,s)) > cutoff2) &
if(abs(H(nOrb+n2h1p+iab,s)) > cutoff2) &
write(*,'(1X,A7,I3,A1,I3,A12,1X,F15.6,1X,F15.6)') &
' (',ia,',',b,') ',H(nBas+n2h1p+iab,s),H(nBas+n2h1p+iab,s)**2
' (',ia,',',b,') ',H(nOrb+n2h1p+iab,s),H(nOrb+n2h1p+iab,s)**2
end do
end do