From 96207d9c58eb43d0a63c6d2174fa796f95e4e48e Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 1 Jun 2020 17:14:42 +0200 Subject: [PATCH] problem minus sign solved --- input/methods | 6 +- src/QuAcK/Bethe_Salpeter.f90 | 2 +- .../Bethe_Salpeter_AB_matrix_dynamic.f90 | 16 ++-- .../Bethe_Salpeter_ZAB_matrix_dynamic.f90 | 20 ++-- src/QuAcK/G0F2.f90 | 92 ++++++++++++++++--- src/QuAcK/QuAcK.f90 | 6 +- src/QuAcK/evGF2.f90 | 11 ++- 7 files changed, 113 insertions(+), 40 deletions(-) diff --git a/input/methods b/input/methods index 63c03d4..cf57016 100644 --- a/input/methods +++ b/input/methods @@ -9,11 +9,11 @@ # CIS CID CISD F F F # RPA RPAx ppRPA - F F F + F T F # G0F2 evGF2 G0F3 evGF3 - F F F F + T F F F # G0W0 evGW qsGW - T F F + F F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/src/QuAcK/Bethe_Salpeter.f90 b/src/QuAcK/Bethe_Salpeter.f90 index 9d19585..2237c46 100644 --- a/src/QuAcK/Bethe_Salpeter.f90 +++ b/src/QuAcK/Bethe_Salpeter.f90 @@ -1,5 +1,5 @@ subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI,eW,eGW,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,EcRPA,EcBSE) + nBas,nC,nO,nV,nR,nS,ERI,eW,eGW,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,EcRPA,EcBSE) ! Compute the Bethe-Salpeter excitation energies diff --git a/src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90 b/src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90 index d370f92..20c1bb4 100644 --- a/src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90 +++ b/src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90 @@ -76,24 +76,24 @@ subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O do kc=1,maxS eps_Ap = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2 - eps_Am = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2 + eps_Am = (- OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2 chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*(+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))/eps_Ap - chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*(+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))/eps_Am + chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*(- OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))/eps_Am eps_Ap = (+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2 - eps_Am = (+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2 + eps_Am = (- OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2 chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*(+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))/eps_Ap - chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*(+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))/eps_Am + chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*(- OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))/eps_Am eps_Bp = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2 - eps_Bm = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2 + eps_Bm = (- OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2 chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*(+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))/eps_Bp - chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*(+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))/eps_Bm + chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*(- OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))/eps_Bm eps_Bp = (+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2 - eps_Bm = (+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2 + eps_Bm = (- OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2 chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*(+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps_Bp - chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*(+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps_Bm + chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*(- OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps_Bm enddo diff --git a/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 b/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 index da0db60..f89999d 100644 --- a/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 +++ b/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 @@ -59,32 +59,32 @@ subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW, do kc=1,maxS eps_Ap = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2 - eps_Am = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2 + eps_Am = (- OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2 chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*((+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))/eps_Ap)**2 - chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*((+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))/eps_Am)**2 + chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*((- OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))/eps_Am)**2 eps_Ap = (+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2 - eps_Am = (+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2 + eps_Am = (- OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2 chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*((+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))/eps_Ap)**2 - chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*((+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))/eps_Am)**2 + chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*((- OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))/eps_Am)**2 eps_Bp = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2 - eps_Bm = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2 + eps_Bm = (- OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2 chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*((+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))/eps_Bp)**2 - chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*((+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))/eps_Bm)**2 + chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*((- OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))/eps_Bm)**2 eps_Bp = (+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2 - eps_Bm = (+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2 + eps_Bm = (- OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2 chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*((+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps_Bp)**2 - chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*((+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps_Bm)**2 + chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*((- OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps_Bm)**2 enddo ZAp(ia,jb) = ZAp(ia,jb) + 2d0*lambda*chi_Ap - ZAm(ia,jb) = ZAm(ia,jb) + 2d0*lambda*chi_Am + ZAm(ia,jb) = ZAm(ia,jb) - 2d0*lambda*chi_Am ZBp(ia,jb) = ZBp(ia,jb) + 2d0*lambda*chi_Bp - ZBm(ia,jb) = ZBm(ia,jb) + 2d0*lambda*chi_Bm + ZBm(ia,jb) = ZBm(ia,jb) - 2d0*lambda*chi_Bm enddo enddo diff --git a/src/QuAcK/G0F2.f90 b/src/QuAcK/G0F2.f90 index 4f03cb6..ce8d8ff 100644 --- a/src/QuAcK/G0F2.f90 +++ b/src/QuAcK/G0F2.f90 @@ -1,4 +1,4 @@ -subroutine G0F2(linearize,nBas,nC,nO,nV,nR,V,e0) +subroutine G0F2(BSE,TDA,singlet_manifold,triplet_manifold,linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Perform a one-shot second-order Green function calculation @@ -7,19 +7,28 @@ subroutine G0F2(linearize,nBas,nC,nO,nV,nR,V,e0) ! Input variables + logical,intent(in) :: BSE + logical,intent(in) :: TDA + logical,intent(in) :: singlet_manifold + logical,intent(in) :: triplet_manifold logical,intent(in) :: linearize + double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nO integer,intent(in) :: nC integer,intent(in) :: nV integer,intent(in) :: nR - double precision,intent(in) :: e0(nBas) - double precision,intent(in) :: V(nBas,nBas,nBas,nBas) + 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) ! Local variables double precision :: eps - double precision :: VV + double precision :: V + double precision :: EcBSE(nspin) double precision,allocatable :: eGF2(:) double precision,allocatable :: Sig(:) double precision,allocatable :: Z(:) @@ -55,10 +64,10 @@ subroutine G0F2(linearize,nBas,nC,nO,nV,nR,V,e0) do j=nC+1,nO do a=nO+1,nBas-nR - eps = e0(p) + e0(a) - e0(i) - e0(j) - VV = (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j) - Sig(p) = Sig(p) + VV/eps - Z(p) = Z(p) + VV/eps**2 + eps = eHF(p) + eHF(a) - eHF(i) - eHF(j) + V = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j) + Sig(p) = Sig(p) + V/eps + Z(p) = Z(p) + V/eps**2 end do end do @@ -70,10 +79,10 @@ subroutine G0F2(linearize,nBas,nC,nO,nV,nR,V,e0) do a=nO+1,nBas-nR do b=nO+1,nBas-nR - eps = e0(p) + e0(i) - e0(a) - e0(b) - VV = (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b) - Sig(p) = Sig(p) + VV/eps - Z(p) = Z(p) + VV/eps**2 + eps = eHF(p) + eHF(i) - eHF(a) - eHF(b) + V = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b) + Sig(p) = Sig(p) + V/eps + Z(p) = Z(p) + V/eps**2 end do end do @@ -84,16 +93,69 @@ subroutine G0F2(linearize,nBas,nC,nO,nV,nR,V,e0) if(linearize) then - eGF2(:) = e0(:) + Z(:)*Sig(:) + eGF2(:) = eHF(:) + Z(:)*Sig(:) else - eGF2(:) = e0(:) + Sig(:) + eGF2(:) = eHF(:) + Sig(:) end if ! Print results - call print_G0F2(nBas,nO,e0,Sig,eGF2,Z) + call print_G0F2(nBas,nO,eHF,Sig,eGF2,Z) +! Perform BSE2 calculation + + if(BSE) then + + call BSE2(TDA,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF2,EcBSE) + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'Tr@BSE2@G0F correlation energy (singlet) =',EcBSE(1) + write(*,'(2X,A50,F20.10)') 'Tr@BSE2@G0F correlation energy (triplet) =',EcBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@BSE2@G0F correlation energy =',EcBSE(1) + EcBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@BSE2@G0F total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2) + 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 ACFDT(exchange_kernel,doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold,eta, & +! nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,Omega,XpY,XmY,rho,EcAC) + +! if(exchange_kernel) then +! +! EcAC(1) = 0.5d0*EcAC(1) +! EcAC(2) = 1.5d0*EcAC(1) +! +! end if + +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy (singlet) =',EcAC(1) +! write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy (triplet) =',EcAC(2) +! write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy =',EcAC(1) + EcAC(2) +! write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 total energy =',ENuc + ERHF + EcAC(1) + EcAC(2) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! end if + + end if end subroutine G0F2 diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index ee79980..62f7723 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -601,7 +601,8 @@ program QuAcK if(doG0F2) then call cpu_time(start_GF2) - call G0F2(linGF,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO,eHF) + call G0F2(BSE,TDA,singlet_manifold,triplet_manifold,linGF, & + eta,nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF) call cpu_time(end_GF2) t_GF2 = end_GF2 - start_GF2 @@ -617,7 +618,8 @@ program QuAcK if(doevGF2) then call cpu_time(start_GF2) - call evGF2(maxSCF_GF,thresh_GF,n_diis_GF,linGF,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO,eHF) + call evGF2(BSE,TDA,maxSCF_GF,thresh_GF,n_diis_GF,singlet_manifold,triplet_manifold,linGF, & + eta,nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF) call cpu_time(end_GF2) t_GF2 = end_GF2 - start_GF2 diff --git a/src/QuAcK/evGF2.f90 b/src/QuAcK/evGF2.f90 index 27539f0..98fb616 100644 --- a/src/QuAcK/evGF2.f90 +++ b/src/QuAcK/evGF2.f90 @@ -1,4 +1,5 @@ -subroutine evGF2(maxSCF,thresh,max_diis,linearize,nBas,nC,nO,nV,nR,V,e0) +subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifold,linearize, & + eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,V,e0) ! Perform eigenvalue self-consistent second-order Green function calculation @@ -7,15 +8,23 @@ subroutine evGF2(maxSCF,thresh,max_diis,linearize,nBas,nC,nO,nV,nR,V,e0) ! Input variables + logical,intent(in) :: BSE + logical,intent(in) :: TDA integer,intent(in) :: maxSCF double precision,intent(in) :: thresh integer,intent(in) :: max_diis + logical,intent(in) :: singlet_manifold + logical,intent(in) :: triplet_manifold logical,intent(in) :: linearize + double precision,intent(in) :: eta 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) :: e0(nBas) double precision,intent(in) :: V(nBas,nBas,nBas,nBas)