From 8469a6b91c696c6f6b786e40409fb5d0b90d339f Mon Sep 17 00:00:00 2001 From: pfloos Date: Tue, 28 Nov 2023 19:15:39 +0100 Subject: [PATCH] ufG0T0pp --- src/GT/ufG0T0pp.f90 | 537 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 537 insertions(+) create mode 100644 src/GT/ufG0T0pp.f90 diff --git a/src/GT/ufG0T0pp.f90 b/src/GT/ufG0T0pp.f90 new file mode 100644 index 0000000..1936d58 --- /dev/null +++ b/src/GT/ufG0T0pp.f90 @@ -0,0 +1,537 @@ +subroutine ufG0T0pp(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) + +! Upfolded G0T0pp equations + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: dotest + + logical,intent(in) :: TDA_W + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + double precision,intent(in) :: ENuc + double precision,intent(in) :: ERHF + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: eHF(nBas) + +! Local variables + + integer :: p + integer :: s + integer :: i,j,k,l + integer :: a,b,c,d + integer :: jb,kc,ia,ja + integer :: klc,kcd,ija,ijb,iab,jab + + logical :: dRPA + integer :: ispin + double precision :: EcRPA + integer :: n2h1p,n2p1h,nH + double precision,external :: Kronecker_delta + double precision,allocatable :: H(:,:) + double precision,allocatable :: eGW(:) + double precision,allocatable :: Z(:) + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) + double precision,allocatable :: Om(:) + double precision,allocatable :: XpY(:,:) + double precision,allocatable :: XmY(:,:) + double precision,allocatable :: rho(:,:,:) + + logical :: verbose = .true. + double precision,parameter :: cutoff1 = 0.01d0 + double precision,parameter :: cutoff2 = 0.01d0 + double precision :: eF + double precision,parameter :: window = 2d0 + + double precision :: start_timing,end_timing,timing + +! Output variables + +! Hello world + + write(*,*) + write(*,*)'******************************************' + write(*,*)'* Restricted Upfolded G0T0pp Calculation *' + write(*,*)'******************************************' + write(*,*) + +! Dimension of the supermatrix + + n2h1p = nO*nO*nV + n2p1h = nV*nV*nO + nH = 1 + n2h1p + n2p1h + +! Memory allocation + + allocate(H(nH,nH),eGW(nH),Z(nH)) + +! Initialization + + dRPA = .true. + EcRPA = 0d0 + + eF = 0.5d0*(eHF(nO+1) + eHF(nO)) + +!-------------------------! +! Main loop over orbitals ! +!-------------------------! + + do p=nO-1,nO + + H(:,:) = 0d0 + + if (TDA_W) then + + ! TDA for W + + write(*,*) 'Tamm-Dancoff approximation actived!' + write(*,*) + + !---------------------------! + ! Compute GW supermatrix ! + !---------------------------! + ! ! + ! | F V2h1p V2p1h | ! + ! | | ! + ! H = | V2h1p C2h1p 0 | ! + ! | | ! + ! | V2p1h 0 C2p1h | ! + ! ! + !---------------------------! + + call wall_time(start_timing) + + !---------! + ! Block F ! + !---------! + +! H(1,1) = eHF(p) +! +! !-------------! +! ! Block V2h1p ! +! !-------------! + +! ija = 0 +! do i=nC+1,nO +! do j=nC+1,nO +! do a=nO+1,nBas-nR +! ija = ija + 1 +! +! H(1 ,1+ija) = sqrt(2d0)*ERI(p,a,i,j) +! H(1+ija,1 ) = sqrt(2d0)*ERI(p,a,i,j) +! +! end do +! end do +! end do + +! !-------------! +! ! Block V2p1h ! +! !-------------! +! +! iab = 0 +! do i=nC+1,nO +! do a=nO+1,nBas-nR +! do b=nO+1,nBas-nR +! iab = iab + 1 +! +! H(1 ,1+n2h1p+iab) = sqrt(2d0)*ERI(p,i,b,a) +! H(1+n2h1p+iab,1 ) = sqrt(2d0)*ERI(p,i,b,a) +! +! end do +! end do +! end do +! +! !-------------! +! ! Block C2h1p ! +! !-------------! +! +! ija = 0 +! do i=nC+1,nO +! do j=nC+1,nO +! do a=nO+1,nBas-nR +! ija = ija + 1 +! +! klc = 0 +! do k=nC+1,nO +! do l=nC+1,nO +! do c=nO+1,nBas-nR +! klc = klc + 1 +! +! H(1+ija,1+klc) & +! = ((eHF(i) + eHF(j) - eHF(a))*Kronecker_delta(j,l)*Kronecker_delta(a,c) & +! - 2d0*ERI(j,c,a,l))*Kronecker_delta(i,k) +! +! end do +! end do +! end do +! +! end do +! end do +! end do +! +! !-------------! +! ! Block C2p1h ! +! !-------------! +! +! iab = 0 +! do i=nC+1,nO +! do a=nO+1,nBas-nR +! do b=nO+1,nBas-nR +! iab = iab + 1 +! +! kcd = 0 +! do k=nC+1,nO +! do c=nO+1,nBas-nR +! do d=nO+1,nBas-nR +! kcd = kcd + 1 +! +! H(1+n2h1p+iab,1+n2h1p+kcd) & +! = ((eHF(a) + eHF(b) - eHF(i))*Kronecker_delta(i,k)*Kronecker_delta(a,c) & +! + 2d0*ERI(a,k,i,c))*Kronecker_delta(b,d) +! +! end do +! end do +! end do +! +! end do +! end do +! end do +! + call wall_time(end_timing) + + timing = end_timing - start_timing + write(*,*) + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for construction of supermatrix = ',timing,' seconds' + write(*,*) + + else + + ! RPA for W + + write(*,*) 'Tamm-Dancoff approximation deactivated!' + write(*,*) + + !---------------------------! + ! Compute GW supermatrix ! + !---------------------------! + ! ! + ! | F W2h1p W2p1h | ! + ! | | ! + ! H = | W2h1p D2h1p 0 | ! + ! | | ! + ! | W2p1h 0 D2p1h | ! + ! ! + !---------------------------! + + ! Memory allocation + + allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & + Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), & + rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs), & + Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & + Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & + rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt)) + + ! alpha-beta block + + ispin = 1 + iblock = 3 + + ! Compute linear response + + allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) + + if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp) + call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eHF,ERI,Cpp) + call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eHF,ERI,Dpp) + + call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin)) + + deallocate(Bpp,Cpp,Dpp) + + call print_excitation_energies('ppRPA@RHF','2p (alpha-beta)',nVVs,Om1s(:)) + call print_excitation_energies('ppRPA@RHF','2h (alpha-beta)',nOOs,Om2s(:)) + + ! alpha-alpha block + + ispin = 2 + iblock = 4 + + ! Compute linear response + + allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) + + if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp) + call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eHF,ERI,Cpp) + call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eHF,ERI,Dpp) + + call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin)) + + deallocate(Bpp,Cpp,Dpp) + + call print_excitation_energies('ppRPA@RHF','2p (alpha-alpha)',nVVt,Om1t) + call print_excitation_energies('ppRPA@RHF','2h (alpha-beta)',nOOt,Om2t) + + !---------------------------------------------- + ! Compute excitation densities + !---------------------------------------------- + + iblock = 3 + call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) + + iblock = 4 + call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) + + call wall_time(start_timing) + + !---------! + ! Block F ! + !---------! + + H(1,1) = eHF(p) + + !-------------! + ! Block D2h1p ! + !-------------! + + ija = 0 + do i=nC+1,nO + do ja=1,nS + ija = ija + 1 + + H(1+ija,1+ija) = eHF(i) - Om(ja) + + end do + end do + + !-------------! + ! Block W2h1p ! + !-------------! + + ija = 0 + do i=nC+1,nO + do ja=1,nS + ija = ija + 1 + + H(1 ,1+ija) = sqrt(2d0)*rho(p,i,ja) + H(1+ija,1 ) = sqrt(2d0)*rho(p,i,ja) + + end do + end do + + !-------------! + ! Block D2p1h ! + !-------------! + + iab = 0 + do ia=1,nS + do b=nO+1,nBas-nR + iab = iab + 1 + + H(1+n2h1p+iab,1+n2h1p+iab) = eHF(b) + Om(ia) + + end do + end do + + !-------------! + ! Block W2p1h ! + !-------------! + + iab = 0 + do ia=1,nS + do b=nO+1,nBas-nR + iab = iab + 1 + + H(1 ,1+n2h1p+iab) = sqrt(2d0)*rho(p,b,ia) + H(1+n2h1p+iab,1 ) = sqrt(2d0)*rho(p,b,ia) + + end do + end do + + ! Memory deallocation + + deallocate(Om,Aph,Bph,XpY,XmY,rho) + + call wall_time(end_timing) + + timing = end_timing - start_timing + write(*,*) + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for construction of supermatrix = ',timing,' seconds' + write(*,*) + + end if + + !-------------------------! + ! Diagonalize supermatrix ! + !-------------------------! + + call wall_time(start_timing) + + call diagonalize_matrix(nH,H,eGW) + + call wall_time(end_timing) + + timing = end_timing - start_timing + write(*,*) + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for diagonalization of supermatrix = ',timing,' seconds' + write(*,*) + + !-----------------! + ! Compute weights ! + !-----------------! + + do s=1,nH + Z(s) = H(1,s)**2 + end do + + !--------------! + ! Dump results ! + !--------------! + + write(*,*)'-------------------------------------------' + write(*,'(1X,A32,I3,A8)')'| G0W0 energies (eV) for orbital',p,' |' + write(*,*)'-------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X)') & + '|','#','|','e_QP','|','Z','|' + write(*,*)'-------------------------------------------' + + do s=1,nH + if(eGW(s) < eF .and. eGW(s) > eF - window) then +! if(Z(s) > cutoff1) then + write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & + '|',s,'|',eGW(s)*HaToeV,'|',Z(s),'|' + end if + end do + + write(*,*)'-------------------------------------------' + write(*,*) + + if(verbose) then + + if(TDA_W) then + + ! TDA printing format + + do s=1,nH + + if(eGW(s) < eF .and. eGW(s) > eF - window) then + + write(*,*)'-------------------------------------------------------------' + write(*,'(1X,A7,1X,I3,A6,I3,A1,1X,A7,F12.6,A13,F6.4,1X)') & + 'Orbital',p,' and #',s,':','e_QP = ',eGW(s)*HaToeV,' eV and Z = ',Z(s) + write(*,*)'-------------------------------------------------------------' + write(*,'(1X,A20,1X,A20,1X,A15,1X)') & + ' Configuration ',' Coefficient ',' Weight ' + write(*,*)'-------------------------------------------------------------' + + if(p <= nO) & + write(*,'(1X,A7,I3,A16,1X,F15.6,1X,F15.6)') & + ' (',p,') ',H(1,s),H(1,s)**2 + if(p > nO) & + write(*,'(1X,A16,I3,A7,1X,F15.6,1X,F15.6)') & + ' (',p,') ',H(1,s),H(1,s)**2 + + ija = 0 + do i=nC+1,nO + do j=nC+1,nO + do a=nO+1,nBas-nR + ija = ija + 1 + + if(abs(H(1+ija,s)) > cutoff2) & + write(*,'(1X,A3,I3,A1,I3,A6,I3,A7,1X,F15.6,1X,F15.6)') & + ' (',i,',',j,') -> (',a,') ',H(1+ija,s),H(1+ija,s)**2 + + end do + end do + end do + + iab = 0 + do i=nC+1,nO + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + iab = iab + 1 + + if(abs(H(1+n2h1p+iab,s)) > cutoff2) & + write(*,'(1X,A7,I3,A6,I3,A1,I3,A3,1X,F15.6,1X,F15.6)') & + ' (',i,') -> (',a,',',b,') ',H(1+n2h1p+iab,s),H(1+n2h1p+iab,s)**2 + + end do + end do + end do + + write(*,*)'-------------------------------------------------------------' + write(*,*) + + end if + + end do + + else + + ! non-TDA printing format + + do s=1,nH + + if(eGW(s) < eF .and. eGW(s) > eF - window) then + + write(*,*)'-------------------------------------------------------------' + write(*,'(1X,A7,1X,I3,A6,I3,A1,1X,A7,F12.6,A13,F6.4,1X)') & + 'Orbital',p,' and #',s,':','e_QP = ',eGW(s)*HaToeV,' eV and Z = ',Z(s) + write(*,*)'-------------------------------------------------------------' + write(*,'(1X,A20,1X,A20,1X,A15,1X)') & + ' Conf. (p,ia) ',' Coefficient ',' Weight ' + write(*,*)'-------------------------------------------------------------' + + if(p <= nO) & + write(*,'(1X,A7,I3,A16,1X,F15.6,1X,F15.6)') & + ' (',p,') ',H(1,s),H(1,s)**2 + if(p > nO) & + write(*,'(1X,A16,I3,A7,1X,F15.6,1X,F15.6)') & + ' (',p,') ',H(1,s),H(1,s)**2 + + ija = 0 + do i=nC+1,nO + do ja=1,nS + ija = ija + 1 + + if(abs(H(1+ija,s)) > cutoff2) & + write(*,'(1X,A7,I3,A1,I3,A12,1X,F15.6,1X,F15.6)') & + ' (',i,',',ja,') ',H(1+ija,s),H(1+ija,s)**2 + + end do + end do + + iab = 0 + do ia=1,nS + do b=nO+1,nBas-nR + iab = iab + 1 + + if(abs(H(1+n2h1p+iab,s)) > cutoff2) & + write(*,'(1X,A7,I3,A1,I3,A12,1X,F15.6,1X,F15.6)') & + ' (',ia,',',b,') ',H(1+n2h1p+iab,s),H(1+n2h1p+iab,s)**2 + + end do + end do + + write(*,*)'-------------------------------------------------------------' + write(*,*) + + end if + + end do + + end if + + end if + + end do + +end subroutine