Skip to content

Commit b177a3d

Browse files
author
lburth
committed
static BSE on MOM ref
1 parent bb3617e commit b177a3d

File tree

6 files changed

+592
-56
lines changed

6 files changed

+592
-56
lines changed

‎src/GW/CVS_UG0W0.f90‎

Lines changed: 18 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -236,21 +236,23 @@ subroutine CVS_UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE
236236
deallocate(Om,XpY,XmY,rho)
237237

238238
! Perform BSE calculation
239-
!
240-
! if(dophBSE) then
241-
!
242-
! call UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,S, &
243-
! ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF,eGW,EcBSE)
244-
!
245-
! write(*,*)
246-
! write(*,*)'-------------------------------------------------------------------------------'
247-
! write(*,'(2X,A60,F20.10,A3)') 'Tr@BSE@G0W0@UHF correlation energy (spin-conserved) = ',EcBSE(1),' au'
248-
! write(*,'(2X,A60,F20.10,A3)') 'Tr@BSE@G0W0@UHF correlation energy (spin-flip) = ',EcBSE(2),' au'
249-
! write(*,'(2X,A60,F20.10,A3)') 'Tr@BSE@G0W0@UHF correlation energy = ',sum(EcBSE),' au'
250-
! write(*,'(2X,A60,F20.10,A3)') 'Tr@BSE@G0W0@UHF total energy = ',ENuc + EUHF + sum(EcBSE),' au'
251-
! write(*,*)'-------------------------------------------------------------------------------'
252-
! write(*,*)
253-
!
239+
240+
if(dophBSE) then
241+
242+
call CVS_UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, &
243+
nCVS,nFC,occupations_fc,virtuals,S, &
244+
ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF,eGW,EcBSE)
245+
246+
write(*,*)
247+
write(*,*)'-------------------------------------------------------------------------------'
248+
write(*,'(2X,A60,F20.10,A3)') 'Tr@BSE@G0W0@UHF correlation energy (spin-conserved) = ',EcBSE(1),' au'
249+
write(*,'(2X,A60,F20.10,A3)') 'Tr@BSE@G0W0@UHF correlation energy (spin-flip) = ',EcBSE(2),' au'
250+
write(*,'(2X,A60,F20.10,A3)') 'Tr@BSE@G0W0@UHF correlation energy = ',sum(EcBSE),' au'
251+
write(*,'(2X,A60,F20.10,A3)') 'Tr@BSE@G0W0@UHF total energy = ',ENuc + EUHF + sum(EcBSE),' au'
252+
write(*,*)'-------------------------------------------------------------------------------'
253+
write(*,*)
254+
255+
! Adiabatic connection not implemented yet
254256
!! Compute the BSE correlation energy via the adiabatic connection
255257
!
256258
! if(doACFDT) then
@@ -269,7 +271,7 @@ subroutine CVS_UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE
269271
!
270272
! end if
271273
!
272-
! end if
274+
end if
273275
!
274276
!! Testing zone
275277
!

‎src/GW/CVS_UGW_phBSE.f90‎

Lines changed: 220 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,220 @@
1+
subroutine CVS_UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, &
2+
nCVS,nFC,occupations,virtuals,S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cW,eW,eGW,EcBSE)
3+
4+
! Compute the Bethe-Salpeter excitation energies
5+
6+
implicit none
7+
include 'parameters.h'
8+
9+
! Input variables
10+
11+
logical,intent(in) :: exchange_kernel
12+
logical,intent(in) :: TDA_W
13+
logical,intent(in) :: TDA
14+
logical,intent(in) :: dBSE
15+
logical,intent(in) :: dTDA
16+
logical,intent(in) :: spin_conserved
17+
logical,intent(in) :: spin_flip
18+
19+
double precision,intent(in) :: eta
20+
integer,intent(in) :: nBas
21+
integer,intent(in) :: nC(nspin)
22+
integer,intent(in) :: nO(nspin)
23+
integer,intent(in) :: nV(nspin)
24+
integer,intent(in) :: nR(nspin)
25+
integer,intent(in) :: nS(nspin)
26+
integer,intent(in) :: nCVS(nspin)
27+
integer,intent(in) :: nFC(nspin)
28+
integer,intent(in) :: occupations(maxval(nO-nFC),nspin)
29+
integer,intent(in) :: virtuals(nBas-minval(nO),nspin)
30+
double precision,intent(in) :: S(nBas,nBas)
31+
double precision,intent(in) :: cW(nBas,nBas,nspin)
32+
double precision,intent(in) :: eW(nBas,nspin)
33+
double precision,intent(in) :: eGW(nBas,nspin)
34+
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
35+
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
36+
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
37+
double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart)
38+
double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart)
39+
40+
! Local variables
41+
42+
integer :: ispin
43+
integer :: isp_W
44+
logical :: dRPA = .false.
45+
logical :: dRPA_W = .true.
46+
47+
double precision,allocatable :: Aph(:,:)
48+
double precision,allocatable :: Bph(:,:)
49+
double precision,allocatable :: KA(:,:)
50+
double precision,allocatable :: KB(:,:)
51+
52+
double precision :: EcRPA
53+
double precision,allocatable :: OmRPA(:)
54+
double precision,allocatable :: XpY_RPA(:,:)
55+
double precision,allocatable :: XmY_RPA(:,:)
56+
double precision,allocatable :: rho_RPA(:,:,:,:)
57+
58+
integer :: nS_aa,nS_bb,nS_sc
59+
integer :: nS_ab,nS_ba,nS_sf
60+
double precision,allocatable :: OmBSE(:)
61+
double precision,allocatable :: XpY_BSE(:,:)
62+
double precision,allocatable :: XmY_BSE(:,:)
63+
64+
! Output variables
65+
66+
double precision,intent(out) :: EcBSE(nspin)
67+
68+
! Memory allocation
69+
70+
nS_aa = (nO(1) - nFC(1))*(nV(1)-nCVS(1))
71+
nS_bb = (nO(2) - nFC(2))*(nV(2)-nCVS(2))
72+
nS_sc = nS_aa + nS_bb
73+
74+
nS_ab = (nO(1) - nFC(1))*(nV(2) - nCVS(2))
75+
nS_ba = (nO(2) - nFC(2))*(nV(1) - nCVS(1))
76+
nS_sf = nS_ab + nS_ba
77+
78+
! TDA
79+
80+
if(TDA) then
81+
write(*,*) 'Tamm-Dancoff approximation activated in phBSE!'
82+
write(*,*)
83+
end if
84+
85+
! Initialization
86+
87+
EcBSE(:) = 0d0
88+
89+
!--------------------------!
90+
! Spin-conserved screening !
91+
!--------------------------!
92+
93+
isp_W = 1
94+
95+
! Compute spin-conserved RPA screening
96+
97+
allocate(Aph(nS_sc,nS_sc),Bph(nS_sc,nS_sc))
98+
allocate(OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin))
99+
100+
call CVS_phULR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nCVS,nFC,occupations,virtuals, &
101+
1d0,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
102+
if(.not.TDA_W) call CVS_phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nCVS,nFC,occupations,virtuals, &
103+
1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
104+
call CVS_phULR(TDA_W,nS_aa,nS_bb,nS_sc,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
105+
call CVS_UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,nCVS,nFC,occupations,virtuals, &
106+
ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
107+
108+
deallocate(Aph,Bph)
109+
110+
!----------------------------!
111+
! Spin-conserved excitations !
112+
!----------------------------!
113+
114+
if(spin_conserved) then
115+
116+
ispin = 1
117+
118+
allocate(Aph(nS_sc,nS_sc),Bph(nS_sc,nS_sc),KA(nS_sc,nS_sc),KB(nS_sc,nS_sc))
119+
allocate(OmBSE(nS_sc),XpY_BSE(nS_sc,nS_sc),XmY_BSE(nS_sc,nS_sc))
120+
121+
! Compute spin-conserved BSE excitation energies
122+
123+
call CVS_phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nCVS,nFC,occupations,virtuals, &
124+
1d0,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
125+
if(.not.TDA) call CVS_phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nCVS,nFC,occupations,virtuals, &
126+
1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
127+
128+
call CVS_UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc, &
129+
nCVS,nFC,occupations,virtuals, &
130+
1d0,OmRPA,rho_RPA,KA)
131+
if(.not.TDA) call CVS_UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc, &
132+
nCVS,nFC,occupations,virtuals, &
133+
1d0,OmRPA,rho_RPA,KB)
134+
135+
Aph(:,:) = Aph(:,:) + KA(:,:)
136+
if(.not.TDA) Bph(:,:) = Bph(:,:) + KB(:,:)
137+
138+
call CVS_phULR(TDA,nS_aa,nS_bb,nS_sc,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
139+
140+
call print_excitation_energies('phBSE@GW@UHF','spin-conserved',nS_sc,OmBSE)
141+
call CVS_phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,nCVS,nFC,occupations,virtuals, &
142+
dipole_int_aa,dipole_int_bb, &
143+
cW,S,OmBSE,XpY_BSE,XmY_BSE)
144+
145+
! !-------------------------------------------------
146+
! ! Compute the dynamical screening at the BSE level
147+
! !-------------------------------------------------
148+
149+
! if(dBSE) &
150+
! call UGW_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,nS_sc, &
151+
! eW,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, &
152+
! OmRPA,rho_RPA,OmBSE,XpY_BSE,XmY_BSE)
153+
154+
deallocate(Aph,Bph,KA,KB)
155+
deallocate(OmBSE,XpY_BSE,XmY_BSE)
156+
157+
end if
158+
159+
!-----------------------!
160+
! Spin-flip excitations !
161+
!-----------------------!
162+
163+
if(spin_flip) then
164+
165+
ispin = 2
166+
167+
allocate(Aph(nS_sf,nS_sf),Bph(nS_sf,nS_sf),KA(nS_sf,nS_sf),KB(nS_sf,nS_sf))
168+
allocate(OmBSE(nS_sf),XpY_BSE(nS_sf,nS_sf),XmY_BSE(nS_sf,nS_sf))
169+
170+
! Compute spin-conserved BSE excitation energies
171+
172+
call CVS_phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,nCVS,nFC,occupations,virtuals, &
173+
1d0,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
174+
if(.not.TDA) call CVS_phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,nCVS,nFC,occupations,virtuals, &
175+
1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
176+
177+
call CVS_UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sf, &
178+
nCVS,nFC,occupations,virtuals, &
179+
1d0,OmRPA,rho_RPA,KA)
180+
if(.not.TDA) call CVS_UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sf, &
181+
nCVS,nFC,occupations,virtuals, &
182+
1d0,OmRPA,rho_RPA,KB)
183+
184+
Aph(:,:) = Aph(:,:) + KA(:,:)
185+
if(.not.TDA) Bph(:,:) = Bph(:,:) + KB(:,:)
186+
187+
call CVS_phULR(TDA,nS_aa,nS_bb,nS_sc,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
188+
189+
call print_excitation_energies('phBSE@GW@UHF','spin-flip',nS_sf,OmBSE)
190+
call CVS_phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,nCVS,nFC,occupations,virtuals,&
191+
dipole_int_aa,dipole_int_bb, &
192+
cW,S,OmBSE,XpY_BSE,XmY_BSE)
193+
194+
! !-------------------------------------------------
195+
! ! Compute the dynamical screening at the BSE level
196+
! !-------------------------------------------------
197+
198+
! if(dBSE) &
199+
! call UGW_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,nS_sc, &
200+
! eW,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, &
201+
! OmRPA,rho_RPA,OmBSE,XpY_BSE,XmY_BSE)
202+
203+
deallocate(Aph,Bph,KA,KB)
204+
deallocate(OmBSE,XpY_BSE,XmY_BSE)
205+
206+
end if
207+
208+
! Scale properly correlation energy if exchange is included in interaction kernel
209+
210+
if(exchange_kernel) then
211+
212+
EcBSE(:) = 0.5d0*EcBSE(:)
213+
214+
else
215+
216+
EcBSE(2) = 0.0d0
217+
218+
end if
219+
220+
end subroutine

0 commit comments

Comments
 (0)