From 93564f350db3bb9bbe56e942af32a0edd54a2eeb Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 17 Aug 2022 14:32:14 +0200 Subject: [PATCH] ppBSE@G0W0 --- input/methods | 2 +- input/options | 4 +- int/ERI.Hu.dat | 2 + int/Kin.Hu.dat | 2 + int/Nuc.Hu.dat | 2 + int/Ov.Hu.dat | 2 + int/x.Hu.dat | 0 int/y.Hu.dat | 0 int/z.Hu.dat | 0 src/GW/Bethe_Salpeter_pp.f90 | 136 ++++++++++++++++++++++++++++++ src/GW/G0W0.f90 | 17 ++++ src/GW/static_screening_WB_pp.f90 | 99 ++++++++++++++++++++++ src/GW/static_screening_WC_pp.f90 | 100 ++++++++++++++++++++++ src/GW/static_screening_WD_pp.f90 | 99 ++++++++++++++++++++++ src/LR/linear_response_BSE.f90 | 2 +- src/LR/linear_response_B_pp.f90 | 56 ++++++------ src/LR/linear_response_C_pp.f90 | 8 +- src/LR/linear_response_D_pp.f90 | 8 +- src/LR/linear_response_pp.f90 | 24 ++---- src/LR/linear_response_pp_BSE.f90 | 121 ++++++++++++++++++++++++++ 20 files changed, 636 insertions(+), 48 deletions(-) create mode 100644 int/ERI.Hu.dat create mode 100644 int/Kin.Hu.dat create mode 100644 int/Nuc.Hu.dat create mode 100644 int/Ov.Hu.dat create mode 100644 int/x.Hu.dat create mode 100644 int/y.Hu.dat create mode 100644 int/z.Hu.dat create mode 100644 src/GW/Bethe_Salpeter_pp.f90 create mode 100644 src/GW/static_screening_WB_pp.f90 create mode 100644 src/GW/static_screening_WC_pp.f90 create mode 100644 src/GW/static_screening_WD_pp.f90 create mode 100644 src/LR/linear_response_pp_BSE.f90 diff --git a/input/methods b/input/methods index 4fb62d7..14799e4 100644 --- a/input/methods +++ b/input/methods @@ -9,7 +9,7 @@ # CIS* CIS(D) CID CISD FCI F F F F F # RPA* RPAx* crRPA ppRPA - F F F F + F F F T # G0F2* evGF2* qsGF2* G0F3 evGF3 F F F F F # G0W0* evGW* qsGW* ufG0W0 ufGW diff --git a/input/options b/input/options index 5919d78..dca8460 100644 --- a/input/options +++ b/input/options @@ -1,5 +1,5 @@ # HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess level_shift stability - 512 0.0000001 T 5 2 1 F 1.0 F + 512 0.0000001 T 5 2 1 F 0.0 F # MP: # CC: maxSCF thresh DIIS n_diis @@ -15,6 +15,6 @@ # ACFDT: AC Kx XBS F F T # BSE: BSE dBSE dTDA evDyn - T T T F + F F T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/int/ERI.Hu.dat b/int/ERI.Hu.dat new file mode 100644 index 0000000..e7b1d13 --- /dev/null +++ b/int/ERI.Hu.dat @@ -0,0 +1,2 @@ +1 1 1 1 1. +2 2 2 2 1. diff --git a/int/Kin.Hu.dat b/int/Kin.Hu.dat new file mode 100644 index 0000000..f31176f --- /dev/null +++ b/int/Kin.Hu.dat @@ -0,0 +1,2 @@ +1 2 -1.0 +2 1 -1.0 diff --git a/int/Nuc.Hu.dat b/int/Nuc.Hu.dat new file mode 100644 index 0000000..9772f82 --- /dev/null +++ b/int/Nuc.Hu.dat @@ -0,0 +1,2 @@ +1 1 0.0 +2 2 0.1 diff --git a/int/Ov.Hu.dat b/int/Ov.Hu.dat new file mode 100644 index 0000000..f52fb30 --- /dev/null +++ b/int/Ov.Hu.dat @@ -0,0 +1,2 @@ +1 1 1.0 +2 2 1.0 diff --git a/int/x.Hu.dat b/int/x.Hu.dat new file mode 100644 index 0000000..e69de29 diff --git a/int/y.Hu.dat b/int/y.Hu.dat new file mode 100644 index 0000000..e69de29 diff --git a/int/z.Hu.dat b/int/z.Hu.dat new file mode 100644 index 0000000..e69de29 diff --git a/src/GW/Bethe_Salpeter_pp.f90 b/src/GW/Bethe_Salpeter_pp.f90 new file mode 100644 index 0000000..32eac69 --- /dev/null +++ b/src/GW/Bethe_Salpeter_pp.f90 @@ -0,0 +1,136 @@ +subroutine Bethe_Salpeter_pp(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE) + +! Compute the Bethe-Salpeter excitation energies at the pp level + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: TDA_W + logical,intent(in) :: TDA + logical,intent(in) :: singlet + logical,intent(in) :: triplet + + 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 + + integer :: ispin + integer :: isp_W + + integer :: nOO + integer :: nVV + + 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 :: Omega1(:) + double precision,allocatable :: X1(:,:) + double precision,allocatable :: Y1(:,:) + + double precision,allocatable :: Omega2(:) + double precision,allocatable :: X2(:,:) + double precision,allocatable :: Y2(:,:) + + double precision,allocatable :: WB(:,:) + double precision,allocatable :: WC(:,:) + double precision,allocatable :: WD(:,:) + +! Output variables + + double precision,intent(out) :: EcBSE(nspin) + +!--------------------------------- +! Compute (singlet) RPA screening +!--------------------------------- + + isp_W = 1 + EcRPA = 0d0 + + ! Memory allocation + + allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS)) + + call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + +!------------------- +! Singlet manifold +!------------------- + + if(singlet) then + + ispin = 1 + EcBSE(ispin) = 0d0 + + nOO = nO*(nO+1)/2 + nVV = nV*(nV+1)/2 + + allocate(Omega1(nVV),X1(nVV,nVV),Y1(nOO,nVV), & + Omega2(nOO),X2(nVV,nOO),Y2(nOO,nOO), & + WB(nVV,nOO),WC(nVV,nVV),WD(nOO,nOO)) + + if(.not.TDA) call static_screening_WB_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,WB) + call static_screening_WC_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,WC) + call static_screening_WD_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,WD) + + ! Compute BSE excitation energies + + call linear_response_pp_BSE(ispin,TDA,.true.,nBas,nC,nO,nV,nR,nOO,nVV,1d0,eGW,ERI,WB,WC,WD, & + Omega1,X1,Y1,Omega2,X2,Y2,EcBSE(ispin)) + + call print_excitation('pp-BSE (N+2)',ispin,nVV,Omega1) + call print_excitation('pp-BSE (N-2)',ispin,nOO,Omega2) + + deallocate(Omega1,X1,Y1,Omega2,X2,Y2,WB,WC,WD) + + end if + +!------------------- +! Triplet manifold +!------------------- + + if(triplet) then + + ispin = 2 + EcBSE(ispin) = 0d0 + + nOO = nO*(nO-1)/2 + nVV = nV*(nV-1)/2 + + allocate(Omega1(nVV),X1(nVV,nVV),Y1(nOO,nVV), & + Omega2(nOO),X2(nVV,nOO),Y2(nOO,nOO), & + WB(nVV,nOO),WC(nVV,nVV),WD(nOO,nOO)) + + if(.not.TDA) call static_screening_WB_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,WB) + call static_screening_WC_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,WC) + call static_screening_WD_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,WD) + + ! Compute BSE excitation energies + + call linear_response_pp_BSE(ispin,TDA,.true.,nBas,nC,nO,nV,nR,nOO,nVV,1d0,eGW,ERI,WB,WC,WD, & + Omega1,X1,Y1,Omega2,X2,Y2,EcBSE(ispin)) + + call print_excitation('pp-BSE (N+2)',ispin,nVV,Omega1) + call print_excitation('pp-BSE (N-2)',ispin,nOO,Omega2) + + deallocate(Omega1,X1,Y1,Omega2,X2,Y2,WB,WC,WD) + + end if + +end subroutine Bethe_Salpeter_pp diff --git a/src/GW/G0W0.f90 b/src/GW/G0W0.f90 index e6f49bc..fa9367e 100644 --- a/src/GW/G0W0.f90 +++ b/src/GW/G0W0.f90 @@ -15,6 +15,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & logical,intent(in) :: doXBS logical,intent(in) :: COHSEX logical,intent(in) :: BSE + logical :: ppBSE = .true. logical,intent(in) :: TDA_W logical,intent(in) :: TDA logical,intent(in) :: dBSE @@ -49,6 +50,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & double precision :: EcRPA double precision :: EcBSE(nspin) double precision :: EcAC(nspin) + double precision :: EcppBSE(nspin) double precision :: EcGM double precision,allocatable :: SigX(:) double precision,allocatable :: SigC(:) @@ -238,4 +240,19 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & end if + if(ppBSE) then + + call Bethe_Salpeter_pp(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eG0W0,EcppBSE) + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@G0W0 correlation energy (singlet) =',EcppBSE(1) + write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@G0W0 correlation energy (triplet) =',3d0*EcppBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@G0W0 correlation energy =',EcppBSE(1) + 3d0*EcppBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@G0W0 total energy =',ENuc + ERHF + EcppBSE(1) + 3d0*EcppBSE(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + + end if + end subroutine G0W0 diff --git a/src/GW/static_screening_WB_pp.f90 b/src/GW/static_screening_WB_pp.f90 new file mode 100644 index 0000000..07047c8 --- /dev/null +++ b/src/GW/static_screening_WB_pp.f90 @@ -0,0 +1,99 @@ +subroutine static_screening_WB_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega,rho,WB) + +! Compute the VVOO block of the static screening W for the pp-BSE + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + integer,intent(in) :: nOO + integer,intent(in) :: nVV + double precision,intent(in) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + double precision :: chi + double precision :: eps + integer :: a,b,i,j,ab,ij,m + +! Output variables + + double precision,intent(out) :: WB(nVV,nOO) + +! Initialization + + WB(:,:) = 0d0 + +!---------------! +! Singlet block ! +!---------------! + + if(ispin == 1) then + + ab = 0 + do a=nO+1,nBas-nR + do b=a,nBas-nR + ab = ab + 1 + ij = 0 + do i=nC+1,nO + do j=i,nO + ij = ij + 1 + + chi = 0d0 + do m=1,nS + eps = Omega(m)**2 + eta**2 + chi = chi + rho(a,j,m)*rho(b,i,m)*Omega(m)/eps + enddo + + WB(ab,ij) = WB(ab,ij) - 4d0*lambda*chi + + end do + end do + end do + end do + + end if + +!---------------! +! Triplet block ! +!---------------! + + if(ispin == 2) then + + ab = 0 + do a=nO+1,nBas-nR + do b=a+1,nBas-nR + ab = ab + 1 + ij = 0 + do i=nC+1,nO + do j=i+1,nO + ij = ij + 1 + + chi = 0d0 + do m=1,nS + eps = Omega(m)**2 + eta**2 + chi = chi + rho(a,j,m)*rho(b,i,m)*Omega(m)/eps + enddo + + WB(ab,ij) = WB(ab,ij) - 4d0*lambda*chi + + end do + end do + end do + end do + + end if + +end subroutine static_screening_WB_pp diff --git a/src/GW/static_screening_WC_pp.f90 b/src/GW/static_screening_WC_pp.f90 new file mode 100644 index 0000000..2c41605 --- /dev/null +++ b/src/GW/static_screening_WC_pp.f90 @@ -0,0 +1,100 @@ +subroutine static_screening_WC_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega,rho,WC) + +! Compute the VVVV block of the static screening W for the pp-BSE + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + integer,intent(in) :: nOO + integer,intent(in) :: nVV + double precision,intent(in) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + double precision :: chi + double precision :: eps + integer :: a,b,c,d,ab,cd,m + +! Output variables + + double precision,intent(out) :: WC(nVV,nVV) + +! Initialization + + WC(:,:) = 0d0 + +!---------------! +! Singlet block ! +!---------------! + + if(ispin == 1) then + + 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 + + chi = 0d0 + do m=1,nS + eps = Omega(m)**2 + eta**2 + chi = chi + rho(a,d,m)*rho(b,c,m)*Omega(m)/eps + enddo + + WC(ab,cd) = WC(ab,cd) - 4d0*lambda*chi + + + end do + end do + end do + end do + + end if + +!---------------! +! Triplet block ! +!---------------! + + if(ispin == 2) then + + ab = 0 + do a=nO+1,nBas-nR + do b=a+1,nBas-nR + ab = ab + 1 + cd = 0 + do c=nO+1,nBas-nR + do d=c+1,nBas-nR + cd = cd + 1 + + chi = 0d0 + do m=1,nS + eps = Omega(m)**2 + eta**2 + chi = chi + rho(a,d,m)*rho(b,c,m)*Omega(m)/eps + enddo + + WC(ab,cd) = WC(ab,cd) - 4d0*lambda*chi + + end do + end do + end do + end do + + end if + +end subroutine static_screening_WC_pp diff --git a/src/GW/static_screening_WD_pp.f90 b/src/GW/static_screening_WD_pp.f90 new file mode 100644 index 0000000..4dfdf9e --- /dev/null +++ b/src/GW/static_screening_WD_pp.f90 @@ -0,0 +1,99 @@ +subroutine static_screening_WD_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega,rho,WD) + +! Compute the OOOO block of the static screening W for the pp-BSE + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + integer,intent(in) :: nOO + integer,intent(in) :: nVV + double precision,intent(in) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + double precision :: chi + double precision :: eps + integer :: i,j,k,l,ij,kl,m + +! Output variables + + double precision,intent(out) :: WD(nOO,nOO) + +! Initialization + + WD(:,:) = 0d0 + +!---------------! +! Singlet block ! +!---------------! + + if(ispin == 1) then + + ij = 0 + do i=nC+1,nO + do j=i,nO + ij = ij + 1 + kl = 0 + do k=nC+1,nO + do l=k,nO + kl = kl + 1 + + chi = 0d0 + do m=1,nS + eps = Omega(m)**2 + eta**2 + chi = chi + rho(i,l,m)*rho(j,k,m)*Omega(m)/eps + enddo + + WD(ij,kl) = WD(ij,kl) - 4d0*lambda*chi + + end do + end do + end do + end do + + end if + +!---------------! +! Triplet block ! +!---------------! + + if(ispin == 2) then + + ij = 0 + do i=nC+1,nO + do j=i+1,nO + ij = ij + 1 + kl = 0 + do k=nC+1,nO + do l=k+1,nO + kl = kl + 1 + + chi = 0d0 + do m=1,nS + eps = Omega(m)**2 + eta**2 + chi = chi + rho(i,l,m)*rho(j,k,m)*Omega(m)/eps + enddo + + WD(ij,kl) = WD(ij,kl) - 4d0*lambda*chi + + end do + end do + end do + end do + + end if + +end subroutine static_screening_WD_pp diff --git a/src/LR/linear_response_BSE.f90 b/src/LR/linear_response_BSE.f90 index cd02e93..46d2b03 100644 --- a/src/LR/linear_response_BSE.f90 +++ b/src/LR/linear_response_BSE.f90 @@ -1,6 +1,6 @@ subroutine linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A_BSE,B_BSE,Ec,Omega,XpY,XmY) -! Compute linear response +! Compute linear response with BSE additional terms implicit none include 'parameters.h' diff --git a/src/LR/linear_response_B_pp.f90 b/src/LR/linear_response_B_pp.f90 index 46c151e..9cc5860 100644 --- a/src/LR/linear_response_B_pp.f90 +++ b/src/LR/linear_response_B_pp.f90 @@ -8,7 +8,13 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,B_pp) ! Input variables integer,intent(in) :: ispin - integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nOO + integer,intent(in) :: nVV double precision,intent(in) :: lambda double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) @@ -27,11 +33,11 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,B_pp) ab = 0 do a=nO+1,nBas-nR - do b=a,nBas-nR + do b=a,nBas-nR ab = ab + 1 ij = 0 do i=nC+1,nO - do j=i,nO + do j=i,nO ij = ij + 1 B_pp(ab,ij) = lambda*(ERI(a,b,i,j) + ERI(a,b,j,i))/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j))) @@ -43,28 +49,6 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,B_pp) end if -! Build the alpha-beta block of the B matrix - - if(ispin == 3) then - - ab = 0 - do a=nO+1,nBas-nR - do b=nO+1,nBas-nR - ab = ab + 1 - ij = 0 - do i=nC+1,nO - do j=nC+1,nO - ij = ij + 1 - - B_pp(ab,ij) = lambda*ERI(a,b,i,j) - - end do - end do - end do - end do - - end if - ! Build the B matrix for the triplet manifold, or alpha-alpha, or in the spin-orbital basis if(ispin == 2 .or. ispin == 4) then @@ -87,4 +71,26 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,B_pp) end if +! Build the alpha-beta block of the B matrix + + if(ispin == 3) then + + ab = 0 + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + ab = ab + 1 + ij = 0 + do i=nC+1,nO + do j=nC+1,nO + ij = ij + 1 + + B_pp(ab,ij) = lambda*ERI(a,b,i,j) + + end do + end do + end do + end do + + end if + end subroutine linear_response_B_pp diff --git a/src/LR/linear_response_C_pp.f90 b/src/LR/linear_response_C_pp.f90 index b9315a7..452975e 100644 --- a/src/LR/linear_response_C_pp.f90 +++ b/src/LR/linear_response_C_pp.f90 @@ -8,7 +8,13 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,C_pp ! Input variables integer,intent(in) :: ispin - integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nOO + integer,intent(in) :: nVV double precision,intent(in) :: lambda double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas) diff --git a/src/LR/linear_response_D_pp.f90 b/src/LR/linear_response_D_pp.f90 index e05ed6c..d6a799b 100644 --- a/src/LR/linear_response_D_pp.f90 +++ b/src/LR/linear_response_D_pp.f90 @@ -8,7 +8,13 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,D_pp ! Input variables integer,intent(in) :: ispin - integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nOO + integer,intent(in) :: nVV double precision,intent(in) :: lambda double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas) diff --git a/src/LR/linear_response_pp.f90 b/src/LR/linear_response_pp.f90 index c457ff8..772b375 100644 --- a/src/LR/linear_response_pp.f90 +++ b/src/LR/linear_response_pp.f90 @@ -63,22 +63,21 @@ subroutine linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,Om call linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,C) call linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,D) -!call matout(nVV,nVV,C) + if(TDA) then X1(:,:) = +C(:,:) Y1(:,:) = 0d0 - call diagonalize_matrix(nVV,X1,Omega1) + if(nVV > 0) call diagonalize_matrix(nVV,X1,Omega1) X2(:,:) = 0d0 Y2(:,:) = -D(:,:) - call diagonalize_matrix(nOO,Y2,Omega2) + if(nOO > 0) call diagonalize_matrix(nOO,Y2,Omega2) else call linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,B) -!call matout(nVV,nOO,B) -!call matout(nOO,nOO,D) + ! Diagonal blocks M( 1:nVV , 1:nVV) = + C(1:nVV,1:nVV) @@ -88,8 +87,8 @@ subroutine linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,Om M( 1:nVV ,nVV+1:nOO+nVV) = - B(1:nVV,1:nOO) M(nVV+1:nOO+nVV, 1:nVV) = + transpose(B(1:nVV,1:nOO)) -!call matout(nOO+nVV,nOO+nVV,M) - ! Diagonalize the p-h matrix + + ! Diagonalize the p-p matrix if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV,M,Omega,Z) @@ -97,16 +96,6 @@ subroutine linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,Om call sort_ppRPA(nOO,nVV,Omega(:),Z(:,:),Omega1(:),X1(:,:),Y1(:,:),Omega2(:),X2(:,:),Y2(:,:)) -! print*,'Omega1' -! call matout(nVV,1,Omega1) -! print*,'X1t.X1 - Y1t.Y1' -! call matout(nVV,nVV,matmul(transpose(X1),X1) - matmul(transpose(Y1),Y1)) - -! print*,'Omega2' -! call matout(nOO,1,Omega2) -! print*,'Y2t.Y2 - X2t.X2' -! call matout(nOO,nOO,matmul(transpose(Y2),Y2) - matmul(transpose(X2),X2)) - end if ! Compute the RPA correlation energy @@ -114,6 +103,7 @@ subroutine linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,Om EcRPA = 0.5d0*( sum(Omega1(:)) - sum(Omega2(:)) - trace_matrix(nVV,C(:,:)) - trace_matrix(nOO,D(:,:)) ) EcRPA1 = +sum(Omega1(:)) - trace_matrix(nVV,C(:,:)) EcRPA2 = -sum(Omega2(:)) - trace_matrix(nOO,D(:,:)) + if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) & print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!' diff --git a/src/LR/linear_response_pp_BSE.f90 b/src/LR/linear_response_pp_BSE.f90 new file mode 100644 index 0000000..1584ba4 --- /dev/null +++ b/src/LR/linear_response_pp_BSE.f90 @@ -0,0 +1,121 @@ +subroutine linear_response_pp_BSE(ispin,TDA,BSE,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,WB,WC,WD,Omega1,X1,Y1,Omega2,X2,Y2,EcBSE) + +! Compute the p-p channel of BSE + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + logical,intent(in) :: TDA + logical,intent(in) :: BSE + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nOO + integer,intent(in) :: nVV + double precision,intent(in) :: lambda + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: WB(nVV,nOO) + double precision,intent(in) :: WC(nVV,nVV) + double precision,intent(in) :: WD(nOO,nOO) + +! Local variables + + integer :: ab,cd,ij,kl + integer :: p,q,r,s + double precision :: trace_matrix + double precision :: EcBSE1 + double precision :: EcBSE2 + double precision,allocatable :: B(:,:) + double precision,allocatable :: C(:,:) + double precision,allocatable :: D(:,:) + double precision,allocatable :: M(:,:) + double precision,allocatable :: Z(:,:) + double precision,allocatable :: Omega(:) + +! Output variables + + double precision,intent(out) :: Omega1(nVV) + double precision,intent(out) :: X1(nVV,nVV) + double precision,intent(out) :: Y1(nOO,nVV) + double precision,intent(out) :: Omega2(nOO) + double precision,intent(out) :: X2(nVV,nOO) + double precision,intent(out) :: Y2(nOO,nOO) + double precision,intent(out) :: EcBSE + +! Memory allocation + + allocate(B(nVV,nOO),C(nVV,nVV),D(nOO,nOO),M(nOO+nVV,nOO+nVV),Z(nOO+nVV,nOO+nVV),Omega(nOO+nVV)) + +!-------------------------------------------------! +! Solve the p-p eigenproblem ! +!-------------------------------------------------! +! ! +! | C B | | X1 X2 | | w1 0 | | X1 X2 | ! +! | | | | = | | | | ! +! | -Bt -D | | Y1 Y2 | | 0 w2 | | Y1 Y2 | ! +! ! +!-------------------------------------------------! + +! Build B, C and D matrices for the pp channel + + call linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,C) + call linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,D) + + if(BSE) then + + C(:,:) = C(:,:) - WC(:,:) + D(:,:) = D(:,:) - WD(:,:) + + end if + + if(TDA) then + + X1(:,:) = +C(:,:) + Y1(:,:) = 0d0 + if(nVV > 0) call diagonalize_matrix(nVV,X1,Omega1) + + X2(:,:) = 0d0 + Y2(:,:) = -D(:,:) + if(nOO > 0) call diagonalize_matrix(nOO,Y2,Omega2) + + else + + call linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,B) + if(BSE) B(:,:) = B(:,:) - WB(:,:) + + ! Diagonal blocks + + M( 1:nVV , 1:nVV) = + C(1:nVV,1:nVV) + M(nVV+1:nVV+nOO,nVV+1:nVV+nOO) = - D(1:nOO,1:nOO) + + ! Off-diagonal blocks + + M( 1:nVV ,nVV+1:nOO+nVV) = - B(1:nVV,1:nOO) + M(nVV+1:nOO+nVV, 1:nVV) = + transpose(B(1:nVV,1:nOO)) + + ! Diagonalize the p-p matrix + + if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV,M,Omega,Z) + + ! Split the various quantities in p-p and h-h parts + + call sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) + + end if + +! Compute the BSE correlation energy + + EcBSE = 0.5d0*( sum(Omega1(:)) - sum(Omega2(:)) - trace_matrix(nVV,C(:,:)) - trace_matrix(nOO,D(:,:)) ) + EcBSE1 = +sum(Omega1(:)) - trace_matrix(nVV,C(:,:)) + EcBSE2 = -sum(Omega2(:)) - trace_matrix(nOO,D(:,:)) + + if(abs(EcBSE - EcBSE1) > 1d-6 .or. abs(EcBSE - EcBSE2) > 1d-6) & + print*,'!!! Issue in pp-BSE linear reponse calculation BSE1 != BSE2 !!!' + +end subroutine linear_response_pp_BSE