diff --git a/input/methods b/input/methods index df12e6b..07cc57d 100644 --- a/input/methods +++ b/input/methods @@ -1,7 +1,7 @@ # RHF UHF KS MOM T F F F # MP2* MP3 MP2-F12 - T F F + F F F # CCD pCCD DCD CCSD CCSD(T) F F F F F # drCCD rCCD crCCD lCCD diff --git a/src/GT/G0T0.f90 b/src/GT/G0T0.f90 index 03db7ea..a472c1c 100644 --- a/src/GT/G0T0.f90 +++ b/src/GT/G0T0.f90 @@ -80,10 +80,10 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing ! Dimensions of the pp-RPA linear reponse matrices -! nOOs = nO*(nO + 1)/2 -! nVVs = nV*(nV + 1)/2 - nOOs = nO*nO - nVVs = nV*nV + nOOs = nO*(nO + 1)/2 + nVVs = nV*(nV + 1)/2 +! nOOs = nO*nO +! nVVs = nV*nV nOOt = nO*(nO - 1)/2 nVVt = nV*(nV - 1)/2 @@ -103,8 +103,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing !---------------------------------------------- ispin = 1 -! iblock = 1 - iblock = 3 + iblock = 1 +! iblock = 3 ! Compute linear response @@ -113,8 +113,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing ! EcRPA(ispin) = 1d0*EcRPA(ispin) -! call print_excitation('pp-RPA (N+2)',iblock,nVVs,Omega1s(:)) -! call print_excitation('pp-RPA (N-2)',iblock,nOOs,Omega2s(:)) + call print_excitation('pp-RPA (N+2)',iblock,nVVs,Omega1s(:)) + call print_excitation('pp-RPA (N-2)',iblock,nOOs,Omega2s(:)) !---------------------------------------------- ! alpha-alpha block @@ -131,8 +131,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing ! EcRPA(ispin) = 2d0*EcRPA(ispin) ! EcRPA(ispin) = 3d0*EcRPA(ispin) -! call print_excitation('pp-RPA (N+2)',iblock,nVVt,Omega1t(:)) -! call print_excitation('pp-RPA (N-2)',iblock,nOOt,Omega2t(:)) + call print_excitation('pp-RPA (N+2)',iblock,nVVt,Omega1t(:)) + call print_excitation('pp-RPA (N-2)',iblock,nOOt,Omega2t(:)) !---------------------------------------------- ! Compute T-matrix version of the self-energy @@ -142,12 +142,12 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing SigT(:) = 0d0 Z(:) = 0d0 -! iblock = 1 - iblock = 3 + iblock = 1 +! iblock = 3 call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO,X1s,Y1s,rho1s,X2s,Y2s,rho2s) - call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,Omega1s,rho1s,Omega2s,rho2s,EcGM,SigT) + call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,sqrt(1d0)*Omega1s,rho1s,sqrt(1d0)*Omega2s,rho2s,EcGM,SigT) call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,Omega1s,rho1s,Omega2s,rho2s,Z) @@ -155,7 +155,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO,X1t,Y1t,rho1t,X2t,Y2t,rho2t) - call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF,Omega1t,rho1t,Omega2t,rho2t,EcGM,SigT) + call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF,Omega1t,sqrt(1.5d0)*rho1t,Omega2t,sqrt(1.5d0)*rho2t,EcGM,SigT) call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF,Omega1t,rho1t,Omega2t,rho2t,Z) @@ -189,8 +189,9 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing ! Compute the ppRPA correlation energy ispin = 1 - iblock = 3 -! iblock = 1 +! iblock = 3 + iblock = 1 + call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eG0T0,ERI_MO, & Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin)) ispin = 2 diff --git a/src/GT/excitation_density_Tmatrix.f90 b/src/GT/excitation_density_Tmatrix.f90 index 33c0ff8..2a366b8 100644 --- a/src/GT/excitation_density_Tmatrix.f90 +++ b/src/GT/excitation_density_Tmatrix.f90 @@ -22,8 +22,8 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r ! Local variables - integer :: k,l - integer :: c,d + integer :: i,j,k,l + integer :: a,b,c,d integer :: p,q integer :: ab,cd,ij,kl double precision,external :: Kronecker_delta @@ -47,14 +47,20 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r do p=nC+1,nBas-nR do q=nC+1,nBas-nR - do ab=1,nVV +! do ab=1,nVV + ab = 0 + do a=nO+1,nBas-nR + do b=a,nBas-nR + ab = ab + 1 cd = 0 do c=nO+1,nBas-nR do d=c,nBas-nR cd = cd + 1 rho1(p,q,ab) = rho1(p,q,ab) & - + (ERI(p,q,c,d) + ERI(p,q,d,c))*X1(cd,ab)/sqrt((1d0 + Kronecker_delta(p,q))*(1d0 + Kronecker_delta(c,d))) + + ERI(p,q,c,d)*X1(cd,ab) +! + ERI(p,q,c,d)*X1(cd,ab)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) +! + (ERI(p,q,c,d) + ERI(p,q,d,c))*X1(cd,ab)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) end do end do @@ -63,20 +69,29 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r do l=k,nO kl = kl + 1 rho1(p,q,ab) = rho1(p,q,ab) & - + (ERI(p,q,k,l) + ERI(p,q,l,k))*Y1(kl,ab)/sqrt((1d0 + Kronecker_delta(p,q))*(1d0 + Kronecker_delta(k,l))) + + ERI(p,q,k,l)*Y1(kl,ab) +! + ERI(p,q,k,l)*Y1(kl,ab)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(k,l))) +! + (ERI(p,q,k,l) + ERI(p,q,l,k))*Y1(kl,ab)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(k,l))) end do end do end do + end do - do ij=1,nOO +! do ij=1,nOO + ij = 0 + do i=nC+1,nO + do j=i,nO + ij = ij + 1 cd = 0 do c=nO+1,nBas-nR do d=c,nBas-nR cd = cd + 1 rho2(p,q,ij) = rho2(p,q,ij) & - + (ERI(p,q,c,d) + ERI(p,q,d,c))*X2(cd,ij)/sqrt((1d0 + Kronecker_delta(p,q))*(1d0 + Kronecker_delta(c,d))) + + ERI(p,q,c,d)*X2(cd,ij) +! + ERI(p,q,c,d)*X2(cd,ij)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(c,d))) +! + (ERI(p,q,c,d) + ERI(p,q,d,c))*X2(cd,ij)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(c,d))) end do end do @@ -85,11 +100,14 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r do l=k,nO kl = kl + 1 rho2(p,q,ij) = rho2(p,q,ij) & - + (ERI(p,q,k,l) + ERI(p,q,l,k))*Y2(kl,ij)/sqrt((1d0 + Kronecker_delta(p,q))*(1d0 + Kronecker_delta(k,l))) + + ERI(p,q,k,l)*Y2(kl,ij) +! + ERI(p,q,k,l)*Y2(kl,ij)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) +! + (ERI(p,q,k,l) + ERI(p,q,l,k))*Y2(kl,ij)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) end do end do end do + end do end do end do diff --git a/src/eDFT/print_UKS.f90 b/src/eDFT/print_UKS.f90 index cb56dda..0d41965 100644 --- a/src/eDFT/print_UKS.f90 +++ b/src/eDFT/print_UKS.f90 @@ -66,7 +66,7 @@ subroutine print_UKS(nBas,nEns,occnum,Ov,wEns,eps,c,ENuc,ET,EV,EH,Ex,Ec,Ew,dipol HOMOb = -huge(0d0) if(iHOMOb > 0) HOMOb = eps(iHOMOb,2) LUMOb = +huge(0d0) - if(iLUMOb <= nBas) LUMOb = eps(iLUMOb,1) + if(iLUMOb <= nBas) LUMOb = eps(iLUMOb,2) HOMO = max(HOMOa,HOMOb) LUMO = min(LUMOa,LUMOb) diff --git a/src/eDFT/print_individual_energy.f90 b/src/eDFT/print_individual_energy.f90 index 483da4e..88ff78b 100644 --- a/src/eDFT/print_individual_energy.f90 +++ b/src/eDFT/print_individual_energy.f90 @@ -149,12 +149,12 @@ subroutine print_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux,LZH,LZx,LZc, write(*,'(A40,F16.10,A3)') ' H Levy-Zahariev shift: ',sum(LZH(:)),' au' write(*,'(A40,F16.10,A3)') ' x Levy-Zahariev shift: ',sum(LZx(:)),' au' write(*,'(A40,F16.10,A3)') ' c Levy-Zahariev shift: ',sum(LZc(:)),' au' - write(*,'(A40,F16.10,A3)') ' Hxc Levy-Zahariev shift: ',sum(LZH(:))+sum(LZx(:))+sum(LZx(:)),' au' + write(*,'(A40,F16.10,A3)') ' Hxc Levy-Zahariev shift: ',sum(LZH(:))+sum(LZx(:))+sum(LZc(:)),' au' write(*,*) write(*,'(A40,F16.10,A3)') ' H Levy-Zahariev shift: ',sum(LZH(:))*HaToeV,' eV' write(*,'(A40,F16.10,A3)') ' x Levy-Zahariev shift: ',sum(LZx(:))*HaToeV,' eV' write(*,'(A40,F16.10,A3)') ' c Levy-Zahariev shift: ',sum(LZc(:))*HaToeV,' eV' - write(*,'(A40,F16.10,A3)') ' Hxc Levy-Zahariev shift: ',(sum(LZH(:))+sum(LZx(:))+sum(LZx(:)))*HaToeV,' eV' + write(*,'(A40,F16.10,A3)') ' Hxc Levy-Zahariev shift: ',(sum(LZH(:))+sum(LZx(:))+sum(LZc(:)))*HaToeV,' eV' write(*,'(A60)') '-------------------------------------------------' write(*,*)