LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-06-26 15:12:17 +02:00
Code Releases Activity

sf-BSE

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-09-23 22:48:54 +02:00
parent 2b1b2096c4
commit 08f3567d20
11 changed files with 275 additions and 238 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
input/options
View File

@ -5,11 +5,11 @@
# CC: maxSCF thresh DIIS n_diis
64 0.0000001 T 5
# spin: singlet triplet spin_conserved spin_flip TDA
T T T F F
T T T T F
# GF: maxSCF thresh DIIS n_diis lin eta renorm
256 0.00001 T 5 T 0.0 3
# GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0
256 0.00001 T 5 F 0.0 F F F F F
256 0.00001 T 5 T 0.0 F F F F F
# ACFDT: AC Kx XBS
F F T
# BSE: BSE dBSE dTDA evDyn

8
src/QuAcK/UG0W0.f90
View File

@ -112,8 +112,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
ispin = 1
call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0, &
eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc)
call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, &
eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sc,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc)
if(print_W) call print_excitation('RPA@UHF',5,nS_sc,Omega_sc)
@ -165,8 +165,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
! Compute the RPA correlation energy
call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0, &
eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc)
call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, &
eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sc,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'

3
src/QuAcK/UHF.f90
View File

@ -232,6 +232,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,EUH
! Compute final UHF energy
call print_UHF(nBas,nO(:),e(:,:),c(:,:,:),ENuc,ET(:),EV(:),EJ(:),Ex(:),EUHF)
call matout(nBas,2,e)
call print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF)
end subroutine UHF

8
src/QuAcK/URPAx.f90
View File

@ -73,8 +73,8 @@ subroutine URPAx(doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO
allocate(Omega_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc))
call unrestricted_linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,e, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_sc,EcRPAx(ispin),Omega_sc,XpY_sc,XmY_sc)
call unrestricted_linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sc,rho_sc,EcRPAx(ispin),Omega_sc,XpY_sc,XmY_sc)
call print_excitation('URPAx ',5,nS_sc,Omega_sc)
! call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
@ -95,8 +95,8 @@ subroutine URPAx(doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO
allocate(Omega_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf))
call unrestricted_linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,1d0,e, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_sf,EcRPAx(ispin),Omega_sf,XpY_sf,XmY_sf)
call unrestricted_linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,nS_sf,1d0,e, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sf,rho_sf,EcRPAx(ispin),Omega_sf,XpY_sf,XmY_sf)
call print_excitation('URPAx ',6,nS_sf,Omega_sf)
! call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))

8
src/QuAcK/UdRPA.f90
View File

@ -73,8 +73,8 @@ subroutine UdRPA(doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO
allocate(Omega_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc))
call unrestricted_linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,e, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc)
call unrestricted_linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sc,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc)
call print_excitation('URPA ',5,nS_sc,Omega_sc)
! call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
@ -95,8 +95,8 @@ subroutine UdRPA(doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO
allocate(Omega_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf))
call unrestricted_linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,1d0,e, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_sf,EcRPA(ispin),Omega_sf,XpY_sf,XmY_sf)
call unrestricted_linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,nS_sf,1d0,e, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sf,rho_sf,EcRPA(ispin),Omega_sf,XpY_sf,XmY_sf)
call print_excitation('URPA ',6,nS_sf,Omega_sf)
! call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))

82
src/QuAcK/unrestricted_Bethe_Salpeter.f90
View File

@ -36,8 +36,8 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,
integer :: ispin
integer :: isp_W
integer :: nS_aa,nS_bb,nS_sc
integer :: nS_ab,nS_ba,nS_sf
double precision,allocatable :: OmRPA_sc(:)
double precision,allocatable :: XpY_RPA_sc(:,:)
double precision,allocatable :: XmY_RPA_sc(:,:)
@ -46,6 +46,11 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,
double precision,allocatable :: XpY_BSE_sc(:,:)
double precision,allocatable :: XmY_BSE_sc(:,:)
integer :: nS_ab,nS_ba,nS_sf
double precision,allocatable :: OmBSE_sf(:)
double precision,allocatable :: XpY_BSE_sf(:,:)
double precision,allocatable :: XmY_BSE_sf(:,:)
! Output variables
double precision,intent(out) :: EcRPA(nspin)
@ -55,38 +60,41 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,
! Spin-conserved excitations !
!----------------------------!
if(spin_conserved) then
isp_W = 1
! Memory allocation
nS_aa = nS(1)
nS_bb = nS(2)
nS_sc = nS_aa + nS_bb
allocate(OmRPA_sc(nS_sc),XpY_RPA_sc(nS_sc,nS_sc),XmY_RPA_sc(nS_sc,nS_sc),rho_RPA_sc(nBas,nBas,nS_sc,nspin))
! Compute spin-conserved RPA screening
call unrestricted_linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, &
eW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA_sc,rho_RPA_sc,EcRPA(isp_W), &
OmRPA_sc,XpY_RPA_sc,XmY_RPA_sc)
! call print_excitation('RPA@UG0W0',5,nS_sc,OmRPA_sc)
call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb, &
XpY_RPA_sc,rho_RPA_sc)
if(spin_conserved) then
ispin = 1
isp_W = 1
EcBSE(ispin) = 0d0
! Memory allocation
nS_aa = nS(1)
nS_bb = nS(2)
nS_sc = nS_aa + nS_bb
allocate(OmRPA_sc(nS_sc),XpY_RPA_sc(nS_sc,nS_sc),XmY_RPA_sc(nS_sc,nS_sc),rho_RPA_sc(nBas,nBas,nS_sc,nspin))
allocate(OmBSE_sc(nS_sc),XpY_BSE_sc(nS_sc,nS_sc),XmY_BSE_sc(nS_sc,nS_sc))
! Compute spin-conserved RPA screening
call unrestricted_linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0, &
eW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_RPA_sc,EcRPA(ispin), &
OmRPA_sc,XpY_RPA_sc,XmY_RPA_sc)
! call print_excitation('RPA@UG0W0',5,nS_sc,OmRPA_sc)
call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb, &
XpY_RPA_sc,rho_RPA_sc)
! Compute spin-conserved BSE excitation energies
OmBSE_sc(:) = OmRPA_sc(:)
call unrestricted_linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0, &
eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_RPA_sc,EcBSE(ispin), &
call unrestricted_linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, &
eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA_sc,rho_RPA_sc,EcBSE(ispin), &
OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc)
call print_excitation('BSE@UG0W0',5,nS_sc,OmBSE_sc)
@ -117,25 +125,25 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,
! Spin-flip excitations !
!-----------------------!
!if(spin_flip) then
if(spin_flip) then
! ispin = 2
! isp_W = 1
! EcBSE(ispin) = 0d0
ispin = 2
! ! Compute (singlet) RPA screening
EcBSE(ispin) = 0d0
! 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))
! Memory allocation
! ! Compute BSE excitation energies
nS_ab = (nO(1) - nC(1))*(nV(2) - nR(2))
nS_ba = (nO(2) - nC(2))*(nV(1) - nR(1))
nS_sf = nS_ab + nS_ba
allocate(OmBSE_sf(nS_sf),XpY_BSE_sf(nS_sf,nS_sf),XmY_BSE_sf(nS_sf,nS_sf))
! OmBSE(:,ispin) = OmRPA(:,ispin)
call unrestricted_linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sc,1d0, &
eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA_sc,rho_RPA_sc,EcBSE(ispin), &
OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf)
! call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
! 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@UG0W0',6,nS_sf,OmBSE_sf)
!-------------------------------------------------
! Compute the dynamical screening at the BSE level
@ -157,6 +165,6 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,
! end if
! end if
end if
end subroutine unrestricted_Bethe_Salpeter

180
src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90
View File

@ -1,4 +1,4 @@
subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, &
subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nSsc,lambda, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega,rho,A_lr)
! Compute the extra term for Bethe-Salpeter equation for linear response in the unrestricted formalism
@ -8,6 +8,7 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt
! Input variables
integer,intent(in) :: ispin
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
@ -16,14 +17,15 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt
integer,intent(in) :: nSa
integer,intent(in) :: nSb
integer,intent(in) :: nSt
integer,intent(in) :: nSsc
double precision,intent(in) :: eta
double precision,intent(in) :: lambda
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Omega(nSt)
double precision,intent(in) :: rho(nBas,nBas,nSt,nspin)
double precision,intent(in) :: Omega(nSsc)
double precision,intent(in) :: rho(nBas,nBas,nSsc,nspin)
! Local variables
@ -35,108 +37,116 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt
double precision,intent(out) :: A_lr(nSt,nSt)
!--------------------------------!
! Build part A of the BSE matrix !
!--------------------------------!
!--------------------------------------------------!
! Build BSE matrix for spin-conserving transitions !
!--------------------------------------------------!
! aaaa block
if(ispin == 1) then
ia = 0
do i=nC(1)+1,nO(1)
do a=nO(1)+1,nBas-nR(1)
ia = ia + 1
jb = 0
do j=nC(1)+1,nO(1)
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
! aaaa block
ia = 0
do i=nC(1)+1,nO(1)
do a=nO(1)+1,nBas-nR(1)
ia = ia + 1
jb = 0
do j=nC(1)+1,nO(1)
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
chi = 0d0
do kc=1,nSsc
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,j,kc,1)*rho(a,b,kc,1)*Omega(kc)/eps
enddo
A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_aaaa(i,b,j,a) + 4d0*lambda*chi
chi = 0d0
do kc=1,nSt
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,j,kc,1)*rho(a,b,kc,1)*Omega(kc)/eps &
+ rho(i,j,kc,1)*rho(a,b,kc,1)*Omega(kc)/eps
enddo
A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_aaaa(i,b,j,a) + 2d0*lambda*chi
enddo
enddo
enddo
enddo
! aabb block
ia = 0
do i=nC(1)+1,nO(1)
do a=nO(1)+1,nBas-nR(1)
ia = ia + 1
jb = 0
do j=nC(2)+1,nO(2)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
! bbbb block
ia = 0
do i=nC(2)+1,nO(2)
do a=nO(2)+1,nBas-nR(2)
ia = ia + 1
jb = 0
do j=nC(2)+1,nO(2)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
chi = 0d0
do kc=1,nSsc
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,j,kc,2)*rho(a,b,kc,2)*Omega(kc)/eps
enddo
A_lr(nSa+ia,nSa+jb) = A_lr(nSa+ia,nSa+jb) - lambda*ERI_bbbb(i,b,j,a) + 4d0*lambda*chi
chi = 0d0
do kc=1,nSt
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,j,kc,1)*rho(a,b,kc,1)*Omega(kc)/eps &
+ rho(i,j,kc,2)*rho(a,b,kc,2)*Omega(kc)/eps
enddo
A_lr(ia,nSa+jb) = A_lr(ia,nSa+jb) - lambda*ERI_aabb(i,b,j,a) + 2d0*lambda*chi
enddo
enddo
enddo
enddo
! bbaa block
end if
ia = 0
do i=nC(2)+1,nO(2)
do a=nO(2)+1,nBas-nR(2)
ia = ia + 1
jb = 0
do j=nC(1)+1,nO(1)
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
chi = 0d0
do kc=1,nSt
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,j,kc,2)*rho(a,b,kc,2)*Omega(kc)/eps &
+ rho(i,j,kc,1)*rho(a,b,kc,1)*Omega(kc)/eps
enddo
!--------------------------------------------!
! Build BSE matrix for spin-flip transitions !
!--------------------------------------------!
A_lr(nSa+ia,jb) = A_lr(nSa+ia,jb) - lambda*ERI_aabb(b,i,a,j) + 2d0*lambda*chi
if(ispin == 2) then
enddo
enddo
enddo
enddo
! abab block
! bbbb block
ia = 0
do i=nC(1)+1,nO(1)
do a=nO(2)+1,nBas-nR(2)
ia = ia + 1
jb = 0
do j=nC(1)+1,nO(1)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
ia = 0
do i=nC(2)+1,nO(2)
do a=nO(2)+1,nBas-nR(2)
ia = ia + 1
jb = 0
do j=nC(2)+1,nO(2)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
chi = 0d0
do kc=1,nSt
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,j,kc,2)*rho(a,b,kc,2)*Omega(kc)/eps &
+ rho(i,j,kc,2)*rho(a,b,kc,2)*Omega(kc)/eps
enddo
chi = 0d0
do kc=1,nSsc
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,j,kc,1)*rho(a,b,kc,2)*Omega(kc)/eps
enddo
A_lr(nSa+ia,nSa+jb) = A_lr(nSa+ia,nSa+jb) - lambda*ERI_bbbb(i,b,j,a) + 2d0*lambda*chi
A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_abab(i,b,j,a) + 4d0*lambda*chi
enddo
enddo
enddo
enddo
end do
end do
end do
end do
! baba block
ia = 0
do i=nC(2)+1,nO(2)
do a=nO(1)+1,nBas-nR(1)
ia = ia + 1
jb = 0
do j=nC(2)+1,nO(2)
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
chi = 0d0
do kc=1,nSsc
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,j,kc,2)*rho(a,b,kc,1)*Omega(kc)/eps
enddo
A_lr(nSa+ia,nSa+jb) = A_lr(nSa+ia,nSa+jb) - lambda*ERI_abab(b,i,a,j) + 4d0*lambda*chi
end do
end do
end do
end do
end if
end subroutine unrestricted_Bethe_Salpeter_A_matrix

185
src/QuAcK/unrestricted_Bethe_Salpeter_B_matrix.f90
View File

@ -1,6 +1,5 @@
subroutine unrestricted_Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab, &
Omega,rho,B_lr)
subroutine unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nSsc,lambda, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega,rho,B_lr)
! Compute the extra term for Bethe-Salpeter equation for linear response
@ -9,6 +8,7 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt
! Input variables
integer,intent(in) :: ispin
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
@ -17,14 +17,15 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt
integer,intent(in) :: nSa
integer,intent(in) :: nSb
integer,intent(in) :: nSt
integer,intent(in) :: nSsc
double precision,intent(in) :: eta
double precision,intent(in) :: lambda
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Omega(nSt)
double precision,intent(in) :: rho(nBas,nBas,nSt,nspin)
double precision,intent(in) :: Omega(nSsc)
double precision,intent(in) :: rho(nBas,nBas,nSsc,nspin)
! Local variables
@ -36,104 +37,118 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt
double precision,intent(out) :: B_lr(nSt,nSt)
! alpha-alpha block
!--------------------------------------------------!
! Build BSE matrix for spin-conserving transitions !
!--------------------------------------------------!
ia = 0
do i=nC(1)+1,nO(1)
do a=nO(1)+1,nBas-nR(1)
ia = ia + 1
jb = 0
do j=nC(1)+1,nO(1)
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
if(ispin == 1) then
! aaaa block
ia = 0
do i=nC(1)+1,nO(1)
do a=nO(1)+1,nBas-nR(1)
ia = ia + 1
jb = 0
do j=nC(1)+1,nO(1)
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
chi = 0d0
do kc=1,nSsc
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,b,kc,1)*rho(a,j,kc,1)*Omega(kc)/eps
enddo
B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI_aaaa(i,j,b,a) + 4d0*lambda*chi
chi = 0d0
do kc=1,nSt
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,b,kc,1)*rho(a,j,kc,1)*Omega(kc)/eps &
+ rho(i,b,kc,1)*rho(a,j,kc,1)*Omega(kc)/eps
enddo
B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI_aaaa(i,j,b,a) + 2d0*lambda*chi
enddo
enddo
enddo
enddo
! alpha-beta block
ia = 0
do i=nC(1)+1,nO(1)
do a=nO(1)+1,nBas-nR(1)
ia = ia + 1
jb = 0
do j=nC(2)+1,nO(2)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
chi = 0d0
do kc=1,nSt
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,b,kc,1)*rho(a,j,kc,1)*Omega(kc)/eps &
+ rho(i,b,kc,2)*rho(a,j,kc,2)*Omega(kc)/eps
! bbbb block
ia = 0
do i=nC(2)+1,nO(2)
do a=nO(2)+1,nBas-nR(2)
ia = ia + 1
jb = 0
do j=nC(2)+1,nO(2)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
chi = 0d0
do kc=1,nSsc
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,b,kc,2)*rho(a,j,kc,2)*Omega(kc)/eps
enddo
B_lr(nSa+ia,nSa+jb) = B_lr(nSa+ia,nSa+jb) - lambda*ERI_bbbb(i,j,b,a) + 4d0*lambda*chi
enddo
B_lr(ia,nSa+jb) = B_lr(ia,nSa+jb) - lambda*ERI_aabb(i,j,b,a) + 2d0*lambda*chi
enddo
enddo
enddo
enddo
! beta-alpha block
end if
ia = 0
do i=nC(2)+1,nO(2)
do a=nO(2)+1,nBas-nR(2)
ia = ia + 1
jb = 0
do j=nC(1)+1,nO(1)
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
chi = 0d0
do kc=1,nSt
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,b,kc,2)*rho(a,j,kc,2)*Omega(kc)/eps &
+ rho(i,b,kc,1)*rho(a,j,kc,1)*Omega(kc)/eps
enddo
B_lr(nSa+ia,jb) = B_lr(nSa+ia,jb) - lambda*ERI_aabb(j,i,a,b) + 2d0*lambda*chi
!--------------------------------------------!
! Build BSE matrix for spin-flip transitions !
!--------------------------------------------!
enddo
enddo
enddo
enddo
if(ispin == 2) then
! beta-beta block
! abba block
ia = 0
do i=nC(2)+1,nO(2)
do a=nO(2)+1,nBas-nR(2)
ia = ia + 1
jb = 0
do j=nC(2)+1,nO(2)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
chi = 0d0
do kc=1,nSt
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,b,kc,2)*rho(a,j,kc,2)*Omega(kc)/eps &
+ rho(i,b,kc,2)*rho(a,j,kc,2)*Omega(kc)/eps
enddo
ia = 0
do i=nC(1)+1,nO(1)
do a=nO(2)+1,nBas-nR(2)
ia = ia + 1
jb = 0
do j=nC(2)+1,nO(2)
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
B_lr(nSa+ia,nSa+jb) = B_lr(nSa+ia,nSa+jb) - lambda*ERI_bbbb(i,j,b,a) + 2d0*lambda*chi
chi = 0d0
do kc=1,nSsc
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,b,kc,1)*rho(a,j,kc,2)*Omega(kc)/eps
enddo
enddo
enddo
enddo
enddo
B_lr(ia,nSa+jb) = B_lr(ia,nSa+jb) - lambda*ERI_abab(i,a,b,j) + 4d0*lambda*chi
end do
end do
end do
end do
! baab block
ia = 0
do i=nC(2)+1,nO(2)
do a=nO(1)+1,nBas-nR(1)
ia = ia + 1
jb = 0
do j=nC(1)+1,nO(1)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
chi = 0d0
do kc=1,nSsc
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,b,kc,2)*rho(a,j,kc,1)*Omega(kc)/eps
enddo
B_lr(nSa+ia,jb) = B_lr(nSa+ia,jb) - lambda*ERI_abab(b,j,i,a) + 4d0*lambda*chi
end do
end do
end do
end do
end if
end subroutine unrestricted_Bethe_Salpeter_B_matrix

17
src/QuAcK/unrestricted_linear_response.f90
View File

@ -1,5 +1,5 @@
subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, &
e,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho,EcRPA,Omega,XpY,XmY)
subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nSsc,lambda, &
e,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_RPA,rho_RPA,EcRPA,Omega,XpY,XmY)
! Compute linear response for unrestricted formalism
@ -21,13 +21,16 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,
integer,intent(in) :: nSa
integer,intent(in) :: nSb
integer,intent(in) :: nSt
integer,intent(in) :: nSsc
double precision,intent(in) :: lambda
double precision,intent(in) :: e(nBas,nspin)
double precision,intent(in) :: rho(nBas,nBas,nSt,nspin)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Omega_RPA(nSsc)
double precision,intent(in) :: rho_RPA(nBas,nBas,nSsc,nspin)
! Local variables
@ -58,8 +61,8 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,A)
if(BSE) &
call unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega,rho,A)
call unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nSsc,lambda, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_RPA,rho_RPA,A)
! Tamm-Dancoff approximation
@ -70,8 +73,8 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,B)
if(BSE) &
call unrestricted_Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega,rho,B)
call unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nSsc,lambda, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_RPA,rho_RPA,B)
end if

4
src/QuAcK/unrestricted_linear_response_A_matrix.f90
View File

@ -143,7 +143,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
jb = jb + 1
A_lr(ia,jb) = (e(a,2) - e(i,1))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
+ lambda*ERI_abab(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI_abab(i,b,j,a)
- (1d0 - delta_dRPA)*lambda*ERI_abab(i,b,j,a)
end do
end do
@ -162,7 +162,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
jb = jb + 1
A_lr(nSa+ia,nSa+jb) = (e(a,1) - e(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
+ lambda*ERI_abab(b,i,j,a) - (1d0 - delta_dRPA)*lambda*ERI_abab(b,i,a,j)
- (1d0 - delta_dRPA)*lambda*ERI_abab(b,i,a,j)
end do
end do

14
src/QuAcK/unrestricted_linear_response_B_matrix.f90
View File

@ -128,7 +128,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
B_lr(:,:) = 0d0
! abab block
! abba block
ia = 0
do i=nC(1)+1,nO(1)
@ -136,28 +136,28 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
ia = ia + 1
jb = 0
do j=nC(2)+1,nO(2)
do b=nO(2)+1,nBas-nR(2)
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
B_lr(ia,jb) = lambda*ERI_abab(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI_abab(i,j,b,a)
B_lr(ia,nSa+jb) = - (1d0 - delta_dRPA)*lambda*ERI_abab(i,a,b,j)
end do
end do
end do
end do
! bbbb block
! baab block
ia = 0
do i=nC(2)+1,nO(2)
do a=nO(1)+1,nBas-nR(1)
ia = ia + 1
jb = 0
do j=nC(2)+1,nO(2)
do b=nO(1)+1,nBas-nR(1)
do j=nC(1)+1,nO(1)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
B_lr(nSa+ia,nSa+jb) = lambda*ERI_abab(j,i,b,a) - (1d0 - delta_dRPA)*lambda*ERI_abab(j,i,a,b)
B_lr(nSa+ia,jb) = - (1d0 - delta_dRPA)*lambda*ERI_abab(b,j,i,a)
end do
end do