mirror of
https://github.com/pfloos/quack
synced 2024-12-22 12:23:50 +01:00
evGT
This commit is contained in:
parent
cb67629aed
commit
90a3570323
@ -15,7 +15,7 @@
|
||||
# G0W0* evGW* qsGW*
|
||||
F F F
|
||||
# G0T0 evGT qsGT
|
||||
T F F
|
||||
T T T
|
||||
# MCMP2
|
||||
F
|
||||
# * unrestricted version available
|
||||
|
||||
@ -13,6 +13,6 @@
|
||||
# ACFDT: AC Kx XBS
|
||||
F F T
|
||||
# BSE: BSE dBSE dTDA evDyn
|
||||
T T T F
|
||||
F T T F
|
||||
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
|
||||
1000000 100000 10 0.3 10000 1234 T
|
||||
|
||||
@ -56,8 +56,6 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,
|
||||
|
||||
integer :: ispin
|
||||
integer :: iblock
|
||||
integer :: dERI
|
||||
integer :: xERI
|
||||
|
||||
double precision :: EcRPA(nspin)
|
||||
double precision,allocatable :: TA(:,:),TB(:,:)
|
||||
@ -84,13 +82,11 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,
|
||||
|
||||
ispin = 1
|
||||
iblock = 3
|
||||
dERI = +1d0
|
||||
xERI = +0d0
|
||||
|
||||
call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eT,ERI, &
|
||||
Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin))
|
||||
|
||||
! call excitation_density_Tmatrix(iblock,dERI,xERI,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
|
||||
! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
|
||||
|
||||
call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,ERI,Omega1s,rho1s,Omega2s,rho2s,TA)
|
||||
if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,ERI,Omega1s,rho1s,Omega2s,rho2s,TB)
|
||||
@ -101,13 +97,11 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,
|
||||
|
||||
ispin = 2
|
||||
iblock = 4
|
||||
dERI = +1d0
|
||||
xERI = -1d0
|
||||
|
||||
call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eT,ERI, &
|
||||
Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin))
|
||||
|
||||
! call excitation_density_Tmatrix(iblock,dERI,xERI,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
|
||||
! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
|
||||
|
||||
call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,ERI,Omega1t,rho1t,Omega2t,rho2t,TA)
|
||||
if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,ERI,Omega1t,rho1t,Omega2t,rho2t,TB)
|
||||
|
||||
@ -44,9 +44,6 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
|
||||
integer :: iblock
|
||||
integer :: nOOs,nOOt
|
||||
integer :: nVVs,nVVt
|
||||
double precision :: dERI
|
||||
double precision :: xERI
|
||||
double precision :: alpha
|
||||
double precision :: EcRPA(nspin)
|
||||
double precision :: EcBSE(nspin)
|
||||
double precision :: EcAC(nspin)
|
||||
@ -140,20 +137,16 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
|
||||
Z(:) = 0d0
|
||||
|
||||
iblock = 3
|
||||
dERI = +1d0
|
||||
xERI = +0d0
|
||||
|
||||
call excitation_density_Tmatrix(iblock,dERI,xERI,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
|
||||
call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
|
||||
|
||||
call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,Omega1s,rho1s,Omega2s,rho2s,SigT)
|
||||
|
||||
call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,Omega1s,rho1s,Omega2s,rho2s,Z)
|
||||
|
||||
iblock = 4
|
||||
dERI = +1d0
|
||||
xERI = -1d0
|
||||
|
||||
call excitation_density_Tmatrix(iblock,dERI,xERI,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
|
||||
call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
|
||||
|
||||
call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF,Omega1t,rho1t,Omega2t,rho2t,SigT)
|
||||
|
||||
|
||||
@ -54,9 +54,6 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
|
||||
integer :: iblock
|
||||
integer :: nOOs,nOOt
|
||||
integer :: nVVs,nVVt
|
||||
double precision :: dERI
|
||||
double precision :: xERI
|
||||
double precision :: alpha
|
||||
double precision :: EcRPA(nspin)
|
||||
double precision :: EcBSE(nspin)
|
||||
double precision :: EcAC(nspin)
|
||||
@ -76,11 +73,6 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
|
||||
double precision,allocatable :: SigT(:)
|
||||
double precision,allocatable :: Z(:)
|
||||
|
||||
double precision,allocatable :: Omega(:,:)
|
||||
double precision,allocatable :: XpY(:,:,:)
|
||||
double precision,allocatable :: XmY(:,:,:)
|
||||
double precision,allocatable :: rho(:,:,:,:)
|
||||
|
||||
! Output variables
|
||||
|
||||
! Hello world
|
||||
@ -140,14 +132,9 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
|
||||
|
||||
! Compute linear response
|
||||
|
||||
call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT(:),ERI_MO(:,:,:,:), &
|
||||
Omega1s(:),X1s(:,:),Y1s(:,:),Omega2s(:),X2s(:,:),Y2s(:,:),EcRPA(ispin))
|
||||
call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,ERI_MO, &
|
||||
Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin))
|
||||
|
||||
! EcRPA(ispin) = 1d0*EcRPA(ispin)
|
||||
|
||||
! call print_excitation('pp-RPA (N+2)',iblock,nVVs,Omega1s(:))
|
||||
! call print_excitation('pp-RPA (N-2)',iblock,nOOs,Omega2s(:))
|
||||
|
||||
!----------------------------------------------
|
||||
! alpha-alpha block
|
||||
!----------------------------------------------
|
||||
@ -157,14 +144,8 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
|
||||
|
||||
! Compute linear response
|
||||
|
||||
call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT(:),ERI_MO(:,:,:,:), &
|
||||
Omega1t(:),X1t(:,:),Y1t(:,:),Omega2t(:),X2t(:,:),Y2t(:,:),EcRPA(ispin))
|
||||
|
||||
! EcRPA(ispin) = 2d0*EcRPA(ispin)
|
||||
! EcRPA(ispin) = 3d0*EcRPA(ispin)
|
||||
|
||||
! call print_excitation('pp-RPA (N+2)',iblock,nVVt,Omega1t(:))
|
||||
! call print_excitation('pp-RPA (N-2)',iblock,nOOt,Omega2t(:))
|
||||
call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,ERI_MO, &
|
||||
Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin))
|
||||
|
||||
!----------------------------------------------
|
||||
! Compute T-matrix version of the self-energy
|
||||
@ -174,32 +155,26 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
|
||||
Z(:) = 0d0
|
||||
|
||||
iblock = 3
|
||||
dERI = +1d0
|
||||
xERI = +0d0
|
||||
alpha = +1d0
|
||||
|
||||
call excitation_density_Tmatrix(iblock,dERI,xERI,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO(:,:,:,:), &
|
||||
X1s(:,:),Y1s(:,:),rho1s(:,:,:),X2s(:,:),Y2s(:,:),rho2s(:,:,:))
|
||||
call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO, &
|
||||
X1s,Y1s,rho1s,X2s,Y2s,rho2s)
|
||||
|
||||
call self_energy_Tmatrix_diag(alpha,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT(:), &
|
||||
Omega1s(:),rho1s(:,:,:),Omega2s(:),rho2s(:,:,:),SigT(:))
|
||||
call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, &
|
||||
Omega1s,rho1s,Omega2s,rho2s,SigT)
|
||||
|
||||
call renormalization_factor_Tmatrix(alpha,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT(:), &
|
||||
Omega1s(:),rho1s(:,:,:),Omega2s(:),rho2s(:,:,:),Z(:))
|
||||
call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, &
|
||||
Omega1s,rho1s,Omega2s,rho2s,Z)
|
||||
|
||||
iblock = 4
|
||||
dERI = +1d0
|
||||
xERI = -1d0
|
||||
alpha = +1d0
|
||||
|
||||
call excitation_density_Tmatrix(iblock,dERI,xERI,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO(:,:,:,:), &
|
||||
X1t(:,:),Y1t(:,:),rho1t(:,:,:),X2t(:,:),Y2t(:,:),rho2t(:,:,:))
|
||||
call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO, &
|
||||
X1t,Y1t,rho1t,X2t,Y2t,rho2t)
|
||||
|
||||
call self_energy_Tmatrix_diag(alpha,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT(:), &
|
||||
Omega1t(:),rho1t(:,:,:),Omega2t(:),rho2t(:,:,:),SigT(:))
|
||||
call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, &
|
||||
Omega1t,rho1t,Omega2t,rho2t,SigT)
|
||||
|
||||
call renormalization_factor_Tmatrix(alpha,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT(:), &
|
||||
Omega1t(:),rho1t(:,:,:),Omega2t(:),rho2t(:,:,:),Z(:))
|
||||
call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, &
|
||||
Omega1t,rho1t,Omega2t,rho2t,Z)
|
||||
|
||||
Z(:) = 1d0/(1d0 - Z(:))
|
||||
|
||||
@ -219,7 +194,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
|
||||
! Dump results
|
||||
!----------------------------------------------
|
||||
|
||||
call print_evGT(nBas,nO,nSCF,Conv,eHF(:),SigT(:),Z(:),eGT(:))
|
||||
call print_evGT(nBas,nO,nSCF,Conv,eHF,SigT,Z,eGT)
|
||||
|
||||
! DIIS extrapolation
|
||||
|
||||
@ -247,12 +222,12 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
|
||||
|
||||
ispin = 1
|
||||
iblock = 3
|
||||
call linear_response_pp(iblock,.false.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT(:),ERI_MO(:,:,:,:), &
|
||||
Omega1s(:),X1s(:,:),Y1s(:,:),Omega2s(:),X2s(:,:),Y2s(:,:),EcRPA(ispin))
|
||||
call linear_response_pp(iblock,.false.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,ERI_MO, &
|
||||
Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin))
|
||||
ispin = 2
|
||||
iblock = 4
|
||||
call linear_response_pp(iblock,.false.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT(:),ERI_MO(:,:,:,:), &
|
||||
Omega1t(:),X1t(:,:),Y1t(:,:),Omega2t(:),X2t(:,:),Y2t(:,:),EcRPA(ispin))
|
||||
call linear_response_pp(iblock,.false.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,ERI_MO, &
|
||||
Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin))
|
||||
EcRPA(1) = EcRPA(1) - EcRPA(2)
|
||||
EcRPA(2) = 3d0*EcRPA(2)
|
||||
|
||||
|
||||
@ -1,4 +1,4 @@
|
||||
subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2)
|
||||
subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2)
|
||||
|
||||
! Compute excitation densities for T-matrix self-energy
|
||||
|
||||
@ -7,8 +7,6 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: ispin
|
||||
double precision,intent(in) :: dERI
|
||||
double precision,intent(in) :: xERI
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
@ -55,8 +53,7 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E
|
||||
do c=nO+1,nBas-nR
|
||||
do d=c,nBas-nR
|
||||
cd = cd + 1
|
||||
rho1(p,q,ab) = rho1(p,q,ab) &
|
||||
+ (dERI*ERI(p,q,c,d) + xERI*ERI(p,q,d,c))*X1(cd,ab)
|
||||
rho1(p,q,ab) = rho1(p,q,ab) + ERI(p,q,c,d)*X1(cd,ab)
|
||||
end do
|
||||
end do
|
||||
|
||||
@ -64,8 +61,7 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E
|
||||
do k=nC+1,nO
|
||||
do l=k,nO
|
||||
kl = kl + 1
|
||||
rho1(p,q,ab) = rho1(p,q,ab) &
|
||||
+ (dERI*ERI(p,q,k,l) + xERI*ERI(p,q,l,k))*Y1(kl,ab)
|
||||
rho1(p,q,ab) = rho1(p,q,ab) + ERI(p,q,k,l)*Y1(kl,ab)
|
||||
end do
|
||||
end do
|
||||
|
||||
@ -77,8 +73,7 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E
|
||||
do c=nO+1,nBas-nR
|
||||
do d=c,nBas-nR
|
||||
cd = cd + 1
|
||||
rho2(p,q,ij) = rho2(p,q,ij) &
|
||||
+ (dERI*ERI(p,q,c,d) + xERI*ERI(p,q,d,c))*X2(cd,ij)
|
||||
rho2(p,q,ij) = rho2(p,q,ij) + ERI(p,q,c,d)*X2(cd,ij)
|
||||
end do
|
||||
end do
|
||||
|
||||
@ -86,8 +81,7 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E
|
||||
do k=nC+1,nO
|
||||
do l=k,nO
|
||||
kl = kl + 1
|
||||
rho2(p,q,ij) = rho2(p,q,ij) &
|
||||
+ (dERI*ERI(p,q,k,l) + xERI*ERI(p,q,l,k))*Y2(kl,ij)
|
||||
rho2(p,q,ij) = rho2(p,q,ij) + ERI(p,q,k,l)*Y2(kl,ij)
|
||||
end do
|
||||
end do
|
||||
|
||||
@ -114,7 +108,7 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E
|
||||
do d=c+1,nBas-nR
|
||||
cd = cd + 1
|
||||
rho1(p,q,ab) = rho1(p,q,ab) &
|
||||
+ (dERI*ERI(p,q,c,d) + xERI*ERI(p,q,d,c))*X1(cd,ab)
|
||||
+ (ERI(p,q,c,d) - ERI(p,q,d,c))*X1(cd,ab)
|
||||
end do
|
||||
end do
|
||||
|
||||
@ -123,7 +117,7 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E
|
||||
do l=k+1,nO
|
||||
kl = kl + 1
|
||||
rho1(p,q,ab) = rho1(p,q,ab) &
|
||||
+ (dERI*ERI(p,q,k,l) + xERI*ERI(p,q,l,k))*Y1(kl,ab)
|
||||
+ (ERI(p,q,k,l) - ERI(p,q,l,k))*Y1(kl,ab)
|
||||
end do
|
||||
end do
|
||||
|
||||
@ -136,7 +130,7 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E
|
||||
do d=c+1,nBas-nR
|
||||
cd = cd + 1
|
||||
rho2(p,q,ij) = rho2(p,q,ij) &
|
||||
+ (dERI*ERI(p,q,c,d) + xERI*ERI(p,q,d,c))*X2(cd,ij)
|
||||
+ (ERI(p,q,c,d) - ERI(p,q,d,c))*X2(cd,ij)
|
||||
end do
|
||||
end do
|
||||
|
||||
@ -145,7 +139,7 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E
|
||||
do l=k+1,nO
|
||||
kl = kl + 1
|
||||
rho2(p,q,ij) = rho2(p,q,ij) &
|
||||
+ (dERI*ERI(p,q,k,l) + xERI*ERI(p,q,l,k))*Y2(kl,ij)
|
||||
+ (ERI(p,q,k,l) - ERI(p,q,l,k))*Y2(kl,ij)
|
||||
end do
|
||||
end do
|
||||
|
||||
@ -171,8 +165,7 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E
|
||||
do c=nO+1,nBas-nR
|
||||
do d=nO+1,nBas-nR
|
||||
cd = cd + 1
|
||||
rho1(p,q,ab) = rho1(p,q,ab) &
|
||||
+ (dERI*ERI(p,q,c,d) + xERI*ERI(p,q,d,c))*X1(cd,ab)
|
||||
rho1(p,q,ab) = rho1(p,q,ab) + ERI(p,q,c,d)*X1(cd,ab)
|
||||
end do
|
||||
end do
|
||||
|
||||
@ -180,8 +173,7 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E
|
||||
do k=nC+1,nO
|
||||
do l=nC+1,nO
|
||||
kl = kl + 1
|
||||
rho1(p,q,ab) = rho1(p,q,ab) &
|
||||
+ (dERI*ERI(p,q,k,l) + xERI*ERI(p,q,l,k))*Y1(kl,ab)
|
||||
rho1(p,q,ab) = rho1(p,q,ab) + ERI(p,q,k,l)*Y1(kl,ab)
|
||||
end do
|
||||
end do
|
||||
|
||||
@ -193,8 +185,7 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E
|
||||
do c=nO+1,nBas-nR
|
||||
do d=nO+1,nBas-nR
|
||||
cd = cd + 1
|
||||
rho2(p,q,ij) = rho2(p,q,ij) &
|
||||
+ (dERI*ERI(p,q,c,d) + xERI*ERI(p,q,d,c))*X2(cd,ij)
|
||||
rho2(p,q,ij) = rho2(p,q,ij) + ERI(p,q,c,d)*X2(cd,ij)
|
||||
end do
|
||||
end do
|
||||
|
||||
@ -202,8 +193,7 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E
|
||||
do k=nC+1,nO
|
||||
do l=nC+1,nO
|
||||
kl = kl + 1
|
||||
rho2(p,q,ij) = rho2(p,q,ij) &
|
||||
+ (dERI*ERI(p,q,k,l) + xERI*ERI(p,q,l,k))*Y2(kl,ij)
|
||||
rho2(p,q,ij) = rho2(p,q,ij) + ERI(p,q,k,l)*Y2(kl,ij)
|
||||
end do
|
||||
end do
|
||||
|
||||
|
||||
120
src/MBPT/print_qsGT.f90
Normal file
120
src/MBPT/print_qsGT.f90
Normal file
@ -0,0 +1,120 @@
|
||||
subroutine print_qsGT(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole)
|
||||
|
||||
! Print one-electron energies and other stuff for qsGT
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nSCF
|
||||
double precision,intent(in) :: ENuc
|
||||
double precision,intent(in) :: ET
|
||||
double precision,intent(in) :: EV
|
||||
double precision,intent(in) :: EJ
|
||||
double precision,intent(in) :: Ex
|
||||
double precision,intent(in) :: EcGM
|
||||
double precision,intent(in) :: EcRPA(nspin)
|
||||
double precision,intent(in) :: Conv
|
||||
double precision,intent(in) :: thresh
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
double precision,intent(in) :: eGT(nBas)
|
||||
double precision,intent(in) :: c(nBas)
|
||||
double precision,intent(in) :: SigC(nBas,nBas)
|
||||
double precision,intent(in) :: Z(nBas)
|
||||
double precision,intent(in) :: dipole(ncart)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: x,ixyz,HOMO,LUMO
|
||||
double precision :: Gap
|
||||
double precision,external :: trace_matrix
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: EqsGT
|
||||
|
||||
! HOMO and LUMO
|
||||
|
||||
HOMO = nO
|
||||
LUMO = HOMO + 1
|
||||
Gap = eGT(LUMO)-eGT(HOMO)
|
||||
|
||||
! Compute energies
|
||||
|
||||
! Dump results
|
||||
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
if(nSCF < 10) then
|
||||
write(*,'(1X,A21,I1,A1,I1,A12)')' Self-consistent qsG',nSCF,'W',nSCF,' calculation'
|
||||
else
|
||||
write(*,'(1X,A21,I2,A1,I2,A12)')' Self-consistent qsG',nSCF,'W',nSCF,' calculation'
|
||||
endif
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
|
||||
'|','#','|','e_HF (eV)','|','Sig_T (eV)','|','Z','|','e_QP (eV)','|'
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
|
||||
do x=1,nBas
|
||||
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
|
||||
'|',x,'|',eHF(x)*HaToeV,'|',SigC(x,x)*HaToeV,'|',Z(x),'|',eGT(x)*HaToeV,'|'
|
||||
enddo
|
||||
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,'(2X,A10,I3)') 'Iteration ',nSCF
|
||||
write(*,'(2X,A14,F15.5)')'Convergence = ',Conv
|
||||
write(*,*)'-------------------------------------------'
|
||||
write(*,'(2X,A30,F15.6,A3)') 'qsGT HOMO energy:',eGT(HOMO)*HaToeV,' eV'
|
||||
write(*,'(2X,A30,F15.6,A3)') 'qsGT LUMO energy:',eGT(LUMO)*HaToeV,' eV'
|
||||
write(*,'(2X,A30,F15.6,A3)') 'qsGT HOMO-LUMO gap :',Gap*HaToeV,' eV'
|
||||
write(*,*)'-------------------------------------------'
|
||||
write(*,'(2X,A30,F15.6,A3)') ' qsGT total energy:',ENuc + EqsGT,' au'
|
||||
write(*,'(2X,A30,F15.6,A3)') ' qsGT exchange energy:',Ex,' au'
|
||||
write(*,'(2X,A30,F15.6,A3)') ' GM@qsGT correlation energy:',EcGM,' au'
|
||||
write(*,'(2X,A30,F15.6,A3)') 'RPA@qsGT correlation energy:',sum(EcRPA(:)),' au'
|
||||
write(*,*)'-------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
! Dump results for final iteration
|
||||
|
||||
if(Conv < thresh) then
|
||||
|
||||
write(*,*)
|
||||
write(*,'(A50)') '---------------------------------------'
|
||||
write(*,'(A32)') ' Summary '
|
||||
write(*,'(A50)') '---------------------------------------'
|
||||
write(*,'(A32,1X,F16.10,A3)') ' One-electron energy: ',ET + EV,' au'
|
||||
write(*,'(A32,1X,F16.10,A3)') ' Kinetic energy: ',ET,' au'
|
||||
write(*,'(A32,1X,F16.10,A3)') ' Potential energy: ',EV,' au'
|
||||
write(*,'(A50)') '---------------------------------------'
|
||||
write(*,'(A32,1X,F16.10,A3)') ' Two-electron energy: ',EJ + Ex,' au'
|
||||
write(*,'(A32,1X,F16.10,A3)') ' Hartree energy: ',EJ,' au'
|
||||
write(*,'(A32,1X,F16.10,A3)') ' Exchange energy: ',Ex,' au'
|
||||
write(*,'(A32,1X,F16.10,A3)') ' Correlation energy: ',EcGM,' au'
|
||||
write(*,'(A50)') '---------------------------------------'
|
||||
write(*,'(A32,1X,F16.10,A3)') ' Electronic energy: ',EqsGT,' au'
|
||||
write(*,'(A32,1X,F16.10,A3)') ' Nuclear repulsion: ',ENuc,' au'
|
||||
write(*,'(A32,1X,F16.10,A3)') ' qsGT energy: ',ENuc + EqsGT,' au'
|
||||
write(*,'(A50)') '---------------------------------------'
|
||||
write(*,'(A35)') ' Dipole moment (Debye) '
|
||||
write(*,'(10X,4A10)') 'X','Y','Z','Tot.'
|
||||
write(*,'(10X,4F10.6)') (dipole(ixyz)*auToD,ixyz=1,ncart),norm2(dipole)*auToD
|
||||
write(*,'(A50)') '-----------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
write(*,'(A50)') '---------------------------------------'
|
||||
write(*,'(A32)') ' qsGT MO coefficients'
|
||||
write(*,'(A50)') '---------------------------------------'
|
||||
call matout(nBas,nBas,c)
|
||||
write(*,*)
|
||||
write(*,'(A50)') '---------------------------------------'
|
||||
write(*,'(A32)') ' qsGT MO energies'
|
||||
write(*,'(A50)') '---------------------------------------'
|
||||
call matout(nBas,1,eGT)
|
||||
write(*,*)
|
||||
|
||||
endif
|
||||
|
||||
end subroutine print_qsGT
|
||||
399
src/MBPT/qsGT.f90
Normal file
399
src/MBPT/qsGT.f90
Normal file
@ -0,0 +1,399 @@
|
||||
subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA, &
|
||||
dBSE,dTDA,evDyn,singlet,triplet,eta,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)
|
||||
|
||||
! Perform a quasiparticle self-consistent GT calculation
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: maxSCF
|
||||
integer,intent(in) :: max_diis
|
||||
double precision,intent(in) :: thresh
|
||||
logical,intent(in) :: doACFDT
|
||||
logical,intent(in) :: exchange_kernel
|
||||
logical,intent(in) :: doXBS
|
||||
logical,intent(in) :: BSE
|
||||
logical,intent(in) :: TDA_T
|
||||
logical,intent(in) :: TDA
|
||||
logical,intent(in) :: dBSE
|
||||
logical,intent(in) :: dTDA
|
||||
logical,intent(in) :: evDyn
|
||||
logical,intent(in) :: singlet
|
||||
logical,intent(in) :: triplet
|
||||
double precision,intent(in) :: eta
|
||||
|
||||
integer,intent(in) :: nNuc
|
||||
double precision,intent(in) :: ZNuc(nNuc)
|
||||
double precision,intent(in) :: rNuc(nNuc,ncart)
|
||||
double precision,intent(in) :: ENuc
|
||||
|
||||
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
|
||||
double precision,intent(in) :: ERHF
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
double precision,intent(in) :: cHF(nBas,nBas)
|
||||
double precision,intent(in) :: PHF(nBas,nBas)
|
||||
double precision,intent(in) :: S(nBas,nBas)
|
||||
double precision,intent(in) :: T(nBas,nBas)
|
||||
double precision,intent(in) :: V(nBas,nBas)
|
||||
double precision,intent(in) :: Hc(nBas,nBas)
|
||||
double precision,intent(in) :: X(nBas,nBas)
|
||||
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(inout):: ERI_MO(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
|
||||
double precision,intent(in) :: dipole_int_MO(nBas,nBas,ncart)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: nSCF
|
||||
integer :: nBasSq
|
||||
integer :: ispin
|
||||
integer :: iblock
|
||||
integer :: n_diis
|
||||
double precision :: ET
|
||||
double precision :: EV
|
||||
double precision :: EJ
|
||||
double precision :: Ex
|
||||
double precision :: EqsGT
|
||||
double precision :: EcRPA(nspin)
|
||||
double precision :: EcBSE(nspin)
|
||||
double precision :: EcAC(nspin)
|
||||
double precision :: EcGM
|
||||
double precision :: Conv
|
||||
double precision :: rcond
|
||||
double precision,external :: trace_matrix
|
||||
double precision :: dipole(ncart)
|
||||
|
||||
integer :: nOOs,nOOt
|
||||
integer :: nVVs,nVVt
|
||||
|
||||
logical :: print_W = .false.
|
||||
double precision,allocatable :: error_diis(:,:)
|
||||
double precision,allocatable :: F_diis(:,:)
|
||||
double precision,allocatable :: c(:,:)
|
||||
double precision,allocatable :: cp(:,:)
|
||||
double precision,allocatable :: eGT(:)
|
||||
double precision,allocatable :: eOld(:)
|
||||
double precision,allocatable :: Omega1s(:),Omega1t(:)
|
||||
double precision,allocatable :: X1s(:,:),X1t(:,:)
|
||||
double precision,allocatable :: Y1s(:,:),Y1t(:,:)
|
||||
double precision,allocatable :: rho1s(:,:,:),rho1t(:,:,:)
|
||||
double precision,allocatable :: Omega2s(:),Omega2t(:)
|
||||
double precision,allocatable :: X2s(:,:),X2t(:,:)
|
||||
double precision,allocatable :: Y2s(:,:),Y2t(:,:)
|
||||
double precision,allocatable :: rho2s(:,:,:),rho2t(:,:,:)
|
||||
double precision,allocatable :: P(:,:)
|
||||
double precision,allocatable :: F(:,:)
|
||||
double precision,allocatable :: Fp(:,:)
|
||||
double precision,allocatable :: J(:,:)
|
||||
double precision,allocatable :: K(:,:)
|
||||
double precision,allocatable :: SigT(:,:)
|
||||
double precision,allocatable :: SigTp(:,:)
|
||||
double precision,allocatable :: SigTm(:,:)
|
||||
double precision,allocatable :: Z(:)
|
||||
double precision,allocatable :: error(:,:)
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'************************************************'
|
||||
write(*,*)'| Self-consistent qsGT calculation |'
|
||||
write(*,*)'************************************************'
|
||||
write(*,*)
|
||||
|
||||
! Dimensions of the pp-RPA linear reponse matrices
|
||||
|
||||
nOOs = nO*nO
|
||||
nVVs = nV*nV
|
||||
|
||||
nOOt = nO*(nO - 1)/2
|
||||
nVVt = nV*(nV - 1)/2
|
||||
|
||||
! Warning
|
||||
|
||||
write(*,*) '!! ERIs in MO basis will be overwritten in qsGT !!'
|
||||
write(*,*)
|
||||
|
||||
! Stuff
|
||||
|
||||
nBasSq = nBas*nBas
|
||||
|
||||
! TDA for T
|
||||
|
||||
if(TDA_T) then
|
||||
write(*,*) 'Tamm-Dancoff approximation for T-matrix!'
|
||||
write(*,*)
|
||||
end if
|
||||
|
||||
! TDA
|
||||
|
||||
if(TDA) then
|
||||
write(*,*) 'Tamm-Dancoff approximation activated!'
|
||||
write(*,*)
|
||||
end if
|
||||
|
||||
! Memory allocation
|
||||
|
||||
allocate(eGT(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
|
||||
J(nBas,nBas),K(nBas,nBas),SigT(nBas,nBas),SigTp(nBas,nBas),SigTm(nBas,nBas),Z(nBas), &
|
||||
error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
|
||||
|
||||
allocate(Omega1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), &
|
||||
Omega2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), &
|
||||
rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs), &
|
||||
Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), &
|
||||
Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), &
|
||||
rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt))
|
||||
|
||||
! Initialization
|
||||
|
||||
nSCF = -1
|
||||
n_diis = 0
|
||||
ispin = 1
|
||||
Conv = 1d0
|
||||
P(:,:) = PHF(:,:)
|
||||
eGT(:) = eHF(:)
|
||||
eOld(:) = eHF(:)
|
||||
c(:,:) = cHF(:,:)
|
||||
F_diis(:,:) = 0d0
|
||||
error_diis(:,:) = 0d0
|
||||
rcond = 1d0
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Main loop
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
do while(Conv > thresh .and. nSCF <= maxSCF)
|
||||
|
||||
! Increment
|
||||
|
||||
nSCF = nSCF + 1
|
||||
|
||||
! Buid Coulomb matrix
|
||||
|
||||
call Coulomb_matrix_AO_basis(nBas,P,ERI_AO,J)
|
||||
|
||||
! Compute exchange part of the self-energy
|
||||
|
||||
call exchange_matrix_AO_basis(nBas,P,ERI_AO,K)
|
||||
|
||||
! AO to MO transformation of two-electron integrals
|
||||
|
||||
call AOtoMO_integral_transform(1,1,1,1,nBas,c,ERI_AO,ERI_MO)
|
||||
|
||||
! Compute linear response
|
||||
|
||||
ispin = 1
|
||||
iblock = 3
|
||||
|
||||
call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,ERI_MO, &
|
||||
Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin))
|
||||
|
||||
ispin = 2
|
||||
iblock = 4
|
||||
|
||||
call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,ERI_MO, &
|
||||
Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin))
|
||||
|
||||
! Compute correlation part of the self-energy
|
||||
|
||||
SigT(:,:) = 0d0
|
||||
Z(:) = 0d0
|
||||
|
||||
iblock = 3
|
||||
|
||||
call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO, &
|
||||
X1s,Y1s,rho1s,X2s,Y2s,rho2s)
|
||||
|
||||
call self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, &
|
||||
Omega1s,rho1s,Omega2s,rho2s,SigT)
|
||||
|
||||
call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, &
|
||||
Omega1s,rho1s,Omega2s,rho2s,Z)
|
||||
|
||||
iblock = 4
|
||||
|
||||
call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO, &
|
||||
X1t,Y1t,rho1t,X2t,Y2t,rho2t)
|
||||
|
||||
call self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, &
|
||||
Omega1t,rho1t,Omega2t,rho2t,SigT)
|
||||
|
||||
call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, &
|
||||
Omega1t,rho1t,Omega2t,rho2t,Z)
|
||||
|
||||
Z(:) = 1d0/(1d0 - Z(:))
|
||||
|
||||
! Make correlation self-energy Hermitian and transform it back to AO basis
|
||||
|
||||
SigTp = 0.5d0*(SigT + transpose(SigT))
|
||||
SigTm = 0.5d0*(SigT - transpose(SigT))
|
||||
|
||||
call MOtoAO_transform(nBas,S,c,SigTp)
|
||||
|
||||
! Solve the quasi-particle equation
|
||||
|
||||
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + SigTp(:,:)
|
||||
|
||||
! Compute commutator and convergence criteria
|
||||
|
||||
error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)
|
||||
|
||||
! DIIS extrapolation
|
||||
|
||||
n_diis = min(n_diis+1,max_diis)
|
||||
if(abs(rcond) > 1d-7) then
|
||||
call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F)
|
||||
else
|
||||
n_diis = 0
|
||||
end if
|
||||
|
||||
! Diagonalize Hamiltonian in AO basis
|
||||
|
||||
Fp = matmul(transpose(X),matmul(F,X))
|
||||
cp(:,:) = Fp(:,:)
|
||||
call diagonalize_matrix(nBas,cp,eGT)
|
||||
c = matmul(X,cp)
|
||||
SigTp = matmul(transpose(c),matmul(SigTp,c))
|
||||
|
||||
! Compute new density matrix in the AO basis
|
||||
|
||||
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
|
||||
|
||||
! Save quasiparticles energy for next cycle
|
||||
|
||||
Conv = maxval(abs(eGT - eOld))
|
||||
eOld(:) = eGT(:)
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Compute total energy
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
! Kinetic energy
|
||||
|
||||
ET = trace_matrix(nBas,matmul(P,T))
|
||||
|
||||
! Potential energy
|
||||
|
||||
EV = trace_matrix(nBas,matmul(P,V))
|
||||
|
||||
! Coulomb energy
|
||||
|
||||
EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
|
||||
|
||||
! Exchange energy
|
||||
|
||||
Ex = 0.25d0*trace_matrix(nBas,matmul(P,K))
|
||||
|
||||
! Total energy
|
||||
|
||||
EqsGT = ET + EV + EJ + Ex
|
||||
|
||||
! Print results
|
||||
|
||||
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
|
||||
call print_qsGT(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigTp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole)
|
||||
|
||||
enddo
|
||||
!------------------------------------------------------------------------
|
||||
! End main loop
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
! Compute the ppRPA correlation energy
|
||||
|
||||
ispin = 1
|
||||
iblock = 3
|
||||
call linear_response_pp(iblock,.false.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,ERI_MO, &
|
||||
Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin))
|
||||
ispin = 2
|
||||
iblock = 4
|
||||
call linear_response_pp(iblock,.false.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,ERI_MO, &
|
||||
Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin))
|
||||
EcRPA(1) = EcRPA(1) - EcRPA(2)
|
||||
EcRPA(2) = 3d0*EcRPA(2)
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@qsGT correlation energy (singlet) =',EcRPA(1)
|
||||
write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@qsGT correlation energy (triplet) =',EcRPA(2)
|
||||
write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@qsGT correlation energy =',EcRPA(1) + EcRPA(2)
|
||||
write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@qsGT total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
! Did it actually converge?
|
||||
|
||||
if(nSCF == maxSCF+1) then
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)' Convergence failed '
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)
|
||||
|
||||
stop
|
||||
|
||||
endif
|
||||
|
||||
! Deallocate memory
|
||||
|
||||
deallocate(c,cp,P,F,Fp,J,K,SigT,SigTp,SigTm,Z,error,error_diis,F_diis)
|
||||
|
||||
! Perform BSE calculation
|
||||
|
||||
if(BSE) then
|
||||
|
||||
call Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
|
||||
Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s,Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,rho1t,rho2t, &
|
||||
ERI_MO,dipole_int_MO,eGT,eGT,EcBSE)
|
||||
|
||||
if(exchange_kernel) then
|
||||
|
||||
EcBSE(1) = 0.5d0*EcBSE(1)
|
||||
EcBSE(2) = 1.5d0*EcBSE(2)
|
||||
|
||||
end if
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGT correlation energy (singlet) =',EcBSE(1)
|
||||
write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGT correlation energy (triplet) =',EcBSE(2)
|
||||
write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGT correlation energy =',EcBSE(1) + EcBSE(2)
|
||||
write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGT total energy =',ENuc + EqsGT + EcBSE(1) + EcBSE(2)
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
! Compute the BSE correlation energy via the adiabatic connection
|
||||
|
||||
if(doACFDT) then
|
||||
|
||||
write(*,*) '------------------------------------------------------'
|
||||
write(*,*) 'Adiabatic connection version of BSE correlation energy'
|
||||
write(*,*) '------------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
if(doXBS) then
|
||||
|
||||
write(*,*) '*** scaled screening version (XBS) ***'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
call ACFDT(exchange_kernel,doXBS,.true.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGT,eGT,EcAC)
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGT correlation energy (singlet) =',EcAC(1)
|
||||
write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGT correlation energy (triplet) =',EcAC(2)
|
||||
write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGT correlation energy =',EcAC(1) + EcAC(2)
|
||||
write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGT total energy =',ENuc + EqsGT + EcAC(1) + EcAC(2)
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
end if
|
||||
|
||||
end subroutine qsGT
|
||||
@ -28,7 +28,7 @@ subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: SigT(nBas)
|
||||
double precision,intent(inout) :: SigT(nBas,nBas)
|
||||
|
||||
!----------------------------------------------
|
||||
! Occupied part of the T-matrix self-energy
|
||||
@ -38,8 +38,8 @@ subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2
|
||||
do q=nC+1,nBas-nR
|
||||
do i=nC+1,nO
|
||||
do cd=1,nVV
|
||||
eps = e(p) + e(i) - Omega1(cd)
|
||||
SigT(p) = SigT(p) + rho1(p,i,cd)*rho1(q,i,cd)*eps/(eps**2 + eta**2)
|
||||
eps = e(p) + e(i) - Omega1(cd)
|
||||
SigT(p,q) = SigT(p,q) + rho1(p,i,cd)*rho1(q,i,cd)*eps/(eps**2 + eta**2)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
@ -51,10 +51,10 @@ subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2
|
||||
|
||||
do p=nC+1,nBas-nR
|
||||
do q=nC+1,nBas-nR
|
||||
do a=1,nV-nR
|
||||
do a=nO+1,nBas-nR
|
||||
do kl=1,nOO
|
||||
eps = e(p) + e(nO+a) - Omega2(kl)
|
||||
SigT(p) = SigT(p) + rho2(p,nO+a,kl)*rho2(q,nO+a,kl)*eps/(eps**2 + eta**2)
|
||||
eps = e(p) + e(a) - Omega2(kl)
|
||||
SigT(p,q) = SigT(p,q) + rho2(p,a,kl)*rho2(q,a,kl)*eps/(eps**2 + eta**2)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
@ -28,7 +28,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,O
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: SigT(nBas)
|
||||
double precision,intent(inout) :: SigT(nBas)
|
||||
|
||||
!----------------------------------------------
|
||||
! Occupied part of the T-matrix self-energy
|
||||
@ -48,10 +48,10 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,O
|
||||
!----------------------------------------------
|
||||
|
||||
do p=nC+1,nBas-nR
|
||||
do a=1,nV-nR
|
||||
do a=nO+1,nBas-nR
|
||||
do kl=1,nOO
|
||||
eps = e(p) + e(nO+a) - Omega2(kl)
|
||||
SigT(p) = SigT(p) + rho2(p,nO+a,kl)**2*eps/(eps**2 + eta**2)
|
||||
eps = e(p) + e(a) - Omega2(kl)
|
||||
SigT(p) = SigT(p) + rho2(p,a,kl)**2*eps/(eps**2 + eta**2)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
@ -1062,6 +1062,26 @@ program QuAcK
|
||||
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Perform qsGT calculation
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
if(doqsGT) then
|
||||
|
||||
call cpu_time(start_qsGT)
|
||||
|
||||
call qsGT(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS, &
|
||||
BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GW,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)
|
||||
|
||||
call cpu_time(end_qsGT)
|
||||
|
||||
t_qsGT = end_qsGT - start_qsGT
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGT = ',t_qsGT,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Information for Monte Carlo calculations
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
Loading…
Reference in New Issue
Block a user