From 4c9841c5cd342d08a73b9212464bad73243cc2e0 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Fri, 27 Oct 2023 13:56:53 +0200 Subject: [PATCH] evGGF2 and evGGW --- input/methods | 2 +- input/options | 2 +- src/GF/G0F2.f90 | 4 +- src/GF/GG0F2.f90 | 32 ++- src/GF/GGF.f90 | 5 +- src/GF/GGF2_phBSE2.f90 | 71 +++++++ src/GF/GGF2_phBSE2_static_kernel_A.f90 | 162 +++++++++++++++ src/GF/GGF2_phBSE2_static_kernel_B.f90 | 161 +++++++++++++++ src/GF/evGF2.f90 | 4 +- src/GF/evGGF2.f90 | 173 ++++++++++++++++ src/GW/GG0W0.f90 | 4 +- src/GW/GGW.f90 | 4 +- src/GW/evGGW.f90 | 267 +++++++++++++++++++++++++ src/HF/print_GHF.f90 | 2 - 14 files changed, 860 insertions(+), 33 deletions(-) create mode 100644 src/GF/GGF2_phBSE2.f90 create mode 100644 src/GF/GGF2_phBSE2_static_kernel_A.f90 create mode 100644 src/GF/GGF2_phBSE2_static_kernel_B.f90 create mode 100644 src/GF/evGGF2.f90 create mode 100644 src/GW/evGGW.f90 diff --git a/input/methods b/input/methods index 7f0c8cf..f0cc512 100644 --- a/input/methods +++ b/input/methods @@ -11,7 +11,7 @@ # phRPA phRPAx crRPA ppRPA F F F F # G0F2 evGF2 qsGF2 G0F3 evGF3 - F F F F F + F T F F F # G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW F F F F F F # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh diff --git a/input/options b/input/options index a899bae..0f3571d 100644 --- a/input/options +++ b/input/options @@ -15,4 +15,4 @@ # ACFDT: AC Kx XBS F F T # BSE: phBSE phBSE2 ppBSE dBSE dTDA - T F F F T + F F F F T diff --git a/src/GF/G0F2.f90 b/src/GF/G0F2.f90 index bfc0c6f..317fc86 100644 --- a/src/GF/G0F2.f90 +++ b/src/GF/G0F2.f90 @@ -82,7 +82,7 @@ subroutine G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regu ! Print results - call MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,eGF,Ec) + call MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec) call print_G0F2(nBas,nO,eHF,SigC,eGF,Z,ENuc,ERHF,Ec) ! Perform BSE2 calculation @@ -96,7 +96,7 @@ subroutine G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regu write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy (singlet) =',EcBSE(1) write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy (triplet) =',EcBSE(2) write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy =',sum(EcBSE(:)) - write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 total energy =',ENuc + EHF + sum(EcBSE(:)) + write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 total energy =',ENuc + ERHF + sum(EcBSE(:)) write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/GF/GG0F2.f90 b/src/GF/GG0F2.f90 index eceb702..4b82ef5 100644 --- a/src/GF/GG0F2.f90 +++ b/src/GF/GG0F2.f90 @@ -31,7 +31,7 @@ subroutine GG0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize, & ! Local variables double precision :: Ec - double precision :: EcBSE(nspin) + double precision :: EcBSE double precision,allocatable :: eGFlin(:) double precision,allocatable :: eGF(:) double precision,allocatable :: SigC(:) @@ -85,33 +85,29 @@ subroutine GG0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize, & ! Perform BSE2 calculation -! if(dophBSE) then -! -! call GF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE) + if(dophBSE) then + + call GGF2_phBSE2(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE) -! write(*,*) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy (singlet) =',EcBSE(1) -! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy (triplet) =',EcBSE(2) -! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy =',sum(EcBSE(:)) -! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 total energy =',ENuc + EHF + sum(EcBSE(:)) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,*) + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@GG0F2 correlation energy =',EcBSE,' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@GG0F2 total energy =',ENuc + EHF + EcBSE,' au' + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) -! end if + end if ! Perform ppBSE2 calculation ! if(doppBSE) then ! -! call GF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE) +! call GGF2_ppBSE2(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE) ! write(*,*) ! write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy (singlet) =',EcBSE(1),' au' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy (triplet) =',3d0*EcBSE(2),' au' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy =',EcBSE(1) + 3d0*EcBSE(2),' au' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 total energy =',ENuc + ERHF + EcBSE(1) + 3d0*EcBSE(2),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@GG0F2 correlation energy =',EcBSE,' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@GG0F2 total energy =',ENuc + ERHF + EcBSE,' au' ! write(*,*)'-------------------------------------------------------------------------------' ! write(*,*) diff --git a/src/GF/GGF.f90 b/src/GF/GGF.f90 index 9c63e33..8570bfa 100644 --- a/src/GF/GGF.f90 +++ b/src/GF/GGF.f90 @@ -78,9 +78,8 @@ subroutine GGF(doG0F2,doevGF2,doqsGF2,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA if(doevGF2) then call wall_time(start_GF) -! call evGGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, & -! linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF, & -! ERI,dipole_int,epsHF) + call evGGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,linearize,eta,regularize, & + nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF) call wall_time(end_GF) t_GF = end_GF - start_GF diff --git a/src/GF/GGF2_phBSE2.f90 b/src/GF/GGF2_phBSE2.f90 new file mode 100644 index 0000000..471d723 --- /dev/null +++ b/src/GF/GGF2_phBSE2.f90 @@ -0,0 +1,71 @@ +subroutine GGF2_phBSE2(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE) + +! Compute the second-order Bethe-Salpeter excitation energies + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + + double precision,intent(in) :: eta + 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) :: eGF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + +! Local variables + + logical :: dRPA = .false. + integer :: ispin + double precision,allocatable :: OmBSE(:) + double precision,allocatable :: XpY(:,:) + double precision,allocatable :: XmY(:,:) + double precision,allocatable :: A_sta(:,:) + double precision,allocatable :: B_sta(:,:) + double precision,allocatable :: KA_sta(:,:) + double precision,allocatable :: KB_sta(:,:) + +! Output variables + + double precision,intent(out) :: EcBSE + +! Memory allocation + + allocate(OmBSE(nS),XpY(nS,nS),XmY(nS,nS),A_sta(nS,nS),KA_sta(nS,nS)) + allocate(B_sta(nS,nS),KB_sta(nS,nS)) + + ispin = 3 + EcBSE = 0d0 + + call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI,A_sta) + if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,B_sta) + + ! Compute static kernel + + call GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KA_sta) + if(.not.TDA) call GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KB_sta) + + A_sta(:,:) = A_sta(:,:) + KA_sta(:,:) + if(.not.TDA) B_sta(:,:) = B_sta(:,:) + KB_sta(:,:) + + ! Compute phBSE2@GF2 excitation energies + + call phLR(TDA,nS,A_sta,B_sta,EcBSE,OmBSE,XpY,XmY) + call print_excitation_energies('phBSE2@GGF2',ispin,nS,OmBSE) + call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY) + + ! Compute dynamic correction for BSE via perturbation theory + +! if(dBSE) & +! call GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY) + +end subroutine diff --git a/src/GF/GGF2_phBSE2_static_kernel_A.f90 b/src/GF/GGF2_phBSE2_static_kernel_A.f90 new file mode 100644 index 0000000..9079f94 --- /dev/null +++ b/src/GF/GGF2_phBSE2_static_kernel_A.f90 @@ -0,0 +1,162 @@ +subroutine GGF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,KA_sta) + +! Compute the resonant part of the static BSE2 matrix + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: nBas,nC,nO,nV,nR,nS + double precision,intent(in) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: eGF(nBas) + +! Local variables + + double precision :: dem,num + integer :: i,j,k,l + integer :: a,b,c,d + integer :: ia,jb + +! Output variables + + double precision,intent(out) :: KA_sta(nS,nS) + +! Initialization + + KA_sta(:,:) = 0d0 + +! Second-order correlation kernel for the block A of the singlet manifold + + if(ispin == 1) then + + jb = 0 +!$omp parallel do default(private) shared(KA_sta,ERI,num,dem,eGF,nO,nBas,eta,nC,nR) +! do j=nC+1,nO +! do b=nO+1,nBas-nR +! jb = (b-nO) + (j-1)*(nBas-nO) +! +! ia = 0 +! do i=nC+1,nO +! do a=nO+1,nBas-nR +! ia = (a-nO) + (i-1)*(nBas-nO) +! +! +! do k=nC+1,nO +! do c=nO+1,nBas-nR +! +! dem = - (eGF(c) - eGF(k)) +! num = 2d0*ERI(j,k,i,c)*ERI(a,c,b,k) - ERI(j,k,i,c)*ERI(a,c,k,b) & +! - ERI(j,k,c,i)*ERI(a,c,b,k) + 2d0*ERI(j,k,c,i)*ERI(a,c,k,b) + +! KA_sta(ia,jb) = KA_sta(ia,jb) - num*dem/(dem**2 + eta**2) +! +! dem = + (eGF(c) - eGF(k)) +! num = 2d0*ERI(j,c,i,k)*ERI(a,k,b,c) - ERI(j,c,i,k)*ERI(a,k,c,b) & +! - ERI(j,c,k,i)*ERI(a,k,b,c) + 2d0*ERI(j,c,k,i)*ERI(a,k,c,b) + +! KA_sta(ia,jb) = KA_sta(ia,jb) + num*dem/(dem**2 + eta**2) +! +! end do +! end do + +! do c=nO+1,nBas-nR +! do d=nO+1,nBas-nR +! +! dem = - (eGF(c) + eGF(d)) +! num = 2d0*ERI(a,j,c,d)*ERI(c,d,i,b) - ERI(a,j,c,d)*ERI(c,d,b,i) & +! - ERI(a,j,d,c)*ERI(c,d,i,b) + 2d0*ERI(a,j,d,c)*ERI(c,d,b,i) + +! KA_sta(ia,jb) = KA_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) +! +! end do +! end do + +! do k=nC+1,nO +! do l=nC+1,nO + +! dem = - (eGF(k) + eGF(l)) +! num = 2d0*ERI(a,j,k,l)*ERI(k,l,i,b) - ERI(a,j,k,l)*ERI(k,l,b,i) & +! - ERI(a,j,l,k)*ERI(k,l,i,b) + 2d0*ERI(a,j,l,k)*ERI(k,l,b,i) + +! KA_sta(ia,jb) = KA_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) +! +! end do +! end do +! +! end do +! end do + +! end do +! end do +!$omp end parallel do + +! end if + +! Second-order correlation kernel for the block A of the triplet manifold + +! if(ispin == 2) then + +! jb = 0 +!$omp parallel do default(private) shared(KA_sta,ERI,num,dem,eGF,nO,nBas,eta,nC,nR) +! do j=nC+1,nO +! do b=nO+1,nBas-nR +! jb = (b-nO) + (j-1)*(nBas-nO) +! +! ia = 0 +! do i=nC+1,nO +! do a=nO+1,nBas-nR +! ia = (a-nO) + (i-1)*(nBas-nO) +! +! do k=nC+1,nO +! do c=nO+1,nBas-nR +! +! dem = - (eGF(c) - eGF(k)) +! num = 2d0*ERI(j,k,i,c)*ERI(a,c,b,k) - ERI(j,k,i,c)*ERI(a,c,k,b) - ERI(j,k,c,i)*ERI(a,c,b,k) + +! KA_sta(ia,jb) = KA_sta(ia,jb) - num*dem/(dem**2 + eta**2) +! +! dem = + (eGF(c) - eGF(k)) +! num = 2d0*ERI(j,c,i,k)*ERI(a,k,b,c) - ERI(j,c,i,k)*ERI(a,k,c,b) - ERI(j,c,k,i)*ERI(a,k,b,c) + +! KA_sta(ia,jb) = KA_sta(ia,jb) + num*dem/(dem**2 + eta**2) +! +! end do +! end do + +! do c=nO+1,nBas-nR +! do d=nO+1,nBas-nR +! +! dem = - (eGF(c) + eGF(d)) +! num = ERI(a,j,c,d)*ERI(c,d,b,i) + ERI(a,j,d,c)*ERI(c,d,i,b) + +! KA_sta(ia,jb) = KA_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) +! +! end do +! end do + +! do k=nC+1,nO +! do l=nC+1,nO + +! dem = - (eGF(k) + eGF(l)) +! num = ERI(a,j,k,l)*ERI(k,l,b,i) + ERI(a,j,l,k)*ERI(k,l,i,b) + +! KA_sta(ia,jb) = KA_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) +! +! end do +! end do +! +! end do +! end do + +! end do +! end do +!$omp end parallel do + + end if + + +end subroutine diff --git a/src/GF/GGF2_phBSE2_static_kernel_B.f90 b/src/GF/GGF2_phBSE2_static_kernel_B.f90 new file mode 100644 index 0000000..c789898 --- /dev/null +++ b/src/GF/GGF2_phBSE2_static_kernel_B.f90 @@ -0,0 +1,161 @@ +subroutine GGF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,KB_sta) + +! Compute the anti-resonant part of the static BSE2 matrix + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: nBas,nC,nO,nV,nR,nS + double precision,intent(in) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: eGF(nBas) + +! Local variables + + double precision :: dem,num + integer :: i,j,k,l + integer :: a,b,c,d + integer :: ia,jb + +! Output variables + + double precision,intent(out) :: KB_sta(nS,nS) + +! Initialization + + KB_sta(:,:) = 0d0 + +! Second-order correlation kernel for the block A of the singlet manifold + + if(ispin == 1) then + + jb = 0 +!$omp parallel do default(private) shared(KB_sta,ERI,num,dem,eGF,nO,nBas,eta,nC,nR) +! do j=nC+1,nO +! do b=nO+1,nBas-nR +! jb = (b-nO) + (j-1)*(nBas-nO) + +! ia = 0 +! do i=nC+1,nO +! do a=nO+1,nBas-nR +! ia = (a-nO) + (i-1)*(nBas-nO) +! +! do k=nC+1,nO +! do c=nO+1,nBas-nR +! +! dem = + eGF(k) - eGF(c) +! num = 2d0*ERI(b,k,i,c)*ERI(a,c,j,k) - ERI(b,k,i,c)*ERI(a,c,k,j) & +! - ERI(b,k,c,i)*ERI(a,c,j,k) + 2d0*ERI(b,k,c,i)*ERI(a,c,k,j) + +! KB_sta(ia,jb) = KB_sta(ia,jb) - num*dem/(dem**2 + eta**2) +! +! dem = - eGF(c) + eGF(k) +! num = 2d0*ERI(b,c,i,k)*ERI(a,k,j,c) - ERI(b,c,i,k)*ERI(a,k,c,j) & +! - ERI(b,c,k,i)*ERI(a,k,j,c) + 2d0*ERI(b,c,k,i)*ERI(a,k,c,j) + +! KB_sta(ia,jb) = KB_sta(ia,jb) - num*dem/(dem**2 + eta**2) +! +! end do +! end do + +! do c=nO+1,nBas-nR +! do d=nO+1,nBas-nR +! +! dem = - eGF(c) - eGF(d) +! num = 2d0*ERI(a,b,c,d)*ERI(c,d,i,j) - ERI(a,b,c,d)*ERI(c,d,j,i) & +! - ERI(a,b,d,c)*ERI(c,d,i,j) + 2d0*ERI(a,b,d,c)*ERI(c,d,j,i) + +! KB_sta(ia,jb) = KB_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) +! +! end do +! end do + +! do k=nC+1,nO +! do l=nC+1,nO + +! dem = + eGF(k) + eGF(l) +! num = 2d0*ERI(a,b,k,l)*ERI(k,l,i,j) - ERI(a,b,k,l)*ERI(k,l,j,i) & +! - ERI(a,b,l,k)*ERI(k,l,i,j) + 2d0*ERI(a,b,l,k)*ERI(k,l,j,i) + +! KB_sta(ia,jb) = KB_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) +! +! end do +! end do +! +! end do +! end do + +! end do +! end do +!$omp end parallel do + +! end if + +! Second-order correlation kernel for the block A of the triplet manifold + +! if(ispin == 2) then + +! jb = 0 +!$omp parallel do default(private) shared(KB_sta,ERI,num,dem,eGF,nO,nBas,eta,nC,nR) +! do j=nC+1,nO +! do b=nO+1,nBas-nR +! jb = (b-nO) + (j-1)*(nBas-nO) +! +! ia = 0 +! do i=nC+1,nO +! do a=nO+1,nBas-nR +! ia = (a-nO) + (i-1)*(nBas-nO) +! +! do k=nC+1,nO +! do c=nO+1,nBas-nR +! +! dem = + eGF(k) - eGF(c) +! num = 2d0*ERI(b,k,i,c)*ERI(a,c,j,k) - ERI(b,k,i,c)*ERI(a,c,k,j) - ERI(b,k,c,i)*ERI(a,c,j,k) + +! KB_sta(ia,jb) = KB_sta(ia,jb) - num*dem/(dem**2 + eta**2) +! +! dem = - eGF(c) + eGF(k) +! num = 2d0*ERI(b,c,i,k)*ERI(a,k,j,c) - ERI(b,c,i,k)*ERI(a,k,c,j) - ERI(b,c,k,i)*ERI(a,k,j,c) + +! KB_sta(ia,jb) = KB_sta(ia,jb) - num*dem/(dem**2 + eta**2) +! +! end do +! end do + +! do c=nO+1,nBas-nR +! do d=nO+1,nBas-nR +! +! dem = - eGF(c) - eGF(d) +! num = ERI(a,b,c,d)*ERI(c,d,j,i) + ERI(a,b,d,c)*ERI(c,d,i,j) + +! KB_sta(ia,jb) = KB_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) +! +! end do +! end do + +! do k=nC+1,nO +! do l=nC+1,nO + +! dem = + eGF(k) + eGF(l) +! num = ERI(a,b,k,l)*ERI(k,l,j,i) + ERI(a,b,l,k)*ERI(k,l,i,j) + +! KB_sta(ia,jb) = KB_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) +! +! end do +! end do +! +! end do +! end do + +! end do +! end do +!$omp end parallel do + + end if + + +end subroutine diff --git a/src/GF/evGF2.f90 b/src/GF/evGF2.f90 index 6984f13..89c4163 100644 --- a/src/GF/evGF2.f90 +++ b/src/GF/evGF2.f90 @@ -107,7 +107,7 @@ subroutine evGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,tr ! Print results - call MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,eGF,Ec) + call MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec) call print_evGF2(nBas,nO,nSCF,Conv,eHF,SigC,Z,eGF,ENuc,ERHF,Ec) ! DIIS extrapolation @@ -153,7 +153,7 @@ subroutine evGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,tr write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGF2 correlation energy (singlet) =',EcBSE(1) write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGF2 correlation energy (triplet) =',EcBSE(2) write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGF2 correlation energy =',sum(EcBSE(:)) - write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGF2 total energy =',ENuc + EHF + sum(EcBSE(:)) + write(*,'(2X,A50,F20.10)') 'Tr@phBSE@evGF2 total energy =',ENuc + ERHF + sum(EcBSE(:)) write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/GF/evGGF2.f90 b/src/GF/evGGF2.f90 new file mode 100644 index 0000000..086c7ff --- /dev/null +++ b/src/GF/evGGF2.f90 @@ -0,0 +1,173 @@ +subroutine evGGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, & + linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) + +! Perform eigenvalue self-consistent second-order Green function calculation + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: dophBSE + logical,intent(in) :: doppBSE + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + integer,intent(in) :: maxSCF + double precision,intent(in) :: thresh + integer,intent(in) :: max_diis + logical,intent(in) :: linearize + double precision,intent(in) :: eta + logical,intent(in) :: regularize + + integer,intent(in) :: nBas + integer,intent(in) :: nO + integer,intent(in) :: nC + 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) :: eHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + +! Local variables + + integer :: nSCF + integer :: n_diis + double precision :: Ec + double precision :: EcBSE + double precision :: Conv + double precision :: rcond + double precision,allocatable :: eGF(:) + double precision,allocatable :: eOld(:) + double precision,allocatable :: SigC(:) + double precision,allocatable :: Z(:) + double precision,allocatable :: error_diis(:,:) + double precision,allocatable :: e_diis(:,:) + +! Hello world + + write(*,*) + write(*,*)'************************************************' + write(*,*)'| Second-order Green function calculation |' + write(*,*)'************************************************' + write(*,*) + +! Memory allocation + + allocate(SigC(nBas),Z(nBas),eGF(nBas),eOld(nBas),error_diis(nBas,max_diis),e_diis(nBas,max_diis)) + +! Initialization + + Conv = 1d0 + nSCF = 0 + n_diis = 0 + e_diis(:,:) = 0d0 + error_diis(:,:) = 0d0 + eGF(:) = eHF(:) + eOld(:) = eHF(:) + rcond = 0d0 + +!------------------------------------------------------------------------ +! Main SCF loop +!------------------------------------------------------------------------ + + do while(Conv > thresh .and. nSCF < maxSCF) + + ! Frequency-dependent second-order contribution + + if(regularize) then + +! call GGF2_reg_self_energy_diag(eta,nBas,nC,nO,nV,nR,eGF,ERI,SigC,Z) + + else + + call GGF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,eGF,ERI,SigC,Z) + + end if + + ! Solve the quasi-particle equation + + if(linearize) then + + else + + write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** ' + write(*,*) + + call GGF2_QP_graph(eta,nBas,nC,nO,nV,nR,eHF,ERI,eOld,eOld,eGF,Z) + + end if + + Conv = maxval(abs(eGF - eOld)) + + ! Print results + + call GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec) + call print_evGF2(nBas,nO,nSCF,Conv,eHF,SigC,Z,eGF,ENuc,ERHF,Ec) + + ! DIIS extrapolation + + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGF-eOld,eGF) + + if(abs(rcond) < 1d-15) n_diis = 0 + + eOld(:) = eGF(:) + + ! Increment + + nSCF = nSCF + 1 + + end do +!------------------------------------------------------------------------ +! End main SCF loop +!------------------------------------------------------------------------ + +! Did it actually converge? + + if(nSCF == maxSCF+1) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + + stop + + end if + +! Perform BSE2 calculation + + if(dophBSE) then + + call GGF2_phBSE2(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE) + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@evGGF2 correlation energy =',EcBSE,' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@evGGF2 total energy =',ENuc + ERHF + EcBSE,' au' + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + + end if + +! Perform ppBSE2 calculation + +! if(doppBSE) then + +! call GGF2_ppBSE2(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE) + +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@evGGF2 correlation energy =',EcBSE),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@evGGF2 total energy =',ENuc + ERHF + EcBSE,' au' +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! end if + +end subroutine diff --git a/src/GW/GG0W0.f90 b/src/GW/GG0W0.f90 index 2d496f9..7498b3b 100644 --- a/src/GW/GG0W0.f90 +++ b/src/GW/GG0W0.f90 @@ -202,7 +202,7 @@ subroutine GG0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,d ! end if -! end if + end if ! if(doppBSE) then @@ -217,6 +217,6 @@ subroutine GG0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,d ! write(*,*)'-------------------------------------------------------------------------------' ! write(*,*) - end if +! end if end subroutine diff --git a/src/GW/GGW.f90 b/src/GW/GGW.f90 index 55cec7f..af1b4f0 100644 --- a/src/GW/GGW.f90 +++ b/src/GW/GGW.f90 @@ -88,8 +88,8 @@ subroutine GGW(doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT, & if(doevGW) then call wall_time(start_GW) -! call evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & -! linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF) + call evGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & + linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF) call wall_time(end_GW) t_GW = end_GW - start_GW diff --git a/src/GW/evGGW.f90 b/src/GW/evGGW.f90 new file mode 100644 index 0000000..923ea1a --- /dev/null +++ b/src/GW/evGGW.f90 @@ -0,0 +1,267 @@ +subroutine evGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & + linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) + +! Perform self-consistent eigenvalue-only GW calculation + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + double precision,intent(in) :: ENuc + double precision,intent(in) :: ERHF + logical,intent(in) :: doACFDT + logical,intent(in) :: exchange_kernel + logical,intent(in) :: doXBS + logical,intent(in) :: dophBSE + logical,intent(in) :: dophBSE2 + logical,intent(in) :: TDA_W + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + logical,intent(in) :: doppBSE + logical,intent(in) :: linearize + double precision,intent(in) :: eta + logical,intent(in) :: regularize + + 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) :: eHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + +! Local variables + + logical :: linear_mixing + logical :: dRPA = .true. + integer :: ispin + integer :: nSCF + integer :: n_diis + double precision :: rcond + double precision :: Conv + double precision :: EcRPA + double precision :: EcBSE + double precision :: EcGM + double precision :: alpha + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) + double precision,allocatable :: error_diis(:,:) + double precision,allocatable :: e_diis(:,:) + double precision,allocatable :: eGW(:) + double precision,allocatable :: eOld(:) + double precision,allocatable :: Z(:) + double precision,allocatable :: SigC(:) + double precision,allocatable :: Om(:) + double precision,allocatable :: XpY(:,:) + double precision,allocatable :: XmY(:,:) + double precision,allocatable :: rho(:,:,:) + +! Hello world + + write(*,*) + write(*,*)'************************************************' + write(*,*)'| Self-consistent evGW calculation |' + write(*,*)'************************************************' + write(*,*) + +! TDA for W + + if(TDA_W) then + write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' + write(*,*) + end if + +! TDA + + if(TDA) then + write(*,*) 'Tamm-Dancoff approximation activated!' + write(*,*) + end if + +! Linear mixing + + linear_mixing = .false. + alpha = 0.2d0 + +! Memory allocation + + allocate(Aph(nS,nS),Bph(nS,nS),eGW(nBas),eOld(nBas),Z(nBas),SigC(nBas), & + Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS),error_diis(nBas,max_diis),e_diis(nBas,max_diis)) + +! Initialization + + nSCF = 0 + ispin = 3 + n_diis = 0 + Conv = 1d0 + e_diis(:,:) = 0d0 + error_diis(:,:) = 0d0 + eGW(:) = eHF(:) + eOld(:) = eGW(:) + Z(:) = 1d0 + rcond = 0d0 + +!------------------------------------------------------------------------ +! Main loop +!------------------------------------------------------------------------ + + do while(Conv > thresh .and. nSCF <= maxSCF) + + ! Compute screening + + call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph) + if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) + + call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) + + ! Compute spectral weights + + call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho) + + ! Compute correlation part of the self-energy + + if(regularize) call GW_regularization(nBas,nC,nO,nV,nR,nS,eGW,Om,rho) + + call GGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) + + ! Solve the quasi-particle equation + + if(linearize) then + + write(*,*) ' *** Quasiparticle energies obtained by linearization *** ' + write(*,*) + + eGW(:) = eHF(:) + SigC(:) + + else + + write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** ' + write(*,*) + + call GGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eOld,eOld,eGW,Z) + + end if + + ! Convergence criteria + + Conv = maxval(abs(eGW - eOld)) + + ! Print results + + call print_evGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM) + + ! Linear mixing or DIIS extrapolation + + if(linear_mixing) then + + eGW(:) = alpha*eGW(:) + (1d0 - alpha)*eOld(:) + + else + + n_diis = min(n_diis+1,max_diis) + if(abs(rcond) > 1d-7) then + call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGW-eOld,eGW) + else + n_diis = 0 + end if + + end if + + ! Save quasiparticles energy for next cycle + + eOld(:) = eGW(:) + + ! Increment + + nSCF = nSCF + 1 + + end do +!------------------------------------------------------------------------ +! End main loop +!------------------------------------------------------------------------ + +! Did it actually converge? + + if(nSCF == maxSCF+1) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + + stop + + end if + +! Deallocate memory + + deallocate(Aph,Bph,eOld,Z,SigC,Om,XpY,XmY,rho,error_diis,e_diis) + +! Perform BSE calculation + + if(dophBSE) then + + call GGW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGW,eGW,EcBSE) + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW correlation energy =',EcBSE + write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW total energy =',ENuc + ERHF + EcBSE + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + +! Compute the BSE correlation energy via the adiabatic connection + +! if(doACFDT) then + +! write(*,*) '------------------------------------------------------' +! write(*,*) 'Adiabatic connection version of BSE correlation energy' +! write(*,*) '------------------------------------------------------' +! write(*,*) + +! if(doXBS) then + +! write(*,*) '*** scaled screening version (XBS) ***' +! write(*,*) + +! end if + +! call GW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcBSE) + +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy (singlet) =',EcBSE(1) +! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy (triplet) =',EcBSE(2) +! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy =',EcBSE(1) + EcBSE(2) +! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! end if + + end if + +! if(doppBSE) then + +! call GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE) + +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy (singlet) =',EcBSE(1) +! write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy (triplet) =',3d0*EcBSE(2) +! write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy =',EcBSE(1) + 3d0*EcBSE(2) +! write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW total energy =',ENuc + ERHF + EcBSE(1) + 3d0*EcBSE(2) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! end if + +end subroutine diff --git a/src/HF/print_GHF.f90 b/src/HF/print_GHF.f90 index 9490602..a5f844b 100644 --- a/src/HF/print_GHF.f90 +++ b/src/HF/print_GHF.f90 @@ -82,8 +82,6 @@ subroutine print_GHF(nBas,nBas2,nO,e,C,P,ENuc,ET,EV,EJ,EK,EHF,dipole) Sz2 = Sz2 + 0.25d0*(Pab(mu,nu)*Pba(nu,mu) - Pba(mu,nu)*Pab(nu,mu)) end do end do - - print*,Sx2,Sy2,Sz2 S2 = Sx2 + Sy2 + Sz2