From c7093e0c2c8cfddeb4a2b05ae85eec1b1c0caf54 Mon Sep 17 00:00:00 2001 From: pfloos Date: Thu, 26 Oct 2023 23:28:53 +0200 Subject: [PATCH] BSE@GHF --- input/methods | 4 +- input/options | 2 +- src/AOtoMO/AOtoMO_transform_GHF.f90 | 22 +++++ src/GW/GG0W0.f90 | 29 +++---- src/GW/GGW.f90 | 10 +-- src/GW/GGW_phBSE.f90 | 123 +++++++++++++++++++++++++++ src/GW/GGW_phBSE_static_kernel_A.f90 | 56 ++++++++++++ src/GW/GGW_phBSE_static_kernel_B.f90 | 56 ++++++++++++ src/HF/print_GHF.f90 | 42 ++++++--- src/RPA/GRPA.f90 | 10 +-- 10 files changed, 308 insertions(+), 46 deletions(-) create mode 100644 src/AOtoMO/AOtoMO_transform_GHF.f90 create mode 100644 src/GW/GGW_phBSE.f90 create mode 100644 src/GW/GGW_phBSE_static_kernel_A.f90 create mode 100644 src/GW/GGW_phBSE_static_kernel_B.f90 diff --git a/input/methods b/input/methods index 0b5c057..a90066f 100644 --- a/input/methods +++ b/input/methods @@ -1,7 +1,7 @@ # RHF UHF GHF ROHF F F T F # MP2 MP3 - T F + F F # CCD pCCD DCD CCSD CCSD(T) F F F F F # drCCD rCCD crCCD lCCD @@ -13,6 +13,6 @@ # G0F2 evGF2 qsGF2 G0F3 evGF3 F F F F F # G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW - F F F F F F + T F F F F F # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh F F F F F F diff --git a/input/options b/input/options index 0f3571d..a899bae 100644 --- a/input/options +++ b/input/options @@ -15,4 +15,4 @@ # ACFDT: AC Kx XBS F F T # BSE: phBSE phBSE2 ppBSE dBSE dTDA - F F F F T + T F F F T diff --git a/src/AOtoMO/AOtoMO_transform_GHF.f90 b/src/AOtoMO/AOtoMO_transform_GHF.f90 new file mode 100644 index 0000000..ac022f7 --- /dev/null +++ b/src/AOtoMO/AOtoMO_transform_GHF.f90 @@ -0,0 +1,22 @@ +subroutine AOtoMO_transform_GHF(nBas,nBas2,Ca,Cb,A,B) + +! Perform AO to MO transformation of a matrix A for given coefficients c + + implicit none + +! Input variables + + integer,intent(in) :: nBas + integer,intent(in) :: nBas2 + double precision,intent(in) :: Ca(nBas,nBas2) + double precision,intent(in) :: Cb(nBas,nBas2) + double precision,intent(inout):: A(nBas,nBas) + +! Output variables + + double precision,intent(inout):: B(nBas2,nBas2) + + B = matmul(transpose(Ca),matmul(A,Ca)) & + + matmul(transpose(Cb),matmul(A,Cb)) + +end subroutine diff --git a/src/GW/GG0W0.f90 b/src/GW/GG0W0.f90 index 6aae026..2d496f9 100644 --- a/src/GW/GG0W0.f90 +++ b/src/GW/GG0W0.f90 @@ -40,7 +40,7 @@ subroutine GG0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,d logical :: dRPA integer :: ispin double precision :: EcRPA - double precision :: EcBSE(nspin) + double precision :: EcBSE double precision :: EcGM double precision,allocatable :: Aph(:,:) double precision,allocatable :: Bph(:,:) @@ -162,25 +162,16 @@ subroutine GG0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,d ! Perform BSE calculation -! if(dophBSE) then + if(dophBSE) then -! call GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE) + call GGW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE) -! if(exchange_kernel) then -! -! EcBSE(1) = 0.5d0*EcBSE(1) -! EcBSE(2) = 1.5d0*EcBSE(2) -! -! end if - -! write(*,*) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0 correlation energy (singlet) =',EcBSE(1),' au' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0 correlation energy (triplet) =',EcBSE(2),' au' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0 correlation energy =',EcBSE(1) + EcBSE(2),' au' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0 total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2),' au' -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,*) + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0 correlation energy =',EcBSE,' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0 total energy =',ENuc + ERHF + EcBSE,' au' + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) ! Compute the BSE correlation energy via the adiabatic connection @@ -226,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 fed9179..55cec7f 100644 --- a/src/GW/GGW.f90 +++ b/src/GW/GGW.f90 @@ -40,11 +40,11 @@ subroutine GGW(doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT, & double precision,intent(in) :: ENuc integer,intent(in) :: nBas - integer,intent(in) :: nC(nspin) - integer,intent(in) :: nO(nspin) - integer,intent(in) :: nV(nspin) - integer,intent(in) :: nR(nspin) - integer,intent(in) :: nS(nspin) + 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 double precision,intent(in) :: epsHF(nBas) diff --git a/src/GW/GGW_phBSE.f90 b/src/GW/GGW_phBSE.f90 new file mode 100644 index 0000000..d11c168 --- /dev/null +++ b/src/GW/GGW_phBSE.f90 @@ -0,0 +1,123 @@ +subroutine GGW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE) + +! Compute the Bethe-Salpeter excitation energies + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: dophBSE2 + logical,intent(in) :: TDA_W + 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) :: eW(nBas) + double precision,intent(in) :: eGW(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + +! Local variables + + logical :: dRPA = .false. + logical :: dRPA_W = .true. + + integer :: ispin + integer :: isp_W + + double precision :: EcRPA + double precision,allocatable :: OmRPA(:) + double precision,allocatable :: XpY_RPA(:,:) + double precision,allocatable :: XmY_RPA(:,:) + double precision,allocatable :: rho_RPA(:,:,:) + + double precision,allocatable :: OmBSE(:) + double precision,allocatable :: XpY_BSE(:,:) + double precision,allocatable :: XmY_BSE(:,:) + + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) + + double precision,allocatable :: KA_sta(:,:) + double precision,allocatable :: KB_sta(:,:) + + double precision,allocatable :: W(:,:,:,:) + +! Output variables + + double precision,intent(out) :: EcBSE + +! Memory allocation + + allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), & + Aph(nS,nS),Bph(nS,nS),KA_sta(nS,nS),KB_sta(nS,nS), & + OmBSE(nS),XpY_BSE(nS,nS),XmY_BSE(nS,nS)) + +!--------------------------------- +! Compute (singlet) RPA screening +!--------------------------------- + + isp_W = 3 + EcRPA = 0d0 + + call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph) + if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) + + call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + + call GGW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA_sta) + call GGW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB_sta) + + ispin = 3 + EcBSE = 0d0 + + ! Compute BSE excitation energies + + call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph) + if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) + + ! Second-order BSE static kernel + +! if(dophBSE2) then + +! allocate(W(nBas,nBas,nBas,nBas)) + +! write(*,*) +! write(*,*) '*** Second-order BSE static kernel activated! ***' +! write(*,*) + +! call GW_phBSE_static_kernel(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,W) +! call GW_phBSE2_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,eW,W,KA_sta) + +! if(.not.TDA) call GW_phBSE2_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,eW,W,KB_sta) + +! deallocate(W) + +! end if + + Aph(:,:) = Aph(:,:) + KA_sta(:,:) + if(.not.TDA) Bph(:,:) = Bph(:,:) + KB_sta(:,:) + + call phLR(TDA,nS,Aph,Bph,EcBSE,OmBSE,XpY_BSE,XmY_BSE) + + call print_excitation_energies('phBSE@GGW',ispin,nS,OmBSE) + call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE) + + !----------------------------------------------------! + ! Compute the dynamical screening at the phBSE level ! + !----------------------------------------------------! + +! if(dBSE) & +! call GW_phBSE_dynamic_perturbation(dophBSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, & +! OmBSE,XpY_BSE,XmY_BSE,KA_sta,KB_sta) + +end subroutine diff --git a/src/GW/GGW_phBSE_static_kernel_A.f90 b/src/GW/GGW_phBSE_static_kernel_A.f90 new file mode 100644 index 0000000..6d38245 --- /dev/null +++ b/src/GW/GGW_phBSE_static_kernel_A.f90 @@ -0,0 +1,56 @@ +subroutine GGW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Om,rho,KA) + +! Compute the BSE static kernel for the resonant block + + implicit none + include 'parameters.h' + +! Input variables + + 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) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: Om(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + double precision :: chi + double precision :: eps + integer :: i,j,a,b,ia,jb,kc + +! Output variables + + double precision,intent(out) :: KA(nS,nS) + +! Compute static kernel + + ia = 0 + do i=nC+1,nO + do a=nO+1,nBas-nR + ia = ia + 1 + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + + chi = 0d0 + do kc=1,nS + eps = Om(kc)**2 + eta**2 + chi = chi + rho(i,j,kc)*rho(a,b,kc)*Om(kc)/eps + enddo + + KA(ia,jb) = 2d0*lambda*chi + + enddo + enddo + enddo + enddo + +end subroutine diff --git a/src/GW/GGW_phBSE_static_kernel_B.f90 b/src/GW/GGW_phBSE_static_kernel_B.f90 new file mode 100644 index 0000000..5afc33f --- /dev/null +++ b/src/GW/GGW_phBSE_static_kernel_B.f90 @@ -0,0 +1,56 @@ +subroutine GGW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Om,rho,KB) + +! Compute the BSE static kernel for the coupling block + + implicit none + include 'parameters.h' + +! Input variables + + 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) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: Om(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + double precision :: chi + double precision :: eps + integer :: i,j,a,b,ia,jb,kc + +! Output variables + + double precision,intent(out) :: KB(nS,nS) + +! Compute static kernel + + ia = 0 + do i=nC+1,nO + do a=nO+1,nBas-nR + ia = ia + 1 + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + + chi = 0d0 + do kc=1,nS + eps = Om(kc)**2 + eta**2 + chi = chi + rho(i,b,kc)*rho(a,j,kc)*Om(kc)/eps + enddo + + KB(ia,jb) = 2d0*lambda*chi + + enddo + enddo + enddo + enddo + +end subroutine diff --git a/src/HF/print_GHF.f90 b/src/HF/print_GHF.f90 index 27a2392..fa4b92f 100644 --- a/src/HF/print_GHF.f90 +++ b/src/HF/print_GHF.f90 @@ -30,6 +30,8 @@ subroutine print_GHF(nBas,nBas2,nO,e,C,P,ENuc,ET,EV,EJ,EK,EHF,dipole) double precision :: Gap double precision :: Sx2,Sy2,Sz2,S2 + double precision,allocatable :: Ca(:,:) + double precision,allocatable :: Cb(:,:) double precision,allocatable :: Paa(:,:) double precision,allocatable :: Pab(:,:) double precision,allocatable :: Pba(:,:) @@ -45,37 +47,49 @@ subroutine print_GHF(nBas,nBas2,nO,e,C,P,ENuc,ET,EV,EJ,EK,EHF,dipole) ! Density matrices - allocate(Paa(nBas,nBas),Pab(nBas,nBas),Pba(nBas,nBas),Pbb(nBas,nBas)) + allocate(Paa(nBas2,nBas2),Pab(nBas2,nBas2),Pba(nBas2,nBas2),Pbb(nBas2,nBas2)) - Paa(:,:) = P( 1:nBas , 1:nBas ) - Pab(:,:) = P( 1:nBas ,nBas+1:nBas2) - Pba(:,:) = P(nBas+1:nBas2, 1:nBas ) - Pbb(:,:) = P(nBas+1:nBas2,nBas+1:nBas2) +! Paa(:,:) = P( 1:nBas , 1:nBas ) +! Pab(:,:) = P( 1:nBas ,nBas+1:nBas2) +! Pba(:,:) = P(nBas+1:nBas2, 1:nBas ) +! Pbb(:,:) = P(nBas+1:nBas2,nBas+1:nBas2) + + allocate(Ca(nBas,nBas2),Cb(nBas,nBas2)) + + Ca(:,:) = C( 1:nBas ,1:nBas2) + Cb(:,:) = C(nBas+1:nBas2,1:nBas2) + + Paa = matmul(transpose(Ca),matmul(P( 1:nBas , 1:nBas ),Ca)) + Pab = matmul(transpose(Ca),matmul(P( 1:nBas ,nBas+1:nBas2),Cb)) + Pba = matmul(transpose(Cb),matmul(P(nBas+1:nBas2, 1:nBas ),Ca)) + Pbb = matmul(transpose(Cb),matmul(P(nBas+1:nBas2,nBas+1:nBas2),Cb)) ! Compute expectation values of S^2 - Sx2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) + 0.25d0*trace_matrix(nBas,Pab+Pba)**2 - do i=1,nBas - do j=1,nBas + Sx2 = 0.25d0*trace_matrix(nBas2,Paa+Pbb) + 0.25d0*trace_matrix(nBas2,Pab+Pba)**2 + do i=1,nBas2 + do j=1,nBas2 Sx2 = Sx2 - 0.5d0*(Paa(i,j)*Pbb(j,i) + Pab(i,j)*Pab(j,i)) end do end do - Sy2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) - 0.25d0*trace_matrix(nBas,Pab+Pba)**2 - do i=1,nBas - do j=1,nBas + Sy2 = 0.25d0*trace_matrix(nBas2,Paa+Pbb) - 0.25d0*trace_matrix(nBas2,Pab+Pba)**2 + do i=1,nBas2 + do j=1,nBas2 Sy2 = Sy2 - 0.5d0*(Paa(i,j)*Pbb(j,i) - Pab(i,j)*Pab(j,i)) end do end do - Sz2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) + 0.25d0*trace_matrix(nBas,Pab-Pba)**2 - do i=1,nBas - do j=1,nBas + Sz2 = 0.25d0*trace_matrix(nBas2,Paa+Pbb) + 0.25d0*trace_matrix(nBas2,Pab-Pba)**2 + do i=1,nBas2 + do j=1,nBas2 Sz2 = Sz2 - 0.25d0*(Paa(i,j)*Pbb(j,i) - Pab(i,j)*Pab(j,i)) Sz2 = Sz2 + 0.25d0*(Pab(i,j)*Pba(j,i) - Pba(i,j)*Pab(j,i)) end do end do + print*,Sx2,Sy2,Sz2 + S2 = Sx2 + Sy2 + Sz2 ! Dump results diff --git a/src/RPA/GRPA.f90 b/src/RPA/GRPA.f90 index 9aae7f6..88fda48 100644 --- a/src/RPA/GRPA.f90 +++ b/src/RPA/GRPA.f90 @@ -15,11 +15,11 @@ subroutine GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas,nC,nO, logical,intent(in) :: doACFDT logical,intent(in) :: exchange_kernel integer,intent(in) :: nBas - integer,intent(in) :: nC(nspin) - integer,intent(in) :: nO(nspin) - integer,intent(in) :: nV(nspin) - integer,intent(in) :: nR(nspin) - integer,intent(in) :: nS(nspin) + 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) :: EHF double precision,intent(in) :: epsHF(nBas)