LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-11-08 15:13:53 +01:00
Code Releases Activity

remove COHSEX and modify phLR in GW routines

Browse Source
This commit is contained in:
Pierre-Francois Loos 2023-07-18 14:59:18 +02:00
parent 092a3f5e6d
commit a2c38aa03b
58 changed files with 374 additions and 1774 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

80
src/AOtoMO/AOtoMO_integral_transform.f90
View File

@ -1,4 +1,4 @@
subroutine AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,c,ERI_AO_basis,ERI_MO_basis)
subroutine AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,c,ERI_AO,ERI_MO)
! AO to MO transformation of two-electron integrals via the semi-direct O(N^5) algorithm
! bra and ket are the spin of (bra1 bra2|ket1 ket2)
@ -11,7 +11,7 @@ subroutine AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,c,ERI_AO_basis,ERI
integer,intent(in) :: bra1,bra2
integer,intent(in) :: ket1,ket2
integer,intent(in) :: nBas
double precision,intent(in) :: ERI_AO_basis(nBas,nBas,nBas,nBas),c(nBas,nBas,nspin)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas),c(nBas,nBas,nspin)
! Local variables
@ -20,7 +20,7 @@ subroutine AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,c,ERI_AO_basis,ERI
! Output variables
double precision,intent(out) :: ERI_MO_basis(nBas,nBas,nBas,nBas)
double precision,intent(out) :: ERI_MO(nBas,nBas,nBas,nBas)
! Memory allocation
@ -28,76 +28,12 @@ subroutine AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,c,ERI_AO_basis,ERI
! Four-index transform via semi-direct O(N^5) algorithm
! scr(:,:,:,:) = 0d0
call dgemm ('T', 'N', nBas**3, nBas, nBas, 1d0, ERI_AO, nBas, c(1,1,bra2), size(c,1), 0d0, scr, nBas**3)
! ! do l=1,nBas
! ! do si=1,nBas
! ! do la=1,nBas
! ! do nu=1,nBas
! ! do mu=1,nBas
! ! scr(mu,nu,la,l) = scr(mu,nu,la,l) + ERI_AO_basis(mu,nu,la,si)*c(si,l,ket2)
! ! enddo
! ! enddo
! ! enddo
! ! enddo
! ! enddo
call dgemm ('T', 'N', nBas**3, nBas, nBas, 1d0, scr, nBas, c(1,1,bra1), size(c,1), 0d0, ERI_MO, nBas**3)
! call dgemm ('N', 'N', nBas*nBas*nBas, nBas, nBas, 1d0, ERI_AO_basis(1,1,1,1), &
! size(ERI_AO_basis,1)*size(ERI_AO_basis,2)*size(ERI_AO_basis,3), c(1,1,ket2), &
! size(c,1), 0d0, scr, size(scr,1)*size(scr,2)*size(scr,3) )
call dgemm ('T', 'N', nBas**3, nBas, nBas, 1d0, ERI_MO, nBas, c(1,1,ket1), size(c,1), 0d0, scr, nBas**3)
! ERI_MO_basis(:,:,:,:) = 0d0
call dgemm ('T', 'N', nBas**3, nBas, nBas, 1d0, scr, nBas, c(1,1,ket2), size(c,1), 0d0, ERI_MO, nBas**3)
! ! do l=1,nBas
! ! do la=1,nBas
! ! do nu=1,nBas
! ! do i=1,nBas
! ! do mu=1,nBas
! ! ERI_MO_basis(i,nu,la,l) = ERI_MO_basis(i,nu,la,l) + c(mu,i,bra1)*scr(mu,nu,la,l)
! ! enddo
! ! enddo
! ! enddo
! ! enddo
! ! enddo
! call dgemm ('T', 'N', nBas, nBas*nBas*nBas, nBas, 1d0, c(1,1,bra1), size(c,1), &
! scr(1,1,1,1), size(scr,1), 0d0, ERI_MO_basis, size(ERI_MO_basis,1))
! scr(:,:,:,:) = 0d0
! do l=1,nBas
! do k=1,nBas
! do la=1,nBas
! do nu=1,nBas
! do i=1,nBas
! scr(i,nu,k,l) = scr(i,nu,k,l) + ERI_MO_basis(i,nu,la,l)*c(la,k,bra2)
! enddo
! enddo
! enddo
! enddo
! enddo
! ERI_MO_basis(:,:,:,:) = 0d0
! do l=1,nBas
! do k=1,nBas
! do j=1,nBas
! do i=1,nBas
! do nu=1,nBas
! ERI_MO_basis(i,j,k,l) = ERI_MO_basis(i,j,k,l) + c(nu,j,ket1)*scr(i,nu,k,l)
! enddo
! ! write(11,'(I5,I5,I5,I5,F16.10)') i,j,k,l,ERI_MO_basis(i,j,k,l)
! enddo
! enddo
! enddo
! enddo
call dgemm ('T', 'N', nBas**3, nBas, nBas, 1d0, ERI_AO_basis, nBas, c(1,1,bra2), size(c,1), 0d0, scr, nBas**3)
call dgemm ('T', 'N', nBas**3, nBas, nBas, 1d0, scr, nBas, c(1,1,bra1), size(c,1), 0d0, ERI_MO_basis, nBas**3)
call dgemm ('T', 'N', nBas**3, nBas, nBas, 1d0, ERI_MO_basis, nBas, c(1,1,ket1), size(c,1), 0d0, scr, nBas**3)
call dgemm ('T', 'N', nBas**3, nBas, nBas, 1d0, scr, nBas, c(1,1,ket2), size(c,1), 0d0, ERI_MO_basis, nBas**3)
end subroutine AOtoMO_integral_transform
end subroutine

2
src/AOtoMO/AOtoMO_oooa.f90
View File

@ -82,4 +82,4 @@ subroutine AOtoMO_oooa(nBas,nO,nA,cO,cA,O,ooOoa)
deallocate(scr1,scr2)
end subroutine AOtoMO_oooa
end subroutine

2
src/AOtoMO/AOtoMO_oooo.f90
View File

@ -82,4 +82,4 @@ subroutine AOtoMO_oooo(nBas,nO,cO,O,ooOoo)
deallocate(scr1,scr2)
end subroutine AOtoMO_oooo
end subroutine

135
src/AOtoMO/AOtoMO_oooooo.f90
View File

@ -1,135 +0,0 @@
subroutine AOtoMO_oooooo(nBas,nO,cO,O,oooOooo)
! AO to MO transformation of three-electron integrals for the block oooooo
! Semi-direct O(N^7) algorithm
implicit none
! Input variables
integer,intent(in) :: nBas,nO
double precision,intent(in) :: cO(nBas,nO),O(nBas,nBas,nBas,nBas,nBas,nBas)
! Local variables
double precision,allocatable :: scr1(:,:,:,:,:,:),scr2(:,:,:,:,:,:)
integer :: mu,nu,la,si,ka,ta,i,j,k,l,m,n
! Output variables
double precision,intent(out) :: oooOooo(nO,nO,nO,nO,nO,nO)
! Memory allocation
allocate(scr1(nBas,nBas,nBas,nBas,nBas,nBas),scr2(nBas,nBas,nBas,nBas,nBas,nBas))
write(*,*)
write(*,'(A42)') '------------------------------------------'
write(*,'(A42)') ' AO to MO transformation for oooooo block '
write(*,'(A42)') '------------------------------------------'
write(*,*)
scr1 = 0d0
do mu=1,nBas
do nu=1,nBas
do la=1,nBas
do si=1,nBas
do ka=1,nBas
do ta=1,nBas
do n=1,nO
scr1(mu,nu,la,si,ka,n) = scr1(mu,nu,la,si,ka,n) + O(mu,nu,la,si,ka,ta)*cO(ta,n)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
scr2 = 0d0
do mu=1,nBas
do nu=1,nBas
do la=1,nBas
do si=1,nBas
do ka=1,nBas
do i=1,nO
do n=1,nO
scr2(i,nu,la,si,ka,n) = scr2(i,nu,la,si,ka,n) + cO(mu,i)*scr1(mu,nu,la,si,ka,n)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
scr1 = 0d0
do nu=1,nBas
do la=1,nBas
do si=1,nBas
do ka=1,nBas
do i=1,nO
do m=1,nO
do n=1,nO
scr1(i,nu,la,si,m,n) = scr1(i,nu,la,si,m,n) + scr2(i,nu,la,si,m,n)*cO(ka,m)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
scr2 = 0d0
do nu=1,nBas
do la=1,nBas
do si=1,nBas
do i=1,nO
do j=1,nO
do m=1,nO
do n=1,nO
scr2(i,j,la,si,m,n) = scr2(i,j,la,si,m,n) + cO(nu,j)*scr1(i,nu,la,si,m,n)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
scr1 = 0d0
do la=1,nBas
do si=1,nBas
do i=1,nO
do j=1,nO
do l=1,nO
do m=1,nO
do n=1,nO
scr1(i,j,la,l,m,n) = scr1(i,j,la,l,m,n) + scr2(i,j,la,si,m,n)*cO(si,l)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
oooOooo = 0d0
do si=1,nBas
do i=1,nO
do j=1,nO
do k=1,nO
do l=1,nO
do m=1,nO
do n=1,nO
oooOooo(i,j,k,l,m,n) = oooOooo(i,j,k,l,m,n) + cO(la,k)*scr1(i,j,la,k,m,n)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
deallocate(scr1,scr2)
end subroutine AOtoMO_oooooo

2
src/AOtoMO/AOtoMO_oovv.f90
View File

@ -74,4 +74,4 @@ subroutine AOtoMO_oovv(nBas,nO,nV,cO,cV,O,ooOvv)
enddo
enddo
end subroutine AOtoMO_oovv
end subroutine

2
src/AOtoMO/AOtoMO_transform.f90
View File

@ -15,4 +15,4 @@ subroutine AOtoMO_transform(nBas,c,A)
A = matmul(transpose(c),matmul(A,c))
end subroutine AOtoMO_transform
end subroutine

2
src/AOtoMO/Coulomb_matrix_AO_basis.f90
View File

@ -31,4 +31,4 @@ subroutine Coulomb_matrix_AO_basis(nBas,P,G,J)
enddo
end subroutine Coulomb_matrix_AO_basis
end subroutine

2
src/AOtoMO/Coulomb_matrix_MO_basis.f90
View File

@ -23,4 +23,4 @@ subroutine Coulomb_matrix_MO_basis(nBas,c,P,G,J)
J = matmul(transpose(c),matmul(J,c))
end subroutine Coulomb_matrix_MO_basis
end subroutine

2
src/AOtoMO/Hartree_matrix_AO_basis.f90
View File

@ -30,4 +30,4 @@ subroutine Hartree_matrix_AO_basis(nBas,P,Hc,G,H)
enddo
enddo
end subroutine Hartree_matrix_AO_basis
end subroutine

2
src/AOtoMO/Hartree_matrix_MO_basis.f90
View File

@ -23,4 +23,4 @@ subroutine Hartree_matrix_MO_basis(nBas,c,P,Hc,G,H)
H = matmul(transpose(c),matmul(H,c))
end subroutine Hartree_matrix_MO_basis
end subroutine

2
src/AOtoMO/MOtoAO_transform.f90
View File

@ -24,4 +24,4 @@ subroutine MOtoAO_transform(nBas,S,c,A)
Sc = matmul(S,c)
A = matmul(Sc,matmul(A,transpose(Sc)))
end subroutine MOtoAO_transform
end subroutine

2
src/AOtoMO/exchange_matrix_AO_basis.f90
View File

@ -30,4 +30,4 @@ subroutine exchange_matrix_AO_basis(nBas,P,ERI,K)
enddo
enddo
end subroutine exchange_matrix_AO_basis
end subroutine

2
src/CI/CID.f90
View File

@ -215,4 +215,4 @@ subroutine CID(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERIi
write(*,*)
endif
end subroutine CID
end subroutine

25
src/CI/CIS.f90
View File

@ -19,10 +19,11 @@ subroutine CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF)
logical :: dump_matrix = .false.
logical :: dump_trans = .false.
logical :: dRPA = .false.
integer :: ispin
integer :: maxS = 10
double precision :: lambda
double precision,allocatable :: A(:,:),Omega(:)
double precision,allocatable :: A(:,:),Om(:)
! Hello world
@ -38,14 +39,14 @@ subroutine CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF)
! Memory allocation
allocate(A(nS,nS),Omega(nS))
allocate(A(nS,nS),Om(nS))
! Compute CIS matrix
if(singlet) then
ispin = 1
call phLR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
if(dump_matrix) then
print*,'CIS matrix (singlet state)'
@ -53,9 +54,9 @@ subroutine CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF)
write(*,*)
endif
call diagonalize_matrix(nS,A,Omega)
call print_excitation('CIS ',ispin,nS,Omega)
call print_transition_vectors_ph(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega,transpose(A),transpose(A))
call diagonalize_matrix(nS,A,Om)
call print_excitation('CIS ',ispin,nS,Om)
call print_transition_vectors_ph(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,transpose(A),transpose(A))
if(dump_trans) then
print*,'Singlet CIS transition vectors'
@ -66,14 +67,14 @@ subroutine CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF)
! Compute CIS(D) correction
maxS = min(maxS,nS)
if(doCIS_D) call D_correction(ispin,nBas,nC,nO,nV,nR,nS,maxS,eHF,ERI,Omega(1:maxS),A(:,1:maxS))
if(doCIS_D) call D_correction(ispin,nBas,nC,nO,nV,nR,nS,maxS,eHF,ERI,Om(1:maxS),A(:,1:maxS))
endif
if(triplet) then
ispin = 2
call phLR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
if(dump_matrix) then
print*,'CIS matrix (triplet state)'
@ -81,9 +82,9 @@ subroutine CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF)
write(*,*)
endif
call diagonalize_matrix(nS,A,Omega)
call print_excitation('CIS ',ispin,nS,Omega)
call print_transition_vectors_ph(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega,transpose(A),transpose(A))
call diagonalize_matrix(nS,A,Om)
call print_excitation('CIS ',ispin,nS,Om)
call print_transition_vectors_ph(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,transpose(A),transpose(A))
if(dump_trans) then
print*,'Triplet CIS transition vectors'
@ -94,7 +95,7 @@ subroutine CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF)
! Compute CIS(D) correction
maxS = min(maxS,nS)
if(doCIS_D) call D_correction(ispin,nBas,nC,nO,nV,nR,nS,maxS,eHF,ERI,Omega(1:maxS),A(:,1:maxS))
if(doCIS_D) call D_correction(ispin,nBas,nC,nO,nV,nR,nS,maxS,eHF,ERI,Om(1:maxS),A(:,1:maxS))
endif

2
src/CI/CISD.f90
View File

@ -307,4 +307,4 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERI
write(*,*)
endif
end subroutine CISD
end subroutine

2
src/CI/D_correction.f90
View File

@ -274,4 +274,4 @@ subroutine D_correction(ispin,nBasin,nCin,nOin,nVin,nRin,nSin,maxS,eHF,ERI,w,X)
!------------------------------------------------------------------------
! End of loop over single excitations
!------------------------------------------------------------------------
end subroutine D_correction
end subroutine

2
src/CI/FCI.f90
View File

@ -28,4 +28,4 @@ subroutine FCI(nBas,nC,nO,nV,nR,ERI,e)
! Diagonalize FCI matrix
end subroutine FCI
end subroutine

2
src/CI/UCIS.f90
View File

@ -134,4 +134,4 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E
endif
end subroutine UCIS
end subroutine

2
src/GF/QP_graph_GF2.f90
View File

@ -73,4 +73,4 @@ subroutine QP_graph_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2lin,ERI,eGF2)
end do
end subroutine QP_graph_GF2
end subroutine

2
src/GF/regularized_self_energy_GF2_diag.f90
View File

@ -83,4 +83,4 @@ subroutine regularized_self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI
Z(:) = 1d0/(1d0 - Z(:))
end subroutine regularized_self_energy_GF2_diag
end subroutine

98
src/GW/Bethe_Salpeter_pp_so.f90
View File

@ -1,98 +0,0 @@
subroutine Bethe_Salpeter_pp_so(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE)
! Compute the Bethe-Salpeter excitation energies at the pp level
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: TDA_W
logical,intent(in) :: TDA
logical,intent(in) :: singlet
logical,intent(in) :: triplet
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: eW(nBas)
double precision,intent(in) :: eGW(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
integer :: ispin
integer :: isp_W
integer :: nOO
integer :: nVV
double precision :: EcRPA
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:)
double precision,allocatable :: Omega1(:)
double precision,allocatable :: X1(:,:)
double precision,allocatable :: Y1(:,:)
double precision,allocatable :: Omega2(:)
double precision,allocatable :: X2(:,:)
double precision,allocatable :: Y2(:,:)
double precision,allocatable :: WB(:,:)
double precision,allocatable :: WC(:,:)
double precision,allocatable :: WD(:,:)
! Output variables
double precision,intent(out) :: EcBSE(nspin)
!---------------------------------
! Compute RPA screening
!---------------------------------
isp_W = 3
EcRPA = 0d0
! Memory allocation
allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS))
call phLR(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
write(*,*) '****************'
write(*,*) '*** SpinOrbs ***'
write(*,*) '****************'
write(*,*)
ispin = 1
EcBSE(:) = 0d0
nOO = nO*(nO-1)/2
nVV = nV*(nV-1)/2
allocate(Omega1(nVV),X1(nVV,nVV),Y1(nOO,nVV), &
Omega2(nOO),X2(nVV,nOO),Y2(nOO,nOO), &
WB(nVV,nOO),WC(nVV,nVV),WD(nOO,nOO))
if(.not.TDA) call static_screening_WB_pp(4,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,WB)
call static_screening_WC_pp(4,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,WC)
call static_screening_WD_pp(4,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,WD)
! Compute BSE excitation energies
call linear_response_pp_BSE(4,TDA,.true.,nBas,nC,nO,nV,nR,nOO,nVV,1d0,eGW,ERI,WB,WC,WD, &
Omega1,X1,Y1,Omega2,X2,Y2,EcBSE(ispin))
call print_excitation('pp-BSE (N+2)',isp_W,nVV,Omega1)
call print_excitation('pp-BSE (N-2)',isp_W,nOO,Omega2)
end subroutine Bethe_Salpeter_pp_so

120
src/GW/G0W0.f90
View File

@ -1,6 +1,6 @@
subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE, &
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, &
ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc)
subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,doppBSE, &
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO, &
dipole_int,PHF,cHF,eHF,Vxc)
! Perform G0W0 calculation
@ -13,10 +13,9 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD
logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel
logical,intent(in) :: doXBS
logical,intent(in) :: COHSEX
logical,intent(in) :: BSE
logical,intent(in) :: BSE2
logical,intent(in) :: ppBSE
logical,intent(in) :: dophBSE
logical,intent(in) :: dophBSE2
logical,intent(in) :: doppBSE
logical,intent(in) :: TDA_W
logical,intent(in) :: TDA
logical,intent(in) :: dBSE
@ -47,30 +46,25 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD
! Local variables
logical :: print_W = .true.
logical :: dRPA
integer :: ispin
double precision :: EcRPA
double precision :: EcBSE(nspin)
double precision :: EcAC(nspin)
double precision :: EcppBSE(nspin)
double precision :: EcGM
double precision,allocatable :: Aph(:,:)
double precision,allocatable :: Bph(:,:)
double precision,allocatable :: SigX(:)
double precision,allocatable :: SigC(:)
double precision,allocatable :: Z(:)
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:)
double precision,allocatable :: Om(:)
double precision,allocatable :: XpY(:,:)
double precision,allocatable :: XmY(:,:)
double precision,allocatable :: rho(:,:,:)
double precision,allocatable :: eGW(:)
double precision,allocatable :: eGWlin(:)
integer :: nBas2
integer :: nC2
integer :: nO2
integer :: nV2
integer :: nR2
integer :: nS2
! Output variables
! Hello world
@ -83,15 +77,9 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD
! Initialization
dRPA = .true.
EcRPA = 0d0
! COHSEX approximation
if(COHSEX) then
write(*,*) 'COHSEX approximation activated!'
write(*,*)
end if
! TDA for W
if(TDA_W) then
@ -112,21 +100,25 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD
! Memory allocation
allocate(SigC(nBas),SigX(nBas),Z(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS),eGW(nBas),eGWlin(nBas))
allocate(Aph(nS,nS),Bph(nS,nS),SigC(nBas),SigX(nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS), &
eGW(nBas),eGWlin(nBas))
!-------------------!
! Compute screening !
!-------------------!
call phLR(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI_MO,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI_MO,Aph)
if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
if(print_W) call print_excitation('RPA@HF ',ispin,nS,OmRPA)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
if(print_W) call print_excitation('RPA@HF ',ispin,nS,Om)
!--------------------------!
! Compute spectral weights !
!--------------------------!
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA)
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY,rho)
!------------------------!
! Compute GW self-energy !
@ -136,12 +128,12 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD
if(regularize) then
call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
call regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
call regularized_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC)
call regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,Z)
else
call GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC,Z)
call GW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z)
end if
@ -165,17 +157,20 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD
write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
write(*,*)
call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,OmRPA,rho_RPA,eGWlin,eGW,regularize)
call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,Om,rho,eGWlin,eGW,regularize)
! Find all the roots of the QP equation if necessary
! call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,OmRPA,rho_RPA,eGWlin)
! call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Om,rho,eGWlin)
end if
! Compute the RPA correlation energy
call phLR(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph)
if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
!--------------!
! Dump results !
@ -185,17 +180,13 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD
! Deallocate memory
! deallocate(SigC,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,eGWlin)
! Plot stuff
! call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,OmRPA,rho_RPA)
deallocate(SigC,Z,Om,XpY,XmY,rho,eGWlin)
! Perform BSE calculation
if(BSE) then
if(dophBSE) then
call GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcBSE)
call GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcBSE)
if(exchange_kernel) then
@ -229,14 +220,14 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD
end if
call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eHF,eGW,EcAC)
call GW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eHF,eGW,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0 correlation energy (singlet) =',EcAC(1),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0 correlation energy (triplet) =',EcAC(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0 correlation energy =',EcAC(1) + EcAC(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0 total energy =',ENuc + ERHF + EcAC(1) + EcAC(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy (singlet) =',EcAC(1),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy (triplet) =',EcAC(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy =',EcAC(1) + EcAC(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 total energy =',ENuc + ERHF + EcAC(1) + EcAC(2),' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -244,40 +235,19 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD
end if
if(ppBSE) then
if(doppBSE) then
call GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcppBSE)
call GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy (singlet) =',EcppBSE(1),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy (triplet) =',3d0*EcppBSE(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy =',EcppBSE(1) + 3d0*EcppBSE(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 total energy =',ENuc + ERHF + EcppBSE(1) + 3d0*EcppBSE(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy (singlet) =',EcBSE(1),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy (triplet) =',3d0*EcBSE(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy =',EcBSE(1) + 3d0*EcBSE(2),' au'
write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 total energy =',ENuc + ERHF + EcBSE(1) + 3d0*EcBSE(2),' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! nBas2 = 2*nBas
! nO2 = 2*nO
! nV2 = 2*nV
! nC2 = 2*nC
! nR2 = 2*nR
! nS2 = nO2*nV2
!
! allocate(seHF(nBas2),seGW(nBas2),sERI(nBas2,nBas2,nBas2,nBas2))
!
! call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF)
! call spatial_to_spin_MO_energy(nBas,eGW,nBas2,seGW)
! call spatial_to_spin_ERI(nBas,ERI_MO,nBas2,sERI)
!
! call GW_ppBSE_so(TDA_W,TDA,singlet,triplet,eta,nBas2,nC2,nO2,nV2,nR2,nS2,sERI,dipole_int,seHF,seGW,EcppBSE)
end if
! if(BSE) call ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eGW)
! if(BSE) call ufXBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,OmRPA,rho_RPA)
if(BSE) call XBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcBSE)
end subroutine

156
src/GW/GW_self_energy.f90
View File

@ -1,4 +1,4 @@
subroutine GW_self_energy(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z)
subroutine GW_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z)
! Compute correlation part of the self-energy and the renormalization factor
@ -7,7 +7,6 @@ subroutine GW_self_energy(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z)
! Input variables
logical,intent(in) :: COHSEX
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
@ -37,113 +36,72 @@ subroutine GW_self_energy(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z)
Sig(:,:) = 0d0
Z(:) = 0d0
!-----------------------------!
! COHSEX static approximation !
!-----------------------------!
if(COHSEX) then
! COHSEX: SEX of the COHSEX correlation self-energy
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do i=nC+1,nO
do jb=1,nS
Sig(p,q) = Sig(p,q) + 4d0*rho(p,i,jb)*rho(q,i,jb)/Omega(jb)
end do
end do
end do
end do
! COHSEX: COH part of the COHSEX correlation self-energy
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do r=nC+1,nBas-nR
do jb=1,nS
Sig(p,q) = Sig(p,q) - 2d0*rho(p,r,jb)*rho(q,r,jb)/Omega(jb)
end do
end do
end do
end do
EcGM = 0d0
do i=nC+1,nO
EcGM = EcGM + 0.5d0*Sig(i,i)
end do
Z(:) = 0d0
else
!----------------!
! GW self-energy !
!----------------!
! Occupied part of the correlation self-energy
! Occupied part of the correlation self-energy
!$OMP PARALLEL &
!$OMP SHARED(Sig,rho,eta,nS,nC,nO,nBas,nR,e,Omega) &
!$OMP PRIVATE(jb,i,q,p,eps,num) &
!$OMP DEFAULT(NONE)
!$OMP DO
do q=nC+1,nBas-nR
do p=nC+1,nBas-nR
do jb=1,nS
do i=nC+1,nO
!$OMP PARALLEL &
!$OMP SHARED(Sig,rho,eta,nS,nC,nO,nBas,nR,e,Omega) &
!$OMP PRIVATE(jb,i,q,p,eps,num) &
!$OMP DEFAULT(NONE)
!$OMP DO
do q=nC+1,nBas-nR
do p=nC+1,nBas-nR
do jb=1,nS
do i=nC+1,nO
eps = e(p) - e(i) + Omega(jb)
num = 2d0*rho(p,i,jb)*rho(q,i,jb)
Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2)
if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do
end do
end do
!$OMP END DO
!$OMP END PARALLEL
! Virtual part of the correlation self-energy
!$OMP PARALLEL &
!$OMP SHARED(Sig,rho,eta,nS,nC,nO,nBas,nR,e,Omega) &
!$OMP PRIVATE(jb,a,q,p,eps,num) &
!$OMP DEFAULT(NONE)
!$OMP DO
do q=nC+1,nBas-nR
do p=nC+1,nBas-nR
do jb=1,nS
do a=nO+1,nBas-nR
eps = e(p) - e(a) - Omega(jb)
num = 2d0*rho(p,a,jb)*rho(q,a,jb)
Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2)
if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do
end do
end do
!$OMP END DO
!$OMP END PARALLEL
! Galitskii-Migdal correlation energy
EcGM = 0d0
do jb=1,nS
do a=nO+1,nBas-nR
do i=nC+1,nO
eps = e(a) - e(i) + Omega(jb)
num = 4d0*rho(a,i,jb)*rho(a,i,jb)
EcGM = EcGM - num*eps/(eps**2 + eta**2)
eps = e(p) - e(i) + Omega(jb)
num = 2d0*rho(p,i,jb)*rho(q,i,jb)
Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2)
if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do
end do
end do
!$OMP END DO
!$OMP END PARALLEL
! Virtual part of the correlation self-energy
!$OMP PARALLEL &
!$OMP SHARED(Sig,rho,eta,nS,nC,nO,nBas,nR,e,Omega) &
!$OMP PRIVATE(jb,a,q,p,eps,num) &
!$OMP DEFAULT(NONE)
!$OMP DO
do q=nC+1,nBas-nR
do p=nC+1,nBas-nR
do jb=1,nS
do a=nO+1,nBas-nR
eps = e(p) - e(a) - Omega(jb)
num = 2d0*rho(p,a,jb)*rho(q,a,jb)
Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2)
if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do
end do
end do
!$OMP END DO
!$OMP END PARALLEL
! Galitskii-Migdal correlation energy
EcGM = 0d0
do jb=1,nS
do a=nO+1,nBas-nR
do i=nC+1,nO
eps = e(a) - e(i) + Omega(jb)
num = 4d0*rho(a,i,jb)*rho(a,i,jb)
EcGM = EcGM - num*eps/(eps**2 + eta**2)
end do
end do
end if
end do
! Compute renormalization factor from derivative

98
src/GW/GW_self_energy_diag.f90
View File

@ -1,4 +1,4 @@
subroutine GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z)
subroutine GW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z)
! Compute diagonal of the correlation part of the self-energy and the renormalization factor
@ -7,7 +7,6 @@ subroutine GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,S
! Input variables
logical,intent(in) :: COHSEX
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
@ -35,91 +34,54 @@ subroutine GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,S
Sig(:) = 0d0
Z(:) = 0d0
!-----------------------------
! COHSEX static self-energy
!-----------------------------
!----------------!
! GW self-energy !
!----------------!
if(COHSEX) then
! Occupied part of the correlation self-energy
! COHSEX: SEX part of the COHSEX correlation self-energy
do p=nC+1,nBas-nR
do i=nC+1,nO
do jb=1,nS
Sig(p) = Sig(p) + 4d0*rho(p,i,jb)**2/Omega(jb)
end do
end do
end do
! COHSEX: COH part of the COHSEX correlation self-energy
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do jb=1,nS
Sig(p) = Sig(p) - 2d0*rho(p,q,jb)**2/Omega(jb)
end do
end do
end do
! GM correlation energy
EcGM = 0d0
do p=nC+1,nBas-nR
do i=nC+1,nO
EcGM = EcGM - Sig(i)
end do
do jb=1,nS
!-----------------------------
! GW self-energy
!-----------------------------
eps = e(p) - e(i) + Omega(jb)
num = 2d0*rho(p,i,jb)**2
Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
else
! Occupied part of the correlation self-energy
do p=nC+1,nBas-nR
do i=nC+1,nO
do jb=1,nS
eps = e(p) - e(i) + Omega(jb)
num = 2d0*rho(p,i,jb)**2
Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
Z(p) = Z(p) - 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
! Virtual part of the correlation self-energy
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
do jb=1,nS
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
do jb=1,nS
eps = e(p) - e(a) - Omega(jb)
num = 2d0*rho(p,a,jb)**2
Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
eps = e(p) - e(a) - Omega(jb)
num = 2d0*rho(p,a,jb)**2
Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do
end do
end do
! GM correlation energy
! Galitskii-Migdal correlation energy
EcGM = 0d0
do i=nC+1,nO
do a=nO+1,nBas-nR
do jb=1,nS
EcGM = 0d0
do i=nC+1,nO
do a=nO+1,nBas-nR
do jb=1,nS
eps = e(a) - e(i) + Omega(jb)
num = 4d0*rho(a,i,jb)**2
EcGM = EcGM - num*eps/(eps**2 + eta**2)
eps = e(a) - e(i) + Omega(jb)
num = 4d0*rho(a,i,jb)**2
EcGM = EcGM - num*eps/(eps**2 + eta**2)
end do
end do
end do
end if
end do
! Compute renormalization factor from derivative

100
src/GW/evGW.f90
View File

@ -1,4 +1,4 @@
subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE, &
subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,doppBSE, &
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF, &
cHF,eHF,Vxc)
@ -17,15 +17,14 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel
logical,intent(in) :: doXBS
logical,intent(in) :: COHSEX
logical,intent(in) :: BSE
logical,intent(in) :: BSE2
logical,intent(in) :: dophBSE
logical,intent(in) :: dophBSE2
logical,intent(in) :: TDA_W
logical,intent(in) :: TDA
logical,intent(in) :: dBSE
logical,intent(in) :: dTDA
logical,intent(in) :: evDyn
logical,intent(in) :: ppBSE
logical,intent(in) :: doppBSE
logical,intent(in) :: singlet
logical,intent(in) :: triplet
logical,intent(in) :: linearize
@ -49,6 +48,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
! Local variables
logical :: linear_mixing
logical :: dRPA = .true.
integer :: ispin
integer :: nSCF
integer :: n_diis
@ -58,10 +58,10 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
double precision :: EcRPA
double precision :: EcBSE(nspin)
double precision :: EcAC(nspin)
double precision :: EcppBSE(nspin)
double precision :: EcGM
double precision :: alpha
double precision :: Dpijb,Dpajb
double precision,allocatable :: Aph(:,:)
double precision,allocatable :: Bph(:,:)
double precision,allocatable :: error_diis(:,:)
double precision,allocatable :: e_diis(:,:)
double precision,allocatable :: eGW(:)
@ -69,21 +69,10 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
double precision,allocatable :: Z(:)
double precision,allocatable :: SigX(:)
double precision,allocatable :: SigC(:)
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:)
double precision,allocatable :: eGWlin(:)
! integer :: nBas2
! integer :: nC2
! integer :: nO2
! integer :: nV2
! integer :: nR2
! integer :: nS2
! double precision,allocatable :: seHF(:),seGW(:),sERI(:,:,:,:)
double precision,allocatable :: Om(:)
double precision,allocatable :: XpY(:,:)
double precision,allocatable :: XmY(:,:)
double precision,allocatable :: rho(:,:,:)
! Hello world
@ -93,13 +82,6 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
write(*,*)'************************************************'
write(*,*)
! COHSEX approximation
if(COHSEX) then
write(*,*) 'COHSEX approximation activated!'
write(*,*)
end if
! TDA for W
if(TDA_W) then
@ -121,8 +103,8 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
! Memory allocation
allocate(eGW(nBas),eOld(nBas),Z(nBas),SigX(nBas),SigC(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), &
rho_RPA(nBas,nBas,nS),error_diis(nBas,max_diis),e_diis(nBas,max_diis),eGWlin(nBas))
allocate(Aph(nS,nS),Bph(nS,nS),eGW(nBas),eOld(nBas),Z(nBas),SigX(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))
! Compute the exchange part of the self-energy
@ -141,8 +123,6 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
Z(:) = 1d0
rcond = 0d0
!------------------------------------------------------------------------
! Main loop
!------------------------------------------------------------------------
@ -151,28 +131,31 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
! Compute screening
call phLR(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph)
if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
! Compute spectral weights
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA)
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY,rho)
! Compute correlation part of the self-energy
if(regularize) then
call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC)
call renormalization_factor_SRG(eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z)
call regularized_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC)
call renormalization_factor_SRG(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,Z)
else
call GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC,Z)
call GW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
end if
! Solve the quasi-particle equation
eGWlin(:) = eHF(:) + SigX(:) + SigC(:) - Vxc(:)
eGW(:) = eHF(:) + SigX(:) + SigC(:) - Vxc(:)
! Linearized or graphical solution?
@ -181,14 +164,14 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
write(*,*)
eGW(:) = eGWlin(:)
eGW(:) = eGW(:)
else
write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
write(*,*)
call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,OmRPA,rho_RPA,eGWlin,eGW,regularize)
call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,Om,rho,eGW,eGW,regularize)
end if
@ -250,13 +233,13 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
! Deallocate memory
deallocate(eOld,Z,SigC,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,error_diis,e_diis)
deallocate(Aph,Bph,eOld,Z,SigC,Om,XpY,XmY,rho,error_diis,e_diis)
! Perform BSE calculation
if(BSE) then
if(dophBSE) then
call GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eGW,eGW,EcBSE)
call GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eGW,eGW,EcBSE)
if(exchange_kernel) then
@ -290,7 +273,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
end if
call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC)
call GW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
@ -305,34 +288,19 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
end if
if(ppBSE) then
if(doppBSE) then
call GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcppBSE)
call GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy (singlet) =',EcppBSE(1)
write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy (triplet) =',3d0*EcppBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy =',EcppBSE(1) + 3d0*EcppBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW total energy =',ENuc + ERHF + EcppBSE(1) + 3d0*EcppBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy (triplet) =',3d0*EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy =',EcBSE(1) + 3d0*EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW total energy =',ENuc + ERHF + EcBSE(1) + 3d0*EcBSE(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! nBas2 = 2*nBas
! nO2 = 2*nO
! nV2 = 2*nV
! nC2 = 2*nC
! nR2 = 2*nR
! nS2 = nO2*nV2
!
! allocate(seHF(nBas2),seGW(nBas2),sERI(nBas2,nBas2,nBas2,nBas2))
!
! call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF)
! call spatial_to_spin_MO_energy(nBas,eGW,nBas2,seGW)
! call spatial_to_spin_ERI(nBas,ERI_MO,nBas2,sERI)
!
! call GW_ppBSE_so(TDA_W,TDA,singlet,triplet,eta,nBas2,nC2,nO2,nV2,nR2,nS2,sERI,dipole_int,seHF,seGW,EcppBSE)
end if
end subroutine

66
src/GW/qsGW.f90
View File

@ -1,4 +1,4 @@
subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn, &
subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,doppBSE, &
singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO, &
ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
@ -15,14 +15,14 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel
logical,intent(in) :: doXBS
logical,intent(in) :: COHSEX
logical,intent(in) :: BSE
logical,intent(in) :: BSE2
logical,intent(in) :: dophBSE
logical,intent(in) :: dophBSE2
logical,intent(in) :: TDA_W
logical,intent(in) :: TDA
logical,intent(in) :: dBSE
logical,intent(in) :: dTDA
logical,intent(in) :: evDyn
logical,intent(in) :: doppBSE
logical,intent(in) :: singlet
logical,intent(in) :: triplet
double precision,intent(in) :: eta
@ -73,13 +73,16 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
double precision,external :: trace_matrix
double precision :: dipole(ncart)
logical :: dRPA = .true.
logical :: print_W = .true.
double precision,allocatable :: error_diis(:,:)
double precision,allocatable :: F_diis(:,:)
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:)
double precision,allocatable :: Aph(:,:)
double precision,allocatable :: Bph(:,:)
double precision,allocatable :: Om(:)
double precision,allocatable :: XpY(:,:)
double precision,allocatable :: XmY(:,:)
double precision,allocatable :: rho(:,:,:)
double precision,allocatable :: c(:,:)
double precision,allocatable :: cp(:,:)
double precision,allocatable :: eGW(:)
@ -112,13 +115,6 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
nBasSq = nBas*nBas
! COHSEX approximation
if(COHSEX) then
write(*,*) 'COHSEX approximation activated!'
write(*,*)
end if
! TDA for W
if(TDA_W) then
@ -137,7 +133,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
allocate(eGW(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),SigCm(nBas,nBas),Z(nBas), &
OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), &
Aph(nS,nS),Bph(nS,nS),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS), &
error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
! Initialization
@ -178,21 +174,24 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
! Compute linear response
call phLR(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
if(print_W) call print_excitation('RPA@qsGW ',ispin,nS,OmRPA)
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph)
if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
if(print_W) call print_excitation('RPA@qsGW ',ispin,nS,Om)
! Compute correlation part of the self-energy
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA)
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY,rho)
if(regularize) then
call regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC)
call regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z)
call regularized_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC)
call regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,Z)
else
call GW_self_energy(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC,Z)
call GW_self_energy(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
endif
@ -291,13 +290,13 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
! Deallocate memory
deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,error,error_diis,F_diis)
deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,Om,XpY,XmY,rho,error,error_diis,F_diis)
! Perform BSE calculation
if(BSE) then
if(dophBSE) then
call GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGW,eGW,EcBSE)
call GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGW,eGW,EcBSE)
if(exchange_kernel) then
@ -331,7 +330,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
end if
call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC)
call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
@ -346,4 +345,19 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,
end if
if(doppBSE) then
call GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eHF,eGW,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@qsGW correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@qsGW correlation energy (triplet) =',3d0*EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@qsGW correlation energy =',EcBSE(1) + 3d0*EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@qsGW total energy =',ENuc + ERHF + EcBSE(1) + 3d0*EcBSE(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
end if
end subroutine

52
src/GW/regularized_renormalization_factor.f90
View File

@ -1,4 +1,4 @@
subroutine regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Z)
subroutine regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Z)
! Compute the regularized version of the GW renormalization factor
@ -7,7 +7,6 @@ subroutine regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,O
! Input variables
logical,intent(in) :: COHSEX
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
@ -41,45 +40,34 @@ subroutine regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,O
kappa = 1d0
! static COHSEX approximation
! Occupied part of the correlation self-energy
if(COHSEX) then
Z(:) = 1d0
return
else
! Occupied part of the correlation self-energy
do p=nC+1,nBas-nR
do i=nC+1,nO
do jb=1,nS
eps = e(p) - e(i) + Omega(jb)
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
dfk = - fk/eps + 4d0*kappa**2*exp(-2d0*eps**2/kappa**2)
Z(p) = Z(p) - 2d0*rho(p,i,jb)**2*dfk
end do
do p=nC+1,nBas-nR
do i=nC+1,nO
do jb=1,nS
eps = e(p) - e(i) + Omega(jb)
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
dfk = - fk/eps + 4d0*kappa**2*exp(-2d0*eps**2/kappa**2)
Z(p) = Z(p) - 2d0*rho(p,i,jb)**2*dfk
end do
end do
end do
! Virtual part of the correlation self-energy
! Virtual part of the correlation self-energy
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
do jb=1,nS
eps = e(p) - e(a) - Omega(jb)
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
dfk = - fk/eps + 4d0*kappa**2*exp(-2d0*eps**2/kappa**2)
Z(p) = Z(p) - 2d0*rho(p,a,jb)**2*dfk
end do
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
do jb=1,nS
eps = e(p) - e(a) - Omega(jb)
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
dfk = - fk/eps + 4d0*kappa**2*exp(-2d0*eps**2/kappa**2)
Z(p) = Z(p) - 2d0*rho(p,a,jb)**2*dfk
end do
end do
end if
end do
! Compute renormalization factor from derivative of SigC
Z(:) = 1d0/(1d0 - Z(:))
end subroutine regularized_renormalization_factor
end subroutine

98
src/GW/regularized_self_energy_correlation.f90
View File

@ -1,4 +1,4 @@
subroutine regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC)
subroutine regularized_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC)
! Compute correlation part of the regularized self-energy
@ -7,7 +7,6 @@ subroutine regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,
! Input variables
logical,intent(in) :: COHSEX
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
@ -44,88 +43,49 @@ subroutine regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,
kappa = 1d0
!-----------------------------!
! COHSEX static approximation !
!-----------------------------!
if(COHSEX) then
! COHSEX: SEX of the COHSEX correlation self-energy
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do i=nC+1,nO
do jb=1,nS
SigC(p,q) = SigC(p,q) + 4d0*rho(p,i,jb)*rho(q,i,jb)/Omega(jb)
end do
end do
end do
end do
! COHSEX: COH part of the COHSEX correlation self-energy
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do r=nC+1,nBas-nR
do jb=1,nS
SigC(p,q) = SigC(p,q) - 2d0*rho(p,r,jb)*rho(q,r,jb)/Omega(jb)
end do
end do
end do
end do
EcGM = 0d0
do i=nC+1,nO
EcGM = EcGM + 0.5d0*SigC(i,i)
end do
else
!----------------!
! GW self-energy !
!----------------!
! Occupied part of the correlation self-energy
! Occupied part of the correlation self-energy
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do i=nC+1,nO
do jb=1,nS
eps = e(p) - e(i) + Omega(jb)
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
SigC(p,q) = SigC(p,q) + 2d0*rho(p,i,jb)*rho(q,i,jb)*fk
end do
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do i=nC+1,nO
do jb=1,nS
eps = e(p) - e(i) + Omega(jb)
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
SigC(p,q) = SigC(p,q) + 2d0*rho(p,i,jb)*rho(q,i,jb)*fk
end do
end do
end do
end do
! Virtual part of the correlation self-energy
! Virtual part of the correlation self-energy
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do a=nO+1,nBas-nR
do jb=1,nS
eps = e(p) - e(a) - Omega(jb)
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
SigC(p,q) = SigC(p,q) + 2d0*rho(p,a,jb)*rho(q,a,jb)*fk
end do
end do
end do
end do
! GM correlation energy
EcGM = 0d0
do i=nC+1,nO
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do a=nO+1,nBas-nR
do jb=1,nS
eps = e(a) - e(i) + Omega(jb)
eps = e(p) - e(a) - Omega(jb)
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
EcGM = EcGM - 4d0*rho(a,i,jb)**2*fk
SigC(p,q) = SigC(p,q) + 2d0*rho(p,a,jb)*rho(q,a,jb)*fk
end do
end do
end do
end do
end if
! GM correlation energy
end subroutine regularized_self_energy_correlation
EcGM = 0d0
do i=nC+1,nO
do a=nO+1,nBas-nR
do jb=1,nS
eps = e(a) - e(i) + Omega(jb)
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
EcGM = EcGM - 4d0*rho(a,i,jb)**2*fk
end do
end do
end do
end subroutine

96
src/GW/regularized_self_energy_correlation_diag.f90
View File

@ -1,4 +1,4 @@
subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC)
subroutine regularized_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC)
! Compute diagonal of the correlation part of the regularized self-energy
@ -7,7 +7,6 @@ subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,
! Input variables
logical,intent(in) :: COHSEX
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
@ -33,78 +32,41 @@ subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,
SigC(:) = 0d0
!-----------------------------
! COHSEX static self-energy
!-----------------------------
if(COHSEX) then
! COHSEX: SEX part of the COHSEX correlation self-energy
do p=nC+1,nBas-nR
do i=nC+1,nO
do jb=1,nS
SigC(p) = SigC(p) + 4d0*rho(p,i,jb)**2/Omega(jb)
end do
end do
end do
! COHSEX: COH part of the COHSEX correlation self-energy
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do jb=1,nS
SigC(p) = SigC(p) - 2d0*rho(p,q,jb)**2/Omega(jb)
end do
end do
end do
! GM correlation energy
EcGM = 0d0
do i=nC+1,nO
EcGM = EcGM - SigC(i)
end do
!-----------------------------
! GW self-energy
!-----------------------------
else
! Occupied part of the correlation self-energy
! Occupied part of the correlation self-energy
do p=nC+1,nBas-nR
do i=nC+1,nO
do jb=1,nS
Dpijb = e(p) - e(i) + Omega(jb)
SigC(p) = SigC(p) + 2d0*rho(p,i,jb)**2*(1d0 - exp(-2d0*eta*Dpijb*Dpijb))/Dpijb
end do
end do
end do
! Virtual part of the correlation self-energy
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
do jb=1,nS
Dpajb = e(p) - e(a) - Omega(jb)
SigC(p) = SigC(p) + 2d0*rho(p,a,jb)**2*(1d0 - exp(-2d0*eta*Dpajb*Dpajb))/Dpajb
end do
end do
end do
! GM correlation energy
EcGM = 0d0
do p=nC+1,nBas-nR
do i=nC+1,nO
do a=nO+1,nBas-nR
do jb=1,nS
EcGM = EcGM - 4d0*rho(a,i,jb)**2
end do
do jb=1,nS
Dpijb = e(p) - e(i) + Omega(jb)
SigC(p) = SigC(p) + 2d0*rho(p,i,jb)**2*(1d0 - exp(-2d0*eta*Dpijb*Dpijb))/Dpijb
end do
end do
end do
end if
! Virtual part of the correlation self-energy
end subroutine regularized_self_energy_correlation_diag
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
do jb=1,nS
Dpajb = e(p) - e(a) - Omega(jb)
SigC(p) = SigC(p) + 2d0*rho(p,a,jb)**2*(1d0 - exp(-2d0*eta*Dpajb*Dpajb))/Dpajb
end do
end do
end do
! Galitskii-Migdal correlation energy
EcGM = 0d0
do i=nC+1,nO
do a=nO+1,nBas-nR
do jb=1,nS
EcGM = EcGM - 4d0*rho(a,i,jb)**2
end do
end do
end do
end subroutine

72
src/GW/renormalization_factor_so.f90
View File

@ -1,72 +0,0 @@
subroutine renormalization_factor_so(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Z)
! Compute renormalization factor for GW
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: COHSEX
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega(nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables
integer :: p,i,a,jb
double precision :: eps
! Output variables
double precision,intent(out) :: Z(nBas)
! Initialize
Z(:) = 0d0
! static COHSEX approximation
if(COHSEX) then
Z(:) = 1d0
return
else
! Occupied part of the correlation self-energy
do p=nC+1,nBas-nR
do i=nC+1,nO
do jb=1,nS
eps = e(p) - e(i) + Omega(jb)
Z(p) = Z(p) - 1d0*rho(p,i,jb)**2*(eps/(eps**2 + eta**2))**2
end do
end do
end do
! Virtual part of the correlation self-energy
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
do jb=1,nS
eps = e(p) - e(a) - Omega(jb)
Z(p) = Z(p) - 1d0*rho(p,a,jb)**2*(eps/(eps**2 + eta**2))**2
end do
end do
end do
end if
! Compute renormalization factor from derivative of SigC
Z(:) = 1d0/(1d0 - Z(:))
end subroutine renormalization_factor_so

111
src/GW/self_energy_correlation_diag_so.f90
View File

@ -1,111 +0,0 @@
subroutine self_energy_correlation_diag_so(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC)
! Compute diagonal of the correlation part of the self-energy
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: COHSEX
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega(nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables
integer :: i,a,p,q,jb
double precision :: eps
! Output variables
double precision,intent(out) :: SigC(nBas)
double precision,intent(out) :: EcGM
! Initialize
SigC(:) = 0d0
!-----------------------------
! COHSEX static self-energy
!-----------------------------
if(COHSEX) then
! COHSEX: SEX part of the COHSEX correlation self-energy
do p=nC+1,nBas-nR
do i=nC+1,nO
do jb=1,nS
SigC(p) = SigC(p) + 2d0*rho(p,i,jb)**2/Omega(jb)
end do
end do
end do
! COHSEX: COH part of the COHSEX correlation self-energy
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do jb=1,nS
SigC(p) = SigC(p) - rho(p,q,jb)**2/Omega(jb)
end do
end do
end do
! GM correlation energy
EcGM = 0d0
do i=nC+1,nO
EcGM = EcGM - SigC(i)
end do
!-----------------------------
! GW self-energy
!-----------------------------
else
! Occupied part of the correlation self-energy
do p=nC+1,nBas-nR
do i=nC+1,nO
do jb=1,nS
eps = e(p) - e(i) + Omega(jb)
SigC(p) = SigC(p) + rho(p,i,jb)**2*eps/(eps**2 + eta**2)
end do
end do
end do
! Virtual part of the correlation self-energy
do p=nC+1,nBas-nR
do a=nO+1,nBas-nR
do jb=1,nS
eps = e(p) - e(a) - Omega(jb)
SigC(p) = SigC(p) + rho(p,a,jb)**2*eps/(eps**2 + eta**2)
end do
end do
end do
! GM correlation energy
EcGM = 0d0
do i=nC+1,nO
do a=nO+1,nBas-nR
do jb=1,nS
eps = e(a) - e(i) + Omega(jb)
EcGM = EcGM - 2d0*rho(a,i,jb)*rho(a,i,jb)*eps/(eps**2 + eta**2)
end do
end do
end do
end if
end subroutine self_energy_correlation_diag_so

211
src/GW/soG0W0.f90
View File

@ -1,211 +0,0 @@
subroutine soG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE, &
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, &
ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eGW)
! Perform G0W0 calculation in the spin-orbital basis
implicit none
include 'parameters.h'
include 'quadrature.h'
! Input variables
logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel
logical,intent(in) :: doXBS
logical,intent(in) :: COHSEX
logical,intent(in) :: BSE
logical,intent(in) :: ppBSE
logical,intent(in) :: TDA_W
logical,intent(in) :: TDA
logical,intent(in) :: dBSE
logical,intent(in) :: dTDA
logical,intent(in) :: evDyn
logical,intent(in) :: singlet
logical,intent(in) :: triplet
logical,intent(in) :: linearize
double precision,intent(in) :: eta
logical,intent(in) :: regularize
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: Vxc(nBas)
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: PHF(nBas,nBas)
! Local variables
logical :: print_W = .true.
integer :: ispin
double precision :: EcRPA
double precision :: EcBSE(nspin)
double precision :: EcAC(nspin)
double precision :: EcppBSE(nspin)
double precision :: EcGM
double precision,allocatable :: SigC(:)
double precision,allocatable :: Z(:)
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:)
double precision,allocatable :: eGWlin(:)
integer :: nBas2
integer :: nC2
integer :: nO2
integer :: nV2
integer :: nR2
integer :: nS2
double precision,allocatable :: seHF(:),seGW(:),sERI(:,:,:,:)
! Output variables
double precision :: eGW(nBas)
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| One-shot soG0W0 calculation |'
write(*,*)'************************************************'
write(*,*)
! Initialization
EcRPA = 0d0
! COHSEX approximation
if(COHSEX) then
write(*,*) 'COHSEX approximation activated!'
write(*,*)
end if
! TDA for W
if(TDA_W) then
write(*,*) 'Tamm-Dancoff approximation for dynamic screening!'
write(*,*)
end if
! TDA
if(TDA) then
write(*,*) 'Tamm-Dancoff approximation activated!'
write(*,*)
end if
! spatial to spin transformation
nBas2 = 2*nBas
nO2 = 2*nO
nV2 = 2*nV
nC2 = 2*nC
nR2 = 2*nR
nS2 = nO2*nV2
allocate(seHF(nBas2),seGW(nBas2),sERI(nBas2,nBas2,nBas2,nBas2))
call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF)
call spatial_to_spin_MO_energy(nBas,eGW,nBas2,seGW)
call spatial_to_spin_ERI(nBas,ERI_MO,nBas2,sERI)
! Spin manifold
ispin = 3
! Memory allocation
allocate(SigC(nBas2),Z(nBas2),OmRPA(nS2),XpY_RPA(nS2,nS2),XmY_RPA(nS2,nS2),rho_RPA(nBas2,nBas2,nS2),eGWlin(nBas2))
!-------------------!
! Compute screening !
!-------------------!
call phLR(ispin,.true.,TDA_W,eta,nBas2,nC2,nO2,nV2,nR2,nS2,1d0,seHF,sERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
if(print_W) call print_excitation('RPA@HF ',ispin,nS2,OmRPA)
!--------------------------!
! Compute spectral weights !
!--------------------------!
call GW_excitation_density(nBas2,nC2,nO2,nR2,nS2,sERI,XpY_RPA,rho_RPA)
!------------------------!
! Compute GW self-energy !
!------------------------!
if(regularize) then
call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
call regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
else
call self_energy_correlation_diag_so(COHSEX,eta,nBas2,nC2,nO2,nV2,nR2,nS2,seHF,OmRPA,rho_RPA,EcGM,SigC)
call renormalization_factor_so(COHSEX,eta,nBas2,nC2,nO2,nV2,nR2,nS2,seHF,OmRPA,rho_RPA,Z)
end if
!-----------------------------------!
! Solve the quasi-particle equation !
!-----------------------------------!
eGWlin(:) = seHF(:) + Z(:)*SigC(:)
! Linearized or graphical solution?
if(linearize) then
write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
write(*,*)
seGW(:) = eGWlin(:)
end if
! Compute the RPA correlation energy
call phLR(ispin,.true.,TDA_W,eta,nBas2,nC2,nO2,nV2,nR2,nS2,1d0,seGW,sERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
!--------------!
! Dump results !
!--------------!
call print_G0W0(nBas2,nO2,seHF,ENuc,ERHF,SigC,Z,seGW,EcRPA,EcGM)
! Deallocate memory
deallocate(SigC,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,eGWlin)
! Perform BSE calculation
if(ppBSE) then
call Bethe_Salpeter_pp_so(TDA_W,TDA,singlet,triplet,eta,nBas2,nC2,nO2,nV2,nR2,nS2,sERI,dipole_int,seHF,seGW,EcppBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@G0W0 correlation energy (singlet) =',EcppBSE(1)
write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@G0W0 correlation energy (triplet) =',3d0*EcppBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@G0W0 correlation energy =',EcppBSE(1) + 3d0*EcppBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@G0W0 total energy =',ENuc + ERHF + EcppBSE(1) + 3d0*EcppBSE(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
end if
end subroutine

2
src/HF/MOM_overlap.f90
View File

@ -49,4 +49,4 @@ subroutine MOM_overlap(nBas,nO,S,cG,c,ON)
end subroutine MOM_overlap
end subroutine

2
src/HF/RHF.f90
View File

@ -204,4 +204,4 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc
call mo_fock_exchange_potential(nBas,c,P,ERI,Vx)
end subroutine RHF
end subroutine

2
src/HF/UHF.f90
View File

@ -258,4 +258,4 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
call mo_fock_exchange_potential(nBas,c(:,:,ispin),P(:,:,ispin),ERI,Vx(:,ispin))
end do
end subroutine UHF
end subroutine

2
src/HF/UHF_stability.f90
View File

@ -171,4 +171,4 @@ subroutine UHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb)
write(*,*)'-------------------------------------------------------------'
write(*,*)
end subroutine UHF_stability
end subroutine

2
src/HF/density.f90
View File

@ -32,4 +32,4 @@ subroutine density(nGrid,nBas,P,AO,rho)
enddo
enddo
end subroutine density
end subroutine

2
src/HF/density_matrix.f90
View File

@ -27,4 +27,4 @@ subroutine density_matrix(nBas,ON,c,P)
enddo
enddo
end subroutine density_matrix
end subroutine

2
src/HF/dipole_moment.f90
View File

@ -50,4 +50,4 @@ subroutine dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole)
enddo
end subroutine dipole_moment
end subroutine

2
src/HF/mix_guess.f90
View File

@ -45,4 +45,4 @@ subroutine mix_guess(nBas,nO,c)
c(:,nO(2),2) = HOMOb(:)
c(:,nO(2)+1,2) = LUMOb(:)
end subroutine mix_guess
end subroutine

2
src/HF/mo_fock_exchange_potential.f90
View File

@ -38,4 +38,4 @@ subroutine mo_fock_exchange_potential(nBas,c,P,ERI,Vx)
deallocate(Fx)
end subroutine mo_fock_exchange_potential
end subroutine

4
src/HF/print_RHF.f90
View File

@ -74,6 +74,4 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF,dipole)
call matout(nBas,1,eHF)
write(*,*)
end subroutine print_RHF
end subroutine

2
src/HF/print_UHF.f90
View File

@ -124,4 +124,4 @@ subroutine print_UHF(nBas,nO,Ov,e,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole)
call matout(nBas,1,e(:,2))
write(*,*)
end subroutine print_UHF
end subroutine

4
src/LR/print_excitation.f90
View File

@ -41,6 +41,4 @@ subroutine print_excitation(method,ispin,nS,Omega)
write(*,*)'-------------------------------------------------------------'
write(*,*)
end subroutine print_excitation
end subroutine

2
src/LR/print_transition_vectors_pp.f90
View File

@ -191,4 +191,4 @@ subroutine print_transition_vectors_pp(spin_allowed,nBas,nC,nO,nV,nR,nOO,nVV,dip
if(nOO > 0) write(*,'(A50,F10.6)') 'Thomas-Reiche-Kuhn sum rule for h-h sector = ',sum(os2(:))
write(*,*)
end subroutine print_transition_vectors_pp
end subroutine

2
src/QuAcK/read_methods.f90
View File

@ -178,4 +178,4 @@ subroutine read_methods(doRHF,doUHF,doRMOM,doUMOM,doKS, &
close(unit=1)
end subroutine read_methods
end subroutine

2
src/QuAcK/read_options.f90
View File

@ -243,4 +243,4 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
close(unit=1)
end subroutine read_options
end subroutine

223
src/RPA/ACFDT_Tmatrix.f90
View File

@ -1,223 +0,0 @@
subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, &
nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, &
Om2t,X2t,Y2t,rho1t,rho2t,ERI,eT,eGT,EcAC)
! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem for the T-matrix
implicit none
include 'parameters.h'
include 'quadrature.h'
! Input variables
logical,intent(in) :: doXBS
logical,intent(in) :: exchange_kernel
logical,intent(in) :: dRPA
logical,intent(in) :: TDA_T
logical,intent(in) :: TDA
logical,intent(in) :: BSE
logical,intent(in) :: singlet
logical,intent(in) :: triplet
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
integer,intent(in) :: nOOs
integer,intent(in) :: nOOt
integer,intent(in) :: nVVs
integer,intent(in) :: nVVt
double precision,intent(in) :: Om1s(nVVs)
double precision,intent(in) :: X1s(nVVs,nVVs)
double precision,intent(in) :: Y1s(nOOs,nVVs)
double precision,intent(in) :: Om2s(nOOs)
double precision,intent(in) :: X2s(nVVs,nOOs)
double precision,intent(in) :: Y2s(nOOs,nOOs)
double precision,intent(in) :: rho1s(nBas,nBas,nVVs)
double precision,intent(in) :: rho2s(nBas,nBas,nOOs)
double precision,intent(in) :: Om1t(nVVt)
double precision,intent(in) :: X1t(nVVt,nVVt)
double precision,intent(in) :: Y1t(nOOt,nVVt)
double precision,intent(in) :: Om2t(nOOt)
double precision,intent(in) :: X2t(nVVt,nOOt)
double precision,intent(in) :: Y2t(nOOt,nOOt)
double precision,intent(in) :: rho1t(nBas,nBas,nVVt)
double precision,intent(in) :: rho2t(nBas,nBas,nOOt)
double precision,intent(in) :: eT(nBas)
double precision,intent(in) :: eGT(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: ispin
integer :: isp_T
integer :: iblock
integer :: iAC
double precision :: lambda
double precision,allocatable :: Ec(:,:)
double precision :: EcRPA(nspin)
double precision,allocatable :: TAs(:,:)
double precision,allocatable :: TBs(:,:)
double precision,allocatable :: TAt(:,:)
double precision,allocatable :: TBt(:,:)
double precision,allocatable :: Om(:,:)
double precision,allocatable :: XpY(:,:,:)
double precision,allocatable :: XmY(:,:,:)
! Output variables
double precision,intent(out) :: EcAC(nspin)
! Memory allocation
allocate(TAs(nS,nS),TBs(nS,nS),TAt(nS,nS),TBt(nS,nS), &
Om(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin))
allocate(Ec(nAC,nspin))
! Antisymmetrized kernel version
if(exchange_kernel) then
write(*,*)
write(*,*) '*** Exchange kernel version ***'
write(*,*)
end if
EcAC(:) = 0d0
Ec(:,:) = 0d0
! Singlet manifold
if(singlet) then
ispin = 1
write(*,*) '--------------'
write(*,*) 'Singlet states'
write(*,*) '--------------'
write(*,*)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(K x P_lambda)'
write(*,*) '-----------------------------------------------------------------------------------'
do iAC=1,nAC
lambda = rAC(iAC)
if(doXBS) then
isp_T = 1
iblock = 3
call ppLR(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,eT,ERI,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T))
call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Om1s,rho1s,Om2s,rho2s,TAs)
if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Om1s,rho1s,Om2s,rho2s,TBs)
isp_T = 2
iblock = 4
call ppLR(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,eT,ERI,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T))
call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Om1t,rho1t,Om2t,rho2t,TAt)
if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Om1t,rho1t,Om2t,rho2t,TBt)
end if
call linear_response_BSE(ispin,.false.,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TAt+TAs,TBt+TBs, &
EcAC(ispin),Om(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)
end do
EcAC(ispin) = 0.5d0*dot_product(wAC,Ec(:,ispin))
if(exchange_kernel) EcAC(ispin) = 0.5d0*EcAC(ispin)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,*)
end if
! Triplet manifold
if(triplet) then
ispin = 2
write(*,*) '--------------'
write(*,*) 'Triplet states'
write(*,*) '--------------'
write(*,*)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(K x P_lambda)'
write(*,*) '-----------------------------------------------------------------------------------'
do iAC=1,nAC
lambda = rAC(iAC)
if(doXBS) then
isp_T = 1
iblock = 3
call ppLR(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,eT,ERI,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T))
call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Om1s,rho1s,Om2s,rho2s,TAs)
if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Om1s,rho1s,Om2s,rho2s,TBs)
isp_T = 2
iblock = 4
call ppLR(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,eT,ERI,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T))
call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Om1t,rho1t,Om2t,rho2t,TAt)
if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Om1t,rho1t,Om2t,rho2t,TBt)
end if
call linear_response_BSE(ispin,.false.,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TAt-TAs,TBt-TBs, &
EcAC(ispin),Om(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)
end do
EcAC(ispin) = 0.5d0*dot_product(wAC,Ec(:,ispin))
if(exchange_kernel) EcAC(ispin) = 1.5d0*EcAC(ispin)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin)
write(*,*) '-----------------------------------------------------------------------------------'
write(*,*)
end if
end subroutine

113
src/RPA/soRPAx.f90
View File

@ -1,113 +0,0 @@
subroutine soRPAx(TDA,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
! Perform random phase approximation calculation with exchange (aka TDHF) in the
! spinorbital basis
implicit none
include 'parameters.h'
include 'quadrature.h'
! Input variables
logical,intent(in) :: TDA
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: ispin
double precision,allocatable :: Omega(:)
double precision,allocatable :: XpY(:,:)
double precision,allocatable :: XmY(:,:)
double precision,allocatable :: X(:,:)
double precision,allocatable :: Y(:,:)
double precision,allocatable :: Xinv(:,:)
double precision,allocatable :: t(:,:,:,:)
double precision :: EcRPAx
integer ::i,a,j,b,k,c,ia,jb,kc
! Hello world
write(*,*)
write(*,*)'***********************************************************'
write(*,*)'| Random phase approximation calculation with exchange |'
write(*,*)'***********************************************************'
write(*,*)
! TDA
if(TDA) then
write(*,*) 'Tamm-Dancoff approximation activated!'
write(*,*) ' => RPAx + TDA = CIS '
write(*,*)
end if
! Initialization
EcRPAx = 0d0
! Memory allocation
allocate(Omega(nS),XpY(nS,nS),XmY(nS,nS))
ispin = 3
call pphLR(ispin,.false.,TDA,0d0,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,EcRPAx,Omega,XpY,XmY)
call print_excitation('soRPAx@HF ',ispin,nS,Omega)
EcRPAx = 0.5d0*EcRPAx
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@RPAx correlation energy =',EcRPAx
write(*,'(2X,A50,F20.10)') 'Tr@RPAx total energy =',ENuc + ERHF + EcRPAx
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! allocate(X(nS,nS),Y(nS,nS),Xinv(nS,nS),t(nO,nO,nV,nV))
! X(:,:) = transpose(0.5d0*(XpY(:,:) + XmY(:,:)))
! Y(:,:) = transpose(0.5d0*(XpY(:,:) - XmY(:,:)))
! call matout(nS,nS,matmul(transpose(X),X)-matmul(transpose(Y),Y))
! call inverse_matrix(nS,X,Xinv)
! t = 0d0
! ia = 0
! do i=1,nO
! do a=1,nV
! ia = ia + 1
! jb = 0
! do j=1,nO
! do b=1,nV
! jb = jb + 1
! kc = 0
! do k=1,nO
! do c=1,nV
! kc = kc + 1
! t(i,j,a,b) = t(i,j,a,b) + Y(ia,kc)*Xinv(kc,jb)
! end do
! end do
! end do
! end do
! end do
! end do
! call matout(nO*nO,nV*nV,t)
end subroutine soRPAx

2
src/RPA/sort_ppRPA.f90
View File

@ -234,4 +234,4 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
! call matout(nVV,nOO,X2)
! call matout(nOO,nOO,Y2)
end subroutine sort_ppRPA
end subroutine

2
src/utils/DIIS_extrapolation.f90
View File

@ -52,4 +52,4 @@ subroutine DIIS_extrapolation(rcond,n_err,n_e,n_diis,error,e,error_in,e_inout)
e_inout(:) = matmul(w(1:n_diis),transpose(e(:,1:n_diis)))
end subroutine DIIS_extrapolation
end subroutine

2
src/utils/chem_to_phys_ERI.f90
View File

@ -30,4 +30,4 @@ subroutine chem_to_phys_ERI(nBas,ERI)
ERI(:,:,:,:) = pERI(:,:,:,:)
end subroutine chem_to_phys_ERI
end subroutine

2
src/utils/level_shifting.f90
View File

@ -34,4 +34,4 @@ subroutine level_shifting(level_shift,nBas,nO,S,c,F)
Sc(:,:) = matmul(S,c)
F(:,:) = matmul(Sc,matmul(F_MO,transpose(Sc)))
end subroutine level_shifting
end subroutine

128
src/utils/read_integrals_sph.f90
View File

@ -1,128 +0,0 @@
subroutine read_integrals_sph(nEl,nBas,S,T,V,Hc,G)
! Read one- and two-electron integrals from files
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
! Local variables
logical :: debug
integer :: nEl(nspin)
integer :: mu,nu,la,si
double precision :: Ov,Kin,Nuc,ERI
double precision :: rs,R,Rinv
! Output variables
double precision,intent(out) :: S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),G(nBas,nBas,nBas,nBas)
! Open file with integrals
debug = .false.
open(unit=1,file='input/sph')
read(1,*)
read(1,*) rs
R = sqrt(dble(sum(nEl(:))))/2d0*rs
Rinv = 1d0/R
print*, 'Scaling integrals by ',R
open(unit=8 ,file='/Users/loos/Integrals/QuAcK_Sph/Ov.dat')
open(unit=9 ,file='/Users/loos/Integrals/QuAcK_Sph/Kin.dat')
open(unit=10,file='/Users/loos/Integrals/QuAcK_Sph/Nuc.dat')
open(unit=11,file='/Users/loos/Integrals/QuAcK_Sph/ERI.dat')
! Read overlap integrals
S(:,:) = 0d0
do
read(8,*,end=8) mu,nu,Ov
S(mu,nu) = Ov
enddo
8 close(unit=8)
! Read kinetic integrals
T(:,:) = 0d0
do
read(9,*,end=9) mu,nu,Kin
T(mu,nu) = Rinv**2*Kin
enddo
9 close(unit=9)
! Read nuclear integrals
V(:,:) = 0d0
do
read(10,*,end=10) mu,nu,Nuc
V(mu,nu) = Nuc
enddo
10 close(unit=10)
! Define core Hamiltonian
Hc(:,:) = T(:,:) + V(:,:)
! Read nuclear integrals
G(:,:,:,:) = 0d0
do
read(11,*,end=11) mu,nu,la,si,ERI
ERI = Rinv*ERI
! <12|34>
G(mu,nu,la,si) = ERI
! <32|14>
G(la,nu,mu,si) = ERI
! <14|32>
G(mu,si,la,nu) = ERI
! <34|12>
G(la,si,mu,nu) = ERI
! <41|23>
G(si,mu,nu,la) = ERI
! <23|41>
G(nu,la,si,mu) = ERI
! <21|43>
G(nu,mu,si,la) = ERI
! <43|21>
G(si,la,nu,mu) = ERI
enddo
11 close(unit=11)
! Print results
if(debug) then
write(*,'(A28)') '----------------------'
write(*,'(A28)') 'Overlap integrals'
write(*,'(A28)') '----------------------'
call matout(nBas,nBas,S)
write(*,*)
write(*,'(A28)') '----------------------'
write(*,'(A28)') 'Kinetic integrals'
write(*,'(A28)') '----------------------'
call matout(nBas,nBas,T)
write(*,*)
write(*,'(A28)') '----------------------'
write(*,'(A28)') 'Nuclear integrals'
write(*,'(A28)') '----------------------'
call matout(nBas,nBas,V)
write(*,*)
write(*,'(A28)') '----------------------'
write(*,'(A28)') 'Electron repulsion integrals'
write(*,'(A28)') '----------------------'
do la=1,nBas
do si=1,nBas
call matout(nBas,nBas,G(1,1,la,si))
enddo
enddo
write(*,*)
endif
end subroutine read_integrals_sph

24
src/utils/rij.f90
View File

@ -1,24 +0,0 @@
subroutine rij(nWalk,r,r12)
! Compute the interelectronic distances
implicit none
! Input variables
integer,intent(in) :: nWalk
double precision,intent(in) :: r(nWalk,1:2,1:3)
! Output variables
double precision,intent(out) :: r12(nWalk)
! Compute
r12(1:nWalk) = (r(1:nWalk,1,1)-r(1:nWalk,2,1))**2 &
+ (r(1:nWalk,1,2)-r(1:nWalk,2,2))**2 &
+ (r(1:nWalk,1,3)-r(1:nWalk,2,3))**2
r12 = sqrt(r12)
end subroutine rij

6
src/utils/sort.f90
View File

@ -8,7 +8,7 @@
call rec_quicksort(x,iorder,isize,1,isize,1)
end subroutine quick_sort
end subroutine
recursive subroutine rec_quicksort(x,iorder,isize,first,last,level)
@ -58,7 +58,7 @@ recursive subroutine rec_quicksort(x,iorder,isize,first,last,level)
call rec_quicksort(x, iorder, isize, j+1, last,level/2)
endif
endif
end subroutine rec_quicksort
end subroutine
subroutine set_order(x,iorder,isize,jsize)
@ -86,4 +86,4 @@ subroutine set_order(x,iorder,isize,jsize)
deallocate(xtmp)
end subroutine set_order
end subroutine

40
src/utils/utils.f90
View File

@ -19,7 +19,7 @@ function Kronecker_delta(i,j) result(delta)
delta = 0d0
endif
end function Kronecker_delta
end function
function KroneckerDelta(i,j) result(delta)
@ -42,7 +42,7 @@ function KroneckerDelta(i,j) result(delta)
delta = 0
endif
end function KroneckerDelta
end function
!------------------------------------------------------------------------
subroutine matout(m,n,A)
@ -73,7 +73,7 @@ subroutine matout(m,n,A)
enddo
enddo
end subroutine matout
end subroutine
!------------------------------------------------------------------------
subroutine trace_vector(n,v,Tr)
@ -101,7 +101,7 @@ subroutine trace_vector(n,v,Tr)
Tr = Tr + v(i)
enddo
end subroutine trace_vector
end subroutine
!------------------------------------------------------------------------
function trace_matrix(n,A) result(Tr)
@ -128,7 +128,7 @@ function trace_matrix(n,A) result(Tr)
Tr = Tr + A(i,i)
enddo
end function trace_matrix
end function
!------------------------------------------------------------------------
subroutine compute_error(nData,Mean,Var,Error)
@ -148,7 +148,7 @@ subroutine compute_error(nData,Mean,Var,Error)
Error = sqrt((Var-Mean**2/nData)/nData/(nData-1d0))
end subroutine compute_error
end subroutine
!------------------------------------------------------------------------
subroutine identity_matrix(N,A)
@ -175,7 +175,7 @@ subroutine identity_matrix(N,A)
A(i,i) = 1d0
enddo
end subroutine identity_matrix
end subroutine
!------------------------------------------------------------------------
subroutine prepend(N,M,A,b)
@ -208,7 +208,7 @@ subroutine prepend(N,M,A,b)
A(i,1) = b(i)
enddo
end subroutine prepend
end subroutine
!------------------------------------------------------------------------
subroutine append(N,M,A,b)
@ -237,7 +237,7 @@ subroutine append(N,M,A,b)
A(i,M) = b(i)
enddo
end subroutine append
end subroutine
!------------------------------------------------------------------------
subroutine AtDA(N,A,D,B)
@ -270,7 +270,7 @@ subroutine AtDA(N,A,D,B)
enddo
enddo
end subroutine AtDA
end subroutine
!------------------------------------------------------------------------
subroutine ADAt(N,A,D,B)
@ -303,7 +303,7 @@ subroutine ADAt(N,A,D,B)
enddo
enddo
end subroutine ADAt
end subroutine
!------------------------------------------------------------------------
subroutine DA(N,D,A)
@ -323,7 +323,7 @@ subroutine DA(N,D,A)
enddo
enddo
end subroutine DA
end subroutine
!------------------------------------------------------------------------
subroutine AD(N,A,D)
@ -344,7 +344,7 @@ subroutine AD(N,A,D)
enddo
enddo
end subroutine AD
end subroutine
!------------------------------------------------------------------------
subroutine print_warning(message)
@ -357,7 +357,7 @@ subroutine print_warning(message)
write(*,*) message
end subroutine print_warning
end subroutine
!------------------------------------------------------------------------
@ -387,7 +387,7 @@ subroutine CalcTrAB(n,A,B,Tr)
enddo
enddo
end subroutine CalcTrAB
end subroutine
!------------------------------------------------------------------------
@ -408,7 +408,7 @@ function EpsilonSwitch(i,j) result(delta)
delta = -1
endif
end function EpsilonSwitch
end function
!------------------------------------------------------------------------
@ -429,7 +429,7 @@ function KappaCross(i,j,k) result(kappa)
kappa = dble(EpsilonSwitch(i,j)*KroneckerDelta(i,k) - EpsilonSwitch(k,i)*KroneckerDelta(i,j))
end function KappaCross
end function
!------------------------------------------------------------------------
@ -472,7 +472,7 @@ subroutine CalcInv3(a,det)
enddo
enddo
end subroutine CalcInv3
end subroutine
!------------------------------------------------------------------------
@ -530,7 +530,7 @@ subroutine CalcInv4(a,det)
enddo
enddo
end subroutine CalcInv4
end subroutine
subroutine wall_time(t)
implicit none
@ -542,5 +542,5 @@ subroutine wall_time(t)
endif
CALL SYSTEM_CLOCK(count=c)
t = dble(c)/dble(rate)
end subroutine wall_time
end subroutine

10
src/utils/wrap_lapack.f90
View File

@ -35,7 +35,7 @@ subroutine diagonalize_general_matrix(N,A,WR,VR)
print*,'Problem in diagonalize_general_matrix (dgeev)!!'
endif
end subroutine diagonalize_general_matrix
end subroutine
subroutine diagonalize_matrix(N,A,e)
@ -72,7 +72,7 @@ subroutine diagonalize_matrix(N,A,e)
print*,'Problem in diagonalize_matrix (dsyev)!!'
endif
end subroutine diagonalize_matrix
end subroutine
subroutine svd(N,A,U,D,Vt)
@ -157,7 +157,7 @@ subroutine inverse_matrix(N,A,B)
deallocate(ipiv,work)
end subroutine inverse_matrix
end subroutine
subroutine linear_solve(N,A,b,x,rcond)
@ -187,7 +187,7 @@ subroutine linear_solve(N,A,b,x,rcond)
! endif
end subroutine linear_solve
end subroutine
subroutine easy_linear_solve(N,A,b,x)
@ -218,5 +218,5 @@ subroutine easy_linear_solve(N,A,b,x)
endif
end subroutine easy_linear_solve
end subroutine