mirror of
https://github.com/pfloos/quack
synced 2025-05-06 23:24:58 +02:00
WIP: Implement fully complex G0W0
This commit is contained in:
Loris Burth
parent
ebbd1df3cc
commit
cdf6f5ecf2
245
src/GW/complex_cRG0W0.f90
Normal file
245
src/GW/complex_cRG0W0.f90
Normal file
@ -0,0 +1,245 @@
|
|||||||
|
subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
|
||||||
|
linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,CAP,dipole_int,eHF)
|
||||||
|
|
||||||
|
! Perform a fully complex G0W0 calculation with CAP
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
include 'parameters.h'
|
||||||
|
include 'quadrature.h'
|
||||||
|
|
||||||
|
! Input variables
|
||||||
|
|
||||||
|
logical,intent(in) :: dotest
|
||||||
|
|
||||||
|
logical,intent(in) :: doACFDT
|
||||||
|
logical,intent(in) :: exchange_kernel
|
||||||
|
logical,intent(in) :: doXBS
|
||||||
|
logical,intent(in) :: dophBSE
|
||||||
|
logical,intent(in) :: dophBSE2
|
||||||
|
logical,intent(in) :: doppBSE
|
||||||
|
logical,intent(in) :: TDA_W
|
||||||
|
logical,intent(in) :: TDA
|
||||||
|
logical,intent(in) :: dBSE
|
||||||
|
logical,intent(in) :: dTDA
|
||||||
|
logical,intent(in) :: singlet
|
||||||
|
logical,intent(in) :: triplet
|
||||||
|
logical,intent(in) :: linearize
|
||||||
|
double precision,intent(in) :: eta
|
||||||
|
logical,intent(in) :: doSRG
|
||||||
|
|
||||||
|
integer,intent(in) :: nBas
|
||||||
|
integer,intent(in) :: nOrb
|
||||||
|
integer,intent(in) :: nC
|
||||||
|
integer,intent(in) :: nO
|
||||||
|
integer,intent(in) :: nV
|
||||||
|
integer,intent(in) :: nR
|
||||||
|
integer,intent(in) :: nS
|
||||||
|
double precision,intent(in) :: ENuc
|
||||||
|
double precision,intent(in) :: ERHF
|
||||||
|
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
|
||||||
|
double precision,intent(in) :: CAP(nOrb,nOrb)
|
||||||
|
double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart)
|
||||||
|
double precision,intent(in) :: eHF(nOrb)
|
||||||
|
|
||||||
|
! Local variables
|
||||||
|
|
||||||
|
logical :: print_W = .false.
|
||||||
|
logical :: plot_self = .false.
|
||||||
|
logical :: dRPA_W
|
||||||
|
integer :: isp_W
|
||||||
|
integer :: p
|
||||||
|
double precision :: flow
|
||||||
|
double precision :: EcRPA
|
||||||
|
double precision :: EcBSE(nspin)
|
||||||
|
double precision :: EcGM
|
||||||
|
double precision,allocatable :: Aph(:,:)
|
||||||
|
double precision,allocatable :: Bph(:,:)
|
||||||
|
double precision,allocatable :: Re_SigC(:)
|
||||||
|
double precision,allocatable :: Im_SigC(:)
|
||||||
|
double precision,allocatable :: Re_Z(:)
|
||||||
|
double precision,allocatable :: Im_Z(:)
|
||||||
|
double precision,allocatable :: Om(:)
|
||||||
|
double precision,allocatable :: XpY(:,:)
|
||||||
|
double precision,allocatable :: XmY(:,:)
|
||||||
|
double precision,allocatable :: rho(:,:,:)
|
||||||
|
|
||||||
|
|
||||||
|
double precision,allocatable :: Re_eGWlin(:)
|
||||||
|
double precision, allocatable :: Im_eGWlin(:)
|
||||||
|
double precision,allocatable :: Re_eGW(:)
|
||||||
|
double precision,allocatable :: Im_eGW(:)
|
||||||
|
double precision, allocatable :: e_cap(:)
|
||||||
|
|
||||||
|
! Hello world
|
||||||
|
|
||||||
|
write(*,*)
|
||||||
|
write(*,*)'***************************************'
|
||||||
|
write(*,*)'* Restricted complex G0W0 Calculation *'
|
||||||
|
write(*,*)'***************************************'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
! Spin manifold and TDA for dynamical screening
|
||||||
|
|
||||||
|
isp_W = 1
|
||||||
|
dRPA_W = .true.
|
||||||
|
|
||||||
|
if(TDA_W) then
|
||||||
|
write(*,*) 'Tamm-Dancoff approximation for dynamical screening!'
|
||||||
|
write(*,*)
|
||||||
|
end if
|
||||||
|
|
||||||
|
! SRG regularization
|
||||||
|
|
||||||
|
flow = 500d0
|
||||||
|
|
||||||
|
if(doSRG) then
|
||||||
|
! Not implemented
|
||||||
|
write(*,*) '*** SRG regularized G0W0 scheme ***'
|
||||||
|
write(*,*) '!!! No SRG with cRG0W0 !!!'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
! Memory allocation
|
||||||
|
|
||||||
|
allocate(Aph(nS,nS),Bph(nS,nS),Re_SigC(nOrb),Im_SigC(nOrb),Re_Z(nOrb),Im_Z(nOrb),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nOrb,nOrb,nS), &
|
||||||
|
Re_eGW(nOrb),Im_eGW(nOrb),Re_eGWlin(nOrb),Im_eGWlin(nOrb),e_cap(nOrb))
|
||||||
|
do p = 1, nOrb
|
||||||
|
e_cap(p) = CAP(p,p)
|
||||||
|
end do
|
||||||
|
!-------------------!
|
||||||
|
! Compute screening !
|
||||||
|
!-------------------!
|
||||||
|
|
||||||
|
call complex_phRLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
|
||||||
|
if(.not.TDA_W) call complex_phRLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
|
||||||
|
|
||||||
|
call complex_phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
|
||||||
|
|
||||||
|
!if(print_W) call print_excitation_energies('phRPA@RHF','singlet',nS,Om)
|
||||||
|
|
||||||
|
!--------------------------!
|
||||||
|
! Compute spectral weights !
|
||||||
|
!--------------------------!
|
||||||
|
|
||||||
|
call complex_RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho)
|
||||||
|
|
||||||
|
!------------------------!
|
||||||
|
! Compute GW self-energy !
|
||||||
|
!------------------------!
|
||||||
|
|
||||||
|
!!!! STOPPED HERE PROCEED WITH IMPLEMENTING COMPLEX_CRGW_SELF_ENERGY_DIAG
|
||||||
|
!!!!
|
||||||
|
call complex_cRGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,Re_SigC,Im_SigC,Re_Z,Im_Z,e_cap)
|
||||||
|
|
||||||
|
!-----------------------------------!
|
||||||
|
! Solve the quasi-particle equation !
|
||||||
|
!-----------------------------------!
|
||||||
|
|
||||||
|
! Linearized or graphical solution?
|
||||||
|
Re_eGWlin(:) = eHF(:) + Re_Z(:)*Re_SigC(:) - Im_Z(:)*Im_SigC(:)
|
||||||
|
Im_eGWlin(:) = e_cap(:) + Re_Z(:)*Im_SigC(:) + Im_Z(:)*Re_SigC(:)
|
||||||
|
|
||||||
|
if(linearize) then
|
||||||
|
|
||||||
|
write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
Re_eGW(:) = Re_eGWlin(:)
|
||||||
|
Im_eGW(:) = Im_eGWlin(:)
|
||||||
|
|
||||||
|
else
|
||||||
|
|
||||||
|
write(*,*) ' *** Quasiparticle energies obtained by root search *** '
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
call cRGW_QP_graph(doSRG,eta,flow,nOrb,nC,nO,nV,nR,nS,eHF,e_cap,Om,rho,Re_eGWlin,Im_eGWlin,eHF,e_cap,Re_eGW,Im_eGW,Re_Z,Im_Z)
|
||||||
|
end if
|
||||||
|
|
||||||
|
! Plot self-energy, renormalization factor, and spectral function
|
||||||
|
|
||||||
|
if(plot_self) call RGW_plot_self_energy(nOrb,eta,nC,nO,nV,nR,nS,eHF,eHF,Om,rho)
|
||||||
|
|
||||||
|
! Cumulant expansion
|
||||||
|
|
||||||
|
! call RGWC(dotest,eta,nOrb,nC,nO,nV,nR,nS,Om,rho,eHF,eHF,eGW,Z)
|
||||||
|
|
||||||
|
! Compute the RPA correlation energy
|
||||||
|
|
||||||
|
call phRLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,Re_eGW,ERI,Aph)
|
||||||
|
if(.not.TDA_W) call phRLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
|
||||||
|
|
||||||
|
call phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
|
||||||
|
|
||||||
|
!--------------!
|
||||||
|
! Dump results !
|
||||||
|
!--------------!
|
||||||
|
|
||||||
|
call print_cRG0W0(nOrb,nO,eHF,ENuc,ERHF,Re_SigC,Im_SigC,Re_Z,Im_Z,Re_eGW,Im_eGW,EcRPA,EcGM,CAP)
|
||||||
|
!---------------------------!
|
||||||
|
! Perform phBSE calculation !
|
||||||
|
!---------------------------!
|
||||||
|
!
|
||||||
|
! if(dophBSE) then
|
||||||
|
!
|
||||||
|
! call RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta, &
|
||||||
|
! nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,Re_eGW,EcBSE)
|
||||||
|
!
|
||||||
|
! write(*,*)
|
||||||
|
! write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF correlation energy (singlet) = ',EcBSE(1),' au'
|
||||||
|
! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF correlation energy (triplet) = ',EcBSE(2),' au'
|
||||||
|
! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF correlation energy = ',sum(EcBSE),' au'
|
||||||
|
! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au'
|
||||||
|
! write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
! write(*,*)
|
||||||
|
!
|
||||||
|
! ! Compute the BSE correlation energy via the adiabatic connection fluctuation dissipation theorem
|
||||||
|
!
|
||||||
|
! if(doACFDT) then
|
||||||
|
!
|
||||||
|
! call RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,eHF,Re_eGW,EcBSE)
|
||||||
|
!
|
||||||
|
! write(*,*)
|
||||||
|
! write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF correlation energy (singlet) = ',EcBSE(1),' au'
|
||||||
|
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF correlation energy (triplet) = ',EcBSE(2),' au'
|
||||||
|
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF correlation energy = ',sum(EcBSE),' au'
|
||||||
|
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au'
|
||||||
|
! write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
! write(*,*)
|
||||||
|
!
|
||||||
|
! end if
|
||||||
|
!
|
||||||
|
! end if
|
||||||
|
!
|
||||||
|
!!---------------------------!
|
||||||
|
!! Perform ppBSE calculation !
|
||||||
|
!!---------------------------!
|
||||||
|
!
|
||||||
|
! if(doppBSE) then
|
||||||
|
!
|
||||||
|
! call RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,Re_eGW,EcBSE)
|
||||||
|
!
|
||||||
|
! write(*,*)
|
||||||
|
! write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF correlation energy (singlet) = ',EcBSE(1),' au'
|
||||||
|
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF correlation energy (triplet) = ',EcBSE(2),' au'
|
||||||
|
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF correlation energy = ',sum(EcBSE),' au'
|
||||||
|
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au'
|
||||||
|
! write(*,*)'-------------------------------------------------------------------------------'
|
||||||
|
! write(*,*)
|
||||||
|
!
|
||||||
|
! end if
|
||||||
|
!
|
||||||
|
!! Testing zone
|
||||||
|
!
|
||||||
|
! if(dotest) then
|
||||||
|
!
|
||||||
|
! call dump_test_value('R','G0W0 correlation energy',EcRPA)
|
||||||
|
! call dump_test_value('R','G0W0 HOMO energy',Re_eGW(nO))
|
||||||
|
! call dump_test_value('R','G0W0 LUMO energy',Re_eGW(nO+1))
|
||||||
|
!
|
||||||
|
! end if
|
||||||
|
!
|
||||||
|
end subroutine
|
||||||
80
src/GW/complex_cRGW_excitation_density.f90
Normal file
80
src/GW/complex_cRGW_excitation_density.f90
Normal file
@ -0,0 +1,80 @@
|
|||||||
|
subroutine RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho)
|
||||||
|
|
||||||
|
! Compute excitation densities
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
! Input variables
|
||||||
|
|
||||||
|
integer,intent(in) :: nOrb
|
||||||
|
integer,intent(in) :: nC
|
||||||
|
integer,intent(in) :: nO
|
||||||
|
integer,intent(in) :: nR
|
||||||
|
integer,intent(in) :: nS
|
||||||
|
complex*16,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
|
||||||
|
complex*16,intent(in) :: XpY(nS,nS)
|
||||||
|
|
||||||
|
! Local variables
|
||||||
|
|
||||||
|
integer :: ia,jb,p,q,j,b
|
||||||
|
complex*16, allocatable :: tmp(:,:,:)
|
||||||
|
|
||||||
|
! Output variables
|
||||||
|
|
||||||
|
complex*16,intent(out) :: rho(nOrb,nOrb,nS)
|
||||||
|
|
||||||
|
if(nOrb .lt. 256) then
|
||||||
|
|
||||||
|
allocate(tmp(nOrb,nOrb,nS))
|
||||||
|
|
||||||
|
!$OMP PARALLEL DEFAULT(NONE) &
|
||||||
|
!$OMP PRIVATE(p, q, j, b, jb) &
|
||||||
|
!$OMP SHARED(nOrb, nC, nO, nR, ERI, tmp)
|
||||||
|
!$OMP DO COLLAPSE(2)
|
||||||
|
do p = 1, nOrb
|
||||||
|
do q = 1, nOrb
|
||||||
|
jb = 0
|
||||||
|
do j = nC+1, nO
|
||||||
|
do b = nO+1, nOrb-nR
|
||||||
|
jb = jb + 1
|
||||||
|
tmp(p,q,jb) = ERI(p,j,q,b)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
!$OMP END DO
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
|
call zgemm("N", "T", nOrb*nOrb, nS, nS, 1.d0, &
|
||||||
|
tmp(1,1,1), nOrb*nOrb, XpY(1,1), nS, &
|
||||||
|
0.d0, rho(1,1,1), nOrb*nOrb)
|
||||||
|
|
||||||
|
deallocate(tmp)
|
||||||
|
|
||||||
|
else
|
||||||
|
|
||||||
|
rho(:,:,:) = 0d0
|
||||||
|
!$OMP PARALLEL &
|
||||||
|
!$OMP SHARED(nC,nOrb,nR,nO,nS,rho,ERI,XpY) &
|
||||||
|
!$OMP PRIVATE(q,p,jb,ia) &
|
||||||
|
!$OMP DEFAULT(NONE)
|
||||||
|
!$OMP DO
|
||||||
|
do q=nC+1,nOrb-nR
|
||||||
|
do p=nC+1,nOrb-nR
|
||||||
|
jb = 0
|
||||||
|
do j=nC+1,nO
|
||||||
|
do b=nO+1,nOrb-nR
|
||||||
|
jb = jb + 1
|
||||||
|
do ia=1,nS
|
||||||
|
rho(p,q,ia) = rho(p,q,ia) + ERI(p,j,q,b)*XpY(ia,jb)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
!$OMP END DO
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
|
endif
|
||||||
|
|
||||||
|
end subroutine
|
||||||
65
src/LR/complex_phRLR.f90
Normal file
65
src/LR/complex_phRLR.f90
Normal file
@ -0,0 +1,65 @@
|
|||||||
|
subroutine phRLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
|
||||||
|
|
||||||
|
! Compute linear response
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
include 'parameters.h'
|
||||||
|
|
||||||
|
! Input variables
|
||||||
|
|
||||||
|
logical,intent(in) :: TDA
|
||||||
|
integer,intent(in) :: nS
|
||||||
|
complex*16,intent(in) :: Aph(nS,nS)
|
||||||
|
complex*16,intent(in) :: Bph(nS,nS)
|
||||||
|
|
||||||
|
! Local variables
|
||||||
|
|
||||||
|
complex*16 :: complex_trace_matrix
|
||||||
|
complex*16 :: t1, t2
|
||||||
|
complex*16,allocatable :: RPA_matrix(:,:)
|
||||||
|
complex*16,allocatable :: Z(:,:)
|
||||||
|
complex*16,allocatable :: tmp(:,:)
|
||||||
|
complex*16,allocatable :: OmOmminus(:)
|
||||||
|
|
||||||
|
! Output variables
|
||||||
|
|
||||||
|
complex*16,intent(out) :: EcRPA
|
||||||
|
complex*16,intent(out) :: Om(nS)
|
||||||
|
complex*16,intent(out) :: XpY(nS,nS)
|
||||||
|
complex*16,intent(out) :: XmY(nS,nS)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
! Tamm-Dancoff approximation
|
||||||
|
|
||||||
|
if(TDA) then
|
||||||
|
|
||||||
|
XpY(:,:) = Aph(:,:)
|
||||||
|
call complex_diagonalize_matrix(nS,XpY,Om)
|
||||||
|
XpY(:,:) = transpose(XpY(:,:))
|
||||||
|
XmY(:,:) = XpY(:,:)
|
||||||
|
|
||||||
|
else
|
||||||
|
|
||||||
|
allocate(RPA_matrix(2*nS,2*nS),OmOmminus(2*nS))
|
||||||
|
|
||||||
|
call complex_diagonalize_matrix(2*nS,RPA_matrix,Om)
|
||||||
|
Om(:) = OmOmminus(1:nS)
|
||||||
|
call complex_vecout(nS,Om)
|
||||||
|
call complex_vecout(nS,Om(nS+1:2*nS))
|
||||||
|
if (maxval(abs(Om(:)+Om(nS+1:2*nS))) > 1e10) &
|
||||||
|
call print_warning('We dont find a Om and -Om structure as solution of the RPA. There might be a problem somewhere.')
|
||||||
|
if(minval(abs(Om(:))) < 0d0) &
|
||||||
|
call print_warning('You may have instabilities in linear response: A-B is not positive definite!!')
|
||||||
|
XpY(:,:) = RPA_matrix(1:nS,1:nS) + RPA_matrix(nS+1:2*nS,nS+1:2*nS)
|
||||||
|
XmY(:,:) = RPA_matrix(1:nS,1:nS) - RPA_matrix(nS+1:2*nS,nS+1:2*nS)
|
||||||
|
|
||||||
|
deallocate(RPA_matrix)
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
! Compute the RPA correlation energy
|
||||||
|
|
||||||
|
EcRPA = 0.5d0*(sum(Om) - complex_trace_matrix(nS,Aph))
|
||||||
|
|
||||||
|
end subroutine
|
||||||
109
src/LR/complex_phRLR_A.f90
Normal file
109
src/LR/complex_phRLR_A.f90
Normal file
@ -0,0 +1,109 @@
|
|||||||
|
subroutine complex_phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
|
||||||
|
|
||||||
|
! Compute resonant block of the ph channel
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
include 'parameters.h'
|
||||||
|
|
||||||
|
! Input variables
|
||||||
|
|
||||||
|
logical,intent(in) :: dRPA
|
||||||
|
integer,intent(in) :: ispin
|
||||||
|
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,intent(in) :: lambda
|
||||||
|
complex*16,intent(in) :: e(nBas)
|
||||||
|
complex*16,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||||
|
|
||||||
|
! Local variables
|
||||||
|
|
||||||
|
double precision :: delta_dRPA
|
||||||
|
|
||||||
|
integer :: i,j,a,b,ia,jb
|
||||||
|
integer :: nn,jb0
|
||||||
|
logical :: i_eq_j
|
||||||
|
complex*16 :: ct1,ct2
|
||||||
|
|
||||||
|
! Output variables
|
||||||
|
|
||||||
|
complex*16,intent(out) :: Aph(nS,nS)
|
||||||
|
|
||||||
|
! Direct RPA
|
||||||
|
|
||||||
|
delta_dRPA = 0d0
|
||||||
|
if(dRPA) delta_dRPA = 1d0
|
||||||
|
|
||||||
|
! Build A matrix for single manifold
|
||||||
|
|
||||||
|
if(ispin == 1) then
|
||||||
|
|
||||||
|
nn = nBas - nR - nO
|
||||||
|
ct1 = 2d0 * lambda
|
||||||
|
ct2 = - (1d0 - delta_dRPA) * lambda
|
||||||
|
!$OMP PARALLEL DEFAULT(NONE) &
|
||||||
|
!$OMP PRIVATE (i, a, j, b, i_eq_j, ia, jb0, jb) &
|
||||||
|
!$OMP SHARED (nC, nO, nR, nBas, nn, ct1, ct2, e, ERI, Aph)
|
||||||
|
!$OMP DO COLLAPSE(2)
|
||||||
|
do i = nC+1, nO
|
||||||
|
do a = nO+1, nBas-nR
|
||||||
|
ia = a - nO + (i - nC - 1) * nn
|
||||||
|
|
||||||
|
do j = nC+1, nO
|
||||||
|
i_eq_j = i == j
|
||||||
|
jb0 = (j - nC - 1) * nn - nO
|
||||||
|
|
||||||
|
do b = nO+1, nBas-nR
|
||||||
|
jb = b + jb0
|
||||||
|
|
||||||
|
Aph(ia,jb) = ct1 * ERI(b,i,j,a) + ct2 * ERI(b,j,a,i)
|
||||||
|
if(i_eq_j) then
|
||||||
|
if(a == b) Aph(ia,jb) = Aph(ia,jb) + e(a) - e(i)
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
!$OMP END DO
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
! Build A matrix for triplet manifold
|
||||||
|
|
||||||
|
if(ispin == 2) then
|
||||||
|
|
||||||
|
nn = nBas - nR - nO
|
||||||
|
ct2 = - (1d0 - delta_dRPA) * lambda
|
||||||
|
!$OMP PARALLEL DEFAULT(NONE) &
|
||||||
|
!$OMP PRIVATE (i, a, j, b, i_eq_j, ia, jb0, jb) &
|
||||||
|
!$OMP SHARED (nC, nO, nR, nBas, nn, ct2, e, ERI, Aph)
|
||||||
|
!$OMP DO COLLAPSE(2)
|
||||||
|
do i = nC+1, nO
|
||||||
|
do a = nO+1, nBas-nR
|
||||||
|
ia = a - nO + (i - nC - 1) * nn
|
||||||
|
|
||||||
|
do j = nC+1, nO
|
||||||
|
i_eq_j = i == j
|
||||||
|
jb0 = (j - nC - 1) * nn - nO
|
||||||
|
|
||||||
|
do b = nO+1, nBas-nR
|
||||||
|
jb = b + jb0
|
||||||
|
|
||||||
|
Aph(ia,jb) = ct2 * ERI(b,j,a,i)
|
||||||
|
if(i_eq_j) then
|
||||||
|
if(a == b) Aph(ia,jb) = Aph(ia,jb) + e(a) - e(i)
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
!$OMP END DO
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
end subroutine
|
||||||
121
src/LR/complex_phRLR_B.f90
Normal file
121
src/LR/complex_phRLR_B.f90
Normal file
@ -0,0 +1,121 @@
|
|||||||
|
subroutine complex_phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
|
||||||
|
|
||||||
|
! Compute the coupling block of the ph channel
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
include 'parameters.h'
|
||||||
|
|
||||||
|
! Input variables
|
||||||
|
|
||||||
|
logical,intent(in) :: dRPA
|
||||||
|
integer,intent(in) :: ispin,nBas,nC,nO,nV,nR,nS
|
||||||
|
double precision,intent(in) :: lambda
|
||||||
|
complex*16,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||||
|
|
||||||
|
! Local variables
|
||||||
|
|
||||||
|
double precision :: delta_dRPA
|
||||||
|
|
||||||
|
integer :: i,j,a,b,ia,jb
|
||||||
|
integer :: nn,jb0
|
||||||
|
complex*16 :: ct1,ct2
|
||||||
|
|
||||||
|
! Output variables
|
||||||
|
|
||||||
|
complex*16,intent(out) :: Bph(nS,nS)
|
||||||
|
|
||||||
|
! Direct RPA
|
||||||
|
|
||||||
|
delta_dRPA = 0d0
|
||||||
|
if(dRPA) delta_dRPA = 1d0
|
||||||
|
|
||||||
|
! Build B matrix for singlet manifold
|
||||||
|
|
||||||
|
if(ispin == 1) then
|
||||||
|
|
||||||
|
nn = nBas - nR - nO
|
||||||
|
ct1 = 2d0 * lambda
|
||||||
|
ct2 = - (1d0 - delta_dRPA) * lambda
|
||||||
|
!$OMP PARALLEL DEFAULT(NONE) &
|
||||||
|
!$OMP PRIVATE (i, a, j, b, ia, jb0, jb) &
|
||||||
|
!$OMP SHARED (nC, nO, nR, nBas, nn, ct1, ct2, ERI, Bph)
|
||||||
|
!$OMP DO COLLAPSE(2)
|
||||||
|
do i = nC+1, nO
|
||||||
|
do a = nO+1, nBas-nR
|
||||||
|
ia = a - nO + (i - nC - 1) * nn
|
||||||
|
|
||||||
|
do j = nC+1, nO
|
||||||
|
jb0 = (j - nC - 1) * nn - nO
|
||||||
|
|
||||||
|
do b = nO+1, nBas-nR
|
||||||
|
jb = b + jb0
|
||||||
|
|
||||||
|
Bph(ia,jb) = ct1 * ERI(b,i,j,a) + ct2 * ERI(b,j,i,a)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
!$OMP END DO
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
|
!ia = 0
|
||||||
|
!do i=nC+1,nO
|
||||||
|
! do a=nO+1,nBas-nR
|
||||||
|
! ia = ia + 1
|
||||||
|
! jb = 0
|
||||||
|
! do j=nC+1,nO
|
||||||
|
! do b=nO+1,nBas-nR
|
||||||
|
! jb = jb + 1
|
||||||
|
! Bph(ia,jb) = 2d0*lambda*ERI(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a)
|
||||||
|
! end do
|
||||||
|
! end do
|
||||||
|
! end do
|
||||||
|
!end do
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
! Build B matrix for triplet manifold
|
||||||
|
|
||||||
|
if(ispin == 2) then
|
||||||
|
|
||||||
|
nn = nBas - nR - nO
|
||||||
|
ct2 = - (1d0 - delta_dRPA) * lambda
|
||||||
|
!$OMP PARALLEL DEFAULT(NONE) &
|
||||||
|
!$OMP PRIVATE (i, a, j, b, ia, jb0, jb) &
|
||||||
|
!$OMP SHARED (nC, nO, nR, nBas, nn, ct2, ERI, Bph)
|
||||||
|
!$OMP DO COLLAPSE(2)
|
||||||
|
do i = nC+1, nO
|
||||||
|
do a = nO+1, nBas-nR
|
||||||
|
ia = a - nO + (i - nC - 1) * nn
|
||||||
|
|
||||||
|
do j = nC+1, nO
|
||||||
|
jb0 = (j - nC - 1) * nn - nO
|
||||||
|
|
||||||
|
do b = nO+1, nBas-nR
|
||||||
|
jb = b + jb0
|
||||||
|
|
||||||
|
Bph(ia,jb) = ct2 * ERI(b,j,i,a)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
!$OMP END DO
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
|
! ia = 0
|
||||||
|
! do i=nC+1,nO
|
||||||
|
! do a=nO+1,nBas-nR
|
||||||
|
! ia = ia + 1
|
||||||
|
! jb = 0
|
||||||
|
! do j=nC+1,nO
|
||||||
|
! do b=nO+1,nBas-nR
|
||||||
|
! jb = jb + 1
|
||||||
|
! Bph(ia,jb) = - (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a)
|
||||||
|
! end do
|
||||||
|
! end do
|
||||||
|
! end do
|
||||||
|
! end do
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
end subroutine
|
||||||
Loading…
x
Reference in New Issue
Block a user