LCPQ
/
quack
mirror of https://github.com/pfloos/quack
4
1
Fork 0
Code Releases Activity

fix bug in UG0W0

This commit is contained in:
Pierre-Francois Loos 2023-10-23 21:36:07 +02:00
parent 120377f611
commit 83e5c9b07c
10 changed files with 26 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
input/methods
View File

@ -1,5 +1,5 @@
# RHF UHF ROHF RMOM UMOM KS
T F F F F F
F T F F F F
# MP2* MP3
F F
# CCD pCCD DCD CCSD CCSD(T)
@ -13,7 +13,7 @@
# G0F2* evGF2* qsGF2* G0F3 evGF3
F F F F F
# G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW
F F F F F F
T F F F F F
# G0T0pp* evGTpp* qsGTpp* G0T0eh evGTeh qsGTeh
T F F F F F
F F F F F F
# * unrestricted version available

2
input/options
View File

@ -9,7 +9,7 @@
# GF: maxSCF thresh DIIS n_diis lin eta renorm reg
256 0.00001 T 5 F 0.0 0 F
# GW: maxSCF thresh DIIS n_diis lin eta TDA_W reg
256 0.00001 T 5 F 0.0 F F
256 0.00001 T 5 T 0.00367493 F F
# GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg
256 0.00001 T 5 F 0.0 F F
# ACFDT: AC Kx XBS

2
src/GT/G0T0pp.f90
View File

@ -66,7 +66,7 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
write(*,*)
write(*,*)'************************************************'
write(*,*)'| One-shot G0T0pp calculation |'
write(*,*)'| One-shot G0T0pp calculation |'
write(*,*)'************************************************'
write(*,*)

2
src/GT/GTeh_plot_self_energy.f90
View File

@ -33,7 +33,7 @@ subroutine GTeh_plot_self_energy(nBas,nC,nO,nV,nR,nS,eHF,eGT,Om,rhoL,rhoR)
! Broadening parameter
eta = 0.1d0
eta = 0.01d0
! Construct grid

2
src/GT/GTpp_plot_self_energy.f90
View File

@ -37,7 +37,7 @@ subroutine GTpp_plot_self_energy(nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,eGT,Om
! Broadening parameter
eta = 0.1d0
eta = 0.01d0
! Construct grid

2
src/GW/G0W0.f90
View File

@ -144,7 +144,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT
end if
call GW_plot_self_energy(nBas,nC,nO,nV,nR,nS,eHF,eHF,Om,rho)
! call GW_plot_self_energy(nBas,nC,nO,nV,nR,nS,eHF,eHF,Om,rho)
! Compute the RPA correlation energy

2
src/GW/GW_plot_self_energy.f90
View File

@ -32,7 +32,7 @@ subroutine GW_plot_self_energy(nBas,nC,nO,nV,nR,nS,eHF,eGW,Om,rho)
! Broadening parameter
eta = 0.1d0
eta = 0.01d0
! Construct grid

8
src/GW/UG0W0.f90
View File

@ -44,6 +44,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons
! Local variables
logical :: print_W = .true.
logical :: dRPA
integer :: is
integer :: ispin
double precision :: EcRPA
@ -76,6 +77,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons
! Initialization
EcRPA = 0d0
dRPA = .true.
! TDA for W
@ -104,11 +106,11 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons
! Compute screening !
!-------------------!
! Spin-conserving transition
! Spin-conserving transitions
ispin = 1
call phULR(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, &
call phULR(ispin,dRPA,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, &
eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,Om,rho,EcRPA,Om,XpY,XmY)
if(print_W) call print_excitation_energies('phRPA@UHF',5,nS_sc,Om)
@ -157,7 +159,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_cons
! Compute RPA correlation energy
call phULR(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, &
call phULR(ispin,dRPA,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, &
eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Om,rho,EcRPA,Om,XpY,XmY)
! Dump results

2
src/GW/UGW_excitation_density.f90
View File

@ -81,7 +81,7 @@ subroutine UGW_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ER
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
rho(p,q,ia,2) = rho(p,q,ia,2) + ERI_aabb(j,p,b,q)*XpY(ia,jb)
rho(p,q,ia,2) = rho(p,q,ia,2) + ERI_aabb(p,j,q,b)*XpY(ia,jb)
enddo
enddo

18
src/LR/phLR_oscillator_strength.f90
View File

@ -21,7 +21,7 @@ subroutine phLR_oscillator_strength(nBas,nC,nO,nV,nR,nS,maxS,dipole_int,Om,XpY,X
! Local variables
integer :: ia,jb,i,j,a,b
integer :: m,jb,i,j,a,b
integer :: ixyz
double precision,allocatable :: f(:,:)
@ -40,21 +40,21 @@ subroutine phLR_oscillator_strength(nBas,nC,nO,nV,nR,nS,maxS,dipole_int,Om,XpY,X
! Compute dipole moments and oscillator strengths
do ia=1,maxS
do m=1,maxS
do ixyz=1,ncart
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
f(ia,ixyz) = f(ia,ixyz) + dipole_int(j,b,ixyz)*XpY(ia,jb)
f(m,ixyz) = f(m,ixyz) + dipole_int(j,b,ixyz)*XpY(m,jb)
end do
end do
end do
end do
f(:,:) = sqrt(2d0)*f(:,:)
do ia=1,maxS
os(ia) = 2d0/3d0*Om(ia)*sum(f(ia,:)**2)
do m=1,maxS
os(m) = 2d0/3d0*Om(m)*sum(f(m,:)**2)
end do
write(*,*) '---------------------------------------------------------------'
@ -62,10 +62,14 @@ subroutine phLR_oscillator_strength(nBas,nC,nO,nV,nR,nS,maxS,dipole_int,Om,XpY,X
write(*,*) '---------------------------------------------------------------'
write(*,'(A3,5A12)') '#','X','Y','Z','dip. str.','osc. str.'
write(*,*) '---------------------------------------------------------------'
do ia=1,maxS
write(*,'(I3,5F12.6)') ia,(f(ia,ixyz),ixyz=1,ncart),sum(f(ia,:)**2),os(ia)
do m=1,maxS
write(*,'(I3,5F12.6)') m,(f(m,ixyz),ixyz=1,ncart),sum(f(m,:)**2),os(m)
end do
write(*,*) '---------------------------------------------------------------'
write(*,*)
! do m=1,maxS
! write(*,'(I3,3F12.6)') m,Om(m),os(m)
! end do
end subroutine