diff --git a/src/GW/GW_ppBSE_dynamic_kernel_B.f90 b/src/GW/GW_ppBSE_dynamic_kernel_B.f90 index 424ded1..c1b5784 100644 --- a/src/GW/GW_ppBSE_dynamic_kernel_B.f90 +++ b/src/GW/GW_ppBSE_dynamic_kernel_B.f90 @@ -25,6 +25,7 @@ subroutine GW_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambd ! Local variables + double precision,external :: Kronecker_delta double precision :: dem,num integer :: m integer :: a,b,i,j @@ -62,7 +63,7 @@ subroutine GW_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambd dem = eGW(j) - Om(m) - eGW(a) num = rho(b,i,m)*rho(a,j,m) - KB_dyn(ab,ij) = KB_dyn(ab,ij) - num*dem/(dem**2 + eta**2) + KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2) dem = eGW(i) - Om(m) - eGW(a) num = rho(a,i,m)*rho(b,j,m) @@ -72,7 +73,9 @@ subroutine GW_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambd dem = eGW(i) - Om(m) - eGW(b) num = rho(b,i,m)*rho(a,j,m) - KB_dyn(ab,ij) = KB_dyn(ab,ij) - num*dem/(dem**2 + eta**2) + KB_dyn(ab,ij) = KB_dyn(ab,ij) + num*dem/(dem**2 + eta**2) + + KB_dyn(ab,ij) = 2d0*KB_dyn(ab,ij)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j))) end do @@ -118,6 +121,8 @@ subroutine GW_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambd KB_dyn(ab,ij) = KB_dyn(ab,ij) - num*dem/(dem**2 + eta**2) + KB_dyn(ab,ij) = 2d0*KB_dyn(ab,ij) + end do end do diff --git a/src/GW/GW_ppBSE_dynamic_kernel_C.f90 b/src/GW/GW_ppBSE_dynamic_kernel_C.f90 index 439f26f..6c6b684 100644 --- a/src/GW/GW_ppBSE_dynamic_kernel_C.f90 +++ b/src/GW/GW_ppBSE_dynamic_kernel_C.f90 @@ -24,6 +24,7 @@ subroutine GW_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,eG ! Local variables + double precision,external :: Kronecker_delta double precision :: dem,num integer :: m integer :: a,b,c,d @@ -64,8 +65,8 @@ subroutine GW_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,eG dem = OmBSE - eGW(c) - Om(m) - eGW(a) num = rho(b,c,m)*rho(a,d,m) - KC_dyn(ab,cd) = KC_dyn(ab,cd) - num*dem/(dem**2 + eta**2) - ZC_dyn(ab,cd) = ZC_dyn(ab,cd) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 + KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) + ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 dem = OmBSE - eGW(d) - Om(m) - eGW(a) num = rho(a,c,m)*rho(b,d,m) @@ -76,8 +77,11 @@ subroutine GW_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,eG dem = OmBSE - eGW(d) - Om(m) - eGW(b) num = rho(b,c,m)*rho(a,d,m) - KC_dyn(ab,cd) = KC_dyn(ab,cd) - num*dem/(dem**2 + eta**2) - ZC_dyn(ab,cd) = ZC_dyn(ab,cd) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 + KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2) + ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 + + KC_dyn(ab,cd) = 2d0*KC_dyn(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) + ZC_dyn(ab,cd) = 2d0*ZC_dyn(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) end do @@ -127,6 +131,9 @@ subroutine GW_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,eG KC_dyn(ab,cd) = KC_dyn(ab,cd) - num*dem/(dem**2 + eta**2) ZC_dyn(ab,cd) = ZC_dyn(ab,cd) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 + KC_dyn(ab,cd) = 2d0*KC_dyn(ab,cd) + ZC_dyn(ab,cd) = 2d0*ZC_dyn(ab,cd) + end do end do diff --git a/src/GW/GW_ppBSE_dynamic_kernel_D.f90 b/src/GW/GW_ppBSE_dynamic_kernel_D.f90 index c8433bf..ed95489 100644 --- a/src/GW/GW_ppBSE_dynamic_kernel_D.f90 +++ b/src/GW/GW_ppBSE_dynamic_kernel_D.f90 @@ -24,6 +24,7 @@ subroutine GW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,eG ! Local variables + double precision,external :: Kronecker_delta double precision :: dem,num integer :: m integer :: i,j,k,l @@ -64,8 +65,8 @@ subroutine GW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,eG dem = - OmBSE + eGW(k) - Om(m) + eGW(i) num = rho(j,k,m)*rho(i,l,m) - KD_dyn(ij,kl) = KD_dyn(ij,kl) - num*dem/(dem**2 + eta**2) - ZD_dyn(ij,kl) = ZD_dyn(ij,kl) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 + KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**2) + ZD_dyn(ij,kl) = ZD_dyn(ij,kl) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 dem = - OmBSE + eGW(l) - Om(m) + eGW(i) num = rho(i,k,m)*rho(j,l,m) @@ -76,8 +77,11 @@ subroutine GW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,eG dem = - OmBSE + eGW(l) - Om(m) + eGW(j) num = rho(j,k,m)*rho(i,l,m) - KD_dyn(ij,kl) = KD_dyn(ij,kl) - num*dem/(dem**2 + eta**2) - ZD_dyn(ij,kl) = ZD_dyn(ij,kl) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 + KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**2) + ZD_dyn(ij,kl) = ZD_dyn(ij,kl) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 + + KD_dyn(ij,kl) = 2d0*KD_dyn(ij,kl)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) + ZD_dyn(ij,kl) = 2d0*ZD_dyn(ij,kl)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) end do @@ -127,6 +131,9 @@ subroutine GW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,eG KD_dyn(ij,kl) = KD_dyn(ij,kl) - num*dem/(dem**2 + eta**2) ZD_dyn(ij,kl) = ZD_dyn(ij,kl) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 + KD_dyn(ij,kl) = 2d0*KD_dyn(ij,kl) + ZD_dyn(ij,kl) = 2d0*ZD_dyn(ij,kl) + end do end do diff --git a/src/GW/GW_ppBSE_dynamic_perturbation.f90 b/src/GW/GW_ppBSE_dynamic_perturbation.f90 index 39afc87..3d45425 100644 --- a/src/GW/GW_ppBSE_dynamic_perturbation.f90 +++ b/src/GW/GW_ppBSE_dynamic_perturbation.f90 @@ -96,7 +96,7 @@ subroutine GW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO, write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)' write(*,*) '---------------------------------------------------------------------------------------------------' - do ij=1,min(nOO,maxOO) + do ij=max(1,nOO+1-maxOO),nOO ! if(.not.dTDA) call GW_ppBSE_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGW,OmRPA,rho_RPA,OmBSE(ab),KB_dyn) call GW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,eGW,OmRPA,rho_RPA,Om2(ij),KD_dyn,ZD_dyn) diff --git a/src/GW/GW_ppBSE_static_kernel_B.f90 b/src/GW/GW_ppBSE_static_kernel_B.f90 index 7da42cb..6110654 100644 --- a/src/GW/GW_ppBSE_static_kernel_B.f90 +++ b/src/GW/GW_ppBSE_static_kernel_B.f90 @@ -56,10 +56,10 @@ subroutine GW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda do m=1,nS eps = Om(m)**2 + eta**2 chi = chi - rho(i,a,m)*rho(j,b,m)*Om(m)/eps & - + rho(i,b,m)*rho(a,j,m)*Om(m)/eps + - rho(i,b,m)*rho(a,j,m)*Om(m)/eps end do - KB(ab,ij) = 2d0*lambda*chi/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j))) + KB(ab,ij) = 4d0*lambda*chi/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j))) end do end do @@ -90,7 +90,7 @@ subroutine GW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda + rho(i,b,m)*rho(a,j,m)*Om(m)/eps end do - KB(ab,ij) = 2d0*lambda*chi + KB(ab,ij) = 4d0*lambda*chi end do end do @@ -121,7 +121,7 @@ subroutine GW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda + rho(i,b,m)*rho(a,j,m)*Om(m)/eps end do - KB(ab,ij) = lambda*chi + KB(ab,ij) = 2d0*lambda*chi end do end do diff --git a/src/GW/GW_ppBSE_static_kernel_C.f90 b/src/GW/GW_ppBSE_static_kernel_C.f90 index 82c6e58..ef21825 100644 --- a/src/GW/GW_ppBSE_static_kernel_C.f90 +++ b/src/GW/GW_ppBSE_static_kernel_C.f90 @@ -55,10 +55,10 @@ subroutine GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ERI do m=1,nS eps = Om(m)**2 + eta**2 chi = chi - rho(a,c,m)*rho(b,d,m)*Om(m)/eps & - + rho(a,d,m)*rho(b,c,m)*Om(m)/eps + - rho(a,d,m)*rho(b,c,m)*Om(m)/eps end do - KC(ab,cd) = 2d0*lambda*chi/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) + KC(ab,cd) = 4d0*lambda*chi/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) end do end do @@ -89,7 +89,7 @@ subroutine GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ERI + rho(a,d,m)*rho(b,c,m)*Om(m)/eps end do - KC(ab,cd) = 2d0*lambda*chi + KC(ab,cd) = 4d0*lambda*chi end do end do @@ -120,7 +120,7 @@ subroutine GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ERI + rho(a,d,m)*rho(b,c,m)*Om(m)/eps end do - KC(ab,cd) = lambda*chi + KC(ab,cd) = 2d0*lambda*chi end do end do diff --git a/src/GW/GW_ppBSE_static_kernel_D.f90 b/src/GW/GW_ppBSE_static_kernel_D.f90 index e0838c7..a54ff90 100644 --- a/src/GW/GW_ppBSE_static_kernel_D.f90 +++ b/src/GW/GW_ppBSE_static_kernel_D.f90 @@ -55,10 +55,10 @@ subroutine GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,ERI do m=1,nS eps = Om(m)**2 + eta**2 chi = chi - rho(i,k,m)*rho(j,l,m)*Om(m)/eps & - + rho(i,l,m)*rho(j,k,m)*Om(m)/eps + - rho(i,l,m)*rho(j,k,m)*Om(m)/eps end do - KD(ij,kl) = 2d0*lambda*chi/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) + KD(ij,kl) = 4d0*lambda*chi/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) end do end do @@ -89,7 +89,7 @@ subroutine GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,ERI + rho(i,l,m)*rho(j,k,m)*Om(m)/eps end do - KD(ij,kl) = 2d0*lambda*chi + KD(ij,kl) = 4d0*lambda*chi end do end do @@ -120,7 +120,7 @@ subroutine GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,ERI + rho(i,l,m)*rho(j,k,m)*Om(m)/eps end do - KD(ij,kl) = lambda*chi + KD(ij,kl) = 2d0*lambda*chi end do end do diff --git a/src/RPA/ppRRPA.f90 b/src/RPA/ppRRPA.f90 index fb5c3c4..65ddf24 100644 --- a/src/RPA/ppRRPA.f90 +++ b/src/RPA/ppRRPA.f90 @@ -76,7 +76,7 @@ subroutine ppRRPA(dotest,TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF, call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA(ispin)) -! call print_transition_vectors_pp(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) + call ppLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) call print_excitation_energies('ppRPA@RHF','2p (singlet)',nVV,Om1) call print_excitation_energies('ppRPA@RHF','2h (singlet)',nOO,Om2) @@ -108,7 +108,7 @@ subroutine ppRRPA(dotest,TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF, call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA(ispin)) -! call print_transition_vectors_pp(.false.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) + call ppLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) call print_excitation_energies('ppRPA@RHF','2p (triplet)',nVV,Om1) call print_excitation_energies('ppRPA@RHF','2h (triplet)',nOO,Om2)