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

big cleaning

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-09-24 11:56:06 +02:00
parent 08f3567d20
commit 8910ead99f
30 changed files with 402 additions and 479 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
examples/molecule.H2
View File

@ -2,4 +2,4 @@
2 1 1 0 0 2 1 1 0 0
# Znuc x y z # Znuc x y z
H 0. 0. 0. H 0. 0. 0.
H 0. 0. 1.4 H 0. 0. 1.399

55
input/basis
View File

@ -1,30 +1,39 @@
1 6 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
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
1 0.3300000 1.0000000
D 1 D 1
1 0.1239000 1.0000000 1 4.0140000 1.0000000
D 1
1 1.0960000 1.0000000
F 1
1 2.5440000 1.0000000

4
input/methods
View File

@ -1,5 +1,5 @@
# RHF UHF MOM # RHF UHF MOM
F T F T F F
# MP2 MP3 MP2-F12 # MP2 MP3 MP2-F12
F F F F F F
# CCD CCSD CCSD(T) # CCD CCSD CCSD(T)
@ -13,7 +13,7 @@
# G0F2 evGF2 G0F3 evGF3 # G0F2 evGF2 G0F3 evGF3
F F F F F F F F
# G0W0 evGW qsGW # G0W0 evGW qsGW
T F F F F T
# G0T0 evGT qsGT # G0T0 evGT qsGT
F F F F F F
# MCMP2 # MCMP2

4
input/molecule
View File

@ -1,4 +1,4 @@
# nAt nEla nElb nCore nRyd # nAt nEla nElb nCore nRyd
1 2 1 0 0 1 5 5 0 0
# Znuc x y z # Znuc x y z
Li 0.0 0.0 0.0 Ne 0.0 0.0 0.0

2
input/molecule.xyz
View File

@ -1,3 +1,3 @@
1 1
Li 0.0000000000 0.0000000000 0.0000000000 Ne 0.0000000000 0.0000000000 0.0000000000

2
input/options
View File

@ -13,6 +13,6 @@
# ACFDT: AC Kx XBS # ACFDT: AC Kx XBS
F F T F F T
# BSE: BSE dBSE dTDA evDyn # BSE: BSE dBSE dTDA evDyn
T F F F T F T T
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
1000000 100000 10 0.3 10000 1234 T 1000000 100000 10 0.3 10000 1234 T

67
src/QuAcK/ACFDT.f90
View File

@ -1,5 +1,4 @@
subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,triplet_manifold,eta, & subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eW,e,EcAC)
nBas,nC,nO,nV,nR,nS,ERI,eW,e,Omega,XpY,XmY,rho,EcAC)
! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem ! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem
@ -15,8 +14,8 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl
logical,intent(in) :: TDA_W logical,intent(in) :: TDA_W
logical,intent(in) :: TDA logical,intent(in) :: TDA
logical,intent(in) :: BSE logical,intent(in) :: BSE
logical,intent(in) :: singlet_manifold logical,intent(in) :: singlet
logical,intent(in) :: triplet_manifold logical,intent(in) :: triplet
double precision,intent(in) :: eta double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR,nS integer,intent(in) :: nBas,nC,nO,nV,nR,nS
@ -24,11 +23,6 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl
double precision,intent(in) :: e(nBas) double precision,intent(in) :: e(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision :: Omega(nS,nspin)
double precision :: XpY(nS,nS,nspin)
double precision :: XmY(nS,nS,nspin)
double precision :: rho(nBas,nBas,nS,nspin)
! Local variables ! Local variables
integer :: ispin integer :: ispin
@ -37,6 +31,16 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl
double precision :: lambda double precision :: lambda
double precision,allocatable :: Ec(:,:) double precision,allocatable :: Ec(:,:)
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 :: Omega(:,:)
double precision,allocatable :: XpY(:,:,:)
double precision,allocatable :: XmY(:,:,:)
! Output variables ! Output variables
double precision,intent(out) :: EcAC(nspin) double precision,intent(out) :: EcAC(nspin)
@ -44,6 +48,8 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl
! Memory allocation ! Memory allocation
allocate(Ec(nAC,nspin)) allocate(Ec(nAC,nspin))
allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS))
allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin))
! Antisymmetrized kernel version ! Antisymmetrized kernel version
@ -58,12 +64,20 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl
EcAC(:) = 0d0 EcAC(:) = 0d0
Ec(:,:) = 0d0 Ec(:,:) = 0d0
! Compute (singlet) RPA screening
isp_W = 1
EcRPA = 0d0
call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,OmRPA, &
rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
! Singlet manifold ! Singlet manifold
if(singlet_manifold) then if(singlet) then
ispin = 1 ispin = 1
isp_W = 1
write(*,*) '--------------' write(*,*) '--------------'
write(*,*) 'Singlet states' write(*,*) 'Singlet states'
@ -80,17 +94,17 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl
if(doXBS) then if(doXBS) then
call linear_response(isp_W,dRPA,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, & call linear_response(isp_W,dRPA,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, &
rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
end if end if
call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, & call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,OmRPA, &
rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) rho_RPA,EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, & call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, &
ERI(:,:,:,:),XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)
@ -98,6 +112,8 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl
EcAC(ispin) = 0.5d0*dot_product(wAC,Ec(:,ispin)) EcAC(ispin) = 0.5d0*dot_product(wAC,Ec(:,ispin))
if(exchange_kernel) EcAC(ispin) = 0.5d0*EcAC(ispin)
write(*,*) '-----------------------------------------------------------------------------------' write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin) write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin)
write(*,*) '-----------------------------------------------------------------------------------' write(*,*) '-----------------------------------------------------------------------------------'
@ -107,7 +123,7 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl
! Triplet manifold ! Triplet manifold
if(triplet_manifold) then if(triplet) then
ispin = 2 ispin = 2
isp_W = 1 isp_W = 1
@ -127,17 +143,16 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl
if(doXBS) then if(doXBS) then
call linear_response(isp_W,dRPA,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, & call linear_response(isp_W,dRPA,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, &
rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
end if end if
call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, & call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,OmRPA, &
rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) rho_RPA,EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, & call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))
ERI(:,:,:,:),XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)
@ -145,6 +160,8 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl
EcAC(ispin) = 0.5d0*dot_product(wAC,Ec(:,ispin)) EcAC(ispin) = 0.5d0*dot_product(wAC,Ec(:,ispin))
if(exchange_kernel) EcAC(ispin) = 1.5d0*EcAC(ispin)
write(*,*) '-----------------------------------------------------------------------------------' write(*,*) '-----------------------------------------------------------------------------------'
write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin) write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin)
write(*,*) '-----------------------------------------------------------------------------------' write(*,*) '-----------------------------------------------------------------------------------'

4
src/QuAcK/BSE2.f90
View File

@ -54,7 +54,7 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, &
! Compute BSE2 excitation energies ! Compute BSE2 excitation energies
call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF(:),ERI(:,:,:,:), & call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF(:),ERI(:,:,:,:), &
rho,EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) OmBSE(:,ispin),rho,EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('BSE2 ',ispin,nS,OmBSE(:,ispin)) call print_excitation('BSE2 ',ispin,nS,OmBSE(:,ispin))
! Compute dynamic correction for BSE via perturbation theory ! Compute dynamic correction for BSE via perturbation theory
@ -88,7 +88,7 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, &
! Compute BSE2 excitation energies ! Compute BSE2 excitation energies
call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF(:),ERI(:,:,:,:), & call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF(:),ERI(:,:,:,:), &
rho,EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) OmBSE(:,ispin),rho,EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('BSE2 ',ispin,nS,OmBSE(:,ispin)) call print_excitation('BSE2 ',ispin,nS,OmBSE(:,ispin))
! Compute dynamic correction for BSE via perturbation theory ! Compute dynamic correction for BSE via perturbation theory

82
src/QuAcK/Bethe_Salpeter.f90
View File

@ -1,5 +1,4 @@
subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, & subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eW,eGW,EcBSE)
nBas,nC,nO,nV,nR,nS,ERI,eW,eGW,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,EcRPA,EcBSE)
! Compute the Bethe-Salpeter excitation energies ! Compute the Bethe-Salpeter excitation energies
@ -13,8 +12,8 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man
logical,intent(in) :: dBSE logical,intent(in) :: dBSE
logical,intent(in) :: dTDA logical,intent(in) :: dTDA
logical,intent(in) :: evDyn logical,intent(in) :: evDyn
logical,intent(in) :: singlet_manifold logical,intent(in) :: singlet
logical,intent(in) :: triplet_manifold logical,intent(in) :: triplet
double precision,intent(in) :: eta double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR,nS integer,intent(in) :: nBas,nC,nO,nV,nR,nS
@ -22,51 +21,55 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man
double precision,intent(in) :: eGW(nBas) double precision,intent(in) :: eGW(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision :: OmRPA(nS,nspin)
double precision :: XpY_RPA(nS,nS,nspin)
double precision :: XmY_RPA(nS,nS,nspin)
double precision :: rho_RPA(nBas,nBas,nS,nspin)
! Local variables ! Local variables
integer :: ispin integer :: ispin
integer :: isp_W integer :: isp_W
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 :: OmBSE(:,:) double precision,allocatable :: OmBSE(:,:)
double precision,allocatable :: XpY_BSE(:,:,:) double precision,allocatable :: XpY_BSE(:,:,:)
double precision,allocatable :: XmY_BSE(:,:,:) double precision,allocatable :: XmY_BSE(:,:,:)
! Output variables ! Output variables
double precision,intent(out) :: EcRPA(nspin)
double precision,intent(out) :: EcBSE(nspin) double precision,intent(out) :: EcBSE(nspin)
! Memory allocation ! Memory allocation
allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS))
allocate(OmBSE(nS,nspin),XpY_BSE(nS,nS,nspin),XmY_BSE(nS,nS,nspin)) allocate(OmBSE(nS,nspin),XpY_BSE(nS,nS,nspin),XmY_BSE(nS,nS,nspin))
!---------------------------------
! Compute (singlet) RPA screening
!---------------------------------
isp_W = 1
EcRPA = 0d0
call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,OmRPA, &
rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
!------------------- !-------------------
! Singlet manifold ! Singlet manifold
!------------------- !-------------------
if(singlet_manifold) then if(singlet) then
ispin = 1 ispin = 1
isp_W = 1
EcBSE(ispin) = 0d0 EcBSE(ispin) = 0d0
! Compute (singlet) RPA screening
call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, &
rho_RPA(:,:,:,ispin),EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA(:,:,ispin),rho_RPA(:,:,:,ispin))
! Compute BSE excitation energies ! Compute BSE excitation energies
OmBSE(:,ispin) = OmRPA(:,ispin) call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,OmRPA, &
rho_RPA,EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & call print_excitation('BSE@GW ',ispin,nS,OmBSE(:,ispin))
rho_RPA(:,:,:,ispin),EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
call print_excitation('BSE ',ispin,nS,OmBSE(:,ispin))
!------------------------------------------------- !-------------------------------------------------
! Compute the dynamical screening at the BSE level ! Compute the dynamical screening at the BSE level
@ -78,12 +81,12 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man
if(evDyn) then if(evDyn) then
call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA, &
XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
else else
call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA, &
XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
end if end if
end if end if
@ -94,25 +97,16 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man
! Triplet manifold ! Triplet manifold
!------------------- !-------------------
if(triplet_manifold) then if(triplet) then
ispin = 2 ispin = 2
isp_W = 1
EcBSE(ispin) = 0d0 EcBSE(ispin) = 0d0
! Compute (singlet) RPA screening
call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, &
rho_RPA(:,:,:,ispin),EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA(:,:,ispin),rho_RPA(:,:,:,ispin))
! Compute BSE excitation energies ! Compute BSE excitation energies
OmBSE(:,ispin) = OmRPA(:,ispin) call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,OmRPA, &
rho_RPA,EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & call print_excitation('BSE@GW ',ispin,nS,OmBSE(:,ispin))
rho_RPA(:,:,:,ispin),EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
call print_excitation('BSE ',ispin,nS,OmBSE(:,ispin))
!------------------------------------------------- !-------------------------------------------------
! Compute the dynamical screening at the BSE level ! Compute the dynamical screening at the BSE level
@ -124,12 +118,12 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man
if(evDyn) then if(evDyn) then
call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA(:,ispin),OmBSE(:,ispin), & call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA, &
XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
else else
call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA(:,ispin),OmBSE(:,ispin), & call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA, &
XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
end if end if
end if end if

25
src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90
View File

@ -1,5 +1,4 @@
subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho, & subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rhO_RPA,OmBSE,Ap,Am,Bp,Bm)
Ap,Am,Bp,Bm)
! Compute the dynamic part of the Bethe-Salpeter equation matrices ! Compute the dynamic part of the Bethe-Salpeter equation matrices
@ -13,8 +12,8 @@ subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O
double precision,intent(in) :: lambda double precision,intent(in) :: lambda
double precision,intent(in) :: eGW(nBas) double precision,intent(in) :: eGW(nBas)
double precision,intent(in) :: OmRPA(nS) double precision,intent(in) :: OmRPA(nS)
double precision,intent(in) :: rho_RPA(nBas,nBas,nS)
double precision,intent(in) :: OmBSE double precision,intent(in) :: OmBSE
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables ! Local variables
@ -60,8 +59,8 @@ subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O
do kc=1,maxS do kc=1,maxS
chi_A = chi_A + rho(i,j,kc)*rho(a,b,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2) chi_A = chi_A + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2)
chi_B = chi_B + rho(i,b,kc)*rho(a,j,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2) chi_B = chi_B + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2)
enddo enddo
@ -80,28 +79,28 @@ subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O
do kc=1,maxS do kc=1,maxS
eps_Ap = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) eps_Ap = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j))
chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*eps_Ap/(eps_Ap**2 + eta**2) chi_Ap = chi_Ap + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps_Ap/(eps_Ap**2 + eta**2)
eps_Ap = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) eps_Ap = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i))
chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*eps_Ap/(eps_Ap**2 + eta**2) chi_Ap = chi_Ap + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps_Ap/(eps_Ap**2 + eta**2)
eps_Am = - OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) eps_Am = - OmBSE - OmRPA(kc) - (eGW(a) - eGW(j))
chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*eps_Am/(eps_Am**2 + eta**2) chi_Am = chi_Am + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps_Am/(eps_Am**2 + eta**2)
eps_Am = - OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) eps_Am = - OmBSE - OmRPA(kc) - (eGW(b) - eGW(i))
chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*eps_Am/(eps_Am**2 + eta**2) chi_Am = chi_Am + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps_Am/(eps_Am**2 + eta**2)
eps_Bp = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)) eps_Bp = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(b))
chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*eps_Bp/(eps_Bp**2 + eta**2) chi_Bp = chi_Bp + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*eps_Bp/(eps_Bp**2 + eta**2)
eps_Bp = + OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)) eps_Bp = + OmBSE - OmRPA(kc) - (eGW(j) - eGW(i))
chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*eps_Bp/(eps_Bp**2 + eta**2) chi_Bp = chi_Bp + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*eps_Bp/(eps_Bp**2 + eta**2)
eps_Bm = - OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)) eps_Bm = - OmBSE - OmRPA(kc) - (eGW(a) - eGW(b))
chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*eps_Bm/(eps_Bm**2 + eta**2) chi_Bm = chi_Bm + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*eps_Bm/(eps_Bm**2 + eta**2)
eps_Bm = - OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)) eps_Bm = - OmBSE - OmRPA(kc) - (eGW(j) - eGW(i))
chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*eps_Bm/(eps_Bm**2 + eta**2) chi_Bm = chi_Bm + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*eps_Bm/(eps_Bm**2 + eta**2)
enddo enddo

10
src/QuAcK/Bethe_Salpeter_A_matrix_dynamic.f90
View File

@ -1,4 +1,4 @@
subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho,A_dyn) subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE,A_dyn)
! Compute the dynamic part of the Bethe-Salpeter equation matrices ! Compute the dynamic part of the Bethe-Salpeter equation matrices
@ -12,8 +12,8 @@ subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,Om
double precision,intent(in) :: lambda double precision,intent(in) :: lambda
double precision,intent(in) :: eGW(nBas) double precision,intent(in) :: eGW(nBas)
double precision,intent(in) :: OmRPA(nS) double precision,intent(in) :: OmRPA(nS)
double precision,intent(in) :: rho_RPA(nBas,nBas,nS)
double precision,intent(in) :: OmBSE double precision,intent(in) :: OmBSE
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables ! Local variables
@ -48,7 +48,7 @@ subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,Om
chi = 0d0 chi = 0d0
do kc=1,maxS do kc=1,maxS
chi = chi + rho(i,j,kc)*rho(a,b,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2) chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2)
enddo enddo
@ -58,10 +58,10 @@ subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,Om
do kc=1,maxS do kc=1,maxS
eps = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) eps = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j))
chi = chi + rho(i,j,kc)*rho(a,b,kc)*eps/(eps**2 + eta**2) chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps/(eps**2 + eta**2)
eps = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) eps = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i))
chi = chi + rho(i,j,kc)*rho(a,b,kc)*eps/(eps**2 + eta**2) chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps/(eps**2 + eta**2)
enddo enddo

10
src/QuAcK/Bethe_Salpeter_B_matrix.f90
View File

@ -1,4 +1,4 @@
subroutine Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,B_lr) subroutine Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,B_lr)
! Compute the extra term for Bethe-Salpeter equation for linear response ! Compute the extra term for Bethe-Salpeter equation for linear response
@ -11,8 +11,8 @@ subroutine Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,
double precision,intent(in) :: eta double precision,intent(in) :: eta
double precision,intent(in) :: lambda double precision,intent(in) :: lambda
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Omega(nS) double precision,intent(in) :: OmRPA(nS)
double precision,intent(in) :: rho(nBas,nBas,nS) double precision,intent(in) :: rho_RPA(nBas,nBas,nS)
! Local variables ! Local variables
@ -35,8 +35,8 @@ subroutine Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,
chi = 0d0 chi = 0d0
do kc=1,nS do kc=1,nS
eps = Omega(kc)**2 + eta**2 eps = OmRPA(kc)**2 + eta**2
chi = chi + rho(i,b,kc)*rho(a,j,kc)*Omega(kc)/eps chi = chi + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*OmRPA(kc)/eps
enddo enddo
B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI(i,j,b,a) + 4d0*lambda*chi B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI(i,j,b,a) + 4d0*lambda*chi

20
src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90
View File

@ -1,4 +1,4 @@
subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho, & subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE, &
ZAp,ZAm,ZBp,ZBm) ZAp,ZAm,ZBp,ZBm)
! Compute the dynamic part of the renormalization for the Bethe-Salpeter equation matrices ! Compute the dynamic part of the renormalization for the Bethe-Salpeter equation matrices
@ -13,8 +13,8 @@ subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,
double precision,intent(in) :: lambda double precision,intent(in) :: lambda
double precision,intent(in) :: eGW(nBas) double precision,intent(in) :: eGW(nBas)
double precision,intent(in) :: OmRPA(nS) double precision,intent(in) :: OmRPA(nS)
double precision,intent(in) :: rho_RPA(nBas,nBas,nS)
double precision,intent(in) :: OmBSE double precision,intent(in) :: OmBSE
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables ! Local variables
@ -63,28 +63,28 @@ subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,
do kc=1,maxS do kc=1,maxS
eps_Ap = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) eps_Ap = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j))
chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*(eps_Ap**2 - eta**2)/(eps_Ap**2 + eta**2)**2 chi_Ap = chi_Ap + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps_Ap**2 - eta**2)/(eps_Ap**2 + eta**2)**2
eps_Ap = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) eps_Ap = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i))
chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*(eps_Ap**2 - eta**2)/(eps_Ap**2 + eta**2)**2 chi_Ap = chi_Ap + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps_Ap**2 - eta**2)/(eps_Ap**2 + eta**2)**2
eps_Am = - OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) eps_Am = - OmBSE - OmRPA(kc) - (eGW(a) - eGW(j))
chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*(eps_Am**2 - eta**2)/(eps_Am**2 + eta**2)**2 chi_Am = chi_Am + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps_Am**2 - eta**2)/(eps_Am**2 + eta**2)**2
eps_Am = - OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) eps_Am = - OmBSE - OmRPA(kc) - (eGW(b) - eGW(i))
chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*(eps_Am**2 - eta**2)/(eps_Am**2 + eta**2)**2 chi_Am = chi_Am + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps_Am**2 - eta**2)/(eps_Am**2 + eta**2)**2
eps_Bp = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)) eps_Bp = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(b))
chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*(eps_Bp**2 - eta**2)/(eps_Bp**2 + eta**2)**2 chi_Bp = chi_Bp + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*(eps_Bp**2 - eta**2)/(eps_Bp**2 + eta**2)**2
eps_Bp = + OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)) eps_Bp = + OmBSE - OmRPA(kc) - (eGW(j) - eGW(i))
chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*(eps_Bp**2 - eta**2)/(eps_Bp**2 + eta**2)**2 chi_Bp = chi_Bp + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*(eps_Bp**2 - eta**2)/(eps_Bp**2 + eta**2)**2
eps_Bm = - OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)) eps_Bm = - OmBSE - OmRPA(kc) - (eGW(a) - eGW(b))
chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*(eps_Bm**2 - eta**2)/(eps_Bm**2 + eta**2)**2 chi_Bm = chi_Bm + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*(eps_Bm**2 - eta**2)/(eps_Bm**2 + eta**2)**2
eps_Bm = - OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)) eps_Bm = - OmBSE - OmRPA(kc) - (eGW(j) - eGW(i))
chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*(eps_Bm**2 - eta**2)/(eps_Bm**2 + eta**2)**2 chi_Bm = chi_Bm + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*(eps_Bm**2 - eta**2)/(eps_Bm**2 + eta**2)**2
enddo enddo

8
src/QuAcK/Bethe_Salpeter_ZA_matrix_dynamic.f90
View File

@ -1,4 +1,4 @@
subroutine Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho,ZA_dyn) subroutine Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE,ZA_dyn)
! Compute the dynamic part of the Bethe-Salpeter equation matrices ! Compute the dynamic part of the Bethe-Salpeter equation matrices
@ -12,8 +12,8 @@ subroutine Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O
double precision,intent(in) :: lambda double precision,intent(in) :: lambda
double precision,intent(in) :: eGW(nBas) double precision,intent(in) :: eGW(nBas)
double precision,intent(in) :: OmRPA(nS) double precision,intent(in) :: OmRPA(nS)
double precision,intent(in) :: rho_RPA(nBas,nBas,nS)
double precision,intent(in) :: OmBSE double precision,intent(in) :: OmBSE
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables ! Local variables
@ -49,10 +49,10 @@ subroutine Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O
do kc=1,maxS do kc=1,maxS
eps = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) eps = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j))
chi = chi + rho(i,j,kc)*rho(a,b,kc)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
eps = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) eps = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i))
chi = chi + rho(i,j,kc)*rho(a,b,kc)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
enddo enddo

12
src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90
View File

@ -1,4 +1,4 @@
subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,OmBSE,XpY,XmY,rho) subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,OmBSE,XpY,XmY)
! Compute dynamical effects via perturbation theory for BSE ! Compute dynamical effects via perturbation theory for BSE
@ -18,10 +18,10 @@ subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,
double precision,intent(in) :: eGW(nBas) double precision,intent(in) :: eGW(nBas)
double precision,intent(in) :: OmRPA(nS) double precision,intent(in) :: OmRPA(nS)
double precision,intent(in) :: rho_RPA(nBas,nBas,nS)
double precision,intent(in) :: OmBSE(nS) double precision,intent(in) :: OmBSE(nS)
double precision,intent(in) :: XpY(nS,nS) double precision,intent(in) :: XpY(nS,nS)
double precision,intent(in) :: XmY(nS,nS) double precision,intent(in) :: XmY(nS,nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables ! Local variables
@ -84,11 +84,11 @@ subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,
! Resonant part of the BSE correction for dynamical TDA ! Resonant part of the BSE correction for dynamical TDA
call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,OmBSE(ia),rho,Ap_dyn) call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia),Ap_dyn)
! Renormalization factor of the resonant parts for dynamical TDA ! Renormalization factor of the resonant parts for dynamical TDA
call Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,OmBSE(ia),rho,ZAp_dyn) call Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia),ZAp_dyn)
ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X))
OmDyn(ia) = dot_product(X,matmul( Ap_dyn,X)) OmDyn(ia) = dot_product(X,matmul( Ap_dyn,X))
@ -97,12 +97,12 @@ subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,
! Resonant and anti-resonant part of the BSE correction ! Resonant and anti-resonant part of the BSE correction
call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,OmBSE(ia),rho, & call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia), &
Ap_dyn,Am_dyn,Bp_dyn,Bm_dyn) Ap_dyn,Am_dyn,Bp_dyn,Bm_dyn)
! Renormalization factor of the resonant and anti-resonant parts ! Renormalization factor of the resonant and anti-resonant parts
call Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,OmBSE(ia),rho, & call Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia), &
ZAp_dyn,ZAm_dyn,ZBp_dyn,ZBm_dyn) ZAp_dyn,ZAm_dyn,ZBp_dyn,ZBm_dyn)
ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) & ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) &

11
src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90
View File

@ -1,4 +1,4 @@
subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,OmBSE,XpY,XmY,rho) subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,OmBSE,XpY,XmY)
! Compute self-consistently the dynamical effects via perturbation theory for BSE ! Compute self-consistently the dynamical effects via perturbation theory for BSE
@ -18,10 +18,10 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,
double precision,intent(in) :: eGW(nBas) double precision,intent(in) :: eGW(nBas)
double precision,intent(in) :: OmRPA(nS) double precision,intent(in) :: OmRPA(nS)
double precision,intent(in) :: rho_RPA(nBas,nBas,nS)
double precision,intent(in) :: OmBSE(nS) double precision,intent(in) :: OmBSE(nS)
double precision,intent(in) :: XpY(nS,nS) double precision,intent(in) :: XpY(nS,nS)
double precision,intent(in) :: XmY(nS,nS) double precision,intent(in) :: XmY(nS,nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables ! Local variables
@ -97,8 +97,7 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,
! Resonant part of the BSE correction ! Resonant part of the BSE correction
call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmOld(ia),rho(:,:,:), & call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmOld(ia),Ap_dyn)
Ap_dyn(:,:))
OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:)))
@ -106,8 +105,8 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,
! Anti-resonant part of the BSE correction ! Anti-resonant part of the BSE correction
call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmOld(ia),rho(:,:,:), & call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmOld(ia), &
Ap_dyn(:,:),Am_dyn(:,:),Bp_dyn(:,:),Bm_dyn(:,:)) Ap_dyn,Am_dyn,Bp_dyn,Bm_dyn)
OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) & OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) &
- dot_product(Y(:),matmul(Am_dyn(:,:),Y(:))) & - dot_product(Y(:),matmul(Am_dyn(:,:),Y(:))) &

121
src/QuAcK/G0W0.f90
View File

@ -1,5 +1,5 @@
subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linearize,eta, & dBSE,dTDA,evDyn,singlet,triplet,linearize,eta, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI,PHF,cHF,eHF,eGW) nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI,PHF,cHF,eHF,eGW)
! Perform G0W0 calculation ! Perform G0W0 calculation
@ -21,8 +21,8 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,
logical,intent(in) :: dBSE logical,intent(in) :: dBSE
logical,intent(in) :: dTDA logical,intent(in) :: dTDA
logical,intent(in) :: evDyn logical,intent(in) :: evDyn
logical,intent(in) :: singlet_manifold logical,intent(in) :: singlet
logical,intent(in) :: triplet_manifold logical,intent(in) :: triplet
logical,intent(in) :: linearize logical,intent(in) :: linearize
double precision,intent(in) :: eta double precision,intent(in) :: eta
@ -39,17 +39,16 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,
logical :: print_W = .true. logical :: print_W = .true.
integer :: ispin integer :: ispin
double precision :: EcRPA(nspin) double precision :: EcRPA
double precision :: EcBSE(nspin) double precision :: EcBSE(nspin)
double precision :: EcAC(nspin) double precision :: EcAC(nspin)
double precision :: EcGM double precision :: EcGM
double precision,allocatable :: SigC(:) double precision,allocatable :: SigC(:)
double precision,allocatable :: Z(:) double precision,allocatable :: Z(:)
double precision,allocatable :: Omega(:,:) double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY(:,:,:) double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY(:,:,:) double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho(:,:,:,:) double precision,allocatable :: rho_RPA(:,:,:)
double precision,allocatable :: rhox(:,:,:,:)
double precision,allocatable :: eGWlin(:) double precision,allocatable :: eGWlin(:)
@ -67,13 +66,13 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,
! Initialization ! Initialization
EcRPA(:) = 0d0 EcRPA = 0d0
! SOSEX correction ! SOSEX correction
if(SOSEX) then if(SOSEX) then
write(*,*) 'SOSEX correction activated!' write(*,*) 'SOSEX correction activated but BUG!'
write(*,*) stop
end if end if
! COHSEX approximation ! COHSEX approximation
@ -103,34 +102,43 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,
! Memory allocation ! Memory allocation
allocate(SigC(nBas),Z(nBas),Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin), & allocate(SigC(nBas),Z(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS),eGWlin(nBas))
rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin),eGWlin(nBas))
! Compute linear response !-------------------!
! Compute screening !
!-------------------!
call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, & call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) eHF,ERI,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
if(print_W) call print_excitation('RPA@HF ',ispin,nS,Omega(:,ispin)) if(print_W) call print_excitation('RPA@HF ',ispin,nS,OmRPA)
! Compute correlation part of the self-energy !--------------------------!
! Compute spectral weights !
!--------------------------!
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
if(SOSEX) call excitation_density_SOSEX(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rhox(:,:,:,ispin)) !------------------------!
! Compute GW self-energy !
!------------------------!
call self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC)
! Compute renormalization factor !--------------------------------!
! Compute renormalization factor !
!--------------------------------!
call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z(:))
! Solve the quasi-particle equation !-----------------------------------!
! Solve the quasi-particle equation !
!-----------------------------------!
eGWlin(:) = eHF(:) + Z(:)*SigC(:) eGWlin(:) = eHF(:) + Z(:)*SigC(:)
! Linearized or graphical solution?
if(linearize) then if(linearize) then
write(*,*) ' *** Quasiparticle energies obtained by linearization *** ' write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
@ -138,51 +146,48 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,
eGW(:) = eGWlin(:) eGW(:) = eGWlin(:)
! Find all the roots of the QP equation if necessary
! call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin)
else else
! Find graphical solution of the QP equation write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
write(*,*)
call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW) call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,OmRPA,rho_RPA,eGWlin,eGW)
! Find all the roots of the QP equation if necessary
! call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin)
end if end if
! Dump results
call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA(ispin),EcGM)
! Compute the RPA correlation energy ! Compute the RPA correlation energy
call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,OmRPA, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
write(*,*) !--------------!
write(*,*)'-------------------------------------------------------------------------------' ! Dump results !
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) call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------' ! Deallocate memory
write(*,*)
deallocate(SigC,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,eGWlin)
! Plot stuff ! Plot stuff
! call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin)) ! call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,OmRPA,rho_RPA)
! Perform BSE calculation ! Perform BSE calculation
if(BSE) then if(BSE) then
call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, & call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,EcBSE)
nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE)
if(exchange_kernel) then if(exchange_kernel) then
EcRPA(1) = 0.5d0*EcRPA(1) EcBSE(1) = 0.5d0*EcBSE(1)
EcRPA(2) = 1.5d0*EcRPA(1) EcBSE(2) = 1.5d0*EcBSE(2)
end if end if
@ -211,15 +216,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,
end if end if
call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet_manifold,triplet_manifold,eta, & call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,EcAC)
nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,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(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'

22
src/QuAcK/RPAx.f90
View File

@ -1,4 +1,4 @@
subroutine RPAx(doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,e) subroutine RPAx(doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
! Perform random phase approximation calculation with exchange (aka TDHF) ! Perform random phase approximation calculation with exchange (aka TDHF)
@ -21,7 +21,7 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,
integer,intent(in) :: nS integer,intent(in) :: nS
double precision,intent(in) :: ENuc double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF double precision,intent(in) :: ERHF
double precision,intent(in) :: e(nBas) double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables ! Local variables
@ -58,12 +58,11 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,
ispin = 1 ispin = 1
call linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & call linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Omega(:,ispin),rho, &
EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('RPAx ',ispin,nS,Omega(:,ispin)) call print_excitation('RPAx@HF ',ispin,nS,Omega(:,ispin))
call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
endif endif
! Triplet manifold ! Triplet manifold
@ -72,9 +71,9 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,
ispin = 2 ispin = 2
call linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & call linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,rho,Omega(:,ispin), &
EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('RPAx ',ispin,nS,Omega(:,ispin)) call print_excitation('RPAx@HF ',ispin,nS,Omega(:,ispin))
call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
endif endif
@ -105,14 +104,7 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,
write(*,*) write(*,*)
call ACFDT(exchange_kernel,.false.,.false.,.false.,.false.,.false.,singlet,triplet,eta, & call ACFDT(exchange_kernel,.false.,.false.,.false.,.false.,.false.,singlet,triplet,eta, &
nBas,nC,nO,nV,nR,nS,ERI,e,e,Omega,XpY,XmY,rho,EcAC) nBas,nC,nO,nV,nR,nS,ERI,eHF,eHF,EcAC)
if(exchange_kernel) then
EcAC(1) = 0.5d0*EcAC(1)
EcAC(2) = 1.5d0*EcAC(2)
end if
write(*,*) write(*,*)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'

34
src/QuAcK/dRPA.f90
View File

@ -1,5 +1,5 @@
subroutine dRPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & subroutine dRPA(doACFDT,exchange_kernel,singlet,triplet,eta, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,e) nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
! Perform a direct random phase approximation calculation ! Perform a direct random phase approximation calculation
@ -11,8 +11,8 @@ subroutine dRPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
logical,intent(in) :: doACFDT logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel logical,intent(in) :: exchange_kernel
logical,intent(in) :: singlet_manifold logical,intent(in) :: singlet
logical,intent(in) :: triplet_manifold logical,intent(in) :: triplet
double precision,intent(in) :: eta double precision,intent(in) :: eta
integer,intent(in) :: nBas integer,intent(in) :: nBas
integer,intent(in) :: nC integer,intent(in) :: nC
@ -22,7 +22,7 @@ subroutine dRPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
integer,intent(in) :: nS integer,intent(in) :: nS
double precision,intent(in) :: ENuc double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF double precision,intent(in) :: ERHF
double precision,intent(in) :: e(nBas) double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables ! Local variables
@ -40,7 +40,7 @@ subroutine dRPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
write(*,*) write(*,*)
write(*,*)'***********************************************' write(*,*)'***********************************************'
write(*,*)'| random-phase approximation calculation |' write(*,*)'| Random-phase approximation calculation |'
write(*,*)'***********************************************' write(*,*)'***********************************************'
write(*,*) write(*,*)
@ -55,32 +55,34 @@ subroutine dRPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
! Singlet manifold ! Singlet manifold
if(singlet_manifold) then if(singlet) then
ispin = 1 ispin = 1
call linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & call linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,rho,Omega(:,ispin), &
EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) call print_excitation('RPA@HF ',ispin,nS,Omega(:,ispin))
call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
endif endif
! Triplet manifold ! Triplet manifold
if(triplet_manifold) then if(triplet) then
ispin = 2 ispin = 2
call linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & call linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,rho,Omega(:,ispin), &
EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) call print_excitation('RPA@HF ',ispin,nS,Omega(:,ispin))
call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
endif endif
if(exchange_kernel) then if(exchange_kernel) then
EcRPA(1) = 0.5d0*EcRPA(1) EcRPA(1) = 0.5d0*EcRPA(1)
EcRPA(2) = 1.5d0*EcRPA(1) EcRPA(2) = 1.5d0*EcRPA(2)
end if end if
@ -103,13 +105,13 @@ subroutine dRPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
write(*,*) '------------------------------------------------------' write(*,*) '------------------------------------------------------'
write(*,*) write(*,*)
call ACFDT(exchange_kernel,.false.,.true.,.false.,.false.,.false.,singlet_manifold,triplet_manifold,eta, & call ACFDT(exchange_kernel,.false.,.true.,.false.,.false.,.false.,singlet,triplet,eta, &
nBas,nC,nO,nV,nR,nS,ERI,e,e,Omega,XpY,XmY,rho,EcAC) nBas,nC,nO,nV,nR,nS,ERI,eHF,eHF,EcAC)
if(exchange_kernel) then if(exchange_kernel) then
EcAC(1) = 0.5d0*EcAC(1) EcAC(1) = 0.5d0*EcAC(1)
EcAC(2) = 1.5d0*EcAC(1) EcAC(2) = 1.5d0*EcAC(2)
end if end if

119
src/QuAcK/evGW.f90
View File

@ -1,5 +1,6 @@
subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,G0W,GW0, & subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eG0W0) G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H, &
ERI,PHF,cHF,eHF,eG0W0)
! Perform self-consistent eigenvalue-only GW calculation ! Perform self-consistent eigenvalue-only GW calculation
@ -26,8 +27,8 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
logical,intent(in) :: evDyn logical,intent(in) :: evDyn
logical,intent(in) :: G0W logical,intent(in) :: G0W
logical,intent(in) :: GW0 logical,intent(in) :: GW0
logical,intent(in) :: singlet_manifold logical,intent(in) :: singlet
logical,intent(in) :: triplet_manifold logical,intent(in) :: triplet
double precision,intent(in) :: eta double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR,nS integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eHF(nBas)
@ -46,7 +47,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
integer :: n_diis integer :: n_diis
double precision :: rcond double precision :: rcond
double precision :: Conv double precision :: Conv
double precision :: EcRPA(nspin) double precision :: EcRPA
double precision :: EcBSE(nspin) double precision :: EcBSE(nspin)
double precision :: EcAC(nspin) double precision :: EcAC(nspin)
double precision :: EcGM double precision :: EcGM
@ -57,11 +58,10 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
double precision,allocatable :: eOld(:) double precision,allocatable :: eOld(:)
double precision,allocatable :: Z(:) double precision,allocatable :: Z(:)
double precision,allocatable :: SigC(:) double precision,allocatable :: SigC(:)
double precision,allocatable :: Omega(:,:) double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY(:,:,:) double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY(:,:,:) double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho(:,:,:,:) double precision,allocatable :: rho_RPA(:,:,:)
double precision,allocatable :: rhox(:,:,:,:)
! Hello world ! Hello world
@ -73,23 +73,45 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
! SOSEX correction ! SOSEX correction
if(SOSEX) write(*,*) 'SOSEX correction activated!' if(SOSEX) then
write(*,*) write(*,*) 'SOSEX correction activated but BUG!'
stop
end if
! COHSEX approximation ! COHSEX approximation
if(COHSEX) write(*,*) 'COHSEX approximation activated!' if(COHSEX) then
write(*,*) write(*,*) 'COHSEX approximation activated!'
write(*,*)
end if
! TDA for W ! TDA for W
if(TDA_W) write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' if(TDA_W) then
write(*,*) write(*,*) 'Tamm-Dancoff approximation for dynamic screening!'
write(*,*)
end if
! TDA ! TDA
if(TDA) write(*,*) 'Tamm-Dancoff approximation activated!' if(TDA) then
write(*,*) write(*,*) 'Tamm-Dancoff approximation activated!'
write(*,*)
end if
! GW0
if(GW0) then
write(*,*) 'GW0 scheme activated!'
write(*,*)
end if
! G0W
if(G0W) then
write(*,*) 'G0W scheme activated!'
write(*,*)
end if
! Linear mixing ! Linear mixing
@ -98,10 +120,8 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
! Memory allocation ! Memory allocation
allocate(eGW(nBas),eOld(nBas),Z(nBas),SigC(nBas),Omega(nS,nspin), & allocate(eGW(nBas),eOld(nBas),Z(nBas),SigC(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), &
XpY(nS,nS,nspin),XmY(nS,nS,nspin), & rho_RPA(nBas,nBas,nS),error_diis(nBas,max_diis),e_diis(nBas,max_diis))
rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin), &
error_diis(nBas,max_diis),e_diis(nBas,max_diis))
! Initialization ! Initialization
@ -121,36 +141,30 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
do while(Conv > thresh .and. nSCF <= maxSCF) do while(Conv > thresh .and. nSCF <= maxSCF)
! Compute linear response ! Compute screening
if(.not. GW0 .or. nSCF == 0) then if(.not. GW0 .or. nSCF == 0) then
call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,OmRPA, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
endif endif
! Compute spectral weights
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
! Compute correlation part of the self-energy ! Compute correlation part of the self-energy
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
if(SOSEX) call excitation_density_SOSEX(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rhox(:,:,:,ispin))
! Correlation self-energy
if(G0W) then if(G0W) then
call self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC) call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, &
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z(:))
else else
call self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eGW, & call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC)
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC) call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z)
call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eGW, &
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z(:))
endif endif
@ -164,8 +178,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_evGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM)
call print_evGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW)
! Linear mixing or DIIS extrapolation ! Linear mixing or DIIS extrapolation
@ -215,23 +228,22 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
endif endif
! Dump the RPA correlation energy ! Deallocate memory
write(*,*) deallocate(eOld,Z,SigC,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,error_diis,e_diis)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@RPA@evGW correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@evGW correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@evGW correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@evGW total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! Perform BSE calculation ! Perform BSE calculation
if(BSE) then if(BSE) then
call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, & call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcBSE)
nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE)
if(exchange_kernel) then
EcBSE(1) = 0.5d0*EcBSE(1)
EcBSE(2) = 1.5d0*EcBSE(2)
end if
write(*,*) write(*,*)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
@ -258,8 +270,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
end if end if
call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet_manifold,triplet_manifold,eta, & call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcAC)
nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,Omega,XpY,XmY,rho,EcAC)
write(*,*) write(*,*)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'

10
src/QuAcK/linear_response.f90
View File

@ -1,4 +1,4 @@
subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,rho,EcRPA,Omega,XpY,XmY) subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Omega_RPA,rho_RPA,EcRPA,Omega,XpY,XmY)
! Compute linear response ! Compute linear response
@ -12,9 +12,11 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E
integer,intent(in) :: ispin,nBas,nC,nO,nV,nR,nS integer,intent(in) :: ispin,nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: lambda double precision,intent(in) :: lambda
double precision,intent(in) :: e(nBas) double precision,intent(in) :: e(nBas)
double precision,intent(in) :: rho(nBas,nBas,nS)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Omega_RPA(nS)
double precision,intent(in) :: rho_RPA(nBas,nBas,nS)
! Local variables ! Local variables
integer :: ia integer :: ia
@ -42,7 +44,7 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E
call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A) call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A)
if(BSE) call Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,A) if(BSE) call Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega_RPA,rho_RPA,A)
! Tamm-Dancoff approximation ! Tamm-Dancoff approximation
@ -51,7 +53,7 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E
call linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B) call linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B)
if(BSE) call Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,B) if(BSE) call Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega_RPA,rho_RPA,B)
end if end if

4
src/QuAcK/plot_GW.f90
View File

@ -1,4 +1,4 @@
subroutine plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega,rho,rhox) subroutine plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega,rho)
! Dump several GW quantities for external plotting ! Dump several GW quantities for external plotting
@ -8,7 +8,7 @@ subroutine plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega,rho,rhox)
! Input variables ! Input variables
integer,intent(in) :: nBas,nC,nO,nV,nR,nS integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: eHF(nBas),eGW(nBas),Omega(nS),rho(nBas,nBas,nS),rhox(nBas,nBas,nS) double precision,intent(in) :: eHF(nBas),eGW(nBas),Omega(nS),rho(nBas,nBas,nS)
! Local variables ! Local variables

4
src/QuAcK/print_G0W0.f90
View File

@ -40,8 +40,8 @@ subroutine print_G0W0(nBas,nO,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA,EcGM)
write(*,'(2X,A30,F15.6)') 'G0W0 LUMO energy (eV):',eGW(LUMO)*HaToeV write(*,'(2X,A30,F15.6)') 'G0W0 LUMO energy (eV):',eGW(LUMO)*HaToeV
write(*,'(2X,A30,F15.6)') 'G0W0 HOMO-LUMO gap (eV):',Gap*HaToeV write(*,'(2X,A30,F15.6)') 'G0W0 HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A30,F15.6)') 'RPA@HF total energy =',ENuc + EHF + EcRPA write(*,'(2X,A30,F15.6)') 'RPA@G0W0 total energy =',ENuc + EHF + EcRPA
write(*,'(2X,A30,F15.6)') 'RPA@HF correlation energy =',EcRPA write(*,'(2X,A30,F15.6)') 'RPA@G0W0 correlation energy =',EcRPA
write(*,'(2X,A30,F15.6)') 'GM@G0W0 total energy =',ENuc + EHF + EcGM write(*,'(2X,A30,F15.6)') 'GM@G0W0 total energy =',ENuc + EHF + EcGM
write(*,'(2X,A30,F15.6)') 'GM@G0W0 correlation energy =',EcGM write(*,'(2X,A30,F15.6)') 'GM@G0W0 correlation energy =',EcGM
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'

26
src/QuAcK/print_evGW.f90
View File

@ -1,4 +1,4 @@
subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW) subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigC,Z,eGW,EcRPA,EcGM)
! Print one-electron energies and other stuff for evGW ! Print one-electron energies and other stuff for evGW
@ -8,9 +8,15 @@ subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW)
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) :: Conv,e(nBas),SigmaC(nBas),Z(nBas),eGW(nBas) double precision,intent(in) :: Conv
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: SigC(nBas)
double precision,intent(in) :: Z(nBas)
double precision,intent(in) :: eGW(nBas)
double precision,intent(in) :: EcRPA
double precision,intent(in) :: EcGM
integer :: x,HOMO,LUMO integer :: p,HOMO,LUMO
double precision :: Gap double precision :: Gap
! HOMO and LUMO ! HOMO and LUMO
@ -32,9 +38,9 @@ subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW)
'|','#','|','e_HF (eV)','|','Sigma_c (eV)','|','Z','|','e_QP (eV)','|' '|','#','|','e_HF (eV)','|','Sigma_c (eV)','|','Z','|','e_QP (eV)','|'
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
do x=1,nBas do p=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)') & 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,'|',e(x)*HaToeV,'|',SigmaC(x)*HaToeV,'|',Z(x),'|',eGW(x)*HaToeV,'|' '|',p,'|',e(p)*HaToeV,'|',SigC(p)*HaToeV,'|',Z(p),'|',eGW(p)*HaToeV,'|'
enddo enddo
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
@ -45,11 +51,11 @@ subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW)
write(*,'(2X,A30,F15.6)') 'evGW LUMO energy (eV):',eGW(LUMO)*HaToeV write(*,'(2X,A30,F15.6)') 'evGW LUMO energy (eV):',eGW(LUMO)*HaToeV
write(*,'(2X,A30,F15.6)') 'evGW HOMO-LUMO gap (eV):',Gap*HaToeV write(*,'(2X,A30,F15.6)') 'evGW HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
! write(*,'(2X,A30,F15.6)') 'RPA@evGW total energy =',ENuc + EHF + EcRPA write(*,'(2X,A30,F15.6)') 'RPA@evGW total energy =',ENuc + EHF + EcRPA
! write(*,'(2X,A30,F15.6)') 'RPA@evGW correlation energy =',EcRPA write(*,'(2X,A30,F15.6)') 'RPA@evGW correlation energy =',EcRPA
! write(*,'(2X,A30,F15.6)') 'GM@evGW total energy =',ENuc + EHF + EcGM write(*,'(2X,A30,F15.6)') 'GM@evGW total energy =',ENuc + EHF + EcGM
! write(*,'(2X,A30,F15.6)') 'GM@evGW correlation energy =',EcGM write(*,'(2X,A30,F15.6)') 'GM@evGW correlation energy =',EcGM
! write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,*) write(*,*)
end subroutine print_evGW end subroutine print_evGW

2
src/QuAcK/print_excitation.f90
View File

@ -26,7 +26,7 @@ subroutine print_excitation(method,ispin,nS,Omega)
write(*,*) write(*,*)
write(*,*)'-------------------------------------------------------------' write(*,*)'-------------------------------------------------------------'
write(*,'(1X,A1,1X,A14,A14,A14,A9,A13)')'|',method,' calculation: ',spin_manifold,' manifold',' |' write(*,'(1X,A14,A14,A14,A9)') method,' calculation: ',spin_manifold,' manifold'
write(*,*)'-------------------------------------------------------------' write(*,*)'-------------------------------------------------------------'
write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') & write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') &
'|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|' '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|'

5
src/QuAcK/print_qsGW.f90
View File

@ -70,8 +70,8 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,Sig
write(*,'(2X,A30,F15.6)') 'qsGW GM total energy =',EqsGW + ENuc + EcGM write(*,'(2X,A30,F15.6)') 'qsGW GM total energy =',EqsGW + ENuc + EcGM
write(*,'(2X,A30,F15.6)') 'qsGW exchange energy =',Ex write(*,'(2X,A30,F15.6)') 'qsGW exchange energy =',Ex
write(*,'(2X,A30,F15.6)') 'qsGW correlation energy =',Ec write(*,'(2X,A30,F15.6)') 'qsGW correlation energy =',Ec
! write(*,'(2X,A30,F15.6)') 'RPA@qsGW correlation energy =',EcRPA write(*,'(2X,A30,F15.6)') 'RPA@qsGW correlation energy =',EcRPA
! write(*,'(2X,A30,F15.6)') 'GM@qsGW correlation energy =',EcGM write(*,'(2X,A30,F15.6)') 'GM@qsGW correlation energy =',EcGM
write(*,*)'-------------------------------------------' write(*,*)'-------------------------------------------'
write(*,*) write(*,*)
@ -95,7 +95,6 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,Sig
write(*,'(A32,1X,F16.10)') ' Electronic energy ',EqsGW write(*,'(A32,1X,F16.10)') ' Electronic energy ',EqsGW
write(*,'(A32,1X,F16.10)') ' Nuclear repulsion ',ENuc write(*,'(A32,1X,F16.10)') ' Nuclear repulsion ',ENuc
write(*,'(A32,1X,F16.10)') ' qsGW energy ',ENuc + EqsGW write(*,'(A32,1X,F16.10)') ' qsGW energy ',ENuc + EqsGW
! write(*,'(A32,1X,F16.10)') ' RPA corr. energy ',EcRPA
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
write(*,*) write(*,*)

91
src/QuAcK/qsGW.f90
View File

@ -1,6 +1,6 @@
subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
COHSEX,SOSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,G0W,GW0,singlet_manifold,triplet_manifold,eta, & G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,S,X,T,V, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,S,X,T,V,Hc,ERI_AO_basis,ERI_MO_basis,PHF,cHF,eHF) Hc,ERI_AO_basis,ERI_MO_basis,PHF,cHF,eHF)
! Perform a quasiparticle self-consistent GW calculation ! Perform a quasiparticle self-consistent GW calculation
@ -25,8 +25,8 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
logical,intent(in) :: evDyn logical,intent(in) :: evDyn
logical,intent(in) :: G0W logical,intent(in) :: G0W
logical,intent(in) :: GW0 logical,intent(in) :: GW0
logical,intent(in) :: singlet_manifold logical,intent(in) :: singlet
logical,intent(in) :: triplet_manifold logical,intent(in) :: triplet
double precision,intent(in) :: eta double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR,nS integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: ENuc double precision,intent(in) :: ENuc
@ -49,7 +49,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
integer :: ispin integer :: ispin
integer :: n_diis integer :: n_diis
double precision :: EqsGW double precision :: EqsGW
double precision :: EcRPA(nspin) double precision :: EcRPA
double precision :: EcBSE(nspin) double precision :: EcBSE(nspin)
double precision :: EcAC(nspin) double precision :: EcAC(nspin)
double precision :: EcGM double precision :: EcGM
@ -58,11 +58,10 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
double precision,external :: trace_matrix double precision,external :: trace_matrix
double precision,allocatable :: error_diis(:,:) double precision,allocatable :: error_diis(:,:)
double precision,allocatable :: F_diis(:,:) double precision,allocatable :: F_diis(:,:)
double precision,allocatable :: Omega(:,:) double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY(:,:,:) double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY(:,:,:) double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho(:,:,:,:) double precision,allocatable :: rho_RPA(:,:,:)
double precision,allocatable :: rhox(:,:,:,:)
double precision,allocatable :: c(:,:) double precision,allocatable :: c(:,:)
double precision,allocatable :: cp(:,:) double precision,allocatable :: cp(:,:)
double precision,allocatable :: eGW(:) double precision,allocatable :: eGW(:)
@ -96,29 +95,37 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
! SOSEX correction ! SOSEX correction
if(SOSEX) write(*,*) 'SOSEX correction activated!' if(SOSEX) then
write(*,*) write(*,*) 'SOSEX correction activated but BUG!'
stop
end if
! COHSEX approximation ! COHSEX approximation
if(COHSEX) write(*,*) 'COHSEX approximation activated!' if(COHSEX) then
write(*,*) write(*,*) 'COHSEX approximation activated!'
write(*,*)
end if
! TDA for W ! TDA for W
if(TDA_W) write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' if(TDA_W) then
write(*,*) write(*,*) 'Tamm-Dancoff approximation for dynamic screening!'
write(*,*)
end if
! TDA ! TDA
if(TDA) write(*,*) 'Tamm-Dancoff approximation activated!' if(TDA) then
write(*,*) write(*,*) 'Tamm-Dancoff approximation activated!'
write(*,*)
end if
! Memory allocation ! Memory allocation
allocate(eGW(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & allocate(eGW(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), & J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),SigCm(nBas,nBas),Z(nBas), &
Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin),rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin), & OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), &
error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
! Initialization ! Initialization
@ -156,29 +163,23 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
if(.not. GW0 .or. nSCF == 0) then if(.not. GW0 .or. nSCF == 0) then
call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO_basis, & call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
endif endif
! Compute correlation part of the self-energy ! Compute correlation part of the self-energy
call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin)) call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY_RPA,rho_RPA)
if(SOSEX) call excitation_density_SOSEX(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rhox(:,:,:,ispin))
if(G0W) then if(G0W) then
call self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & call self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC) call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, &
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z)
else else
call self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eGW, & call self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC)
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC) call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z)
call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eGW, &
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z)
endif endif
@ -221,7 +222,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
! Print results ! Print results
! call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) ! call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,SigCp,Z,EcRPA(ispin),EcGM,EqsGW) call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,SigCp,Z,EcRPA,EcGM,EqsGW)
! Increment ! Increment
@ -258,23 +259,22 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
endif endif
! Dump RPA correlation energy ! Deallocate memory
write(*,*) deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,error,error_diis,F_diis)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@RPA@qsGW correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@qsGW correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@qsGW correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@qsGW total energy =',ENuc + EqsGW + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! Perform BSE calculation ! Perform BSE calculation
if(BSE) then if(BSE) then
call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, & call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,eGW,EcBSE)
nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE)
if(exchange_kernel) then
EcBSE(1) = 0.5d0*EcBSE(1)
EcBSE(2) = 1.5d0*EcBSE(2)
end if
write(*,*) write(*,*)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
@ -301,8 +301,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
end if end if
call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet_manifold,triplet_manifold,eta, & call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,eGW,EcAC)
nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,eGW,Omega,XpY,XmY,rho,EcAC)
write(*,*) write(*,*)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'

38
src/QuAcK/renormalization_factor.f90
View File

@ -1,4 +1,4 @@
subroutine renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,Z) subroutine renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Z)
! Compute renormalization factor for GW ! Compute renormalization factor for GW
@ -8,13 +8,11 @@ subroutine renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,r
! Input variables ! Input variables
logical,intent(in) :: COHSEX logical,intent(in) :: COHSEX
logical,intent(in) :: SOSEX
double precision,intent(in) :: eta double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR,nS integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: e(nBas) double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega(nS) double precision,intent(in) :: Omega(nS)
double precision,intent(in) :: rho(nBas,nBas,nS) double precision,intent(in) :: rho(nBas,nBas,nS)
double precision,intent(in) :: rhox(nBas,nBas,nS)
! Local variables ! Local variables
@ -36,40 +34,6 @@ subroutine renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,r
Z(:) = 1d0 Z(:) = 1d0
return return
! SOSEX correction
elseif(SOSEX) then
! Occupied part of the correlation self-energy
do x=nC+1,nBas-nR
do i=nC+1,nO
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
eps = e(x) - e(i) + Omega(jb)
Z(x) = Z(x) - (rho(x,i,jb)/eps)*(rhox(x,i,jb)/eps)
end do
end do
end do
end do
! Virtual part of the correlation self-energy
do x=nC+1,nBas-nR
do a=nO+1,nBas-nR
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
eps = e(x) - e(a) - Omega(jb)
Z(x) = Z(x) - (rho(x,a,jb)/eps)*(rhox(x,a,jb)/eps)
end do
end do
end do
end do
else else
! Occupied part of the correlation self-energy ! Occupied part of the correlation self-energy

42
src/QuAcK/self_energy_correlation.f90
View File

@ -1,4 +1,4 @@
subroutine self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,EcGM,SigC) subroutine self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC)
! Compute correlation part of the self-energy ! Compute correlation part of the self-energy
@ -8,13 +8,11 @@ subroutine self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,
! Input variables ! Input variables
logical,intent(in) :: COHSEX logical,intent(in) :: COHSEX
logical,intent(in) :: SOSEX
double precision,intent(in) :: eta double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR,nS integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: e(nBas) double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega(nS) double precision,intent(in) :: Omega(nS)
double precision,intent(in) :: rho(nBas,nBas,nS) double precision,intent(in) :: rho(nBas,nBas,nS)
double precision,intent(in) :: rhox(nBas,nBas,nS)
! Local variables ! Local variables
@ -30,7 +28,9 @@ subroutine self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,
SigC = 0d0 SigC = 0d0
! COHSEX static approximation !-----------------------------!
! COHSEX static approximation !
!-----------------------------!
if(COHSEX) then if(COHSEX) then
@ -65,6 +65,10 @@ subroutine self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,
else else
!----------------!
! GW self-energy !
!----------------!
! Occupied part of the correlation self-energy ! Occupied part of the correlation self-energy
do x=nC+1,nBas-nR do x=nC+1,nBas-nR
@ -91,36 +95,6 @@ subroutine self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,
enddo enddo
enddo enddo
if(SOSEX) then
! SOSEX: occupied part of the correlation self-energy
do x=nC+1,nBas-nR
do y=nC+1,nBas-nR
do i=nC+1,nO
do jb=1,nS
eps = e(x) - e(i) + Omega(jb)
SigC(x,y) = SigC(x,y) - rho(x,i,jb)*rhox(y,i,jb)*eps/(eps**2 + eta**2)
enddo
enddo
enddo
enddo
! SOSEX: virtual part of the correlation self-energy
do x=nC+1,nBas-nR
do y=nC+1,nBas-nR
do a=nO+1,nBas-nR
do jb=1,nS
eps = e(x) - e(a) - Omega(jb)
SigC(x,y) = SigC(x,y) - rho(x,a,jb)*rhox(y,a,jb)*eps/(eps**2 + eta**2)
enddo
enddo
enddo
enddo
endif
endif endif
end subroutine self_energy_correlation end subroutine self_energy_correlation

43
src/QuAcK/self_energy_correlation_diag.f90
View File

@ -1,4 +1,4 @@
subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,EcGM,SigC) subroutine self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC)
! Compute diagonal of the correlation part of the self-energy ! Compute diagonal of the correlation part of the self-energy
@ -8,7 +8,6 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,O
! Input variables ! Input variables
logical,intent(in) :: COHSEX logical,intent(in) :: COHSEX
logical,intent(in) :: SOSEX
double precision,intent(in) :: eta double precision,intent(in) :: eta
integer,intent(in) :: nBas integer,intent(in) :: nBas
integer,intent(in) :: nC integer,intent(in) :: nC
@ -19,7 +18,6 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,O
double precision,intent(in) :: e(nBas) double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega(nS) double precision,intent(in) :: Omega(nS)
double precision,intent(in) :: rho(nBas,nBas,nS) double precision,intent(in) :: rho(nBas,nBas,nS)
double precision,intent(in) :: rhox(nBas,nBas,nS)
! Local variables ! Local variables
@ -69,45 +67,6 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,O
EcGM = EcGM - SigC(i) EcGM = EcGM - SigC(i)
end do end do
!-----------------------------
! SOSEX self-energy *BUG*
!-----------------------------
elseif(SOSEX) then
! SOSEX: 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)*rhox(p,i,jb)*eps/(eps**2 + eta**2)
end do
end do
end do
! SOSEX: 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)*rhox(p,a,jb)*eps/(eps**2 + eta**2)
end do
end do
end do
! GM correlation energy
do i=nC+1,nO
do a=nO+1,nBas-nR
do jb=1,nS
eps = e(a) - e(i) + Omega(jb)
EcGM = EcGM + rho(a,i,jb)*rhox(a,i,jb)*eps/(eps**2 + eta**2)
end do
end do
end do
!----------------------------- !-----------------------------
! GW self-energy ! GW self-energy
!----------------------------- !-----------------------------