mirror of
https://github.com/pfloos/quack
synced 2025-01-03 18:16:03 +01:00
UG0W0 OK
This commit is contained in:
parent
1d3e156838
commit
282cbcb517
@ -13,6 +13,6 @@
|
|||||||
# ACFDT: AC Kx XBS
|
# ACFDT: AC Kx XBS
|
||||||
F F T
|
F F T
|
||||||
# BSE: BSE dBSE dTDA evDyn
|
# BSE: BSE dBSE dTDA evDyn
|
||||||
T T F F
|
F T F F
|
||||||
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
|
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
|
||||||
1000000 100000 10 0.3 10000 1234 T
|
1000000 100000 10 0.3 10000 1234 T
|
||||||
|
|||||||
@ -1,6 +1,6 @@
|
|||||||
subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
|
subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
|
||||||
dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linearize,eta, &
|
dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linearize,eta, &
|
||||||
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eGW)
|
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI,PHF,cHF,eHF,eGW)
|
||||||
|
|
||||||
! Perform G0W0 calculation
|
! Perform G0W0 calculation
|
||||||
|
|
||||||
@ -33,7 +33,6 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,
|
|||||||
double precision,intent(in) :: cHF(nBas,nBas)
|
double precision,intent(in) :: cHF(nBas,nBas)
|
||||||
double precision,intent(in) :: PHF(nBas,nBas)
|
double precision,intent(in) :: PHF(nBas,nBas)
|
||||||
double precision,intent(in) :: Hc(nBas,nBas)
|
double precision,intent(in) :: Hc(nBas,nBas)
|
||||||
double precision,intent(in) :: H(nBas,nBas)
|
|
||||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||||
|
|
||||||
! Local variables
|
! Local variables
|
||||||
|
|||||||
@ -4,7 +4,7 @@ program QuAcK
|
|||||||
include 'parameters.h'
|
include 'parameters.h'
|
||||||
|
|
||||||
logical :: doSph
|
logical :: doSph
|
||||||
logical :: unrestricted
|
logical :: unrestricted = .false.
|
||||||
logical :: doRHF,doUHF,doMOM
|
logical :: doRHF,doUHF,doMOM
|
||||||
logical :: doMP2,doMP3,doMP2F12
|
logical :: doMP2,doMP3,doMP2F12
|
||||||
logical :: doCCD,doCCSD,doCCSDT
|
logical :: doCCD,doCCSD,doCCSDT
|
||||||
|
|||||||
@ -51,16 +51,14 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
|
|||||||
integer :: nSa
|
integer :: nSa
|
||||||
integer :: nSb
|
integer :: nSb
|
||||||
integer :: nSt
|
integer :: nSt
|
||||||
double precision :: EcRPA(nspin)
|
double precision :: EcRPA
|
||||||
double precision :: EcBSE(nspin)
|
double precision :: EcBSE
|
||||||
double precision :: EcAC(nspin)
|
double precision :: EcAC(nspin)
|
||||||
double precision,allocatable :: SigC(:,:)
|
double precision,allocatable :: SigC(:,:)
|
||||||
double precision,allocatable :: Z(:,:)
|
double precision,allocatable :: Z(:,:)
|
||||||
double precision,allocatable :: Omega(:)
|
double precision,allocatable :: Omega(:)
|
||||||
double precision,allocatable :: XpY_a(:,:)
|
double precision,allocatable :: XpY(:,:)
|
||||||
double precision,allocatable :: XpY_b(:,:)
|
double precision,allocatable :: XmY(:,:)
|
||||||
double precision,allocatable :: XmY_a(:,:)
|
|
||||||
double precision,allocatable :: XmY_b(:,:)
|
|
||||||
double precision,allocatable :: rho(:,:,:,:)
|
double precision,allocatable :: rho(:,:,:,:)
|
||||||
|
|
||||||
double precision,allocatable :: eGWlin(:,:)
|
double precision,allocatable :: eGWlin(:,:)
|
||||||
@ -80,7 +78,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
|
|||||||
|
|
||||||
! Initialization
|
! Initialization
|
||||||
|
|
||||||
EcRPA(:) = 0d0
|
EcRPA = 0d0
|
||||||
|
|
||||||
! COHSEX approximation
|
! COHSEX approximation
|
||||||
|
|
||||||
@ -103,53 +101,45 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
|
|||||||
nSb = nS(2)
|
nSb = nS(2)
|
||||||
nSt = nSa + nSb
|
nSt = nSa + nSb
|
||||||
|
|
||||||
allocate(SigC(nBas,nspin),Z(nBas,nspin),Omega(nSt),XpY_a(nSa,nSa),XpY_b(nSb,nSb),XmY_a(nSa,nSa),XmY_b(nSb,nSb), &
|
allocate(SigC(nBas,nspin),Z(nBas,nspin),Omega(nSt),XpY(nSt,nSt),XmY(nSt,nSt), &
|
||||||
rho(nBas,nBas,nSt,nspin),eGWlin(nBas,nspin))
|
rho(nBas,nBas,nSt,nspin),eGWlin(nBas,nspin))
|
||||||
|
|
||||||
! Compute linear response
|
!-------------------!
|
||||||
|
! Compute screening !
|
||||||
|
!-------------------!
|
||||||
|
|
||||||
!----------------------------------------------
|
! Spin-conserving transition
|
||||||
! alpha-alpha block
|
|
||||||
!----------------------------------------------
|
|
||||||
|
|
||||||
ispin = 1
|
ispin = 1
|
||||||
iblock = 3
|
|
||||||
|
|
||||||
call linear_response(iblock,.true.,TDA_W,.false.,eta,nBas,nC(ispin),nO(ispin),nV(ispin),nR(ispin),nSa,1d0, &
|
call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0, &
|
||||||
eHF(:,ispin),ERI_aa,rho(:,:,1:nSa,ispin),EcRPA(ispin),Omega(1:nSa),XpY_a,XmY_a)
|
eHF,ERI_aa,ERI_ab,ERI_bb,rho(:,:,:,ispin),EcRPA,Omega,XpY,XmY)
|
||||||
|
|
||||||
if(print_W) call print_excitation('RPA@HF alpha',iblock,nSa,Omega(1:nSa))
|
if(print_W) call print_excitation('RPA@UHF',3,nSt,Omega)
|
||||||
|
|
||||||
!----------------------------------------------
|
!----------------------!
|
||||||
! alpha-beta block
|
! Excitation densities !
|
||||||
!----------------------------------------------
|
!----------------------!
|
||||||
|
|
||||||
ispin = 2
|
|
||||||
iblock = 3
|
|
||||||
|
|
||||||
call linear_response(iblock,.true.,TDA_W,.false.,eta,nBas,nC(ispin),nO(ispin),nV(ispin),nR(ispin),nSb,1d0, &
|
|
||||||
eHF(:,ispin),ERI_bb,rho(:,:,nSa+1:nSt,ispin),EcRPA(ispin),Omega(nSa+1:nSt),XpY_b,XmY_b)
|
|
||||||
|
|
||||||
if(print_W) call print_excitation('RPA@HF beta ',iblock,nSb,Omega(nSa+1:nSt))
|
|
||||||
|
|
||||||
!----------------------------------------------
|
|
||||||
! Excitation densities for alpha self-energy
|
|
||||||
!----------------------------------------------
|
|
||||||
|
|
||||||
call unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_ab,ERI_bb,XpY_a,XpY_b,rho)
|
call unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_ab,ERI_bb,XpY,rho)
|
||||||
|
|
||||||
!----------------------
|
!---------------------!
|
||||||
! Compute self-energy
|
! Compute self-energy !
|
||||||
!----------------------
|
!---------------------!
|
||||||
|
|
||||||
call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,eHF,Omega,rho,SigC)
|
call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,eHF,Omega,rho,SigC)
|
||||||
|
|
||||||
! Compute renormalization factor
|
!--------------------------------!
|
||||||
|
! Compute renormalization factor !
|
||||||
|
!--------------------------------!
|
||||||
|
|
||||||
! call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, &
|
! call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, &
|
||||||
! Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z(:))
|
! Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z(:))
|
||||||
|
|
||||||
! Solve the quasi-particle equation
|
!-----------------------------------!
|
||||||
|
! Solve the quasi-particle equation !
|
||||||
|
!-----------------------------------!
|
||||||
|
|
||||||
Z(:,:) = 1d0
|
Z(:,:) = 1d0
|
||||||
eGWlin(:,:) = eHF(:,:) + Z(:,:)*SigC(:,:)
|
eGWlin(:,:) = eHF(:,:) + Z(:,:)*SigC(:,:)
|
||||||
|
|
||||||
@ -174,7 +164,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
|
|||||||
! Dump results
|
! Dump results
|
||||||
|
|
||||||
do ispin=1,nspin
|
do ispin=1,nspin
|
||||||
call print_G0W0(nBas,nO(ispin),eHF(:,ispin),ENuc,EUHF,SigC(:,ispin),Z(:,ispin),eGW(:,ispin),EcRPA(ispin))
|
call print_G0W0(nBas,nO(ispin),eHF(:,ispin),ENuc,EUHF,SigC(:,ispin),Z(:,ispin),eGW(:,ispin),EcRPA)
|
||||||
end do
|
end do
|
||||||
|
|
||||||
! Compute the RPA correlation energy
|
! Compute the RPA correlation energy
|
||||||
@ -184,10 +174,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
|
|||||||
|
|
||||||
write(*,*)
|
write(*,*)
|
||||||
write(*,*)'-------------------------------------------------------------------------------'
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy (singlet) =',EcRPA(1)
|
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy =',EcRPA
|
||||||
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy (triplet) =',EcRPA(2)
|
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 total energy =',ENuc + EUHF + EcRPA
|
||||||
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy =',EcRPA(1) + EcRPA(2)
|
|
||||||
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 total energy =',ENuc + EUHF + EcRPA(1) + EcRPA(2)
|
|
||||||
write(*,*)'-------------------------------------------------------------------------------'
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
write(*,*)
|
write(*,*)
|
||||||
|
|
||||||
|
|||||||
@ -40,11 +40,11 @@ subroutine print_G0W0(nBas,nO,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA,EcGM)
|
|||||||
write(*,'(2X,A30,F15.6)') 'G0W0 LUMO energy (eV):',eGW(LUMO)*HaToeV
|
write(*,'(2X,A30,F15.6)') 'G0W0 LUMO energy (eV):',eGW(LUMO)*HaToeV
|
||||||
write(*,'(2X,A30,F15.6)') 'G0W0 HOMO-LUMO gap (eV):',Gap*HaToeV
|
write(*,'(2X,A30,F15.6)') 'G0W0 HOMO-LUMO gap (eV):',Gap*HaToeV
|
||||||
write(*,*)'-------------------------------------------------------------------------------'
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
! write(*,'(2X,A30,F15.6)') 'RPA@G0W0 total energy =',ENuc + EHF + EcRPA
|
write(*,'(2X,A30,F15.6)') 'RPA@HF total energy =',ENuc + EHF + EcRPA
|
||||||
! write(*,'(2X,A30,F15.6)') 'RPA@G0W0 correlation energy =',EcRPA
|
write(*,'(2X,A30,F15.6)') 'RPA@HF correlation energy =',EcRPA
|
||||||
! write(*,'(2X,A30,F15.6)') 'GM@G0W0 total energy =',ENuc + EHF + EcGM
|
write(*,'(2X,A30,F15.6)') 'GM@G0W0 total energy =',ENuc + EHF + EcGM
|
||||||
! write(*,'(2X,A30,F15.6)') 'GM@G0W0 correlation energy =',EcGM
|
write(*,'(2X,A30,F15.6)') 'GM@G0W0 correlation energy =',EcGM
|
||||||
! write(*,*)'-------------------------------------------------------------------------------'
|
write(*,*)'-------------------------------------------------------------------------------'
|
||||||
write(*,*)
|
write(*,*)
|
||||||
|
|
||||||
end subroutine print_G0W0
|
end subroutine print_G0W0
|
||||||
|
|||||||
@ -66,7 +66,7 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,O
|
|||||||
|
|
||||||
EcGM = 0d0
|
EcGM = 0d0
|
||||||
do i=nC+1,nO
|
do i=nC+1,nO
|
||||||
EcGM = EcGM + 0.5d0*SigC(i)
|
EcGM = EcGM - SigC(i)
|
||||||
end do
|
end do
|
||||||
|
|
||||||
!-----------------------------
|
!-----------------------------
|
||||||
@ -143,7 +143,7 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,O
|
|||||||
do a=nO+1,nBas-nR
|
do a=nO+1,nBas-nR
|
||||||
do jb=1,nS
|
do jb=1,nS
|
||||||
eps = e(a) - e(i) + Omega(jb)
|
eps = e(a) - e(i) + Omega(jb)
|
||||||
EcGM = EcGM - 2d0*rho(a,i,jb)*rho(a,i,jb)*eps/(eps**2 + eta**2)
|
EcGM = EcGM - 4d0*rho(a,i,jb)*rho(a,i,jb)*eps/(eps**2 + eta**2)
|
||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
|
|||||||
@ -1,4 +1,4 @@
|
|||||||
subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_ab,ERI_bb,XpY_a,XpY_b,rho)
|
subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_ab,ERI_bb,XpY,rho)
|
||||||
|
|
||||||
! Compute excitation densities for unrestricted reference
|
! Compute excitation densities for unrestricted reference
|
||||||
|
|
||||||
@ -17,8 +17,7 @@ subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_
|
|||||||
double precision,intent(in) :: ERI_aa(nBas,nBas,nBas,nBas)
|
double precision,intent(in) :: ERI_aa(nBas,nBas,nBas,nBas)
|
||||||
double precision,intent(in) :: ERI_ab(nBas,nBas,nBas,nBas)
|
double precision,intent(in) :: ERI_ab(nBas,nBas,nBas,nBas)
|
||||||
double precision,intent(in) :: ERI_bb(nBas,nBas,nBas,nBas)
|
double precision,intent(in) :: ERI_bb(nBas,nBas,nBas,nBas)
|
||||||
double precision,intent(in) :: XpY_a(nSa,nSa)
|
double precision,intent(in) :: XpY(nSt,nSt)
|
||||||
double precision,intent(in) :: XpY_b(nSb,nSb)
|
|
||||||
|
|
||||||
! Local variables
|
! Local variables
|
||||||
|
|
||||||
@ -32,34 +31,34 @@ subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_
|
|||||||
|
|
||||||
rho(:,:,:,:) = 0d0
|
rho(:,:,:,:) = 0d0
|
||||||
|
|
||||||
!-------------
|
!-------------!
|
||||||
! alpha block
|
! alpha block !
|
||||||
!-------------
|
!-------------!
|
||||||
|
|
||||||
do p=nC(1)+1,nBas-nR(1)
|
do p=nC(1)+1,nBas-nR(1)
|
||||||
do q=nC(1)+1,nBas-nR(1)
|
do q=nC(1)+1,nBas-nR(1)
|
||||||
|
|
||||||
! Same-spin contribution
|
! Same-spin contribution
|
||||||
do ia=1,nSa
|
do ia=1,nSt
|
||||||
jb = 0
|
jb = 0
|
||||||
do j=nC(1)+1,nO(1)
|
do j=nC(1)+1,nO(1)
|
||||||
do b=nO(1)+1,nBas-nR(1)
|
do b=nO(1)+1,nBas-nR(1)
|
||||||
jb = jb + 1
|
jb = jb + 1
|
||||||
|
|
||||||
rho(p,q,ia,1) = rho(p,q,ia,1) + ERI_aa(p,j,q,b)*XpY_a(ia,jb)
|
rho(p,q,ia,1) = rho(p,q,ia,1) + ERI_aa(p,j,q,b)*XpY(ia,jb)
|
||||||
|
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
! Opposite-spin contribution
|
! Opposite-spin contribution
|
||||||
do ia=1,nSb
|
do ia=1,nSt
|
||||||
jb = 0
|
jb = nSa
|
||||||
do j=nC(2)+1,nO(2)
|
do j=nC(2)+1,nO(2)
|
||||||
do b=nO(2)+1,nBas-nR(2)
|
do b=nO(2)+1,nBas-nR(2)
|
||||||
jb = jb + 1
|
jb = jb + 1
|
||||||
|
|
||||||
rho(p,q,nSa+ia,1) = rho(p,q,nSa+ia,1) + ERI_ab(p,j,q,b)*XpY_b(ia,jb)
|
rho(p,q,ia,1) = rho(p,q,ia,1) + ERI_ab(p,j,q,b)*XpY(ia,jb)
|
||||||
|
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
@ -68,34 +67,34 @@ subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_
|
|||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
!------------
|
!------------!
|
||||||
! Beta block
|
! Beta block !
|
||||||
!------------
|
!------------!
|
||||||
|
|
||||||
do p=nC(2)+1,nBas-nR(2)
|
do p=nC(2)+1,nBas-nR(2)
|
||||||
do q=nC(2)+1,nBas-nR(2)
|
do q=nC(2)+1,nBas-nR(2)
|
||||||
|
|
||||||
! Same-spin contribution
|
|
||||||
do ia=1,nSb
|
|
||||||
jb = 0
|
|
||||||
do j=nC(2)+1,nO(2)
|
|
||||||
do b=nO(2)+1,nBas-nR(2)
|
|
||||||
jb = jb + 1
|
|
||||||
|
|
||||||
rho(p,q,ia,2) = rho(p,q,ia,2) + ERI_bb(p,j,q,b)*XpY_b(ia,jb)
|
|
||||||
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
! Opposite-spin contribution
|
! Opposite-spin contribution
|
||||||
do ia=1,nSa
|
do ia=1,nSt
|
||||||
jb = 0
|
jb = 0
|
||||||
do j=nC(1)+1,nO(1)
|
do j=nC(1)+1,nO(1)
|
||||||
do b=nO(1)+1,nBas-nR(1)
|
do b=nO(1)+1,nBas-nR(1)
|
||||||
jb = jb + 1
|
jb = jb + 1
|
||||||
|
|
||||||
rho(p,q,nSb+ia,2) = rho(p,q,nSb+ia,2) + ERI_ab(j,p,b,q)*XpY_a(ia,jb)
|
rho(p,q,ia,2) = rho(p,q,ia,2) + ERI_ab(j,p,b,q)*XpY(ia,jb)
|
||||||
|
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
! Same-spin contribution
|
||||||
|
do ia=1,nSt
|
||||||
|
jb = nSa
|
||||||
|
do j=nC(2)+1,nO(2)
|
||||||
|
do b=nO(2)+1,nBas-nR(2)
|
||||||
|
jb = jb + 1
|
||||||
|
|
||||||
|
rho(p,q,ia,2) = rho(p,q,ia,2) + ERI_bb(p,j,q,b)*XpY(ia,jb)
|
||||||
|
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|||||||
113
src/QuAcK/unrestricted_linear_response.f90
Normal file
113
src/QuAcK/unrestricted_linear_response.f90
Normal file
@ -0,0 +1,113 @@
|
|||||||
|
subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, &
|
||||||
|
e,ERI_aa,ERI_ab,ERI_bb,rho,EcRPA,Omega,XpY,XmY)
|
||||||
|
|
||||||
|
! Compute linear response for unrestricted formalism
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
include 'parameters.h'
|
||||||
|
|
||||||
|
! Input variables
|
||||||
|
|
||||||
|
integer,intent(in) :: ispin
|
||||||
|
logical,intent(in) :: dRPA
|
||||||
|
logical,intent(in) :: TDA
|
||||||
|
logical,intent(in) :: BSE
|
||||||
|
double precision,intent(in) :: eta
|
||||||
|
integer,intent(in) :: nBas
|
||||||
|
integer,intent(in) :: nC(nspin)
|
||||||
|
integer,intent(in) :: nO(nspin)
|
||||||
|
integer,intent(in) :: nV(nspin)
|
||||||
|
integer,intent(in) :: nR(nspin)
|
||||||
|
integer,intent(in) :: nSa
|
||||||
|
integer,intent(in) :: nSb
|
||||||
|
integer,intent(in) :: nSt
|
||||||
|
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_aa(nBas,nBas,nBas,nBas)
|
||||||
|
double precision,intent(in) :: ERI_ab(nBas,nBas,nBas,nBas)
|
||||||
|
double precision,intent(in) :: ERI_bb(nBas,nBas,nBas,nBas)
|
||||||
|
|
||||||
|
! Local variables
|
||||||
|
|
||||||
|
integer :: ia
|
||||||
|
double precision :: trace_matrix
|
||||||
|
double precision,allocatable :: A(:,:)
|
||||||
|
double precision,allocatable :: B(:,:)
|
||||||
|
double precision,allocatable :: ApB(:,:)
|
||||||
|
double precision,allocatable :: AmB(:,:)
|
||||||
|
double precision,allocatable :: AmBSq(:,:)
|
||||||
|
double precision,allocatable :: AmBIv(:,:)
|
||||||
|
double precision,allocatable :: Z(:,:)
|
||||||
|
|
||||||
|
! Output variables
|
||||||
|
|
||||||
|
double precision,intent(out) :: EcRPA
|
||||||
|
double precision,intent(out) :: Omega(nSt)
|
||||||
|
double precision,intent(out) :: XpY(nSt,nSt)
|
||||||
|
double precision,intent(out) :: XmY(nSt,nSt)
|
||||||
|
|
||||||
|
! Memory allocation
|
||||||
|
|
||||||
|
allocate(A(nSt,nSt),B(nSt,nSt),ApB(nSt,nSt),AmB(nSt,nSt),AmBSq(nSt,nSt),AmBIv(nSt,nSt),Z(nSt,nSt))
|
||||||
|
|
||||||
|
! Build A and B matrices
|
||||||
|
|
||||||
|
call unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,e,ERI_aa,ERI_ab,ERI_bb,A)
|
||||||
|
|
||||||
|
! if(BSE) call Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,A)
|
||||||
|
|
||||||
|
! Tamm-Dancoff approximation
|
||||||
|
|
||||||
|
B = 0d0
|
||||||
|
if(.not. TDA) then
|
||||||
|
|
||||||
|
call unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aa,ERI_ab,ERI_bb,B)
|
||||||
|
|
||||||
|
! if(BSE) call Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,B)
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
! Build A + B and A - B matrices
|
||||||
|
|
||||||
|
ApB = A + B
|
||||||
|
AmB = A - B
|
||||||
|
|
||||||
|
! Diagonalize linear response matrix
|
||||||
|
|
||||||
|
call diagonalize_matrix(nSt,AmB,Omega)
|
||||||
|
|
||||||
|
if(minval(Omega) < 0d0) &
|
||||||
|
call print_warning('You may have instabilities in linear response: A-B is not positive definite!!')
|
||||||
|
|
||||||
|
do ia=1,nSt
|
||||||
|
if(Omega(ia) < 0d0) Omega(ia) = 0d0
|
||||||
|
end do
|
||||||
|
|
||||||
|
call ADAt(nSt,AmB,1d0*sqrt(Omega),AmBSq)
|
||||||
|
call ADAt(nSt,AmB,1d0/sqrt(Omega),AmBIv)
|
||||||
|
|
||||||
|
Z = matmul(AmBSq,matmul(ApB,AmBSq))
|
||||||
|
|
||||||
|
call diagonalize_matrix(nSt,Z,Omega)
|
||||||
|
|
||||||
|
if(minval(Omega) < 0d0) &
|
||||||
|
call print_warning('You may have instabilities in linear response: negative excitations!!')
|
||||||
|
|
||||||
|
do ia=1,nSt
|
||||||
|
if(Omega(ia) < 0d0) Omega(ia) = 0d0
|
||||||
|
end do
|
||||||
|
|
||||||
|
Omega = sqrt(Omega)
|
||||||
|
|
||||||
|
XpY = matmul(transpose(Z),AmBSq)
|
||||||
|
call DA(nSt,1d0/sqrt(Omega),XpY)
|
||||||
|
|
||||||
|
XmY = matmul(transpose(Z),AmBIv)
|
||||||
|
call DA(nSt,1d0*sqrt(Omega),XmY)
|
||||||
|
|
||||||
|
! Compute the RPA correlation energy
|
||||||
|
|
||||||
|
EcRPA = 0.5d0*(sum(Omega) - trace_matrix(nSt,A))
|
||||||
|
|
||||||
|
end subroutine unrestricted_linear_response
|
||||||
126
src/QuAcK/unrestricted_linear_response_A_matrix.f90
Normal file
126
src/QuAcK/unrestricted_linear_response_A_matrix.f90
Normal file
@ -0,0 +1,126 @@
|
|||||||
|
subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, &
|
||||||
|
e,ERI_aa,ERI_ab,ERI_bb,A_lr)
|
||||||
|
|
||||||
|
! Compute linear response
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
include 'parameters.h'
|
||||||
|
|
||||||
|
! Input variables
|
||||||
|
|
||||||
|
logical,intent(in) :: dRPA
|
||||||
|
integer,intent(in) :: ispin
|
||||||
|
integer,intent(in) :: nBas
|
||||||
|
integer,intent(in) :: nC(nspin)
|
||||||
|
integer,intent(in) :: nO(nspin)
|
||||||
|
integer,intent(in) :: nV(nspin)
|
||||||
|
integer,intent(in) :: nR(nspin)
|
||||||
|
integer,intent(in) :: nSa
|
||||||
|
integer,intent(in) :: nSb
|
||||||
|
integer,intent(in) :: nSt
|
||||||
|
double precision,intent(in) :: lambda
|
||||||
|
double precision,intent(in) :: e(nBas,nspin)
|
||||||
|
double precision,intent(in) :: ERI_aa(nBas,nBas,nBas,nBas)
|
||||||
|
double precision,intent(in) :: ERI_ab(nBas,nBas,nBas,nBas)
|
||||||
|
double precision,intent(in) :: ERI_bb(nBas,nBas,nBas,nBas)
|
||||||
|
|
||||||
|
! Local variables
|
||||||
|
|
||||||
|
double precision :: delta_dRPA
|
||||||
|
double precision,external :: Kronecker_delta
|
||||||
|
|
||||||
|
integer :: i,j,a,b,ia,jb
|
||||||
|
|
||||||
|
! Output variables
|
||||||
|
|
||||||
|
double precision,intent(out) :: A_lr(nSt,nSt)
|
||||||
|
|
||||||
|
! Direct RPA
|
||||||
|
|
||||||
|
delta_dRPA = 0d0
|
||||||
|
if(dRPA) delta_dRPA = 1d0
|
||||||
|
|
||||||
|
!-----------------------------------------------
|
||||||
|
! Build A matrix for spin-conserving transitions
|
||||||
|
!-----------------------------------------------
|
||||||
|
|
||||||
|
if(ispin == 1) then
|
||||||
|
|
||||||
|
! alpha-alpha 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
|
||||||
|
|
||||||
|
A_lr(ia,jb) = (e(a,1) - e(i,1))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
|
||||||
|
+ lambda*ERI_aa(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI_aa(i,b,j,a)
|
||||||
|
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
|
||||||
|
! 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
|
||||||
|
|
||||||
|
A_lr(ia,nSa+jb) = lambda*ERI_ab(i,b,a,j)
|
||||||
|
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
|
||||||
|
! beta-alpha 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(1)+1,nO(1)
|
||||||
|
do b=nO(1)+1,nBas-nR(1)
|
||||||
|
jb = jb + 1
|
||||||
|
|
||||||
|
A_lr(nSa+ia,jb) = lambda*ERI_ab(b,i,j,a)
|
||||||
|
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
|
||||||
|
! beta-beta 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
|
||||||
|
|
||||||
|
A_lr(nSa+ia,nSa+jb) = (e(a,2) - e(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
|
||||||
|
+ lambda*ERI_bb(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI_bb(i,b,j,a)
|
||||||
|
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
|
||||||
|
end subroutine unrestricted_linear_response_A_matrix
|
||||||
123
src/QuAcK/unrestricted_linear_response_B_matrix.f90
Normal file
123
src/QuAcK/unrestricted_linear_response_B_matrix.f90
Normal file
@ -0,0 +1,123 @@
|
|||||||
|
subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, &
|
||||||
|
ERI_aa,ERI_ab,ERI_bb,B_lr)
|
||||||
|
|
||||||
|
! Compute linear response
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
include 'parameters.h'
|
||||||
|
|
||||||
|
! Input variables
|
||||||
|
|
||||||
|
logical,intent(in) :: dRPA
|
||||||
|
integer,intent(in) :: ispin
|
||||||
|
integer,intent(in) :: nBas
|
||||||
|
integer,intent(in) :: nC(nspin)
|
||||||
|
integer,intent(in) :: nO(nspin)
|
||||||
|
integer,intent(in) :: nV(nspin)
|
||||||
|
integer,intent(in) :: nR(nspin)
|
||||||
|
integer,intent(in) :: nSa
|
||||||
|
integer,intent(in) :: nSb
|
||||||
|
integer,intent(in) :: nSt
|
||||||
|
double precision,intent(in) :: lambda
|
||||||
|
double precision,intent(in) :: ERI_aa(nBas,nBas,nBas,nBas)
|
||||||
|
double precision,intent(in) :: ERI_ab(nBas,nBas,nBas,nBas)
|
||||||
|
double precision,intent(in) :: ERI_bb(nBas,nBas,nBas,nBas)
|
||||||
|
|
||||||
|
! Local variables
|
||||||
|
|
||||||
|
double precision :: delta_dRPA
|
||||||
|
double precision,external :: Kronecker_delta
|
||||||
|
|
||||||
|
integer :: i,j,a,b,ia,jb
|
||||||
|
|
||||||
|
! Output variables
|
||||||
|
|
||||||
|
double precision,intent(out) :: B_lr(nSt,nSt)
|
||||||
|
|
||||||
|
! Direct RPA
|
||||||
|
|
||||||
|
delta_dRPA = 0d0
|
||||||
|
if(dRPA) delta_dRPA = 1d0
|
||||||
|
|
||||||
|
!-----------------------------------------------
|
||||||
|
! Build A matrix for spin-conserving transitions
|
||||||
|
!-----------------------------------------------
|
||||||
|
|
||||||
|
if(ispin == 1) then
|
||||||
|
|
||||||
|
! alpha-alpha 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
|
||||||
|
|
||||||
|
B_lr(ia,jb) = lambda*ERI_aa(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI_aa(i,j,b,a)
|
||||||
|
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
|
||||||
|
! 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
|
||||||
|
|
||||||
|
B_lr(ia,nSa+jb) = lambda*ERI_ab(i,j,a,b)
|
||||||
|
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
|
||||||
|
! beta-alpha 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(1)+1,nO(1)
|
||||||
|
do b=nO(1)+1,nBas-nR(1)
|
||||||
|
jb = jb + 1
|
||||||
|
|
||||||
|
B_lr(nSa+ia,jb) = lambda*ERI_ab(j,i,b,a)
|
||||||
|
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
|
||||||
|
! beta-beta 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
|
||||||
|
|
||||||
|
B_lr(nSa+ia,nSa+jb) = lambda*ERI_bb(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI_bb(i,j,b,a)
|
||||||
|
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
|
||||||
|
end subroutine unrestricted_linear_response_B_matrix
|
||||||
@ -33,9 +33,9 @@ subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSa,nS
|
|||||||
|
|
||||||
SigC(:,:) = 0d0
|
SigC(:,:) = 0d0
|
||||||
|
|
||||||
!--------------
|
!--------------!
|
||||||
! Spin-up part
|
! Spin-up part !
|
||||||
!--------------
|
!--------------!
|
||||||
|
|
||||||
! Occupied part of the correlation self-energy
|
! Occupied part of the correlation self-energy
|
||||||
|
|
||||||
@ -59,9 +59,9 @@ subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSa,nS
|
|||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
|
|
||||||
!----------------
|
!----------------!
|
||||||
! Spin-down part
|
! Spin-down part !
|
||||||
!----------------
|
!----------------!
|
||||||
|
|
||||||
! Occupied part of the correlation self-energy
|
! Occupied part of the correlation self-energy
|
||||||
|
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user