From ba01dd71faf2ed3776d04a4e51636cd70b065bf0 Mon Sep 17 00:00:00 2001 From: pfloos Date: Tue, 15 Oct 2024 02:14:51 +0800 Subject: [PATCH] looking for a bug in UGW_phBSE --- src/GW/UG0W0.f90 | 8 ++-- src/GW/UGW_phACFDT.f90 | 15 +++---- src/GW/UGW_phBSE.f90 | 67 +++++++++++++++++++--------- src/GW/UGW_phBSE_static_kernel_A.f90 | 35 +++++++-------- src/GW/UGW_phBSE_static_kernel_B.f90 | 34 +++++++------- src/LR/phULR.f90 | 10 ----- src/LR/phULR_transition_vectors.f90 | 14 +++--- src/LR/print_excitation_energies.f90 | 2 +- 8 files changed, 95 insertions(+), 90 deletions(-) diff --git a/src/GW/UG0W0.f90 b/src/GW/UG0W0.f90 index 77dd2a4..44773a2 100644 --- a/src/GW/UG0W0.f90 +++ b/src/GW/UG0W0.f90 @@ -115,8 +115,8 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD isp_W = 1 - call phULR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) - if(.not.TDA) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) + call phULR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) + if(.not.TDA_W) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) call phULR(TDA_W,nSa,nSb,nSt,Aph,Bph,EcRPA,Om,XpY,XmY) @@ -172,8 +172,8 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD ! Compute RPA correlation energy - call phULR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) - if(.not.TDA) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) + call phULR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) + if(.not.TDA_W) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) call phULR(TDA_W,nSa,nSb,nSt,Aph,Bph,EcRPA,Om,XpY,XmY) diff --git a/src/GW/UGW_phACFDT.f90 b/src/GW/UGW_phACFDT.f90 index 91bbc99..9801c0c 100644 --- a/src/GW/UGW_phACFDT.f90 +++ b/src/GW/UGW_phACFDT.f90 @@ -118,10 +118,8 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,spin_conserved,spin_flip, call phULR(TDA_W,nS_aa,nS_bb,nS_sc,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA) call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) - call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,KA) - call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,KB) + call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,OmRPA,rho_RPA,KA) + call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,OmRPA,rho_RPA,KB) ! Spin-conserved manifold @@ -152,10 +150,8 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,spin_conserved,spin_flip, call phULR(TDA_W,nS_aa,nS_bb,nS_sc,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA) call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) - call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,KA) - call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda, & - ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,KB) + call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,OmRPA,rho_RPA,KA) + call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,OmRPA,rho_RPA,KB) end if @@ -218,6 +214,9 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,spin_conserved,spin_flip, call phULR(TDA_W,nS_aa,nS_bb,nS_sc,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA) call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) + call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,OmRPA,rho_RPA,KA) + call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,OmRPA,rho_RPA,KB) + end if call phULR_A(isp_W,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,lambda,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) diff --git a/src/GW/UGW_phBSE.f90 b/src/GW/UGW_phBSE.f90 index d2430aa..1bf4cdc 100644 --- a/src/GW/UGW_phBSE.f90 +++ b/src/GW/UGW_phBSE.f90 @@ -37,6 +37,13 @@ subroutine UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_fli integer :: ispin integer :: isp_W + logical :: dRPA = .false. + logical :: dRPA_W = .true. + + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) + double precision,allocatable :: KA(:,:) + double precision,allocatable :: KB(:,:) double precision :: EcRPA double precision,allocatable :: OmRPA(:) @@ -64,30 +71,35 @@ subroutine UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_fli 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)) - -!-----! -! TDA ! -!-----! +! TDA if(TDA) then write(*,*) 'Tamm-Dancoff approximation activated in phBSE!' write(*,*) end if +! Initialization + + EcBSE(:) = 0d0 + !--------------------------! ! Spin-conserved screening ! !--------------------------! isp_W = 1 - EcRPA = 0d0 ! Compute spin-conserved RPA screening - call phULR(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) + allocate(Aph(nS_sc,nS_sc),Bph(nS_sc,nS_sc)) + allocate(OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin)) + call phULR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) + if(.not.TDA_W) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) + + call phULR(TDA_W,nS_aa,nS_bb,nS_sc,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA) call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) + + deallocate(Aph,Bph) !----------------------------! ! Spin-conserved excitations ! @@ -96,14 +108,23 @@ subroutine UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_fli if(spin_conserved) then ispin = 1 - EcBSE(ispin) = 0d0 + allocate(Aph(nS_sc,nS_sc),Bph(nS_sc,nS_sc),KA(nS_sc,nS_sc),KB(nS_sc,nS_sc)) allocate(OmBSE(nS_sc),XpY_BSE(nS_sc,nS_sc),XmY_BSE(nS_sc,nS_sc)) ! Compute spin-conserved BSE excitation energies - call phULR(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,OmRPA,rho_RPA,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE) + call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) + if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) + + call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,OmRPA,rho_RPA,KA) + if(.not.TDA) call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,OmRPA,rho_RPA,KB) + + Aph(:,:) = Aph(:,:) + KA(:,:) +! if(.not.TDA) Bph(:,:) = Bph(:,:) + KB(:,:) + + call phULR(TDA,nS_aa,nS_bb,nS_sc,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE) + call print_excitation_energies('phBSE@GW@UHF','spin-conserved',nS_sc,OmBSE) call phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & cW,S,OmBSE,XpY_BSE,XmY_BSE) @@ -117,6 +138,7 @@ subroutine UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_fli eW,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, & OmRPA,rho_RPA,OmBSE,XpY_BSE,XmY_BSE) + deallocate(Aph,Bph,KA,KB) deallocate(OmBSE,XpY_BSE,XmY_BSE) end if @@ -128,17 +150,21 @@ subroutine UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_fli if(spin_flip) then ispin = 2 - EcBSE(ispin) = 0d0 - ! Memory allocation - + allocate(Aph(nS_sf,nS_sf),Bph(nS_sf,nS_sf),KA(nS_sf,nS_sf),KB(nS_sf,nS_sf)) allocate(OmBSE(nS_sf),XpY_BSE(nS_sf,nS_sf),XmY_BSE(nS_sf,nS_sf)) - ! Compute spin-flip BSE excitation energies + ! Compute spin-conserved BSE excitation energies + call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,1d0,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) + if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) - call phULR(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,OmRPA,rho_RPA,EcBSE(ispin), & - OmBSE,XpY_BSE,XmY_BSE) + call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sf,1d0,OmRPA,rho_RPA,KA) + if(.not.TDA) call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sf,1d0,OmRPA,rho_RPA,KB) + + Aph(:,:) = Aph(:,:) + KA(:,:) + if(.not.TDA) Bph(:,:) = Bph(:,:) + KB(:,:) + + call phULR(TDA,nS_aa,nS_bb,nS_sc,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE) call print_excitation_energies('phBSE@GW@UHF','spin-flip',nS_sf,OmBSE) call phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, & @@ -153,17 +179,16 @@ subroutine UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_fli eW,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, & OmRPA,rho_RPA,OmBSE,XpY_BSE,XmY_BSE) + deallocate(Aph,Bph,KA,KB) deallocate(OmBSE,XpY_BSE,XmY_BSE) end if - ! Scale properly correlation energy if exchange is included in interaction kernel if(exchange_kernel) then - EcBSE(1) = 0.5d0*EcBSE(1) - EcBSE(2) = 0.5d0*EcBSE(2) + EcBSE(:) = 0.5d0*EcBSE(:) else diff --git a/src/GW/UGW_phBSE_static_kernel_A.f90 b/src/GW/UGW_phBSE_static_kernel_A.f90 index fbd04b1..ef08ef3 100644 --- a/src/GW/UGW_phBSE_static_kernel_A.f90 +++ b/src/GW/UGW_phBSE_static_kernel_A.f90 @@ -1,5 +1,4 @@ -subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,eGW, & - ERI_aaaa,ERI_aabb,ERI_bbbb,Omega,rho,A_lr) +subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,Om,rho,KA) ! Compute the extra term for Bethe-Salpeter equation for linear response in the unrestricted formalism @@ -20,11 +19,7 @@ subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s integer,intent(in) :: nS_sc double precision,intent(in) :: eta double precision,intent(in) :: lambda - 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,intent(in) :: Omega(nS_sc) + double precision,intent(in) :: Om(nS_sc) double precision,intent(in) :: rho(nBas,nBas,nS_sc,nspin) ! Local variables @@ -35,7 +30,7 @@ subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s ! Output variables - double precision,intent(out) :: A_lr(nSt,nSt) + double precision,intent(out) :: KA(nSt,nSt) !--------------------------------------------------! ! Build BSE matrix for spin-conserving transitions ! @@ -56,11 +51,11 @@ subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s chi = 0d0 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 + eps = Om(kc)**2 + eta**2 + chi = chi + rho(i,j,kc,1)*rho(a,b,kc,1)*Om(kc)/eps end do - A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_aaaa(i,b,j,a) + 2d0*lambda*chi + KA(ia,jb) = KA(ia,jb) + 2d0*lambda*chi end do end do @@ -80,11 +75,11 @@ subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s chi = 0d0 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 + eps = Om(kc)**2 + eta**2 + chi = chi + rho(i,j,kc,2)*rho(a,b,kc,2)*Om(kc)/eps end do - A_lr(nSa+ia,nSa+jb) = A_lr(nSa+ia,nSa+jb) - lambda*ERI_bbbb(i,b,j,a) + 2d0*lambda*chi + KA(nSa+ia,nSa+jb) = KA(nSa+ia,nSa+jb) + 2d0*lambda*chi end do end do @@ -112,11 +107,11 @@ subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s chi = 0d0 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 + eps = Om(kc)**2 + eta**2 + chi = chi + rho(i,j,kc,1)*rho(a,b,kc,2)*Om(kc)/eps end do - A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_aabb(i,b,j,a) + 2d0*lambda*chi + KA(ia,jb) = KA(ia,jb) + 2d0*lambda*chi end do end do @@ -136,11 +131,11 @@ subroutine UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s chi = 0d0 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 + eps = Om(kc)**2 + eta**2 + chi = chi + rho(i,j,kc,2)*rho(a,b,kc,1)*Om(kc)/eps end do - A_lr(nSa+ia,nSa+jb) = A_lr(nSa+ia,nSa+jb) - lambda*ERI_aabb(b,i,a,j) + 2d0*lambda*chi + KA(nSa+ia,nSa+jb) = KA(nSa+ia,nSa+jb) + 2d0*lambda*chi end do end do diff --git a/src/GW/UGW_phBSE_static_kernel_B.f90 b/src/GW/UGW_phBSE_static_kernel_B.f90 index c357ecc..c374bf2 100644 --- a/src/GW/UGW_phBSE_static_kernel_B.f90 +++ b/src/GW/UGW_phBSE_static_kernel_B.f90 @@ -1,5 +1,4 @@ -subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, & - ERI_aaaa,ERI_aabb,ERI_bbbb,Omega,rho,B_lr) +subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,Om,rho,KB) ! Compute the extra term for Bethe-Salpeter equation for linear response @@ -20,10 +19,7 @@ subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s 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) :: Omega(nS_sc) + double precision,intent(in) :: Om(nS_sc) double precision,intent(in) :: rho(nBas,nBas,nS_sc,nspin) ! Local variables @@ -34,7 +30,7 @@ subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s ! Output variables - double precision,intent(out) :: B_lr(nSt,nSt) + double precision,intent(out) :: KB(nSt,nSt) !--------------------------------------------------! ! Build BSE matrix for spin-conserving transitions ! @@ -55,11 +51,11 @@ subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s chi = 0d0 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 + eps = Om(kc)**2 + eta**2 + chi = chi + rho(i,b,kc,1)*rho(a,j,kc,1)*Om(kc)/eps end do - B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI_aaaa(i,j,b,a) + 2d0*lambda*chi + KB(ia,jb) = KB(ia,jb) + 2d0*lambda*chi end do end do @@ -80,11 +76,11 @@ subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s chi = 0d0 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 + eps = Om(kc)**2 + eta**2 + chi = chi + rho(i,b,kc,2)*rho(a,j,kc,2)*Om(kc)/eps end do - B_lr(nSa+ia,nSa+jb) = B_lr(nSa+ia,nSa+jb) - lambda*ERI_bbbb(i,j,b,a) + 2d0*lambda*chi + KB(nSa+ia,nSa+jb) = KB(nSa+ia,nSa+jb) + 2d0*lambda*chi end do end do @@ -113,11 +109,11 @@ subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s chi = 0d0 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 + eps = Om(kc)**2 + eta**2 + chi = chi + rho(i,b,kc,1)*rho(a,j,kc,2)*Om(kc)/eps end do - B_lr(ia,nSa+jb) = B_lr(ia,nSa+jb) - lambda*ERI_aabb(i,j,b,a) + 2d0*lambda*chi + KB(ia,nSa+jb) = KB(ia,nSa+jb) + 2d0*lambda*chi end do end do @@ -137,11 +133,11 @@ subroutine UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_s chi = 0d0 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 + eps = Om(kc)**2 + eta**2 + chi = chi + rho(i,b,kc,2)*rho(a,j,kc,1)*Om(kc)/eps end do - B_lr(nSa+ia,jb) = B_lr(nSa+ia,jb) - lambda*ERI_aabb(j,i,a,b) + 2d0*lambda*chi + KB(nSa+ia,jb) = KB(nSa+ia,jb) + 2d0*lambda*chi end do end do diff --git a/src/LR/phULR.f90 b/src/LR/phULR.f90 index 8f9c640..58a6ff5 100644 --- a/src/LR/phULR.f90 +++ b/src/LR/phULR.f90 @@ -34,12 +34,6 @@ subroutine phULR(TDA,nSa,nSb,nSt,Aph,Bph,EcRPA,Om,XpY,XmY) allocate(ApB(nSt,nSt),AmB(nSt,nSt),AmBSq(nSt,nSt),AmBIv(nSt,nSt),Z(nSt,nSt)) -! Build A and B matrices - -! if(BSE) & -! call UGW_phBSE_static_kernel_A(ispin,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,e, & -! ERI_aaaa,ERI_aabb,ERI_bbbb,Om,rho,Aph) - ! Tamm-Dancoff approximation if(TDA) then @@ -51,10 +45,6 @@ subroutine phULR(TDA,nSa,nSb,nSt,Aph,Bph,EcRPA,Om,XpY,XmY) else -! if(BSE) & -! call UGW_phBSE_static_kernel_B(ispin,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, & -! ERI_aaaa,ERI_aabb,ERI_bbbb,Om,rho,Bph) - ApB(:,:) = Aph(:,:) + Bph(:,:) AmB(:,:) = Aph(:,:) - Bph(:,:) diff --git a/src/LR/phULR_transition_vectors.f90 b/src/LR/phULR_transition_vectors.f90 index 2e0f128..b01ec2c 100644 --- a/src/LR/phULR_transition_vectors.f90 +++ b/src/LR/phULR_transition_vectors.f90 @@ -1,4 +1,4 @@ -subroutine phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole_int_aa,dipole_int_bb,c,S,Omega,XpY,XmY) +subroutine phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole_int_aa,dipole_int_bb,c,S,Om,XpY,XmY) ! Print transition vectors for linear response calculation @@ -21,14 +21,14 @@ subroutine phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole double precision :: dipole_int_bb(nBas,nBas,ncart) double precision,intent(in) :: c(nBas,nBas,nspin) double precision,intent(in) :: S(nBas,nBas) - double precision,intent(in) :: Omega(nSt) + double precision,intent(in) :: Om(nSt) double precision,intent(in) :: XpY(nSt,nSt) double precision,intent(in) :: XmY(nSt,nSt) ! Local variables integer :: ia,jb,j,b - integer :: maxS = 20 + integer :: maxS = 10 double precision,parameter :: thres_vec = 0.1d0 double precision,allocatable :: X(:) double precision,allocatable :: Y(:) @@ -44,11 +44,11 @@ subroutine phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole os(:) = 0d0 if(ispin == 1) call phULR_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS, & - dipole_int_aa,dipole_int_bb,Omega,XpY,XmY,os) + dipole_int_aa,dipole_int_bb,Om,XpY,XmY,os) ! Compute - call S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S,Omega,XpY,XmY,S2) + call S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S,Om,XpY,XmY,S2) ! Print details about spin-conserved excitations @@ -61,7 +61,7 @@ subroutine phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole print*,'-------------------------------------------------------------' write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') & - ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV',' f = ',os(ia),' = ',S2(ia) + ' Excitation n. ',ia,': ',Om(ia)*HaToeV,' eV',' f = ',os(ia),' = ',S2(ia) print*,'-------------------------------------------------------------' ! Spin-up transitions @@ -117,7 +117,7 @@ subroutine phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole print*,'-------------------------------------------------------------' write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') & - ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV',' f = ',os(ia),' = ',S2(ia) + ' Excitation n. ',ia,': ',Om(ia)*HaToeV,' eV',' f = ',os(ia),' = ',S2(ia) print*,'-------------------------------------------------------------' ! Spin-up transitions diff --git a/src/LR/print_excitation_energies.f90 b/src/LR/print_excitation_energies.f90 index 4c67600..8d1448c 100644 --- a/src/LR/print_excitation_energies.f90 +++ b/src/LR/print_excitation_energies.f90 @@ -14,7 +14,7 @@ subroutine print_excitation_energies(method,manifold,nS,Om) ! Local variables - integer,parameter :: maxS = 20 + integer,parameter :: maxS = 10 integer :: m write(*,*)