symmetrized upfolded phBSE scheme

src/GT/RGTpp_phBSE_dynamic_perturbation.f90

@@ -47,7 +47,6 @@ subroutine RGTpp_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,n
   integer                       :: ia
 
   integer                       :: maxS = 10
-  double precision              :: gapGT
 
   double precision,allocatable  :: OmDyn(:)
   double precision,allocatable  :: ZDyn(:)
@@ -111,13 +110,9 @@ subroutine RGTpp_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,n
 ! Dump results !
 !--------------!
 
-  gapGT = eGT(nO+1) - eGT(nO) 
-
   write(*,*) '---------------------------------------------------------------------------------------------------'
   write(*,*) ' First-order dynamical correction to static Bethe-Salpeter excitation energies                     '
   write(*,*) '---------------------------------------------------------------------------------------------------'
-  write(*,'(A57,F10.6,A3)') ' BSE neutral excitation must be lower than the GT gap = ',gapGT*HaToeV,' eV'
-  write(*,*) '---------------------------------------------------------------------------------------------------'
   write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)'
   write(*,*) '---------------------------------------------------------------------------------------------------'
 

src/GW/RGW_phBSE.f90

@@ -146,6 +146,8 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple
     ! Upfolded phBSE !
     !----------------!
 
+    call RGW_phBSE_upfolded_sym(ispin,nOrb,nOrb,nC,nO,nV,nR,nS,ERI,rho_RPA,OmRPA,eW)
+
     call RGW_phBSE_upfolded(ispin,nOrb,nOrb,nC,nO,nV,nR,nS,ERI,rho_RPA,OmRPA,eGW)
 
   end if

src/GW/RGW_phBSE_dynamic_perturbation.f90

@@ -35,7 +35,6 @@ subroutine RGW_phBSE_dynamic_perturbation(dophBSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,
   integer                       :: ia
 
   integer                       :: maxS = 10
-  double precision              :: gapGW
 
   double precision,allocatable  :: Om_dyn(:)
   double precision,allocatable  :: Z_dyn(:)
@@ -73,13 +72,9 @@ subroutine RGW_phBSE_dynamic_perturbation(dophBSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,
 
   end if
 
-  gapGW = eGW(nO+1) - eGW(nO) 
-
   write(*,*) '---------------------------------------------------------------------------------------------------'
   write(*,*) ' First-order dynamical correction to static Bethe-Salpeter excitation energies                     '
   write(*,*) '---------------------------------------------------------------------------------------------------'
-  write(*,'(A57,F10.6,A3)') ' BSE neutral excitation must be lower than the GW gap = ',gapGW*HaToeV,' eV'
-  write(*,*) '---------------------------------------------------------------------------------------------------'
   write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)'
   write(*,*) '---------------------------------------------------------------------------------------------------'
 

src/GW/RGW_phBSE_upfolded.f90

@@ -27,7 +27,7 @@ subroutine RGW_phBSE_upfolded(ispin,nBas,nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eGW)
   integer                       :: a,b,c,d
   integer                       :: m,ia,jb,iam,kcm
   integer,parameter             :: maxH = 20
-  double precision              :: tmp
+  double precision              :: Jph,Kph,C2h2p
 
   integer                       :: n1h,n1p,n1h1p,n2h2p,nH
   double precision,external     :: Kronecker_delta
@@ -153,15 +153,15 @@ subroutine RGW_phBSE_upfolded(ispin,nBas,nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eGW)
 
             ! Jph
 
-            tmp = - sqrt(2d0)*Kronecker_delta(a,c)*rho(i,k,m)
+            Jph = - sqrt(2d0)*Kronecker_delta(a,c)*rho(i,k,m)
-            H(ia             ,n1h1p+kcm) = tmp
+            H(ia             ,n1h1p+kcm) = Jph
-            H(n1h1p+n2h2p+kcm,ia)        = tmp
+            H(n1h1p+n2h2p+kcm,ia)        = Jph
 
             ! Kph
            
-            tmp = sqrt(2d0)*Kronecker_delta(i,k)*rho(a,c,m)
+            Kph = + sqrt(2d0)*Kronecker_delta(i,k)*rho(a,c,m)
-            H(ia       ,n1h1p+n2h2p+kcm) = tmp
+            H(ia       ,n1h1p+n2h2p+kcm) = Kph
-            H(n1h1p+kcm,ia)              = tmp
+            H(n1h1p+kcm,ia)              = Kph
 
           end do
         end do
@@ -180,9 +180,9 @@ subroutine RGW_phBSE_upfolded(ispin,nBas,nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eGW)
       do m=1,nS
         iam = iam + 1
 
-        tmp = Om(m) + eGW(a) - eGW(i) 
+        C2h2p = Om(m) + eGW(a) - eGW(i) 
-        H(n1h1p      +iam,n1h1p      +iam) = tmp
+        H(n1h1p      +iam,n1h1p      +iam) = C2h2p
-        H(n1h1p+n2h2p+iam,n1h1p+n2h2p+iam) = tmp
+        H(n1h1p+n2h2p+iam,n1h1p+n2h2p+iam) = C2h2p
 
       end do
     end do

src/GW/RGW_phBSE_upfolded_sym.f90

@@ -0,0 +1,259 @@
+subroutine RGW_phBSE_upfolded_sym(ispin,nBas,nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eHF)
+
+! Upfolded phBSE@GW (TDA only!)
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: ispin
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nOrb
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+  double precision,intent(in)   :: ERI(nOrb,nOrb,nOrb,nOrb)
+  double precision,intent(in)   :: rho(nOrb,nOrb,nS)
+  double precision,intent(in)   :: Om(nS)
+  double precision,intent(in)   :: eHF(nOrb)
+
+! Local variables
+
+  integer                       :: s
+  integer                       :: i,j,k,l
+  integer                       :: a,b,c,d
+  integer                       :: m,ia,jb,iam,kcm
+  integer,parameter             :: maxH = 20
+  double precision              :: Jph,Kph,C2h2p
+
+  integer                       :: n1h,n1p,n1h1p,n2h2p,nH
+  double precision,external     :: Kronecker_delta
+
+  integer,allocatable           :: order(:)
+  double precision,allocatable  :: H(:,:)
+  double precision,allocatable  :: OmBSE(:)
+  double precision,allocatable  :: Z(:)
+
+! Output variables
+
+! Hello world
+
+  write(*,*)
+  write(*,*)'*********************************************'
+  write(*,*)'* Symmetrized Upfolded phBSE@GW Calculation *'
+  write(*,*)'*********************************************'
+  write(*,*)
+
+! TDA for W
+
+  write(*,*) 'Tamm-Dancoff approximation by default!'
+  write(*,*)
+
+! Dimension of the supermatrix
+
+  n1h = nO
+  n1p = nV
+
+  n1h1p = n1h*n1p
+  n2h2p = n1h1p*n1h1p
+     nH = n1h1p + n2h2p 
+
+! Memory allocation
+
+  allocate(order(nH),H(nH,nH),OmBSE(nH),Z(nH))
+
+! Initialization
+
+  H(:,:) = 0d0
+
+!-------------------------------!
+!  Compute BSE supermatrix      !
+!-------------------------------!
+!                               !
+!     |     A     Jph + Kph |   ! 
+! H = |                     |   ! 
+!     | Jph + Kph   C2h2p   |   ! 
+!                               !
+!-------------------------------!
+
+  !----------------------!
+  ! Block A for singlets !
+  !----------------------!
+
+  if(ispin == 1) then
+
+    ia = 0
+    do i=nC+1,nO
+      do a=nO+1,nOrb-nR
+        ia = ia + 1
+       
+        jb = 0
+        do j=nC+1,nO
+          do b=nO+1,nOrb-nR
+          jb = jb + 1
+       
+            H(ia,jb) = (eHF(a) - eHF(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
+                     + 2d0*ERI(i,b,a,j) - ERI(i,b,j,a)
+       
+            do k=nC+1,nO
+              do m=1,nS
+
+                H(ia,jb) = H(ia,jb) & 
+                         + Kronecker_delta(i,j)*rho(a,k,m)*rho(b,k,m)/(eHF(a) - eHF(k) + Om(m)) &
+                         + Kronecker_delta(i,j)*rho(a,k,m)*rho(b,k,m)/(eHF(b) - eHF(k) + Om(m))
+
+              end do
+            end do
+       
+            do c=nO+1,nOrb-nR
+              do m=1,nS
+
+                H(ia,jb) = H(ia,jb) & 
+                         - Kronecker_delta(a,b)*rho(i,c,m)*rho(j,c,m)/(eHF(i) - eHF(c) - Om(m)) &
+                         - Kronecker_delta(a,b)*rho(i,c,m)*rho(j,c,m)/(eHF(j) - eHF(c) - Om(m))
+
+              end do
+            end do
+       
+          end do
+        end do
+ 
+      end do
+    end do
+
+  end if
+
+  !----------------------!
+  ! Block A for triplets !
+  !----------------------!
+
+  if(ispin == 2) then
+
+    ia = 0
+    do i=nC+1,nO
+      do a=nO+1,nOrb-nR
+        ia = ia + 1
+ 
+        jb = 0
+        do j=nC+1,nO
+          do b=nO+1,nOrb-nR
+            jb = jb + 1
+       
+            H(ia,jb) = (eHF(a) - eHF(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
+                     - ERI(i,b,j,a)
+
+            do k=nC+1,nO
+              do m=1,nS
+
+                H(ia,jb) = H(ia,jb) & 
+                         + Kronecker_delta(i,j)*rho(a,k,m)*rho(b,k,m)/(eHF(a) - eHF(k) + Om(m)) &
+                         + Kronecker_delta(i,j)*rho(a,k,m)*rho(b,k,m)/(eHF(b) - eHF(k) + Om(m))
+
+              end do
+            end do
+       
+            do c=nO+1,nOrb-nR
+              do m=1,nS
+
+                H(ia,jb) = H(ia,jb) & 
+                         - Kronecker_delta(a,b)*rho(i,c,m)*rho(j,c,m)/(eHF(i) - eHF(c) - Om(m)) &
+                         - Kronecker_delta(a,b)*rho(i,c,m)*rho(j,c,m)/(eHF(j) - eHF(c) - Om(m))
+
+              end do
+            end do
+       
+          end do
+        end do
+       
+      end do
+    end do
+
+  end if
+
+  !------------------!
+  ! Blocks Jph & Kph !
+  !------------------!
+
+  ia = 0
+  do i=nC+1,nO
+    do a=nO+1,nOrb-nR
+      ia = ia + 1
+
+      kcm = 0
+      do k=nC+1,nO
+        do c=nO+1,nOrb-nR
+          do m=1,nS
+            kcm = kcm + 1
+
+            ! Jph + Kph
+
+            Jph = - sqrt(2d0)*Kronecker_delta(a,c)*rho(i,k,m)
+            Kph = + sqrt(2d0)*Kronecker_delta(i,k)*rho(a,c,m)
+            H(ia       ,n1h1p+kcm) = Jph + Kph
+            H(n1h1p+kcm,ia)        = Jph + Kph
+
+          end do
+        end do
+      end do
+
+    end do
+  end do
+
+  !-------------!
+  ! Block 2h2p !
+  !-------------!
+
+  iam = 0
+  do i=nC+1,nO
+    do a=nO+1,nOrb-nR
+      do m=1,nS
+        iam = iam + 1
+
+        C2h2p = Om(m) + eHF(a) - eHF(i) 
+        H(n1h1p+iam,n1h1p+iam) = C2h2p
+
+      end do
+    end do
+  end do
+
+!-------------------------!
+! Diagonalize supermatrix !
+!-------------------------!
+
+  call diagonalize_matrix(nH,H,OmBSE)
+
+!-----------------!
+! Compute weights !
+!-----------------!
+
+  Z(:) = 0d0
+  do s=1,nH
+    do ia=1,n1h1p
+      Z(s) = Z(s) + H(ia,s)**2
+    end do
+  end do
+
+!--------------!
+! Dump results !
+!--------------!
+
+  write(*,*)'-------------------------------------------'
+  write(*,*)'  Upfolded phBSE excitation energies (eV)  '
+  write(*,*)'-------------------------------------------'
+  write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X)') &
+            '|','#','|','OmBSE (eV)','|','Z','|'
+  write(*,*)'-------------------------------------------'
+
+  do s=1,min(nH,maxH)
+    if(Z(s) > 1d-7) &
+      write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
+      '|',s,'|',OmBSE(s)*HaToeV,'|',Z(s),'|'
+  end do
+
+  write(*,*)'-------------------------------------------'
+  write(*,*)
+
+end subroutine 

src/GW/UGW_phBSE_dynamic_perturbation.f90

@@ -41,7 +41,6 @@ subroutine UGW_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nSa
   integer                       :: ia
 
   integer,parameter             :: maxS = 10
-  double precision              :: gapGW
 
   double precision,allocatable  :: OmDyn(:)
   double precision,allocatable  :: ZDyn(:)
@@ -63,13 +62,9 @@ subroutine UGW_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nSa
     write(*,*)
   end if
 
-  gapGW = min(eGW(nO(1)+1,1),eGW(nO(2)+1,2)) - max(eGW(nO(1),1),eGW(nO(2),2))
-
   write(*,*) '---------------------------------------------------------------------------------------------------'
   write(*,*) ' First-order dynamical correction to static Bethe-Salpeter excitation energies                     '
   write(*,*) '---------------------------------------------------------------------------------------------------'
-  write(*,'(A57,F10.6,A3)') ' BSE neutral excitation must be lower than the GW gap = ',gapGW*HaToeV,' eV'
-  write(*,*) '---------------------------------------------------------------------------------------------------'
   write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)'
   write(*,*) '---------------------------------------------------------------------------------------------------'