From 83e5c9b07c5bea380eae0073d00b7ec50a60b956 Mon Sep 17 00:00:00 2001 From: pfloos Date: Mon, 23 Oct 2023 21:36:07 +0200 Subject: [PATCH] fix bug in UG0W0 --- input/methods | 6 +++--- input/options | 2 +- src/GT/G0T0pp.f90 | 2 +- src/GT/GTeh_plot_self_energy.f90 | 2 +- src/GT/GTpp_plot_self_energy.f90 | 2 +- src/GW/G0W0.f90 | 2 +- src/GW/GW_plot_self_energy.f90 | 2 +- src/GW/UG0W0.f90 | 8 +++++--- src/GW/UGW_excitation_density.f90 | 2 +- src/LR/phLR_oscillator_strength.f90 | 18 +++++++++++------- 10 files changed, 26 insertions(+), 20 deletions(-) diff --git a/input/methods b/input/methods index cc87fea..4d84f78 100644 --- a/input/methods +++ b/input/methods @@ -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 diff --git a/input/options b/input/options index bfb8d30..bcb52d7 100644 --- a/input/options +++ b/input/options @@ -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 diff --git a/src/GT/G0T0pp.f90 b/src/GT/G0T0pp.f90 index 5c28f4d..7a1e94c 100644 --- a/src/GT/G0T0pp.f90 +++ b/src/GT/G0T0pp.f90 @@ -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(*,*) diff --git a/src/GT/GTeh_plot_self_energy.f90 b/src/GT/GTeh_plot_self_energy.f90 index da9f9c6..3245a06 100644 --- a/src/GT/GTeh_plot_self_energy.f90 +++ b/src/GT/GTeh_plot_self_energy.f90 @@ -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 diff --git a/src/GT/GTpp_plot_self_energy.f90 b/src/GT/GTpp_plot_self_energy.f90 index 87c4039..3d3b8d5 100644 --- a/src/GT/GTpp_plot_self_energy.f90 +++ b/src/GT/GTpp_plot_self_energy.f90 @@ -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 diff --git a/src/GW/G0W0.f90 b/src/GW/G0W0.f90 index e08f3c8..94f6c20 100644 --- a/src/GW/G0W0.f90 +++ b/src/GW/G0W0.f90 @@ -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 diff --git a/src/GW/GW_plot_self_energy.f90 b/src/GW/GW_plot_self_energy.f90 index 7a35911..44b0dea 100644 --- a/src/GW/GW_plot_self_energy.f90 +++ b/src/GW/GW_plot_self_energy.f90 @@ -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 diff --git a/src/GW/UG0W0.f90 b/src/GW/UG0W0.f90 index 1cc8b86..1d3c9e6 100644 --- a/src/GW/UG0W0.f90 +++ b/src/GW/UG0W0.f90 @@ -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 diff --git a/src/GW/UGW_excitation_density.f90 b/src/GW/UGW_excitation_density.f90 index 571051b..ebbf439 100644 --- a/src/GW/UGW_excitation_density.f90 +++ b/src/GW/UGW_excitation_density.f90 @@ -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 diff --git a/src/LR/phLR_oscillator_strength.f90 b/src/LR/phLR_oscillator_strength.f90 index 5085ce1..f8cb55d 100644 --- a/src/LR/phLR_oscillator_strength.f90 +++ b/src/LR/phLR_oscillator_strength.f90 @@ -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