PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-11-04 13:13:51 +01:00
Code Releases Activity

UCIS

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-09-24 14:39:37 +02:00
parent 8910ead99f
commit ff58cd17c6
16 changed files with 251 additions and 135 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

55
input/basis
View File

@ -1,39 +1,30 @@
1 10
1 6
S 8
1 24350.0000000 0.0005020
2 3650.0000000 0.0038810
3 829.6000000 0.0199970
4 234.0000000 0.0784180
5 75.6100000 0.2296760
6 26.7300000 0.4327220
7 9.9270000 0.3506420
8 1.1020000 -0.0076450
1 1469.0000000 0.0007660
2 220.5000000 0.0058920
3 50.2600000 0.0296710
4 14.2400000 0.1091800
5 4.5810000 0.2827890
6 1.5800000 0.4531230
7 0.5640000 0.2747740
8 0.0734500 0.0097510
S 8
1 24350.0000000 -0.0001180
2 3650.0000000 -0.0009150
3 829.6000000 -0.0047370
4 234.0000000 -0.0192330
5 75.6100000 -0.0603690
6 26.7300000 -0.1425080
7 9.9270000 -0.1777100
8 1.1020000 0.6058360
1 1469.0000000 -0.0001200
2 220.5000000 -0.0009230
3 50.2600000 -0.0046890
4 14.2400000 -0.0176820
5 4.5810000 -0.0489020
6 1.5800000 -0.0960090
7 0.5640000 -0.1363800
8 0.0734500 0.5751020
S 1
1 2.8360000 1.0000000
S 1
1 0.3782000 1.0000000
1 0.0280500 1.0000000
P 3
1 54.7000000 0.0171510
2 12.4300000 0.1076560
3 3.6790000 0.3216810
1 1.5340000 0.0227840
2 0.2749000 0.1391070
3 0.0736200 0.5003750
P 1
1 1.1430000 1.0000000
P 1
1 0.3300000 1.0000000
1 0.0240300 1.0000000
D 1
1 4.0140000 1.0000000
D 1
1 1.0960000 1.0000000
F 1
1 2.5440000 1.0000000
1 0.1239000 1.0000000

4
input/methods
View File

@ -1,5 +1,5 @@
# RHF UHF MOM
T F F
F T F
# MP2 MP3 MP2-F12
F F F
# CCD CCSD CCSD(T)
@ -13,7 +13,7 @@
# G0F2 evGF2 G0F3 evGF3
F F F F
# G0W0 evGW qsGW
F F T
T F F
# G0T0 evGT qsGT
F F F
# MCMP2

4
input/molecule
View File

@ -1,4 +1,4 @@
# nAt nEla nElb nCore nRyd
1 5 5 0 0
1 2 1 0 0
# Znuc x y z
Ne 0.0 0.0 0.0
Li 0.0 0.0 0.0

2
input/molecule.xyz
View File

@ -1,3 +1,3 @@
1
Ne 0.0000000000 0.0000000000 0.0000000000
Li 0.0000000000 0.0000000000 0.0000000000

8
input/options
View File

@ -5,14 +5,14 @@
# CC: maxSCF thresh DIIS n_diis
64 0.0000001 T 5
# spin: singlet triplet spin_conserved spin_flip TDA
T T T T F
T T T T T
# 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 T 0.0 F F F F F
256 0.00001 T 5 F 0.0 F F F F F
# ACFDT: AC Kx XBS
F F T
T F T
# BSE: BSE dBSE dTDA evDyn
T F T T
T T T T
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
1000000 100000 10 0.3 10000 1234 T

127
src/QuAcK/UCIS.f90 Normal file
View File

@ -0,0 +1,127 @@
subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,eHF)
! Perform configuration interaction single calculation`
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: spin_conserved
logical,intent(in) :: spin_flip
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) :: nS(nspin)
double precision,intent(in) :: eHF(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,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas)
! Local variables
logical :: dump_matrix = .false.
logical :: dump_trans = .false.
integer :: ispin
double precision :: lambda
integer :: nS_aa,nS_bb,nS_sc
double precision,allocatable :: A_sc(:,:)
double precision,allocatable :: Omega_sc(:)
integer :: nS_ab,nS_ba,nS_sf
double precision,allocatable :: A_sf(:,:)
double precision,allocatable :: Omega_sf(:)
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Configuration Interaction Singles |'
write(*,*)'************************************************'
write(*,*)
! Adiabatic connection scaling
lambda = 1d0
!----------------------------!
! Spin-conserved transitions !
!----------------------------!
if(spin_conserved) then
ispin = 1
! Memory allocation
nS_aa = nS(1)
nS_bb = nS(2)
nS_sc = nS_aa + nS_bb
allocate(A_sc(nS_sc,nS_sc),Omega_sc(nS_sc))
call unrestricted_linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eHF, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,A_sc)
if(dump_matrix) then
print*,'CIS matrix (spin-conserved transitions)'
call matout(nS_sc,nS_sc,A_sc)
write(*,*)
endif
call diagonalize_matrix(nS_sc,A_sc,Omega_sc)
call print_excitation('UCIS ',5,nS_sc,Omega_sc)
if(dump_trans) then
print*,'Spin-conserved CIS transition vectors'
call matout(nS_sc,nS_sc,A_sc)
write(*,*)
endif
deallocate(A_sc,Omega_sc)
endif
!-----------------------!
! Spin-flip transitions !
!-----------------------!
if(spin_flip) then
ispin = 2
! Memory allocation
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(A_sf(nS_sf,nS_sf),Omega_sf(nS_sf))
call unrestricted_linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,nS_sf,lambda,eHF, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,A_sf)
if(dump_matrix) then
print*,'CIS matrix (spin-conserved transitions)'
call matout(nS_sf,nS_sf,A_sf)
write(*,*)
endif
call diagonalize_matrix(nS_sf,A_sf,Omega_sf)
call print_excitation('UCIS ',6,nS_sf,Omega_sf)
if(dump_trans) then
print*,'Spin-flip CIS transition vectors'
call matout(nS_sf,nS_sf,A_sf)
write(*,*)
endif
deallocate(A_sf,Omega_sf)
endif
end subroutine UCIS

46
src/QuAcK/UG0W0.f90
View File

@ -47,13 +47,16 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
logical :: print_W = .true.
integer :: is
integer :: ispin
double precision :: EcRPA(nspin)
double precision :: EcRPA
double precision :: EcBSE(nspin)
double precision :: EcAC(nspin)
double precision,allocatable :: SigC(:,:)
double precision,allocatable :: Z(:,:)
integer :: nS_aa,nS_bb,nS_sc
double precision,allocatable :: Omega_sc(:),XpY_sc(:,:),XmY_sc(:,:),rho_sc(:,:,:,:)
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:,:)
double precision,allocatable :: eGWlin(:,:)
@ -101,8 +104,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
nS_bb = nS(2)
nS_sc = nS_aa + nS_bb
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))
allocate(SigC(nBas,nspin),Z(nBas,nspin),OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc), &
rho_RPA(nBas,nBas,nS_sc,nspin),eGWlin(nBas,nspin))
!-------------------!
! Compute screening !
@ -113,27 +116,27 @@ 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,nS_sc,1d0, &
eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sc,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc)
eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
if(print_W) call print_excitation('RPA@UHF',5,nS_sc,Omega_sc)
if(print_W) call print_excitation('RPA@UHF',5,nS_sc,OmRPA)
!----------------------!
! Excitation densities !
!----------------------!
call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_sc,rho_sc)
call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
!---------------------!
! Compute self-energy !
!---------------------!
call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,Omega_sc,rho_sc,SigC)
call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC)
!--------------------------------!
! Compute renormalization factor !
!--------------------------------!
call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,Omega_sc,rho_sc,Z)
call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,Z)
!-----------------------------------!
! Solve the quasi-particle equation !
@ -153,31 +156,24 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
! Find graphical solution of the QP equation
do is=1,nspin
call unrestricted_QP_graph(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eta,eHF(:,is),Omega_sc, &
rho_sc,eGWlin(:,is),eGW(:,is))
call unrestricted_QP_graph(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eta,eHF(:,is),OmRPA, &
rho_RPA(:,:,:,is),eGWlin(:,is),eGW(:,is))
end do
end if
! Compute RPA correlation energy
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,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
! Dump results
call print_UG0W0(nBas,nO,eHF,ENuc,EUHF,SigC,Z,eGW,EcRPA)
! 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,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(*,*)'-------------------------------------------------------------------------------'
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)
deallocate(OmRPA,XpY_RPA,XmY_RPA,rho_RPA)
! Perform BSE calculation
@ -185,7 +181,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta, &
nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab, &
eHF,eGW,EcRPA,EcBSE)
eHF,eGW,EcBSE)
! if(exchange_kernel) then
!

4
src/QuAcK/UHF.f90
View File

@ -166,7 +166,8 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,EUH
n_diis = min(n_diis+1,max_diis)
do ispin=1,nspin
call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis,err_diis(:,:,ispin),F_diis(:,:,ispin),err(:,:,ispin),F(:,:,ispin))
if(nO(ispin) > 1) call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis,err_diis(:,:,ispin),F_diis(:,:,ispin), &
err(:,:,ispin),F(:,:,ispin))
end do
! Reset DIIS if required
@ -232,7 +233,6 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,EUH
! Compute final UHF energy
call matout(nBas,2,e)
call print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF)
end subroutine UHF

2
src/QuAcK/URPAx.f90
View File

@ -78,6 +78,7 @@ subroutine URPAx(doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO
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))
deallocate(Omega_sc,XpY_sc,XmY_sc)
endif
@ -100,6 +101,7 @@ subroutine URPAx(doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO
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))
deallocate(Omega_sf,XpY_sf,XmY_sf)
endif

6
src/QuAcK/USigmaC.f90
View File

@ -18,7 +18,7 @@ double precision function USigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho)
integer,intent(in) :: nS
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega(nS)
double precision,intent(in) :: rho(nBas,nBas,nS,nspin)
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables
@ -34,14 +34,14 @@ double precision function USigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho)
do i=nC+1,nO
do jb=1,nS
eps = w - e(i) + Omega(jb)
USigmaC = uSigmaC + rho(p,i,jb,1)**2*eps/(eps**2 + eta**2)
USigmaC = uSigmaC + rho(p,i,jb)**2*eps/(eps**2 + eta**2)
end do
end do
do a=nO+1,nBas-nR
do jb=1,nS
eps = w - e(a) - Omega(jb)
USigmaC = USigmaC + rho(p,a,jb,1)**2*eps/(eps**2 + eta**2)
USigmaC = USigmaC + rho(p,a,jb)**2*eps/(eps**2 + eta**2)
end do
end do

6
src/QuAcK/dUSigmaC.f90
View File

@ -18,7 +18,7 @@ double precision function dUSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho)
integer,intent(in) :: nS
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega(nS)
double precision,intent(in) :: rho(nBas,nBas,nS,nspin)
double precision,intent(in) :: rho(nBas,nBas,nS)
! Local variables
@ -34,7 +34,7 @@ double precision function dUSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho)
do i=nC+1,nO
do jb=1,nS
eps = w - e(i) + Omega(jb)
dUSigmaC = dUSigmaC + rho(p,i,jb,1)**2*(eps/(eps**2 + eta**2))**2
dUSigmaC = dUSigmaC + rho(p,i,jb)**2*(eps/(eps**2 + eta**2))**2
end do
end do
@ -43,7 +43,7 @@ double precision function dUSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho)
do a=nO+1,nBas-nR
do jb=1,nS
eps = w - e(a) - Omega(jb)
dUSigmaC = dUSigmaC + rho(p,a,jb,1)**2*(eps/(eps**2 + eta**2))**2
dUSigmaC = dUSigmaC + rho(p,a,jb)**2*(eps/(eps**2 + eta**2))**2
end do
end do

9
src/QuAcK/print_UG0W0.f90
View File

@ -33,9 +33,10 @@ subroutine print_UG0W0(nBas,nO,e,ENuc,EHF,SigC,Z,eGW,EcRPA)
write(*,*)' Unrestricted one-shot G0W0 calculation (eV)'
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,'(A1,A3,A1,A30,A1,A30,A1,A30,A1,A30,A1)') &
'|',' ','|','e_HF ','|','Sig_c ','|','Z ','|','e_QP ','|'
write(*,'(A1,A3,A1,2A15,A1,2A15,A1,2A15,A1,2A15,A1)') &
'|','#','|','e_HF up','e_HF dw','|','Sig_c up','Sig_c dw','|', &
'Z up','Z dw','|','e_QP up','e_QP dw','|'
'|','#','|','up ','dw ','|','up ','dw ','|','up ','dw ','|','up ','dw ','|'
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
@ -52,8 +53,8 @@ subroutine print_UG0W0(nBas,nO,e,ENuc,EHF,SigC,Z,eGW,EcRPA)
write(*,'(2X,A30,F15.6)') 'G0W0 HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,'(2X,A30,F15.6)') 'RPA@HF total energy =',ENuc + EHF + EcRPA
write(*,'(2X,A30,F15.6)') 'RPA@HF correlation energy =',EcRPA
write(*,'(2X,A30,F15.6)') 'RPA@G0W0 total energy =',ENuc + EHF + EcRPA
write(*,'(2X,A30,F15.6)') 'RPA@G0W0 correlation energy =',EcRPA
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,*)

47
src/QuAcK/unrestricted_Bethe_Salpeter.f90
View File

@ -1,6 +1,6 @@
subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta, &
nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab, &
eW,eGW,EcRPA,EcBSE)
eW,eGW,EcBSE)
! Compute the Bethe-Salpeter excitation energies
@ -38,10 +38,12 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,
integer :: isp_W
integer :: nS_aa,nS_bb,nS_sc
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 :: EcRPA
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:,:)
double precision,allocatable :: OmBSE_sc(:)
double precision,allocatable :: XpY_BSE_sc(:,:)
double precision,allocatable :: XmY_BSE_sc(:,:)
@ -53,48 +55,46 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,
! Output variables
double precision,intent(out) :: EcRPA(nspin)
double precision,intent(out) :: EcBSE(nspin)
!----------------------------!
! Spin-conserved excitations !
!----------------------------!
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))
allocate(OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin))
!--------------------------!
! Spin-conserved screening !
!--------------------------!
isp_W = 1
EcRPA = 0d0
! 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)
eW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
! 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,rho_RPA)
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)
!----------------------------!
! Spin-conserved excitations !
!----------------------------!
if(spin_conserved) then
ispin = 1
EcBSE(ispin) = 0d0
allocate(OmBSE_sc(nS_sc),XpY_BSE_sc(nS_sc,nS_sc),XmY_BSE_sc(nS_sc,nS_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,nS_sc,1d0, &
eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA_sc,rho_RPA_sc,EcBSE(ispin), &
eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcBSE(ispin), &
OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc)
call print_excitation('BSE@UG0W0',5,nS_sc,OmBSE_sc)
@ -128,7 +128,6 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,
if(spin_flip) then
ispin = 2
EcBSE(ispin) = 0d0
! Memory allocation
@ -140,7 +139,7 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,
allocate(OmBSE_sf(nS_sf),XpY_BSE_sf(nS_sf,nS_sf),XmY_BSE_sf(nS_sf,nS_sf))
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), &
eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcBSE(ispin), &
OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf)
call print_excitation('BSE@UG0W0',6,nS_sf,OmBSE_sf)

24
src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90
View File

@ -1,4 +1,4 @@
subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nSsc,lambda, &
subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,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
@ -17,15 +17,15 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n
integer,intent(in) :: nSa
integer,intent(in) :: nSb
integer,intent(in) :: nSt
integer,intent(in) :: nSsc
integer,intent(in) :: nS_sc
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(nSsc)
double precision,intent(in) :: rho(nBas,nBas,nSsc,nspin)
double precision,intent(in) :: Omega(nS_sc)
double precision,intent(in) :: rho(nBas,nBas,nS_sc,nspin)
! Local variables
@ -55,12 +55,12 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n
jb = jb + 1
chi = 0d0
do kc=1,nSsc
do kc=1,nS_sc
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
A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_aaaa(i,b,j,a) + 2d0*lambda*chi
enddo
enddo
@ -79,12 +79,12 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n
jb = jb + 1
chi = 0d0
do kc=1,nSsc
do kc=1,nS_sc
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
A_lr(nSa+ia,nSa+jb) = A_lr(nSa+ia,nSa+jb) - lambda*ERI_bbbb(i,b,j,a) + 2d0*lambda*chi
enddo
enddo
@ -111,12 +111,12 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n
jb = jb + 1
chi = 0d0
do kc=1,nSsc
do kc=1,nS_sc
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,j,kc,1)*rho(a,b,kc,2)*Omega(kc)/eps
enddo
A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_abab(i,b,j,a) + 4d0*lambda*chi
A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_abab(i,b,j,a) + 2d0*lambda*chi
end do
end do
@ -135,12 +135,12 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n
jb = jb + 1
chi = 0d0
do kc=1,nSsc
do kc=1,nS_sc
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
A_lr(nSa+ia,nSa+jb) = A_lr(nSa+ia,nSa+jb) - lambda*ERI_abab(b,i,a,j) + 2d0*lambda*chi
end do
end do

24
src/QuAcK/unrestricted_Bethe_Salpeter_B_matrix.f90
View File

@ -1,4 +1,4 @@
subroutine unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nSsc,lambda, &
subroutine unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega,rho,B_lr)
! Compute the extra term for Bethe-Salpeter equation for linear response
@ -17,15 +17,15 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n
integer,intent(in) :: nSa
integer,intent(in) :: nSb
integer,intent(in) :: nSt
integer,intent(in) :: nSsc
integer,intent(in) :: nS_sc
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(nSsc)
double precision,intent(in) :: rho(nBas,nBas,nSsc,nspin)
double precision,intent(in) :: Omega(nS_sc)
double precision,intent(in) :: rho(nBas,nBas,nS_sc,nspin)
! Local variables
@ -55,12 +55,12 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n
jb = jb + 1
chi = 0d0
do kc=1,nSsc
do kc=1,nS_sc
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
B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI_aaaa(i,j,b,a) + 2d0*lambda*chi
enddo
enddo
@ -80,12 +80,12 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n
jb = jb + 1
chi = 0d0
do kc=1,nSsc
do kc=1,nS_sc
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
B_lr(nSa+ia,nSa+jb) = B_lr(nSa+ia,nSa+jb) - lambda*ERI_bbbb(i,j,b,a) + 2d0*lambda*chi
enddo
enddo
@ -113,12 +113,12 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n
jb = jb + 1
chi = 0d0
do kc=1,nSsc
do kc=1,nS_sc
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,b,kc,1)*rho(a,j,kc,2)*Omega(kc)/eps
enddo
B_lr(ia,nSa+jb) = B_lr(ia,nSa+jb) - lambda*ERI_abab(i,a,b,j) + 4d0*lambda*chi
B_lr(ia,nSa+jb) = B_lr(ia,nSa+jb) - lambda*ERI_abab(i,a,b,j) + 2d0*lambda*chi
end do
end do
@ -137,12 +137,12 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,n
jb = jb + 1
chi = 0d0
do kc=1,nSsc
do kc=1,nS_sc
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
B_lr(nSa+ia,jb) = B_lr(nSa+ia,jb) - lambda*ERI_abab(b,j,i,a) + 2d0*lambda*chi
end do
end do

18
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,nSsc,lambda, &
e,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_RPA,rho_RPA,EcRPA,Omega,XpY,XmY)
subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, &
e,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcRPA,Omega,XpY,XmY)
! Compute linear response for unrestricted formalism
@ -21,7 +21,7 @@ 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
integer,intent(in) :: nS_sc
double precision,intent(in) :: lambda
double precision,intent(in) :: e(nBas,nspin)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
@ -29,8 +29,8 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,
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)
double precision,intent(in) :: OmRPA(nS_sc)
double precision,intent(in) :: rho_RPA(nBas,nBas,nS_sc,nspin)
! Local variables
@ -61,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(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)
call unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,A)
! Tamm-Dancoff approximation
@ -73,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(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)
call unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,B)
end if