PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-23 04:43:42 +01:00
Code Releases Activity

Done with UACFDT

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-10-07 22:51:30 +02:00
parent efbe27068e
commit 1abe13c86d
14 changed files with 108 additions and 125 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
examples/molecule.He
View File

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

6
input/methods
View File

@ -1,7 +1,7 @@
# RHF UHF MOM # RHF UHF MOM
F T F F T F
# MP2* MP3 MP2-F12 # MP2* MP3 MP2-F12
F F F T F F
# CCD CCSD CCSD(T) # CCD CCSD CCSD(T)
F F F F F F
# drCCD rCCD lCCD pCCD # drCCD rCCD lCCD pCCD
@ -9,11 +9,11 @@
# CIS* CIS(D) CID CISD # CIS* CIS(D) CID CISD
F F F F F F F F
# RPA* RPAx* ppRPA # RPA* RPAx* ppRPA
T F F F F F
# G0F2 evGF2 G0F3 evGF3 # G0F2 evGF2 G0F3 evGF3
F F F F F F F F
# G0W0* evGW* qsGW # G0W0* evGW* qsGW
F F F T F F
# G0T0 evGT qsGT # G0T0 evGT qsGT
F F F F F F
# MCMP2 # MCMP2

4
input/options
View File

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

14
src/MBPT/G0W0.f90
View File

@ -202,9 +202,9 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
if(doACFDT) then if(doACFDT) then
write(*,*) '------------------------------------------------------' write(*,*) '--------------------------------------------------------------'
write(*,*) 'Adiabatic connection version of BSE correlation energy' write(*,*) ' Adiabatic connection version of BSE@UG0W0 correlation energy '
write(*,*) '------------------------------------------------------' write(*,*) '--------------------------------------------------------------'
write(*,*) write(*,*)
if(doXBS) then if(doXBS) then
@ -218,10 +218,10 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
write(*,*) write(*,*)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy (singlet) =',EcAC(1) write(*,'(2X,A50,F20.10)') 'AC@BSE@UG0W0 correlation energy (singlet) =',EcAC(1)
write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy (triplet) =',EcAC(2) write(*,'(2X,A50,F20.10)') 'AC@BSE@UG0W0 correlation energy (triplet) =',EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy =',EcAC(1) + EcAC(2) write(*,'(2X,A50,F20.10)') 'AC@BSE@UG0W0 correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 total energy =',ENuc + ERHF + EcAC(1) + EcAC(2) write(*,'(2X,A50,F20.10)') 'AC@BSE@UG0W0 total energy =',ENuc + ERHF + EcAC(1) + EcAC(2)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,*) write(*,*)

27
src/MBPT/UG0W0.f90
View File

@ -192,10 +192,10 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
write(*,*) write(*,*)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@BSE@G0W0 correlation energy (spin-conserved) =',EcBSE(1) write(*,'(2X,A50,F20.10)') 'Tr@BSE@UG0W0 correlation energy (spin-conserved) =',EcBSE(1)
write(*,'(2X,A50,F20.10)') 'Tr@BSE@G0W0 correlation energy (spin-flip) =',EcBSE(2) write(*,'(2X,A50,F20.10)') 'Tr@BSE@UG0W0 correlation energy (spin-flip) =',EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@BSE@G0W0 correlation energy =',EcBSE(1) + EcBSE(2) write(*,'(2X,A50,F20.10)') 'Tr@BSE@UG0W0 correlation energy =',EcBSE(1) + EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@BSE@G0W0 total energy =',ENuc + EUHF + EcBSE(1) + EcBSE(2) write(*,'(2X,A50,F20.10)') 'Tr@BSE@UG0W0 total energy =',ENuc + EUHF + EcBSE(1) + EcBSE(2)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,*) write(*,*)
@ -215,22 +215,15 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
end if end if
call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,spin_conserved,spin_flip,eta, & call unrestricted_ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,spin_conserved,spin_flip,eta, &
nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGW,EcAC) nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGW,EcAC)
if(exchange_kernel) then
EcAC(1) = 0.5d0*EcAC(1)
EcAC(2) = 1.5d0*EcAC(1)
end if
write(*,*) write(*,*)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy (spin-conserved) =',EcAC(1) write(*,'(2X,A50,F20.10)') 'AC@BSE@UG0W0 correlation energy (spin-conserved) =',EcAC(1)
write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy (spin-flip) =',EcAC(2) write(*,'(2X,A50,F20.10)') 'AC@BSE@UG0W0 correlation energy (spin-flip) =',EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy =',EcAC(1) + EcAC(2) write(*,'(2X,A50,F20.10)') 'AC@BSE@UG0W0 correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 total energy =',ENuc + EUHF + EcAC(1) + EcAC(2) write(*,'(2X,A50,F20.10)') 'AC@BSE@UG0W0 total energy =',ENuc + EUHF + EcAC(1) + EcAC(2)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,*) write(*,*)

70
src/MBPT/evUGW.f90
View File

@ -1,6 +1,6 @@
subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, &
G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc, & G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc, &
ERHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF,eG0W0) EUHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF,eG0W0)
! Perform self-consistent eigenvalue-only GW calculation ! Perform self-consistent eigenvalue-only GW calculation
@ -13,7 +13,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE
integer,intent(in) :: max_diis integer,intent(in) :: max_diis
double precision,intent(in) :: thresh double precision,intent(in) :: thresh
double precision,intent(in) :: ENuc double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF double precision,intent(in) :: EUHF
logical,intent(in) :: doACFDT logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel logical,intent(in) :: exchange_kernel
logical,intent(in) :: doXBS logical,intent(in) :: doXBS
@ -59,7 +59,6 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE
double precision :: EcRPA double precision :: EcRPA
double precision :: EcBSE(nspin) double precision :: EcBSE(nspin)
double precision :: EcAC(nspin) double precision :: EcAC(nspin)
double precision :: EcGM
double precision :: alpha double precision :: alpha
double precision,allocatable :: error_diis(:,:,:) double precision,allocatable :: error_diis(:,:,:)
double precision,allocatable :: e_diis(:,:,:) double precision,allocatable :: e_diis(:,:,:)
@ -192,7 +191,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE
! Print results ! Print results
call print_evUGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA) call print_evUGW(nBas,nO,nSCF,Conv,eHF,ENuc,EUHF,SigC,Z,eGW,EcRPA)
! Linear mixing or DIIS extrapolation ! Linear mixing or DIIS extrapolation
@ -256,50 +255,51 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE
call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, & call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, &
S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eGW,eGW,EcBSE) S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eGW,eGW,EcBSE)
! if(exchange_kernel) then if(exchange_kernel) then
! EcBSE(1) = 0.5d0*EcBSE(1) EcBSE(1) = 0.5d0*EcBSE(1)
! EcBSE(2) = 1.5d0*EcBSE(2) EcBSE(2) = 1.5d0*EcBSE(2)
! end if end if
! write(*,*) write(*,*)
! write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
! write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW correlation energy (singlet) =',EcBSE(1) write(*,'(2X,A50,F20.10)') 'Tr@BSE@evUGW correlation energy (spin-conserved) =',EcBSE(1)
! write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW correlation energy (triplet) =',EcBSE(2) write(*,'(2X,A50,F20.10)') 'Tr@BSE@evUGW correlation energy (spin-flip) =',EcBSE(2)
! write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW correlation energy =',EcBSE(1) + EcBSE(2) write(*,'(2X,A50,F20.10)') 'Tr@BSE@evUGW correlation energy =',EcBSE(1) + EcBSE(2)
! write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2) write(*,'(2X,A50,F20.10)') 'Tr@BSE@evUGW total energy =',ENuc + EUHF + EcBSE(1) + EcBSE(2)
! write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
! write(*,*) write(*,*)
! Compute the BSE correlation energy via the adiabatic connection ! Compute the BSE correlation energy via the adiabatic connection
! if(doACFDT) then if(doACFDT) then
! write(*,*) '------------------------------------------------------' write(*,*) '--------------------------------------------------------------'
! write(*,*) 'Adiabatic connection version of BSE correlation energy' write(*,*) ' Adiabatic connection version of BSE@evUGW correlation energy '
! write(*,*) '------------------------------------------------------' write(*,*) '--------------------------------------------------------------'
! write(*,*) write(*,*)
! if(doXBS) then if(doXBS) then
! write(*,*) '*** scaled screening version (XBS) ***' write(*,*) '*** scaled screening version (XBS) ***'
! write(*,*) write(*,*)
! end if end if
! call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcAC) call unrestricted_ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,spin_conserved,spin_flip, &
eta,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eGW,eGW,EcAC)
! write(*,*) write(*,*)
! write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy (singlet) =',EcAC(1) write(*,'(2X,A50,F20.10)') 'AC@BSE@evUGW correlation energy (spin-conserved) =',EcAC(1)
! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy (triplet) =',EcAC(2) write(*,'(2X,A50,F20.10)') 'AC@BSE@evUGW correlation energy (spin-flip) =',EcAC(2)
! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy =',EcAC(1) + EcAC(2) write(*,'(2X,A50,F20.10)') 'AC@BSE@evUGW correlation energy =',EcAC(1) + EcAC(2)
! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW total energy =',ENuc + ERHF + EcAC(1) + EcAC(2) write(*,'(2X,A50,F20.10)') 'AC@BSE@evUGW total energy =',ENuc + EUHF + EcAC(1) + EcAC(2)
! write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
! write(*,*) write(*,*)
! end if end if
endif endif

8
src/MBPT/unrestricted_Bethe_Salpeter.f90
View File

@ -65,6 +65,10 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,
nS_bb = nS(2) nS_bb = nS(2)
nS_sc = nS_aa + nS_bb nS_sc = nS_aa + nS_bb
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(OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(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))
!--------------------------! !--------------------------!
@ -126,10 +130,6 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,
! Memory allocation ! 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(OmBSE_sf(nS_sf),XpY_BSE_sf(nS_sf,nS_sf),XmY_BSE_sf(nS_sf,nS_sf)) allocate(OmBSE_sf(nS_sf),XpY_BSE_sf(nS_sf,nS_sf),XmY_BSE_sf(nS_sf,nS_sf))
! Compute spin-flip BSE excitation energies ! Compute spin-flip BSE excitation energies

4
src/QuAcK/QuAcK.f90
View File

@ -652,12 +652,12 @@ program QuAcK
call cpu_time(start_RPA) call cpu_time(start_RPA)
if(unrestricted) then if(unrestricted) then
call UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & call URPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, &
ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa,dipole_int_bb,eHF,cHF,S) ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa,dipole_int_bb,eHF,cHF,S)
else else
call dRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF) call RPA(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF)
end if end if
call cpu_time(end_RPA) call cpu_time(end_RPA)

1
src/RPA/ACFDT.f90
View File

@ -126,7 +126,6 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nB
if(triplet) then if(triplet) then
ispin = 2 ispin = 2
isp_W = 1
write(*,*) '--------------' write(*,*) '--------------'
write(*,*) 'Triplet states' write(*,*) 'Triplet states'

4
src/RPA/dRPA.f90 → src/RPA/RPA.f90
View File

@ -1,4 +1,4 @@
subroutine dRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) subroutine RPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
! Perform a direct random phase approximation calculation ! Perform a direct random phase approximation calculation
@ -135,4 +135,4 @@ subroutine dRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR
end if end if
end subroutine dRPA end subroutine RPA

32
src/RPA/UdRPA.f90 → src/RPA/URPA.f90
View File

@ -1,4 +1,4 @@
subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & subroutine URPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, &
ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,e,c,S) ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,e,c,S)
! Perform random phase approximation calculation with exchange (aka TDHF) in the unrestricted formalism ! Perform random phase approximation calculation with exchange (aka TDHF) in the unrestricted formalism
@ -68,7 +68,7 @@ subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n
! Initialization ! Initialization
EcRPA(:) = 0d0 EcRPA(:) = 0d0
EcAC(:) = 0d0 EcAC(:) = 0d0
! Spin-conserved transitions ! Spin-conserved transitions
@ -90,6 +90,7 @@ subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n
call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, &
c,S,Omega_sc,XpY_sc,XmY_sc) c,S,Omega_sc,XpY_sc,XmY_sc)
deallocate(Omega_sc,XpY_sc,XmY_sc)
endif endif
@ -113,6 +114,8 @@ subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n
call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, & call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, &
c,S,Omega_sf,XpY_sf,XmY_sf) c,S,Omega_sf,XpY_sf,XmY_sf)
deallocate(Omega_sf,XpY_sf,XmY_sf)
endif endif
if(exchange_kernel) then if(exchange_kernel) then
@ -135,30 +138,23 @@ subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n
if(doACFDT) then if(doACFDT) then
write(*,*) '-------------------------------------------------------' write(*,*) '---------------------------------------------------------'
write(*,*) 'Adiabatic connection version of RPA correlation energy' write(*,*) ' Adiabatic connection version of URPA correlation energy '
write(*,*) '-------------------------------------------------------' write(*,*) '---------------------------------------------------------'
write(*,*) write(*,*)
call unrestricted_ACFDT(exchange_kernel,.false.,.false.,.false.,.false.,.false.,spin_conserved,spin_flip,eta, & call unrestricted_ACFDT(exchange_kernel,.false.,.true.,.false.,TDA,.false.,spin_conserved,spin_flip,eta, &
nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,e,e,EcAC) nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,e,e,EcAC)
if(exchange_kernel) then
EcAC(1) = 0.5d0*EcAC(1)
EcAC(2) = 1.5d0*EcAC(2)
end if
write(*,*) write(*,*)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'AC@RPA correlation energy (spin-conserved) =',EcAC(1) write(*,'(2X,A50,F20.10)') 'AC@URPA correlation energy (spin-conserved) =',EcAC(1)
write(*,'(2X,A50,F20.10)') 'AC@RPA correlation energy (spin-flip) =',EcAC(2) write(*,'(2X,A50,F20.10)') 'AC@URPA correlation energy (spin-flip) =',EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@RPA correlation energy =',EcAC(1) + EcAC(2) write(*,'(2X,A50,F20.10)') 'AC@URPA correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@RPA total energy =',ENuc + EUHF + EcAC(1) + EcAC(2) write(*,'(2X,A50,F20.10)') 'AC@URPA total energy =',ENuc + EUHF + EcAC(1) + EcAC(2)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,*) write(*,*)
end if end if
end subroutine UdRPA end subroutine URPA

23
src/RPA/URPAx.f90
View File

@ -139,27 +139,20 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n
if(doACFDT) then if(doACFDT) then
write(*,*) '--------------------------------------------------------' write(*,*) '----------------------------------------------------------'
write(*,*) 'Adiabatic connection version of URPAx correlation energy' write(*,*) ' Adiabatic connection version of URPAx correlation energy '
write(*,*) '--------------------------------------------------------' write(*,*) '----------------------------------------------------------'
write(*,*) write(*,*)
call unrestricted_ACFDT(exchange_kernel,.false.,.false.,.false.,.false.,.false.,spin_conserved,spin_flip,eta, & call unrestricted_ACFDT(exchange_kernel,.false.,.false.,.false.,TDA,.false.,spin_conserved,spin_flip,eta, &
nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,e,e,EcAC) nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,e,e,EcAC)
if(exchange_kernel) then
EcAC(1) = 0.5d0*EcAC(1)
EcAC(2) = 1.5d0*EcAC(2)
end if
write(*,*) write(*,*)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy (spin-conserved) =',EcAC(1) write(*,'(2X,A50,F20.10)') 'AC@URPAx correlation energy (spin-conserved) =',EcAC(1)
write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy (spin-flip) =',EcAC(2) write(*,'(2X,A50,F20.10)') 'AC@URPAx correlation energy (spin-flip) =',EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy =',EcAC(1) + EcAC(2) write(*,'(2X,A50,F20.10)') 'AC@URPAx correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@RPAx total energy =',ENuc + EUHF + EcAC(1) + EcAC(2) write(*,'(2X,A50,F20.10)') 'AC@URPAx total energy =',ENuc + EUHF + EcAC(1) + EcAC(2)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,*) write(*,*)

26
src/RPA/unrestricted_ACFDT.f90
View File

@ -43,7 +43,7 @@ subroutine unrestricted_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_cons
double precision,allocatable :: OmRPA(:) double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:) double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:) double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:) double precision,allocatable :: rho_RPA(:,:,:,:)
integer :: nS_aa,nS_bb,nS_sc integer :: nS_aa,nS_bb,nS_sc
double precision,allocatable :: Omega_sc(:) double precision,allocatable :: Omega_sc(:)
@ -62,7 +62,6 @@ subroutine unrestricted_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_cons
! Memory allocation ! Memory allocation
allocate(Ec(nAC,nspin)) allocate(Ec(nAC,nspin))
allocate(OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc))
! Antisymmetrized kernel version ! Antisymmetrized kernel version
@ -88,6 +87,12 @@ subroutine unrestricted_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_cons
nS_bb = nS(2) nS_bb = nS(2)
nS_sc = nS_aa + nS_bb nS_sc = nS_aa + nS_bb
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(OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin))
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, & 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,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
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,rho_RPA)
@ -122,10 +127,10 @@ subroutine unrestricted_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_cons
end if end if
call unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,e, & call unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,e, &
ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,Omega_sc,XpY_sc,XmY_sc) ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcAC(ispin),Omega_sc,XpY_sc,XmY_sc)
call unrestricted_ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf, & call unrestricted_ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc, &
ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_sc,XmY_sc,Ec(iAC,ispin)) ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_sc,XmY_sc,Ec(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)
@ -149,14 +154,9 @@ subroutine unrestricted_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_cons
if(spin_flip) then if(spin_flip) then
ispin = 2 ispin = 2
isp_W = 1
! Memory allocation ! 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(Omega_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf)) allocate(Omega_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf))
write(*,*) '--------------------' write(*,*) '--------------------'
@ -180,10 +180,10 @@ subroutine unrestricted_ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_cons
end if end if
call unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,nS_sc,lambda,e, & call unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sc,lambda,e, &
ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,Omega_sf,XpY_sf,XmY_sf) ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcAC(ispin),Omega_sf,XpY_sf,XmY_sf)
call unrestricted_ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc, & call unrestricted_ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf, &
ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_sf,XmY_sf,Ec(iAC,ispin)) ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_sf,XmY_sf,Ec(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)

12
src/RPA/unrestricted_ACFDT_correlation_energy.f90
View File

@ -28,7 +28,7 @@ subroutine unrestricted_ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,n
! Local variables ! Local variables
integer :: i,j,a,b integer :: i,j,a,b
integer :: ia,jb,kc integer :: ia,jb
double precision :: delta_Kx double precision :: delta_Kx
double precision,allocatable :: Ap(:,:) double precision,allocatable :: Ap(:,:)
double precision,allocatable :: Bp(:,:) double precision,allocatable :: Bp(:,:)
@ -183,7 +183,9 @@ subroutine unrestricted_ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,n
end do end do
end do end do
ia = 0 ! abba block
ia = 0
do i=nC(1)+1,nO(1) do i=nC(1)+1,nO(1)
do a=nO(2)+1,nBas-nR(2) do a=nO(2)+1,nBas-nR(2)
ia = ia + 1 ia = ia + 1
@ -224,8 +226,8 @@ subroutine unrestricted_ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,n
X(:,:) = 0.5d0*(XpY(:,:) + XmY(:,:)) X(:,:) = 0.5d0*(XpY(:,:) + XmY(:,:))
Y(:,:) = 0.5d0*(XpY(:,:) - XmY(:,:)) Y(:,:) = 0.5d0*(XpY(:,:) - XmY(:,:))
EcAC = trace_matrix(nS,matmul(X,matmul(Bp,transpose(Y))) + matmul(Y,matmul(Bp,transpose(X)))) & EcAC = trace_matrix(nSt,matmul(X,matmul(Bp,transpose(Y))) + matmul(Y,matmul(Bp,transpose(X)))) &
+ trace_matrix(nS,matmul(X,matmul(Ap,transpose(X))) + matmul(Y,matmul(Ap,transpose(Y)))) & + trace_matrix(nSt,matmul(X,matmul(Ap,transpose(X))) + matmul(Y,matmul(Ap,transpose(Y)))) &
- trace_matrix(nS,Ap) - trace_matrix(nSt,Ap)
end subroutine unrestricted_ACFDT_correlation_energy end subroutine unrestricted_ACFDT_correlation_energy