4
1
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:35:36 +01:00

spin conserved and spin flip

This commit is contained in:
Pierre-Francois Loos 2020-09-22 23:08:47 +02:00
parent 8f5b1779de
commit b8bf488a9a
7 changed files with 143 additions and 112 deletions

View File

@ -4,8 +4,8 @@
# CC: maxSCF thresh DIIS n_diis
64 0.0000001 T 5
# spin: singlet triplet TDA
T T F
# spin: singlet triplet spin_conserved spinf_flip TDA
T T T F 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

View File

@ -1,15 +1,15 @@
subroutine AOtoMO_integral_transform(bra,ket,nBas,c,ERI_AO_basis,ERI_MO_basis)
subroutine AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,c,ERI_AO_basis,ERI_MO_basis)
! AO to MO transformation of two-electron integrals via the semi-direct O(N^5) algorithm
! bra and ket are the spin of (bra|ket)
! bra and ket are the spin of (bra1 bra2|ket1 ket2)
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: bra
integer,intent(in) :: ket
integer,intent(in) :: bra1,bra2
integer,intent(in) :: ket1,ket2
integer,intent(in) :: nBas
double precision,intent(in) :: ERI_AO_basis(nBas,nBas,nBas,nBas),c(nBas,nBas,nspin)
@ -35,7 +35,7 @@ subroutine AOtoMO_integral_transform(bra,ket,nBas,c,ERI_AO_basis,ERI_MO_basis)
do la=1,nBas
do nu=1,nBas
do mu=1,nBas
scr(mu,nu,la,l) = scr(mu,nu,la,l) + ERI_AO_basis(mu,nu,la,si)*c(si,l,ket)
scr(mu,nu,la,l) = scr(mu,nu,la,l) + ERI_AO_basis(mu,nu,la,si)*c(si,l,ket2)
enddo
enddo
enddo
@ -49,7 +49,7 @@ subroutine AOtoMO_integral_transform(bra,ket,nBas,c,ERI_AO_basis,ERI_MO_basis)
do nu=1,nBas
do i=1,nBas
do mu=1,nBas
ERI_MO_basis(i,nu,la,l) = ERI_MO_basis(i,nu,la,l) + c(mu,i,bra)*scr(mu,nu,la,l)
ERI_MO_basis(i,nu,la,l) = ERI_MO_basis(i,nu,la,l) + c(mu,i,bra1)*scr(mu,nu,la,l)
enddo
enddo
enddo
@ -63,7 +63,7 @@ subroutine AOtoMO_integral_transform(bra,ket,nBas,c,ERI_AO_basis,ERI_MO_basis)
do la=1,nBas
do nu=1,nBas
do i=1,nBas
scr(i,nu,k,l) = scr(i,nu,k,l) + ERI_MO_basis(i,nu,la,l)*c(la,k,bra)
scr(i,nu,k,l) = scr(i,nu,k,l) + ERI_MO_basis(i,nu,la,l)*c(la,k,bra2)
enddo
enddo
enddo
@ -77,7 +77,7 @@ subroutine AOtoMO_integral_transform(bra,ket,nBas,c,ERI_AO_basis,ERI_MO_basis)
do j=1,nBas
do i=1,nBas
do nu=1,nBas
ERI_MO_basis(i,j,k,l) = ERI_MO_basis(i,j,k,l) + c(nu,j,ket)*scr(i,nu,k,l)
ERI_MO_basis(i,j,k,l) = ERI_MO_basis(i,j,k,l) + c(nu,j,ket1)*scr(i,nu,k,l)
enddo
! print*,i,k,j,l,ERI_MO_basis(i,j,k,l)
enddo

View File

@ -55,11 +55,12 @@ program QuAcK
double precision,allocatable :: S(:,:),T(:,:),V(:,:),Hc(:,:),H(:,:),X(:,:)
double precision,allocatable :: ERI_AO(:,:,:,:)
double precision,allocatable :: ERI_MO(:,:,:,:)
integer :: bra
integer :: ket
integer :: bra1,bra2
integer :: ket1,ket2
double precision,allocatable :: ERI_MO_aaaa(:,:,:,:)
double precision,allocatable :: ERI_MO_aabb(:,:,:,:)
double precision,allocatable :: ERI_MO_bbbb(:,:,:,:)
double precision,allocatable :: ERI_MO_abab(:,:,:,:)
double precision,allocatable :: ERI_ERF_AO(:,:,:,:)
double precision,allocatable :: ERI_ERF_MO(:,:,:,:)
double precision,allocatable :: F12(:,:,:,:),Yuk(:,:,:,:),FC(:,:,:,:,:,:)
@ -101,8 +102,10 @@ program QuAcK
double precision :: thresh_CC
logical :: DIIS_CC
logical :: singlet_manifold
logical :: triplet_manifold
logical :: singlet
logical :: triplet
logical :: spin_conserved
logical :: spin_flip
logical :: TDA
integer :: maxSCF_GF,n_diis_GF,renormGF
@ -156,7 +159,7 @@ program QuAcK
call read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type, &
maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, &
singlet_manifold,triplet_manifold,TDA, &
singlet,triplet,spin_conserved,spin_flip,TDA, &
maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF, &
maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW, &
COHSEX,SOSEX,TDA_W,G0W,GW0, &
@ -334,21 +337,42 @@ program QuAcK
! 4-index transform for (aa|aa) block
bra = 1
ket = 1
call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_aaaa)
bra1 = 1
bra2 = 1
ket1 = 1
ket2 = 1
call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_aaaa)
! 4-index transform for (aa|bb) block
bra = 1
ket = 2
call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_aabb)
bra1 = 1
bra2 = 1
ket1 = 2
ket2 = 2
call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_aabb)
! 4-index transform for (bb|bb) block
bra = 2
ket = 2
call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_bbbb)
bra1 = 2
bra2 = 2
ket1 = 2
ket2 = 2
call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_bbbb)
if(spin_flip) then
allocate(ERI_MO_abab(nBas,nBas,nBas,nBas))
! 4-index transform for (ab|ab) block
bra1 = 1
bra2 = 2
ket1 = 1
ket2 = 2
call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_abab)
end if
else
@ -358,9 +382,11 @@ program QuAcK
! 4-index transform
bra = 1
ket = 1
call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO)
bra1 = 1
bra2 = 1
ket1 = 1
ket2 = 1
call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO)
end if
@ -560,7 +586,7 @@ program QuAcK
if(doCIS) then
call cpu_time(start_CIS)
call CIS(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ERI_MO,eHF)
call CIS(singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI_MO,eHF)
call cpu_time(end_CIS)
t_CIS = end_CIS - start_CIS
@ -576,7 +602,7 @@ program QuAcK
if(doCID) then
call cpu_time(start_CID)
! call CID(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ERI_MO,eHF)
! call CID(singlet,triplet,nBas,nC,nO,nV,nR,ERI_MO,eHF)
call cpu_time(end_CID)
t_CID = end_CID - start_CID
@ -592,7 +618,7 @@ program QuAcK
if(doCISD) then
call cpu_time(start_CISD)
call CISD(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ERI_MO,eHF)
call CISD(singlet,triplet,nBas,nC,nO,nV,nR,ERI_MO,eHF)
call cpu_time(end_CISD)
t_CISD = end_CISD - start_CISD
@ -608,7 +634,7 @@ program QuAcK
if(doRPA) then
call cpu_time(start_RPA)
call RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,0d0, &
call RPA(doACFDT,exchange_kernel,singlet,triplet,0d0, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call cpu_time(end_RPA)
@ -625,7 +651,7 @@ program QuAcK
if(doRPAx) then
call cpu_time(start_RPAx)
call RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,0d0, &
call RPAx(doACFDT,exchange_kernel,singlet,triplet,0d0, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call cpu_time(end_RPAx)
@ -642,7 +668,7 @@ program QuAcK
if(doppRPA) then
call cpu_time(start_ppRPA)
call ppRPA(singlet_manifold,triplet_manifold, &
call ppRPA(singlet,triplet, &
nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF)
call cpu_time(end_ppRPA)
@ -659,7 +685,7 @@ program QuAcK
! if(doADC) then
! call cpu_time(start_ADC)
! call ADC(singlet_manifold,triplet_manifold,maxSCF_GF,thresh_GF,n_diis_GF, &
! call ADC(singlet,triplet,maxSCF_GF,thresh_GF,n_diis_GF, &
! nBas,nC,nO,nV,nR,eHF,ERI_MO)
! call cpu_time(end_ADC)
@ -676,7 +702,7 @@ program QuAcK
if(doG0F2) then
call cpu_time(start_GF2)
call G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linGF, &
call G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linGF, &
eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call cpu_time(end_GF2)
@ -694,7 +720,7 @@ program QuAcK
call cpu_time(start_GF2)
call evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF_GF,thresh_GF,n_diis_GF, &
singlet_manifold,triplet_manifold,linGF, &
singlet,triplet,linGF, &
eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call cpu_time(end_GF2)
@ -748,12 +774,12 @@ program QuAcK
if(unrestricted) then
call UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn, &
singlet_manifold,triplet_manifold,linGW,eta_GW,nBas,nC,nO,nV,nR,nS, &
spin_conserved,spin_flip,linGW,eta_GW,nBas,nC,nO,nV,nR,nS, &
ENuc,EUHF,Hc,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,PHF,cHF,eHF,eG0W0)
else
call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linGW,eta_GW, &
dBSE,dTDA,evDyn,singlet,triplet,linGW,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_MO,PHF,cHF,eHF,eG0W0)
end if
@ -774,7 +800,7 @@ program QuAcK
call cpu_time(start_evGW)
call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, &
BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta_GW, &
BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0)
call cpu_time(end_evGW)
@ -792,7 +818,7 @@ program QuAcK
call cpu_time(start_qsGW)
call qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, &
BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta_GW, &
BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,PHF,cHF,eHF)
call cpu_time(end_qsGW)
@ -812,7 +838,7 @@ program QuAcK
call cpu_time(start_G0T0)
call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA, &
dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linGW,eta_GW, &
dBSE,dTDA,evDyn,singlet,triplet,linGW,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0T0)
call cpu_time(end_G0T0)
@ -830,7 +856,7 @@ program QuAcK
call cpu_time(start_evGT)
call evGT(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS, &
BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta_GW, &
BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0T0)
call cpu_time(end_evGT)
@ -947,7 +973,7 @@ program QuAcK
call cpu_time(start_G0W0)
call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linGW,eta_GW, &
dBSE,dTDA,evDyn,singlet,triplet,linGW,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_ERF_MO,PHF,cHF,eHF,eG0W0)
call cpu_time(end_G0W0)
@ -961,7 +987,7 @@ program QuAcK
call cpu_time(start_G0T0)
call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn, &
singlet_manifold,triplet_manifold,linGW,eta_GW, &
singlet,triplet,linGW,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_ERF_MO,eHF,eG0T0)
call cpu_time(end_G0T0)

View File

@ -45,21 +45,13 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
logical :: print_W = .true.
integer :: ispin
integer :: iblock
integer :: bra
integer :: ket
integer :: nSa
integer :: nSb
integer :: nSt
double precision :: EcRPA
double precision :: EcBSE
double precision :: EcRPA(nspin)
double precision :: EcBSE(nspin)
double precision :: EcAC(nspin)
double precision,allocatable :: SigC(:,:)
double precision,allocatable :: Z(:,:)
double precision,allocatable :: Omega(:)
double precision,allocatable :: XpY(:,:)
double precision,allocatable :: XmY(:,:)
double precision,allocatable :: rho(:,:,:,:)
integer :: nS_aa,nS_bb,nS_sc
double precision,allocatable :: Omega_sc(:),XpY_sc(:,:),XmY_sc(:,:),rho_sc(:,:,:,:)
double precision,allocatable :: eGWlin(:,:)
@ -103,12 +95,12 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
! Memory allocation
nSa = nS(1)
nSb = nS(2)
nSt = nSa + nSb
nS_aa = nS(1)
nS_bb = nS(2)
nS_sc = nS_aa + nS_bb
allocate(SigC(nBas,nspin),Z(nBas,nspin),Omega(nSt),XpY(nSt,nSt),XmY(nSt,nSt), &
rho(nBas,nBas,nSt,nspin),eGWlin(nBas,nspin))
allocate(SigC(nBas,nspin),Z(nBas,nspin),Omega_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc), &
rho_sc(nBas,nBas,nS_sc,nspin),eGWlin(nBas,nspin))
!-------------------!
! Compute screening !
@ -118,28 +110,28 @@ 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,nSa,nSb,nSt,1d0, &
eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,rho,EcRPA,Omega,XpY,XmY)
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,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc)
if(print_W) call print_excitation('RPA@UHF',5,nSt,Omega)
if(print_W) call print_excitation('RPA@UHF',5,nS_sc,Omega_sc)
!----------------------!
! Excitation densities !
!----------------------!
call unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,rho)
call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_sc,rho_sc)
!---------------------!
! 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,nS_aa,nS_bb,nS_sc,eHF,Omega_sc,rho_sc,SigC)
!--------------------------------!
! Compute renormalization factor !
!--------------------------------!
call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,eHF,Omega,rho,Z)
call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,eHF,Omega_sc,rho_sc,Z)
!-----------------------------------!
! Solve the quasi-particle equation !
@ -171,23 +163,26 @@ 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,nSa,nSb,nSt,1d0, &
eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,rho,EcRPA,Omega,XpY,XmY)
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,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy =',EcRPA
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 total energy =',ENuc + EUHF + EcRPA
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy =',EcRPA(ispin)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 total energy =',ENuc + EUHF + EcRPA(ispin)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! Free memory
deallocate(Omega_sc,XpY_sc,XmY_sc,rho_sc)
! Perform BSE calculation
if(BSE) then
call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta, &
nBas,nC,nO,nV,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ERI_bbbb, &
eHF,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE)
nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGW,EcRPA,EcBSE)
! if(exchange_kernel) then
!

View File

@ -149,7 +149,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
! AO to MO transformation of two-electron integrals
call AOtoMO_integral_transform(nBas,c,ERI_AO_basis,ERI_MO_basis)
call AOtoMO_integral_transform(1,1,1,1,nBas,c,ERI_AO_basis,ERI_MO_basis)
! Compute linear response

View File

@ -1,6 +1,6 @@
subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type, &
maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, &
singlet_manifold,triplet_manifold,TDA, &
singlet,triplet,spin_conserved,spin_flip,TDA, &
maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF, &
maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW, &
COHSEX,SOSEX,TDA_W,G0W,GW0, &
@ -26,8 +26,10 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
logical,intent(out) :: DIIS_CC
integer,intent(out) :: n_diis_CC
logical,intent(out) :: singlet_manifold
logical,intent(out) :: triplet_manifold
logical,intent(out) :: singlet
logical,intent(out) :: triplet
logical,intent(out) :: spin_conserved
logical,intent(out) :: spin_flip
logical,intent(out) :: TDA
integer,intent(out) :: maxSCF_GF
@ -113,16 +115,20 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
! Read excited state options
singlet_manifold = .false.
triplet_manifold = .false.
TDA = .false.
singlet = .false.
triplet = .false.
spin_conserved = .false.
spin_flip = .false.
TDA = .false.
read(1,*)
read(1,*) answer1,answer2,answer3
read(1,*) answer1,answer2,answer3,answer4,answer5
if(answer1 == 'T') singlet_manifold = .true.
if(answer2 == 'T') triplet_manifold = .true.
if(answer3 == 'T') TDA = .true.
if(answer1 == 'T') singlet = .true.
if(answer2 == 'T') triplet = .true.
if(answer3 == 'T') spin_conserved = .true.
if(answer4 == 'T') spin_flip = .true.
if(answer5 == 'T') TDA = .true.
! Read Green function options

View File

@ -1,6 +1,5 @@
subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta, &
nBas,nC,nO,nV,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ERI_bbbb, &
eW,eGW,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,EcRPA,EcBSE)
nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eW,eGW,EcRPA,EcBSE)
! Compute the Bethe-Salpeter excitation energies
@ -23,36 +22,32 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,
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
integer,intent(in) :: nS(nspin)
double precision,intent(in) :: eW(nBas,nspin)
double precision,intent(in) :: eGW(nBas,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 :: OmRPA(nSt)
double precision :: XpY_RPA(nSt,nSt)
double precision :: XmY_RPA(nSt,nSt)
double precision :: rho_RPA(nBas,nBas,nSt,nspin)
! Local variables
integer :: ispin
integer :: isp_W
double precision,allocatable :: OmBSE(:)
double precision,allocatable :: XpY_BSE(:,:)
double precision,allocatable :: XmY_BSE(:,:)
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(:,:)
double precision,allocatable :: rho_RPA_sc(:,:,:,:)
double precision,allocatable :: OmBSE_sc(:)
double precision,allocatable :: XpY_BSE_sc(:,:)
double precision,allocatable :: XmY_BSE_sc(:,:)
! Output variables
double precision,intent(out) :: EcRPA
double precision,intent(out) :: EcBSE
! Memory allocation
allocate(OmBSE(nSt),XpY_BSE(nSt,nSt),XmY_BSE(nSt,nSt))
double precision,intent(out) :: EcRPA(nspin)
double precision,intent(out) :: EcBSE(nspin)
!----------------------------!
! Spin-conserved excitations !
@ -62,23 +57,32 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,
ispin = 1
isp_W = 1
EcBSE = 0d0
EcBSE(ispin) = 0d0
! Compute spin-conserved RPA screening
! Memory allocation
call unrestricted_linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0, &
eW,ERI_aaaa,ERI_aabb,ERI_bbbb,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
nS_aa = nS(1)
nS_bb = nS(2)
nS_sc = nS_aa + nS_bb
call unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
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 BSE excitation energies
! Compute spin-conserved RPA screening
OmBSE(:) = OmRPA(:)
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,rho_RPA_sc,EcRPA(ispin),OmRPA_sc,XpY_RPA_sc,XmY_RPA_sc)
call unrestricted_linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0, &
eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,rho_RPA,EcBSE,OmBSE,XpY_BSE,XmY_BSE)
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)
call print_excitation('BSE@UG0W0',5,nSt,OmBSE)
! 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,rho_RPA_sc,EcBSE(ispin),OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc)
call print_excitation('BSE@UG0W0',5,nS_sc,OmBSE_sc)
!-------------------------------------------------
! Compute the dynamical screening at the BSE level