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

dynamic screening in BSE

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-04-23 23:13:15 +02:00
parent e855727d8a
commit 01e9f785c2
15 changed files with 141 additions and 245 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
examples/molecule.CO
View File

@ -1,5 +1,5 @@
# nAt nEla nElb nCore nRyd # nAt nEla nElb nCore nRyd
2 7 7 0 0 2 7 7 0 0
# Znuc x y z # Znuc x y z
C 0. 0. 0. C 0. 0. -1.24942055
O 0. 0. 2.8 O 0. 0. 0.89266692

6
examples/molecule.H2O
View File

@ -1,6 +1,6 @@
# nAt nEla nElb nCore nRyd # nAt nEla nElb nCore nRyd
3 5 5 0 0 3 5 5 0 0
# Znuc x y z # Znuc x y z
O. 0. 0. 0. O 0. 0. -0.13209669
H -1.430428706 0. -1.107156965 H 0. 1.43152878 0.97970006
H 1.430428706 0. -1.107156965 H 0. -1.43152878 0.97970006

4
examples/molecule.HCl
View File

@ -1,5 +1,5 @@
# nAt nEla nElb nCore nRyd # nAt nEla nElb nCore nRyd
2 9 9 0 0 2 9 9 0 0
# Znuc x y z # Znuc x y z
H 0. 0. 0. H 0. 0. -0.02489783
Cl 0. 0. 3.2 Cl 0. 0. 2.38483140

4
examples/molecule.N2
View File

@ -1,5 +1,5 @@
# nAt nEla nElb nCore nRyd # nAt nEla nElb nCore nRyd
2 7 7 0 0 2 7 7 0 0
# Znuc x y z # Znuc x y z
N 0. 0. 0. N 0. 0. -1.04008632
N 0. 0. 2.5 N 0. 0. +1.04008632

86
input/basis
View File

@ -1,83 +1,9 @@
1 5 1 3
S 3 S 3
1 13.0100000 0.0196850 1 38.3600000 0.0238090
2 1.9620000 0.1379770 2 5.7700000 0.1548910
3 0.4446000 0.4781480 3 1.2400000 0.4699870
S 1 S 1
1 0.1220000 1.0000000 1 0.2976000 1.0000000
S 1
1 0.0297400 1.0000000
P 1 P 1
1 0.7270000 1.0000000 1 1.2750000 1.0000000
P 1
1 0.1410000 1.0000000
2 9
S 8
1 9046.0000000 0.0007000
2 1357.0000000 0.0053890
3 309.3000000 0.0274060
4 87.7300000 0.1032070
5 28.5600000 0.2787230
6 10.2100000 0.4485400
7 3.8380000 0.2782380
8 0.7466000 0.0154400
S 8
1 9046.0000000 -0.0001530
2 1357.0000000 -0.0012080
3 309.3000000 -0.0059920
4 87.7300000 -0.0245440
5 28.5600000 -0.0674590
6 10.2100000 -0.1580780
7 3.8380000 -0.1218310
8 0.7466000 0.5490030
S 1
1 0.2248000 1.0000000
S 1
1 0.0612400 1.0000000
P 3
1 13.5500000 0.0399190
2 2.9170000 0.2171690
3 0.7973000 0.5103190
P 1
1 0.2185000 1.0000000
P 1
1 0.0561100 1.0000000
D 1
1 0.8170000 1.0000000
D 1
1 0.2300000 1.0000000
3 9
S 8
1 11720.0000000 0.0007100
2 1759.0000000 0.0054700
3 400.8000000 0.0278370
4 113.7000000 0.1048000
5 37.0300000 0.2830620
6 13.2700000 0.4487190
7 5.0250000 0.2709520
8 1.0130000 0.0154580
S 8
1 11720.0000000 -0.0001600
2 1759.0000000 -0.0012630
3 400.8000000 -0.0062670
4 113.7000000 -0.0257160
5 37.0300000 -0.0709240
6 13.2700000 -0.1654110
7 5.0250000 -0.1169550
8 1.0130000 0.5573680
S 1
1 0.3023000 1.0000000
S 1
1 0.0789600 1.0000000
P 3
1 17.7000000 0.0430180
2 3.8540000 0.2289130
3 1.0460000 0.5087280
P 1
1 0.2753000 1.0000000
P 1
1 0.0685600 1.0000000
D 1
1 1.1850000 1.0000000
D 1
1 0.3320000 1.0000000

6
input/molecule
View File

@ -1,6 +1,4 @@
# nAt nEla nElb nCore nRyd # nAt nEla nElb nCore nRyd
3 8 8 0 0 1 1 1 0 0
# Znuc x y z # Znuc x y z
H 1.18163475 0.00000000 -1.17386890 He 0.0 0.0 0.0
N -0.44776863 0.00000000 -0.03589263
O 0.21099695 0.00000000 2.15462460

6
input/molecule.xyz
View File

@ -1,5 +1,3 @@
3 1
H 0.6252942263 0.0000000000 -0.6211847152 He 0.0000000000 0.0000000000 0.0000000000
N -0.2369489718 0.0000000000 -0.0189935632
O 0.1116547855 0.0000000000 1.1401783185

2
input/options
View File

@ -9,7 +9,7 @@
# GF: maxSCF thresh DIIS n_diis lin renorm # GF: maxSCF thresh DIIS n_diis lin renorm
256 0.00001 T 5 T 3 256 0.00001 T 5 T 3
# GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta # GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta
256 0.00001 T 5 F F T F F F T 0.000 256 0.00001 T 5 F F T F F F T 0.0
# ACFDT: AC Kx XBS # ACFDT: AC Kx XBS
F F T F F T
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift

86
input/weight
View File

@ -1,83 +1,9 @@
1 5 1 3
S 3 S 3
1 13.0100000 0.0196850 1 38.3600000 0.0238090
2 1.9620000 0.1379770 2 5.7700000 0.1548910
3 0.4446000 0.4781480 3 1.2400000 0.4699870
S 1 S 1
1 0.1220000 1.0000000 1 0.2976000 1.0000000
S 1
1 0.0297400 1.0000000
P 1 P 1
1 0.7270000 1.0000000 1 1.2750000 1.0000000
P 1
1 0.1410000 1.0000000
2 9
S 8
1 9046.0000000 0.0007000
2 1357.0000000 0.0053890
3 309.3000000 0.0274060
4 87.7300000 0.1032070
5 28.5600000 0.2787230
6 10.2100000 0.4485400
7 3.8380000 0.2782380
8 0.7466000 0.0154400
S 8
1 9046.0000000 -0.0001530
2 1357.0000000 -0.0012080
3 309.3000000 -0.0059920
4 87.7300000 -0.0245440
5 28.5600000 -0.0674590
6 10.2100000 -0.1580780
7 3.8380000 -0.1218310
8 0.7466000 0.5490030
S 1
1 0.2248000 1.0000000
S 1
1 0.0612400 1.0000000
P 3
1 13.5500000 0.0399190
2 2.9170000 0.2171690
3 0.7973000 0.5103190
P 1
1 0.2185000 1.0000000
P 1
1 0.0561100 1.0000000
D 1
1 0.8170000 1.0000000
D 1
1 0.2300000 1.0000000
3 9
S 8
1 11720.0000000 0.0007100
2 1759.0000000 0.0054700
3 400.8000000 0.0278370
4 113.7000000 0.1048000
5 37.0300000 0.2830620
6 13.2700000 0.4487190
7 5.0250000 0.2709520
8 1.0130000 0.0154580
S 8
1 11720.0000000 -0.0001600
2 1759.0000000 -0.0012630
3 400.8000000 -0.0062670
4 113.7000000 -0.0257160
5 37.0300000 -0.0709240
6 13.2700000 -0.1654110
7 5.0250000 -0.1169550
8 1.0130000 0.5573680
S 1
1 0.3023000 1.0000000
S 1
1 0.0789600 1.0000000
P 3
1 17.7000000 0.0430180
2 3.8540000 0.2289130
3 1.0460000 0.5087280
P 1
1 0.2753000 1.0000000
P 1
1 0.0685600 1.0000000
D 1
1 1.1850000 1.0000000
D 1
1 0.3320000 1.0000000

45
src/QuAcK/Bethe_Salpeter.f90
View File

@ -1,5 +1,5 @@
subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, & subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eW,eGW,OmRPA,XpY,XmY,rho,EcRPA,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
@ -19,14 +19,17 @@ subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, &
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 :: OmRPA(nS,nspin)
double precision :: XpY(nS,nS,nspin) double precision :: XpY_RPA(nS,nS,nspin)
double precision :: XmY(nS,nS,nspin) double precision :: XmY_RPA(nS,nS,nspin)
double precision :: rho(nBas,nBas,nS,nspin) double precision :: rho_RPA(nBas,nBas,nS,nspin)
! Local variables ! Local variables
integer :: ispin integer :: ispin
double precision,allocatable :: OmBSE(:,:) double precision,allocatable :: OmBSE(:,:)
double precision,allocatable :: XpY_BSE(:,:,:)
double precision,allocatable :: XmY_BSE(:,:,:)
double precision,allocatable :: rho_BSE(:,:,:,:)
! Output variables ! Output variables
@ -35,7 +38,7 @@ subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, &
! Memory allocation ! Memory allocation
allocate(OmBSE(nS,nspin)) allocate(OmBSE(nS,nspin),XpY_BSE(nS,nS,nspin),XmY_BSE(nS,nS,nspin),rho_BSE(nBas,nBas,nS,nspin))
!------------------- !-------------------
! Singlet manifold ! Singlet manifold
@ -46,19 +49,25 @@ subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, &
ispin = 1 ispin = 1
EcBSE(ispin) = 0d0 EcBSE(ispin) = 0d0
! Compute RPA screening
call linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & call linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),OmRPA(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) rho_RPA(:,:,:,ispin),EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA(:,:,ispin),rho_RPA(:,:,:,ispin))
! Compute BSE excitation energies
OmBSE(:,ispin) = OmRPA(:,ispin) OmBSE(:,ispin) = OmRPA(:,ispin)
call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
rho(:,:,:,ispin),EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) rho_RPA(:,:,:,ispin),EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
call print_excitation('BSE ',ispin,nS,OmBSE(:,ispin)) call print_excitation('BSE ',ispin,nS,OmBSE(:,ispin))
! call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_BSE(:,:,ispin),rho_BSE(:,:,:,ispin))
! Compute dynamic correction for BSE via perturbation theory ! Compute dynamic correction for BSE via perturbation theory
call Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,OmRPA(:,ispin),OmBSE(:,ispin), & call Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), &
XpY(:,:,ispin),XmY(:,:,ispin),rho(:,:,:,ispin)) XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin))
end if end if
@ -71,19 +80,25 @@ subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, &
ispin = 2 ispin = 2
EcBSE(ispin) = 0d0 EcBSE(ispin) = 0d0
! Compute RPA screening
call linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & call linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),OmRPA(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) rho_RPA(:,:,:,ispin),EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA(:,:,ispin),rho_RPA(:,:,:,ispin))
! Compute BSE excitation energies
OmBSE(:,ispin) = OmRPA(:,ispin) OmBSE(:,ispin) = OmRPA(:,ispin)
call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
rho(:,:,:,ispin),EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) rho_RPA(:,:,:,ispin),EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
call print_excitation('BSE ',ispin,nS,OmBSE(:,ispin)) call print_excitation('BSE ',ispin,nS,OmBSE(:,ispin))
! call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_BSE(:,:,ispin),rho_BSE(:,:,:,ispin))
! Compute dynamic correction for BSE via perturbation theory ! Compute dynamic correction for BSE via perturbation theory
call Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,OmRPA(:,ispin),OmBSE(:,ispin), & call Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), &
XpY(:,:,ispin),XmY(:,:,ispin),rho(:,:,:,ispin)) XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin))
end if end if

3
src/QuAcK/Bethe_Salpeter_A_matrix.f90
View File

@ -46,7 +46,4 @@ subroutine Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,
enddo enddo
enddo enddo
! print*,'BSE A'
! call matout(nS,nS,A_lr)
end subroutine Bethe_Salpeter_A_matrix end subroutine Bethe_Salpeter_A_matrix

32
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,OmRPA,OmBSE,rho,A_lr) subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho,A_dyn)
! Compute the dynamic part of the Bethe-Salpeter equation matrices ! Compute the dynamic part of the Bethe-Salpeter equation matrices
@ -10,21 +10,31 @@ subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,OmRPA,
integer,intent(in) :: nBas,nC,nO,nV,nR,nS integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: eta double precision,intent(in) :: eta
double precision,intent(in) :: lambda double precision,intent(in) :: lambda
double precision,intent(in) :: eGW(nBas)
double precision,intent(in) :: OmRPA(nS) double precision,intent(in) :: OmRPA(nS)
double precision,intent(in) :: OmBSE double precision,intent(in) :: OmBSE
double precision,intent(in) :: rho(nBas,nBas,nS) double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables ! Local variables
integer :: maxS
double precision :: chi double precision :: chi
double precision :: eps double precision :: eps
integer :: i,j,a,b,ia,jb,kc integer :: i,j,a,b,ia,jb,kc
! Output variables ! Output variables
double precision,intent(out) :: A_lr(nS,nS) double precision,intent(out) :: A_dyn(nS,nS)
A_lr(:,:) = 0d0 ! Initialization
A_dyn(:,:) = 0d0
! Number of poles taken into account
maxS = nS
! Build dynamic A matrix
ia = 0 ia = 0
do i=nC+1,nO do i=nC+1,nO
@ -36,27 +46,27 @@ subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,OmRPA,
jb = jb + 1 jb = jb + 1
chi = 0d0 chi = 0d0
do kc=1,nS do kc=1,maxS
eps = OmRPA(kc)**2 + eta**2 eps = OmRPA(kc)**2 + eta**2
chi = chi + rho(i,j,kc)*rho(a,b,kc)*OmRPA(kc)/eps chi = chi + rho(i,j,kc)*rho(a,b,kc)*OmRPA(kc)/eps
enddo enddo
A_lr(ia,jb) = A_lr(ia,jb) - 4d0*lambda*chi A_dyn(ia,jb) = A_dyn(ia,jb) - 4d0*lambda*chi
chi = 0d0 chi = 0d0
do kc=1,nS do kc=1,maxS
eps = (OmBSE - OmRPA(kc))**2 + eta**2 eps = (OmBSE - OmRPA(kc) - (eGW(a) - eGW(i)))**2 + eta**2
chi = chi + rho(i,j,kc)*rho(a,b,kc)*(OmBSE - OmRPA(kc))/eps chi = chi + rho(i,j,kc)*rho(a,b,kc)*(OmBSE - OmRPA(kc) - (eGW(a) - eGW(i)))/eps
eps = (OmBSE + OmRPA(kc))**2 + eta**2 eps = (OmBSE - OmRPA(kc) - (eGW(b) - eGW(j)))**2 + eta**2
chi = chi - rho(i,j,kc)*rho(a,b,kc)*(OmBSE + OmRPA(kc))/eps chi = chi + rho(i,j,kc)*rho(a,b,kc)*(OmBSE - OmRPA(kc) - (eGW(b) - eGW(j)))/eps
enddo enddo
A_lr(ia,jb) = A_lr(ia,jb) - 2d0*lambda*chi A_dyn(ia,jb) = A_dyn(ia,jb) - 2d0*lambda*chi
enddo enddo
enddo enddo

3
src/QuAcK/Bethe_Salpeter_B_matrix.f90
View File

@ -46,7 +46,4 @@ subroutine Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,
enddo enddo
enddo enddo
! print*,'BSE B'
! call matout(nS,nS,B_lr)
end subroutine Bethe_Salpeter_B_matrix end subroutine Bethe_Salpeter_B_matrix

33
src/QuAcK/Bethe_Salpeter_B_matrix_dynamic.f90
View File

@ -1,4 +1,4 @@
subroutine Bethe_Salpeter_B_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,OmRPA,OmBSE,rho,B_lr) subroutine Bethe_Salpeter_B_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho,B_dyn)
! Compute the dynamic part of the Bethe-Salpeter equation matrices ! Compute the dynamic part of the Bethe-Salpeter equation matrices
@ -10,21 +10,31 @@ subroutine Bethe_Salpeter_B_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,OmRPA,
integer,intent(in) :: nBas,nC,nO,nV,nR,nS integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: eta double precision,intent(in) :: eta
double precision,intent(in) :: lambda double precision,intent(in) :: lambda
double precision,intent(in) :: eGW(nBas)
double precision,intent(in) :: OmRPA(nS) double precision,intent(in) :: OmRPA(nS)
double precision,intent(in) :: OmBSE double precision,intent(in) :: OmBSE
double precision,intent(in) :: rho(nBas,nBas,nS) double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables ! Local variables
integer :: maxS
double precision :: chi double precision :: chi
double precision :: eps double precision :: eps
integer :: i,j,a,b,ia,jb,kc integer :: i,j,a,b,ia,jb,kc
! Output variables ! Output variables
double precision,intent(out) :: B_lr(nS,nS) double precision,intent(out) :: B_dyn(nS,nS)
B_lr(:,:) = 0d0 ! Initialization
B_dyn(:,:) = 0d0
! Number of poles taken into account
maxS = nS
! Build dynamic A matrix
ia = 0 ia = 0
do i=nC+1,nO do i=nC+1,nO
@ -36,28 +46,27 @@ subroutine Bethe_Salpeter_B_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,OmRPA,
jb = jb + 1 jb = jb + 1
chi = 0d0 chi = 0d0
do kc=1,nS do kc=1,maxS
eps = OmRPA(kc)**2 + eta**2 eps = OmRPA(kc)**2 + eta**2
chi = chi + rho(i,b,kc)*rho(a,j,kc)*OmRPA(kc)/eps chi = chi + rho(i,b,kc)*rho(a,j,kc)*OmRPA(kc)/eps
enddo enddo
B_lr(ia,jb) = B_lr(ia,jb) - 4d0*lambda*chi B_dyn(ia,jb) = B_dyn(ia,jb) - 4d0*lambda*chi
chi = 0d0 chi = 0d0
do kc=1,nS do kc=1,maxS
eps = (OmBSE - OmRPA(kc))**2 + eta**2 eps = (OmBSE - OmRPA(kc) - (eGW(a) - eGW(i)))**2 + eta**2
chi = chi + rho(i,b,kc)*rho(a,j,kc)*(OmBSE - OmRPA(kc))/eps chi = chi + rho(i,b,kc)*rho(a,j,kc)*(OmBSE - OmRPA(kc) - (eGW(a) - eGW(i)))/eps
eps = (OmBSE + OmRPA(kc))**2 + eta**2
chi = chi - rho(i,b,kc)*rho(a,j,kc)*(OmBSE + OmRPA(kc))/eps
eps = (OmBSE - OmRPA(kc) + (eGW(b) - eGW(j)))**2 + eta**2
chi = chi + rho(i,b,kc)*rho(a,j,kc)*(OmBSE - OmRPA(kc) + (eGW(b) - eGW(j)))/eps
enddo enddo
B_lr(ia,jb) = B_lr(ia,jb) - 2d0*lambda*chi B_dyn(ia,jb) = B_dyn(ia,jb) - 2d0*lambda*chi
enddo enddo
enddo enddo

66
src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90
View File

@ -1,4 +1,4 @@
subroutine Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,OmRPA,OmBSE,XpY,XmY,rho) subroutine Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,OmBSE,XpY,XmY,rho)
! Compute dynamical effects via perturbation theory for BSE ! Compute dynamical effects via perturbation theory for BSE
@ -9,13 +9,19 @@ subroutine Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,OmRPA
logical,intent(in) :: TDA logical,intent(in) :: TDA
double precision,intent(in) :: eta double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR,nS 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 :: OmRPA(nS) double precision,intent(in) :: eGW(nBas)
double precision :: OmBSE(nS) double precision,intent(in) :: OmRPA(nS)
double precision :: XpY(nS,nS) double precision,intent(in) :: OmBSE(nS)
double precision :: XmY(nS,nS) double precision,intent(in) :: XpY(nS,nS)
double precision :: rho(nBas,nBas,nS) double precision,intent(in) :: XmY(nS,nS)
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables ! Local variables
@ -23,44 +29,58 @@ subroutine Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,OmRPA
integer,parameter :: maxS = 10 integer,parameter :: maxS = 10
double precision,allocatable :: OmDyn(:) double precision,allocatable :: OmDyn(:)
double precision,allocatable :: ZDyn(:)
double precision,allocatable :: X(:) double precision,allocatable :: X(:)
double precision,allocatable :: Y(:) double precision,allocatable :: Y(:)
double precision,allocatable :: A_dyn(:,:) double precision,allocatable :: A_dyn(:,:)
double precision,allocatable :: B_dyn(:,:) double precision,allocatable :: B_dyn(:,:)
double precision,allocatable :: Z_dyn(:,:)
! Memory allocation ! Memory allocation
allocate(OmDyn(nS),X(nS),Y(nS),A_dyn(nS,nS),B_dyn(nS,nS)) allocate(OmDyn(nS),ZDyn(nS),X(nS),Y(nS),A_dyn(nS,nS),Z_dyn(nS,nS))
write(*,*) '-----------------------------------------------------------------------------------------------------------' write(*,*) '---------------------------------------------------------------------------------------------------'
write(*,'(2X,A5,1X,A30,1X,A30,1X,A30)') '#','Static excitation (eV)','Dynamic correction (eV)','Dynamic excitation (eV)' write(*,*) ' First-order dynamical correction to static Bethe-Salpeter excitation energies '
write(*,*) '-----------------------------------------------------------------------------------------------------------' write(*,*) '---------------------------------------------------------------------------------------------------'
write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)'
write(*,*) '---------------------------------------------------------------------------------------------------'
do ia=1,min(nS,maxS) do ia=1,min(nS,maxS)
X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:)) X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:))
Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:)) Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:))
call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,OmRPA(:),OmBSE(ia),rho(:,:,:),A_dyn(:,:)) call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:),A_dyn(:,:))
call Bethe_Salpeter_Z_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:),Z_dyn(:,:))
if(TDA) then ! First-order correction
OmDyn(ia) = dot_product(X(:),matmul(A_dyn(:,:),X(:))) if(.true.) then
else ZDyn(ia) = dot_product(X(:),matmul(Z_dyn(:,:),X(:)))
call Bethe_Salpeter_B_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,OmRPA(:),OmBSE(ia),rho(:,:,:),B_dyn(:,:)) else
OmDyn(ia) = dot_product(X(:),matmul(A_dyn(:,:),X(:))) & allocate(B_dyn(nS,nS))
- dot_product(Y(:),matmul(A_dyn(:,:),Y(:))) & call Bethe_Salpeter_B_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:),B_dyn(:,:))
+ dot_product(X(:),matmul(B_dyn(:,:),Y(:))) &
- dot_product(Y(:),matmul(B_dyn(:,:),X(:)))
end if OmDyn(ia) = dot_product(X(:),matmul(A_dyn(:,:),X(:))) &
- dot_product(Y(:),matmul(A_dyn(:,:),Y(:))) &
+ dot_product(X(:),matmul(B_dyn(:,:),Y(:))) &
- dot_product(Y(:),matmul(B_dyn(:,:),X(:)))
write(*,'(2X,I5,15X,F15.6,15X,F15.6,15X,F15.6)') ia,OmBSE(ia)*HaToeV,OmDyn(ia)*HaToeV,(OmBSE(ia)+OmDyn(ia))*HaToeV end if
! Renormalization factor
ZDyn(ia) = 1d0/(1d0 - ZDyn(ia))
OmDyn(ia) = ZDyn(ia)*dot_product(X(:),matmul(A_dyn(:,:),X(:)))
write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6,5X,F15.6)') &
ia,OmBSE(ia)*HaToeV,(OmBSE(ia)+OmDyn(ia))*HaToeV,OmDyn(ia)*HaToeV,ZDyn(ia)
end do end do
write(*,*) '-----------------------------------------------------------------------------------------------------------' write(*,*) '---------------------------------------------------------------------------------------------------'
write(*,*) write(*,*)
end subroutine Bethe_Salpeter_dynamic_perturbation end subroutine Bethe_Salpeter_dynamic_perturbation