mirror of
https://github.com/pfloos/quack
synced 2024-12-22 20:34:46 +01:00
evGT
This commit is contained in:
parent
26aa96f1e1
commit
2f835304f7
68
input/basis
68
input/basis
@ -1,49 +1,39 @@
|
|||||||
1 9
|
1 10
|
||||||
S 8
|
S 8
|
||||||
1 1469.0000000 0.0007660
|
1 24350.0000000 0.0005020
|
||||||
2 220.5000000 0.0058920
|
2 3650.0000000 0.0038810
|
||||||
3 50.2600000 0.0296710
|
3 829.6000000 0.0199970
|
||||||
4 14.2400000 0.1091800
|
4 234.0000000 0.0784180
|
||||||
5 4.5810000 0.2827890
|
5 75.6100000 0.2296760
|
||||||
6 1.5800000 0.4531230
|
6 26.7300000 0.4327220
|
||||||
7 0.5640000 0.2747740
|
7 9.9270000 0.3506420
|
||||||
8 0.0734500 0.0097510
|
8 1.1020000 -0.0076450
|
||||||
S 8
|
S 8
|
||||||
1 1469.0000000 -0.0001200
|
1 24350.0000000 -0.0001180
|
||||||
2 220.5000000 -0.0009230
|
2 3650.0000000 -0.0009150
|
||||||
3 50.2600000 -0.0046890
|
3 829.6000000 -0.0047370
|
||||||
4 14.2400000 -0.0176820
|
4 234.0000000 -0.0192330
|
||||||
5 4.5810000 -0.0489020
|
5 75.6100000 -0.0603690
|
||||||
6 1.5800000 -0.0960090
|
6 26.7300000 -0.1425080
|
||||||
7 0.5640000 -0.1363800
|
7 9.9270000 -0.1777100
|
||||||
8 0.0734500 0.5751020
|
8 1.1020000 0.6058360
|
||||||
S 1
|
S 1
|
||||||
1 0.0280500 1.0000000
|
1 2.8360000 1.0000000
|
||||||
S 1
|
S 1
|
||||||
1 0.0086400 1.0000000
|
1 0.3782000 1.0000000
|
||||||
P 3
|
P 3
|
||||||
1 1.5340000 0.0227840
|
1 54.7000000 0.0171510
|
||||||
2 0.2749000 0.1391070
|
2 12.4300000 0.1076560
|
||||||
3 0.0736200 0.5003750
|
3 3.6790000 0.3216810
|
||||||
P 1
|
P 1
|
||||||
1 0.0240300 1.0000000
|
1 1.1430000 1.0000000
|
||||||
P 1
|
P 1
|
||||||
1 0.0057900 1.0000000
|
1 0.3300000 1.0000000
|
||||||
D 1
|
D 1
|
||||||
1 0.1239000 1.0000000
|
1 4.0140000 1.0000000
|
||||||
D 1
|
D 1
|
||||||
1 0.0725000 1.0000000
|
1 1.0960000 1.0000000
|
||||||
2 5
|
F 1
|
||||||
S 3
|
1 2.5440000 1.0000000
|
||||||
1 13.0100000 0.0196850
|
|
||||||
2 1.9620000 0.1379770
|
|
||||||
3 0.4446000 0.4781480
|
|
||||||
S 1
|
|
||||||
1 0.1220000 1.0000000
|
|
||||||
S 1
|
|
||||||
1 0.0297400 1.0000000
|
|
||||||
P 1
|
|
||||||
1 0.7270000 1.0000000
|
|
||||||
P 1
|
|
||||||
1 0.1410000 1.0000000
|
|
||||||
|
|
||||||
|
|||||||
@ -13,6 +13,6 @@
|
|||||||
# G0W0 evGW qsGW
|
# G0W0 evGW qsGW
|
||||||
F F F
|
F F F
|
||||||
# G0T0 evGT qsGT
|
# G0T0 evGT qsGT
|
||||||
T F F
|
T T F
|
||||||
# MCMP2
|
# MCMP2
|
||||||
F
|
F
|
||||||
|
|||||||
@ -1,5 +1,4 @@
|
|||||||
# nAt nEla nElb nCore nRyd
|
# nAt nEla nElb nCore nRyd
|
||||||
2 2 2 0 0
|
1 5 5 0 0
|
||||||
# Znuc x y z
|
# Znuc x y z
|
||||||
Li 0. 0. 0.
|
Ne 0.0 0.0 0.0
|
||||||
H 0. 0. 3.099
|
|
||||||
|
|||||||
@ -1,4 +1,3 @@
|
|||||||
2
|
1
|
||||||
|
|
||||||
Li 0.0000000000 0.0000000000 0.0000000000
|
Ne 0.0000000000 0.0000000000 0.0000000000
|
||||||
H 0.0000000000 0.0000000000 1.6399202947
|
|
||||||
|
|||||||
68
input/weight
68
input/weight
@ -1,49 +1,39 @@
|
|||||||
1 9
|
1 10
|
||||||
S 8
|
S 8
|
||||||
1 1469.0000000 0.0007660
|
1 24350.0000000 0.0005020
|
||||||
2 220.5000000 0.0058920
|
2 3650.0000000 0.0038810
|
||||||
3 50.2600000 0.0296710
|
3 829.6000000 0.0199970
|
||||||
4 14.2400000 0.1091800
|
4 234.0000000 0.0784180
|
||||||
5 4.5810000 0.2827890
|
5 75.6100000 0.2296760
|
||||||
6 1.5800000 0.4531230
|
6 26.7300000 0.4327220
|
||||||
7 0.5640000 0.2747740
|
7 9.9270000 0.3506420
|
||||||
8 0.0734500 0.0097510
|
8 1.1020000 -0.0076450
|
||||||
S 8
|
S 8
|
||||||
1 1469.0000000 -0.0001200
|
1 24350.0000000 -0.0001180
|
||||||
2 220.5000000 -0.0009230
|
2 3650.0000000 -0.0009150
|
||||||
3 50.2600000 -0.0046890
|
3 829.6000000 -0.0047370
|
||||||
4 14.2400000 -0.0176820
|
4 234.0000000 -0.0192330
|
||||||
5 4.5810000 -0.0489020
|
5 75.6100000 -0.0603690
|
||||||
6 1.5800000 -0.0960090
|
6 26.7300000 -0.1425080
|
||||||
7 0.5640000 -0.1363800
|
7 9.9270000 -0.1777100
|
||||||
8 0.0734500 0.5751020
|
8 1.1020000 0.6058360
|
||||||
S 1
|
S 1
|
||||||
1 0.0280500 1.0000000
|
1 2.8360000 1.0000000
|
||||||
S 1
|
S 1
|
||||||
1 0.0086400 1.0000000
|
1 0.3782000 1.0000000
|
||||||
P 3
|
P 3
|
||||||
1 1.5340000 0.0227840
|
1 54.7000000 0.0171510
|
||||||
2 0.2749000 0.1391070
|
2 12.4300000 0.1076560
|
||||||
3 0.0736200 0.5003750
|
3 3.6790000 0.3216810
|
||||||
P 1
|
P 1
|
||||||
1 0.0240300 1.0000000
|
1 1.1430000 1.0000000
|
||||||
P 1
|
P 1
|
||||||
1 0.0057900 1.0000000
|
1 0.3300000 1.0000000
|
||||||
D 1
|
D 1
|
||||||
1 0.1239000 1.0000000
|
1 4.0140000 1.0000000
|
||||||
D 1
|
D 1
|
||||||
1 0.0725000 1.0000000
|
1 1.0960000 1.0000000
|
||||||
2 5
|
F 1
|
||||||
S 3
|
1 2.5440000 1.0000000
|
||||||
1 13.0100000 0.0196850
|
|
||||||
2 1.9620000 0.1379770
|
|
||||||
3 0.4446000 0.4781480
|
|
||||||
S 1
|
|
||||||
1 0.1220000 1.0000000
|
|
||||||
S 1
|
|
||||||
1 0.0297400 1.0000000
|
|
||||||
P 1
|
|
||||||
1 0.7270000 1.0000000
|
|
||||||
P 1
|
|
||||||
1 0.1410000 1.0000000
|
|
||||||
|
|
||||||
|
|||||||
@ -182,6 +182,28 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m
|
|||||||
|
|
||||||
call print_G0T0(nBas,nO,eHF(:),ENuc,ERHF,SigT(:),Z(:),eG0T0(:),EcRPA(:))
|
call print_G0T0(nBas,nO,eHF(:),ENuc,ERHF,SigT(:),Z(:),eG0T0(:),EcRPA(:))
|
||||||
|
|
||||||
|
! Compute the ppRPA correlation energy
|
||||||
|
|
||||||
|
ispin = 1
|
||||||
|
iblock = 3
|
||||||
|
call linear_response_pp(iblock,.false.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eG0T0(:),ERI(:,:,:,:), &
|
||||||
|
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,eG0T0(:),ERI(:,:,:,:), &
|
||||||
|
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@RPA@G0W0 correlation energy (singlet) =',EcRPA(1)
|
||||||
|
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy (triplet) =',EcRPA(2)
|
||||||
|
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy =',EcRPA(1) + EcRPA(2)
|
||||||
|
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
|
||||||
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
! Perform BSE calculation
|
! Perform BSE calculation
|
||||||
|
|
||||||
if(BSE) then
|
if(BSE) then
|
||||||
|
|||||||
@ -685,7 +685,6 @@ program QuAcK
|
|||||||
if(doG0T0) then
|
if(doG0T0) then
|
||||||
|
|
||||||
call cpu_time(start_G0T0)
|
call cpu_time(start_G0T0)
|
||||||
! call soG0T0(eta,nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF)
|
|
||||||
call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA, &
|
call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA, &
|
||||||
singlet_manifold,triplet_manifold,linGW,eta, &
|
singlet_manifold,triplet_manifold,linGW,eta, &
|
||||||
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO_basis,eHF)
|
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO_basis,eHF)
|
||||||
@ -698,6 +697,24 @@ program QuAcK
|
|||||||
|
|
||||||
end if
|
end if
|
||||||
|
|
||||||
|
!------------------------------------------------------------------------
|
||||||
|
! Perform evGT calculatiom
|
||||||
|
!------------------------------------------------------------------------
|
||||||
|
|
||||||
|
if(doevGT) then
|
||||||
|
|
||||||
|
call cpu_time(start_evGT)
|
||||||
|
call evGT(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_manifold, &
|
||||||
|
eta,nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO_basis,eHF,eG0T0)
|
||||||
|
call cpu_time(end_evGT)
|
||||||
|
|
||||||
|
t_evGT = end_evGT - start_evGT
|
||||||
|
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGT = ',t_evGT,' seconds'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
!------------------------------------------------------------------------
|
||||||
! Information for Monte Carlo calculations
|
! Information for Monte Carlo calculations
|
||||||
!------------------------------------------------------------------------
|
!------------------------------------------------------------------------
|
||||||
|
|||||||
318
src/QuAcK/evGT.f90
Normal file
318
src/QuAcK/evGT.f90
Normal file
@ -0,0 +1,318 @@
|
|||||||
|
subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_manifold, &
|
||||||
|
eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eG0T0)
|
||||||
|
|
||||||
|
! Perform eigenvalue self-consistent calculation with a T-matrix self-energy (evGT)
|
||||||
|
|
||||||
|
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
|
||||||
|
logical,intent(in) :: singlet_manifold
|
||||||
|
logical,intent(in) :: triplet_manifold
|
||||||
|
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) :: ENuc
|
||||||
|
double precision,intent(in) :: ERHF
|
||||||
|
double precision,intent(in) :: eHF(nBas)
|
||||||
|
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||||
|
double precision,intent(in) :: eG0T0(nBas)
|
||||||
|
|
||||||
|
|
||||||
|
! Local variables
|
||||||
|
|
||||||
|
logical :: linear_mixing
|
||||||
|
integer :: nSCF
|
||||||
|
integer :: n_diis
|
||||||
|
double precision :: rcond
|
||||||
|
double precision :: Conv
|
||||||
|
integer :: ispin
|
||||||
|
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)
|
||||||
|
double precision,allocatable :: error_diis(:,:)
|
||||||
|
double precision,allocatable :: e_diis(:,:)
|
||||||
|
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 :: 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
|
||||||
|
|
||||||
|
write(*,*)
|
||||||
|
write(*,*)'************************************************'
|
||||||
|
write(*,*)'| Self-consistent evGT 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
|
||||||
|
|
||||||
|
! Memory allocation
|
||||||
|
|
||||||
|
allocate(Omega1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), &
|
||||||
|
Omega2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), &
|
||||||
|
rho1s(nBas,nO,nVVs),rho2s(nBas,nV,nOOs), &
|
||||||
|
Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), &
|
||||||
|
Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), &
|
||||||
|
rho1t(nBas,nO,nVVt),rho2t(nBas,nV,nOOt), &
|
||||||
|
eGT(nBas),eOld(nBas),Z(nBas),SigT(nBas), &
|
||||||
|
error_diis(nBas,max_diis),e_diis(nBas,max_diis))
|
||||||
|
|
||||||
|
! Initialization
|
||||||
|
|
||||||
|
nSCF = 0
|
||||||
|
n_diis = 0
|
||||||
|
Conv = 1d0
|
||||||
|
e_diis(:,:) = 0d0
|
||||||
|
error_diis(:,:) = 0d0
|
||||||
|
eGT(:) = eG0T0(:)
|
||||||
|
eOld(:) = eGT(:)
|
||||||
|
Z(:) = 1d0
|
||||||
|
|
||||||
|
!------------------------------------------------------------------------
|
||||||
|
! Main loop
|
||||||
|
!------------------------------------------------------------------------
|
||||||
|
|
||||||
|
do while(Conv > thresh .and. nSCF <= maxSCF)
|
||||||
|
|
||||||
|
!----------------------------------------------
|
||||||
|
! alpha-beta block
|
||||||
|
!----------------------------------------------
|
||||||
|
|
||||||
|
ispin = 1
|
||||||
|
iblock = 3
|
||||||
|
|
||||||
|
! Compute linear response
|
||||||
|
|
||||||
|
call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF(:),ERI(:,:,:,:), &
|
||||||
|
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
|
||||||
|
!----------------------------------------------
|
||||||
|
|
||||||
|
ispin = 2
|
||||||
|
iblock = 4
|
||||||
|
|
||||||
|
! Compute linear response
|
||||||
|
|
||||||
|
call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF(:),ERI(:,:,:,:), &
|
||||||
|
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(:))
|
||||||
|
|
||||||
|
!----------------------------------------------
|
||||||
|
! Compute T-matrix version of the self-energy
|
||||||
|
!----------------------------------------------
|
||||||
|
|
||||||
|
SigT(:) = 0d0
|
||||||
|
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(:,:,:,:), &
|
||||||
|
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 renormalization_factor_Tmatrix(alpha,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(:,:,:,:), &
|
||||||
|
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 renormalization_factor_Tmatrix(alpha,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT(:), &
|
||||||
|
Omega1t(:),rho1t(:,:,:),Omega2t(:),rho2t(:,:,:),Z(:))
|
||||||
|
|
||||||
|
Z(:) = 1d0/(1d0 - Z(:))
|
||||||
|
|
||||||
|
! Solve the quasi-particle equation
|
||||||
|
|
||||||
|
!----------------------------------------------
|
||||||
|
! Solve the quasi-particle equation
|
||||||
|
!----------------------------------------------
|
||||||
|
|
||||||
|
eGT(:) = eHF(:) + SigT(:)
|
||||||
|
|
||||||
|
! Convergence criteria
|
||||||
|
|
||||||
|
Conv = maxval(abs(eGT(:) - eOld(:)))
|
||||||
|
|
||||||
|
!----------------------------------------------
|
||||||
|
! Dump results
|
||||||
|
!----------------------------------------------
|
||||||
|
|
||||||
|
call print_evGT(nBas,nO,nSCF,Conv,eHF(:),SigT(:),Z(:),eGT(:))
|
||||||
|
|
||||||
|
! DIIS extrapolation
|
||||||
|
|
||||||
|
n_diis = min(n_diis+1,max_diis)
|
||||||
|
call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGT(:)-eOld(:),eGT(:))
|
||||||
|
|
||||||
|
! Reset DIIS if required
|
||||||
|
|
||||||
|
if(abs(rcond) < 1d-15) n_diis = 0
|
||||||
|
|
||||||
|
! Save quasiparticles energy for next cycle
|
||||||
|
|
||||||
|
eOld(:) = eGT(:)
|
||||||
|
|
||||||
|
! Increment
|
||||||
|
|
||||||
|
nSCF = nSCF + 1
|
||||||
|
|
||||||
|
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(:,:,:,:), &
|
||||||
|
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(:,:,:,:), &
|
||||||
|
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@evGT correlation energy (singlet) =',EcRPA(1)
|
||||||
|
write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@evGT correlation energy (triplet) =',EcRPA(2)
|
||||||
|
write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@evGT correlation energy =',EcRPA(1) + EcRPA(2)
|
||||||
|
write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@evGT total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
|
||||||
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
|
||||||
|
! Perform BSE calculation
|
||||||
|
|
||||||
|
if(BSE) then
|
||||||
|
|
||||||
|
allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin),rho(nBas,nBas,nS,nspin))
|
||||||
|
|
||||||
|
call Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, &
|
||||||
|
nBas,nC,nO,nV,nR,nS,ERI,eGT,eGT,Omega,XpY,XmY,rho,EcRPA,EcBSE)
|
||||||
|
|
||||||
|
if(exchange_kernel) then
|
||||||
|
|
||||||
|
EcRPA(1) = 0.5d0*EcRPA(1)
|
||||||
|
EcRPA(2) = 1.5d0*EcRPA(1)
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
write(*,*)
|
||||||
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGT correlation energy (singlet) =',EcBSE(1)
|
||||||
|
write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGT correlation energy (triplet) =',EcBSE(2)
|
||||||
|
write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGT correlation energy =',EcBSE(1) + EcBSE(2)
|
||||||
|
write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGT total energy =',ENuc + ERHF + 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,BSE,singlet_manifold,triplet_manifold,eta, &
|
||||||
|
nBas,nC,nO,nV,nR,nS,ERI,eGT,eGT,Omega,XpY,XmY,rho,EcAC)
|
||||||
|
|
||||||
|
if(exchange_kernel) then
|
||||||
|
|
||||||
|
EcAC(1) = 0.5d0*EcAC(1)
|
||||||
|
EcAC(2) = 1.5d0*EcAC(1)
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
write(*,*)
|
||||||
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
write(*,'(2X,A50,F20.10)') 'AC@BSE@evGT correlation energy (singlet) =',EcAC(1)
|
||||||
|
write(*,'(2X,A50,F20.10)') 'AC@BSE@evGT correlation energy (triplet) =',EcAC(2)
|
||||||
|
write(*,'(2X,A50,F20.10)') 'AC@BSE@evGT correlation energy =',EcAC(1) + EcAC(2)
|
||||||
|
write(*,'(2X,A50,F20.10)') 'AC@BSE@evGT total energy =',ENuc + ERHF + EcAC(1) + EcAC(2)
|
||||||
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
end subroutine evGT
|
||||||
@ -78,8 +78,6 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
|
|||||||
if(COHSEX) write(*,*) 'COHSEX approximation activated!'
|
if(COHSEX) write(*,*) 'COHSEX approximation activated!'
|
||||||
write(*,*)
|
write(*,*)
|
||||||
|
|
||||||
! Switch off exchange for G0W0
|
|
||||||
|
|
||||||
! Linear mixing
|
! Linear mixing
|
||||||
|
|
||||||
linear_mixing = .false.
|
linear_mixing = .false.
|
||||||
@ -154,7 +152,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
|
|||||||
! Print results
|
! Print results
|
||||||
|
|
||||||
! call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
|
! call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
|
||||||
call print_evGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA(ispin),EcGM)
|
call print_evGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW)
|
||||||
|
|
||||||
! Linear mixing or DIIS extrapolation
|
! Linear mixing or DIIS extrapolation
|
||||||
|
|
||||||
|
|||||||
@ -51,11 +51,11 @@ subroutine print_G0T0(nBas,nO,e,ENuc,ERHF,SigT,Z,eGW,EcRPA)
|
|||||||
write(*,'(2X,A40,F15.6)') 'G0T0 LUMO energy (eV) :',eGW(LUMO)*HaToeV
|
write(*,'(2X,A40,F15.6)') 'G0T0 LUMO energy (eV) :',eGW(LUMO)*HaToeV
|
||||||
write(*,'(2X,A40,F15.6)') 'G0T0 HOMO-LUMO gap (eV) :',Gap*HaToeV
|
write(*,'(2X,A40,F15.6)') 'G0T0 HOMO-LUMO gap (eV) :',Gap*HaToeV
|
||||||
write(*,*)'-------------------------------------------------------------------------------'
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
write(*,'(2X,A40,F15.6)') 'RPA@G0T0 correlation energy (singlet) =',EcRPA(1)
|
! write(*,'(2X,A40,F15.6)') 'RPA@G0T0 correlation energy (singlet) =',EcRPA(1)
|
||||||
write(*,'(2X,A40,F15.6)') 'RPA@G0T0 correlation energy (triplet) =',EcRPA(2)
|
! write(*,'(2X,A40,F15.6)') 'RPA@G0T0 correlation energy (triplet) =',EcRPA(2)
|
||||||
write(*,'(2X,A40,F15.6)') 'RPA@G0T0 correlation energy =',EcRPA(1) + EcRPA(2)
|
! write(*,'(2X,A40,F15.6)') 'RPA@G0T0 correlation energy =',EcRPA(1) + EcRPA(2)
|
||||||
write(*,'(2X,A40,F15.6)') 'RPA@G0T0 total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
|
! write(*,'(2X,A40,F15.6)') 'RPA@G0T0 total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
|
||||||
write(*,*)'-------------------------------------------------------------------------------'
|
! write(*,*)'-------------------------------------------------------------------------------'
|
||||||
write(*,*)
|
write(*,*)
|
||||||
|
|
||||||
end subroutine print_G0T0
|
end subroutine print_G0T0
|
||||||
|
|||||||
54
src/QuAcK/print_evGT.f90
Normal file
54
src/QuAcK/print_evGT.f90
Normal file
@ -0,0 +1,54 @@
|
|||||||
|
subroutine print_evGT(nBas,nO,nSCF,Conv,eHF,SigT,Z,eGT)
|
||||||
|
|
||||||
|
! Print one-electron energies and other stuff for evGT
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
include 'parameters.h'
|
||||||
|
|
||||||
|
integer,intent(in) :: nBas
|
||||||
|
integer,intent(in) :: nO
|
||||||
|
integer,intent(in) :: nSCF
|
||||||
|
double precision,intent(in) :: Conv
|
||||||
|
double precision,intent(in) :: eHF(nBas)
|
||||||
|
double precision,intent(in) :: SigT(nBas)
|
||||||
|
double precision,intent(in) :: Z(nBas)
|
||||||
|
double precision,intent(in) :: eGT(nBas)
|
||||||
|
|
||||||
|
integer :: x,HOMO,LUMO
|
||||||
|
double precision :: Gap
|
||||||
|
|
||||||
|
! HOMO and LUMO
|
||||||
|
|
||||||
|
HOMO = nO
|
||||||
|
LUMO = HOMO + 1
|
||||||
|
Gap = eGT(LUMO)-eGT(HOMO)
|
||||||
|
|
||||||
|
! Dump results
|
||||||
|
|
||||||
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
if(nSCF < 10) then
|
||||||
|
write(*,'(1X,A21,I1,A1,I1,A12)')' Self-consistent evG',nSCF,'T',nSCF,' calculation'
|
||||||
|
else
|
||||||
|
write(*,'(1X,A21,I2,A1,I2,A12)')' Self-consistent evG',nSCF,'T',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)','|','Sigma_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,'|',SigT(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)') 'evGT HOMO energy (eV):',eGT(HOMO)*HaToeV
|
||||||
|
write(*,'(2X,A30,F15.6)') 'evGT LUMO energy (eV):',eGT(LUMO)*HaToeV
|
||||||
|
write(*,'(2X,A30,F15.6)') 'evGT HOMO-LUMO gap (eV):',Gap*HaToeV
|
||||||
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
end subroutine print_evGT
|
||||||
@ -1,4 +1,4 @@
|
|||||||
subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA,EcGM)
|
subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW)
|
||||||
|
|
||||||
! Print one-electron energies and other stuff for evGW
|
! Print one-electron energies and other stuff for evGW
|
||||||
|
|
||||||
@ -8,8 +8,6 @@ subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA,EcGM)
|
|||||||
integer,intent(in) :: nBas,nO,nSCF
|
integer,intent(in) :: nBas,nO,nSCF
|
||||||
double precision,intent(in) :: ENuc
|
double precision,intent(in) :: ENuc
|
||||||
double precision,intent(in) :: EHF
|
double precision,intent(in) :: EHF
|
||||||
double precision,intent(in) :: EcRPA
|
|
||||||
double precision,intent(in) :: EcGM
|
|
||||||
double precision,intent(in) :: Conv,e(nBas),SigmaC(nBas),Z(nBas),eGW(nBas)
|
double precision,intent(in) :: Conv,e(nBas),SigmaC(nBas),Z(nBas),eGW(nBas)
|
||||||
|
|
||||||
integer :: x,HOMO,LUMO
|
integer :: x,HOMO,LUMO
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user