Skip to content

Commit 7a5608a

Browse files
author
lburth
committed
dynamic perturbation BSE for CVS
1 parent 7f47101 commit 7a5608a

File tree

3 files changed

+387
-14
lines changed

3 files changed

+387
-14
lines changed

‎src/GW/CVS_UGW_phBSE.f90‎

Lines changed: 16 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -142,14 +142,15 @@ subroutine CVS_UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin
142142
dipole_int_aa,dipole_int_bb, &
143143
cW,S,OmBSE,XpY_BSE,XmY_BSE)
144144

145-
! !-------------------------------------------------
146-
! ! Compute the dynamical screening at the BSE level
147-
! !-------------------------------------------------
145+
!-------------------------------------------------
146+
! Compute the dynamical screening at the BSE level
147+
!-------------------------------------------------
148148

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)
149+
if(dBSE) &
150+
call CVS_UGW_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,nS_sc, &
151+
nCVS,nFC,occupations,virtuals, &
152+
eW,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, &
153+
OmRPA,rho_RPA,OmBSE,XpY_BSE,XmY_BSE)
153154

154155
deallocate(Aph,Bph,KA,KB)
155156
deallocate(OmBSE,XpY_BSE,XmY_BSE)
@@ -191,14 +192,15 @@ subroutine CVS_UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin
191192
dipole_int_aa,dipole_int_bb, &
192193
cW,S,OmBSE,XpY_BSE,XmY_BSE)
193194

194-
! !-------------------------------------------------
195-
! ! Compute the dynamical screening at the BSE level
196-
! !-------------------------------------------------
195+
!-------------------------------------------------
196+
! Compute the dynamical screening at the BSE level
197+
!-------------------------------------------------
197198

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)
199+
if(dBSE) &
200+
call CVS_UGW_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,nS_sc, &
201+
nCVS,nFC,occupations,virtuals, &
202+
eW,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, &
203+
OmRPA,rho_RPA,OmBSE,XpY_BSE,XmY_BSE)
202204

203205
deallocate(Aph,Bph,KA,KB)
204206
deallocate(OmBSE,XpY_BSE,XmY_BSE)

‎src/GW/CVS_UGW_phBSE_dynamic_kernel_A.f90‎

Lines changed: 265 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,265 @@
1+
subroutine CVS_UGW_phBSE_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,eGW, &
2+
nCVS,nFC,occupations,virtuals, &
3+
ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,OmBSE,A_dyn,ZA_dyn)
4+
5+
! Compute the extra term for dynamical Bethe-Salpeter equation for linear response in the unrestricted formalism
6+
7+
implicit none
8+
include 'parameters.h'
9+
10+
! Input variables
11+
12+
integer,intent(in) :: ispin
13+
integer,intent(in) :: nBas
14+
integer,intent(in) :: nC(nspin)
15+
integer,intent(in) :: nO(nspin)
16+
integer,intent(in) :: nV(nspin)
17+
integer,intent(in) :: nR(nspin)
18+
integer,intent(in) :: nSa
19+
integer,intent(in) :: nSb
20+
integer,intent(in) :: nSt
21+
integer,intent(in) :: nS_sc
22+
integer,intent(in) :: nCVS(nspin)
23+
integer,intent(in) :: nFC(nspin)
24+
integer,intent(in) :: occupations(maxval(nO-nFC),nspin)
25+
integer,intent(in) :: virtuals(nBas-minval(nO),nspin)
26+
double precision,intent(in) :: eta
27+
double precision,intent(in) :: lambda
28+
double precision,intent(in) :: eGW(nBas,nspin)
29+
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
30+
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
31+
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
32+
double precision,intent(in) :: OmRPA(nS_sc)
33+
double precision,intent(in) :: rho_RPA(nBas,nBas,nS_sc,nspin)
34+
double precision,intent(in) :: OmBSE
35+
36+
! Local variables
37+
38+
double precision :: chi
39+
double precision :: eps
40+
integer :: i,j,a,b,ia,jb,kc
41+
42+
! Output variables
43+
44+
double precision,intent(out) :: A_dyn(nSt,nSt)
45+
double precision,intent(out) :: ZA_dyn(nSt,nSt)
46+
47+
!--------------------------------------------------!
48+
! Build BSE matrix for spin-conserving transitions !
49+
!--------------------------------------------------!
50+
51+
A_dyn(:,:) = 0d0
52+
53+
if(ispin == 1) then
54+
55+
! aaaa block
56+
57+
ia = 0
58+
do i=1,nO(1)-nFC(1)
59+
do a=nCVS(1)+1,nBas-nO(1)
60+
ia = ia + 1
61+
jb = 0
62+
do j=1,nO(1)-nFC(1)
63+
do b=nCVS(1)+1,nBas-nO(1)
64+
jb = jb + 1
65+
66+
chi = 0d0
67+
do kc=1,nS_sc
68+
chi = chi + rho_RPA(occupations(i,1),occupations(j,1),kc,1)&
69+
* rho_RPA(virtuals(a,1),virtuals(b,1),kc,1) &
70+
* OmRPA(kc)/(OmRPA(kc)**2 + eta**2)
71+
end do
72+
73+
A_dyn(ia,jb) = A_dyn(ia,jb) - 2d0*lambda*chi
74+
75+
chi = 0d0
76+
do kc=1,nS_sc
77+
78+
eps = + OmBSE - OmRPA(kc) - (eGW(virtuals(a,1),1) - eGW(occupations(j,1),1))
79+
chi = chi + rho_RPA(occupations(i,1),occupations(j,1),kc,1)*rho_RPA(virtuals(a,1),virtuals(b,1),kc,1)*eps/(eps**2 + eta**2)
80+
81+
eps = + OmBSE - OmRPA(kc) - (eGW(virtuals(b,1),1) - eGW(occupations(i,1),1))
82+
chi = chi + rho_RPA(occupations(i,1),occupations(j,1),kc,1)*rho_RPA(virtuals(a,1),virtuals(b,1),kc,1)*eps/(eps**2 + eta**2)
83+
84+
end do
85+
86+
A_dyn(ia,jb) = A_dyn(ia,jb) - lambda*chi
87+
88+
chi = 0d0
89+
do kc=1,nS_sc
90+
91+
eps = + OmBSE - OmRPA(kc) - (eGW(virtuals(a,1),1) - eGW(occupations(j,1),1))
92+
chi = chi + rho_RPA(occupations(i,1),occupations(j,1),kc,1)*rho_RPA(virtuals(a,1),virtuals(b,1),kc,1)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
93+
94+
eps = + OmBSE - OmRPA(kc) - (eGW(virtuals(b,1),1) - eGW(occupations(i,1),1))
95+
chi = chi + rho_RPA(i,j,kc,1)*rho_RPA(a,b,kc,1)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
96+
97+
end do
98+
99+
ZA_dyn(ia,jb) = ZA_dyn(ia,jb) + lambda*chi
100+
101+
end do
102+
end do
103+
end do
104+
end do
105+
106+
! bbbb block
107+
108+
ia = 0
109+
do i=1,nO(2)-nFC(2)
110+
do a=nCVS(2)+1,nBas-nO(2)
111+
ia = ia + 1
112+
jb = 0
113+
do j=1,nO(2)-nFC(2)
114+
do b=nCVS(2)+1,nBas-nO(2)
115+
jb = jb + 1
116+
117+
chi = 0d0
118+
do kc=1,nS_sc
119+
chi = chi + rho_RPA(occupations(i,2),occupations(j,2),kc,2)&
120+
* rho_RPA(virtuals(a,2),virtuals(b,2),kc,2) &
121+
* OmRPA(kc)/(OmRPA(kc)**2 + eta**2)
122+
end do
123+
124+
A_dyn(nSa+ia,nSa+jb) = A_dyn(nSa+ia,nSa+jb) - 2d0*lambda*chi
125+
126+
chi = 0d0
127+
do kc=1,nS_sc
128+
129+
eps = + OmBSE - OmRPA(kc) - (eGW(virtuals(a,2),2) - eGW(occupations(j,2),2))
130+
chi = chi + rho_RPA(occupations(i,2),occupations(j,2),kc,2)*rho_RPA(virtuals(a,2),virtuals(b,2),kc,2)*eps/(eps**2 + eta**2)
131+
132+
eps = + OmBSE - OmRPA(kc) - (eGW(virtuals(b,2),2) - eGW(occupations(i,2),2))
133+
chi = chi + rho_RPA(occupations(i,2),occupations(j,2),kc,2)*rho_RPA(virtuals(a,2),virtuals(b,2),kc,2)*eps/(eps**2 + eta**2)
134+
135+
end do
136+
137+
A_dyn(nSa+ia,nSa+jb) = A_dyn(nSa+ia,nSa+jb) - lambda*chi
138+
139+
chi = 0d0
140+
do kc=1,nS_sc
141+
142+
eps = + OmBSE - OmRPA(kc) - (eGW(virtuals(a,2),2) - eGW(occupations(j,2),2))
143+
chi = chi + rho_RPA(occupations(i,2),occupations(j,2),kc,2)*rho_RPA(virtuals(a,2),virtuals(b,2),kc,2)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
144+
145+
eps = + OmBSE - OmRPA(kc) - (eGW(virtuals(b,2),2) - eGW(occupations(i,2),2))
146+
chi = chi + rho_RPA(occupations(i,2),occupations(j,2),kc,2)*rho_RPA(virtuals(a,2),virtuals(b,2),kc,2)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
147+
148+
end do
149+
150+
ZA_dyn(nSa+ia,nSa+jb) = ZA_dyn(nSa+ia,nSa+jb) + lambda*chi
151+
152+
end do
153+
end do
154+
end do
155+
end do
156+
157+
end if
158+
159+
!--------------------------------------------!
160+
! Build BSE matrix for spin-flip transitions !
161+
!--------------------------------------------!
162+
163+
if(ispin == 2) then
164+
165+
! abab block
166+
167+
ia = 0
168+
do i=1,nO(1)-nFC(1)
169+
do a=nCVS(2)+1,nBas-nO(2)
170+
ia = ia + 1
171+
jb = 0
172+
do j=1,nO(1)-nFC(1)
173+
do b=nCVS(2)+1,nBas-nO(2)
174+
jb = jb + 1
175+
176+
chi = 0d0
177+
do kc=1,nS_sc
178+
chi = chi + rho_RPA(occupations(i,1),occupations(j,1),kc,1)*rho_RPA(virtuals(a,2),virtuals(b,2),kc,2)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2)
179+
end do
180+
181+
A_dyn(ia,jb) = A_dyn(ia,jb) - 2d0*lambda*chi
182+
183+
chi = 0d0
184+
do kc=1,nS_sc
185+
186+
eps = + OmBSE - OmRPA(kc) - (eGW(virtuals(a,2),2) - eGW(occupations(j,1),1))
187+
chi = chi + rho_RPA(occupations(i,1),occupations(j,1),kc,1)*rho_RPA(virtuals(a,2),virtuals(b,2),kc,2)*eps/(eps**2 + eta**2)
188+
189+
eps = + OmBSE - OmRPA(kc) - (eGW(virtuals(b,2),2) - eGW(occupations(i,1),1))
190+
chi = chi + rho_RPA(occupations(i,1),occupations(j,1),kc,1)*rho_RPA(virtuals(a,2),virtuals(b,2),kc,2)*eps/(eps**2 + eta**2)
191+
192+
end do
193+
194+
A_dyn(ia,jb) = A_dyn(ia,jb) - lambda*chi
195+
196+
chi = 0d0
197+
do kc=1,nS_sc
198+
199+
eps = + OmBSE - OmRPA(kc) - (eGW(virtuals(a,2),2) - eGW(occupations(j,1),1))
200+
chi = chi + rho_RPA(occupations(i,1),occupations(j,1),kc,1)*rho_RPA(virtuals(a,2),virtuals(b,2),kc,2)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
201+
202+
eps = + OmBSE - OmRPA(kc) - (eGW(virtuals(b,2),2) - eGW(occupations(i,1),1))
203+
chi = chi + rho_RPA(occupations(i,1),occupations(j,1),kc,1)*rho_RPA(virtuals(a,2),virtuals(b,2),kc,2)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
204+
205+
end do
206+
207+
ZA_dyn(ia,jb) = ZA_dyn(ia,jb) + lambda*chi
208+
209+
end do
210+
end do
211+
end do
212+
end do
213+
214+
! baba block
215+
216+
ia = 0
217+
do i=1,nO(2)-nFC(2)
218+
do a=nCVS(1)+1,nBas-nO(1)
219+
ia = ia + 1
220+
jb = 0
221+
do j=1,nO(2)-nFC(2)
222+
do b=nCVS(1)+1,nBas-nO(1)
223+
jb = jb + 1
224+
225+
chi = 0d0
226+
do kc=1,nS_sc
227+
chi = chi + rho_RPA(occupations(i,2),occupations(j,2),kc,2)*rho_RPA(virtuals(a,1),virtuals(b,1),kc,1)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2)
228+
end do
229+
230+
A_dyn(nSa+ia,nSa+jb) = A_dyn(nSa+ia,nSa+jb) - 2d0*lambda*chi
231+
232+
chi = 0d0
233+
do kc=1,nS_sc
234+
235+
eps = + OmBSE - OmRPA(kc) - (eGW(virtuals(a,1),1) - eGW(occupations(j,2),2))
236+
chi = chi + rho_RPA(occupations(i,2),occupations(j,2),kc,2)*rho_RPA(virtuals(a,1),virtuals(b,1),kc,1)*eps/(eps**2 + eta**2)
237+
238+
eps = + OmBSE - OmRPA(kc) - (eGW(virtuals(b,1),1) - eGW(occupations(i,2),2))
239+
chi = chi + rho_RPA(occupations(i,2),occupations(j,2),kc,2)*rho_RPA(virtuals(a,1),virtuals(b,1),kc,1)*eps/(eps**2 + eta**2)
240+
241+
end do
242+
243+
A_dyn(nSa+ia,nSa+jb) = A_dyn(nSa+ia,nSa+jb) - lambda*chi
244+
245+
chi = 0d0
246+
do kc=1,nS_sc
247+
248+
eps = + OmBSE - OmRPA(kc) - (eGW(virtuals(a,1),1) - eGW(occupations(j,2),2))
249+
chi = chi + rho_RPA(occupations(i,2),occupations(j,2),kc,2)*rho_RPA(virtuals(a,1),virtuals(b,1),kc,1)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
250+
251+
eps = + OmBSE - OmRPA(kc) - (eGW(virtuals(b,1),1) - eGW(occupations(i,2),2))
252+
chi = chi + rho_RPA(occupations(i,2),occupations(j,2),kc,2)*rho_RPA(virtuals(a,1),virtuals(b,1),kc,1)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
253+
254+
end do
255+
256+
ZA_dyn(nSa+ia,nSa+jb) = ZA_dyn(nSa+ia,nSa+jb) + lambda*chi
257+
258+
end do
259+
end do
260+
end do
261+
end do
262+
263+
end if
264+
265+
end subroutine

0 commit comments

Comments
 (0)