fixing UGW_phBSE and UGW_phACFDT

mol/H2.xyz

@@ -1,4 +1,4 @@
 2
 
 H     0.00000000      0.0 0.0
-H     0.00000000      0.0 1.135
+H     0.00000000      0.0 1.235

src/GW/RGW_phBSE.f90

@@ -65,9 +65,7 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple
            Aph(nS,nS),Bph(nS,nS),KA_sta(nS,nS),KB_sta(nS,nS), &
            OmBSE(nS),XpY_BSE(nS,nS),XmY_BSE(nS,nS))
 
-!-----!
+! TDA 
-! TDA !
-!-----!
 
   if(TDA) then
     write(*,*) 'Tamm-Dancoff approximation activated in phBSE!'

src/GW/UG0W0.f90

@@ -46,9 +46,9 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD
 ! Local variables
 
   logical                       :: print_W = .true.
-  logical                       :: dRPA
+  logical                       :: dRPA_W 
   integer                       :: is
-  integer                       :: ispin
+  integer                       :: isp_W
   double precision              :: flow
   double precision              :: EcRPA
   double precision              :: EcGM(nspin)
@@ -78,7 +78,7 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD
 ! Initialization
 
   EcRPA = 0d0
-  dRPA = .true.
+  dRPA_W = .true.
 
 ! TDA for W
 
@@ -113,10 +113,10 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD
 
 ! Spin-conserving transitions
 
-  ispin = 1
+  isp_W = 1
 
-  call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
+               call phULR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
-  if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
+  if(.not.TDA) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
 
   call phULR(TDA_W,nSa,nSb,nSt,Aph,Bph,EcRPA,Om,XpY,XmY)
   
@@ -172,8 +172,8 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD
 
 ! Compute RPA correlation energy
 
-  call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
+               call phULR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
-  if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
+  if(.not.TDA) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
     
   call phULR(TDA_W,nSa,nSb,nSt,Aph,Bph,EcRPA,Om,XpY,XmY)
 
@@ -205,7 +205,7 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD
 
     if(doACFDT) then
 
-      call UGW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,dophBSE,spin_conserved,spin_flip,eta, & 
+      call UGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,spin_conserved,spin_flip,eta, & 
                        nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGW,EcBSE)
 
       write(*,*)

src/GW/UGW_phACFDT.f90

@@ -1,5 +1,5 @@
-subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,spin_flip,eta, &
+subroutine UGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,spin_conserved,spin_flip,eta, &
-                       nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eW,e,EcAC)
+                       nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eW,eGW,EcAC)
 
 ! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem
 
@@ -10,10 +10,8 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s
 ! Input variables
 
   logical,intent(in)            :: doXBS
-  logical,intent(in)            :: dRPA
   logical,intent(in)            :: TDA_W
   logical,intent(in)            :: TDA
-  logical,intent(in)            :: BSE
   logical,intent(in)            :: exchange_kernel
   logical,intent(in)            :: spin_conserved
   logical,intent(in)            :: spin_flip
@@ -26,19 +24,27 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s
   integer,intent(in)            :: nR(nspin)
   integer,intent(in)            :: nS(nspin)
   double precision,intent(in)   :: eW(nBas,nspin)
-  double precision,intent(in)   :: e(nBas,nspin)
+  double precision,intent(in)   :: eGW(nBas,nspin)
   double precision,intent(in)   :: ERI_aaaa(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: ERI_aabb(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: ERI_bbbb(nBas,nBas,nBas,nBas)
 
 ! Local variables
 
+  logical                       :: dRPA   = .false.
+  logical                       :: dRPA_W = .true.
+
   integer                       :: ispin
   integer                       :: isp_W
   integer                       :: iAC
   double precision              :: lambda
   double precision,allocatable  :: Ec(:,:)
 
+  double precision,allocatable  :: Aph(:,:)
+  double precision,allocatable  :: Bph(:,:)
+  double precision,allocatable  :: KA(:,:)
+  double precision,allocatable  :: KB(:,:)
+
   double precision              :: EcRPA
   double precision,allocatable  :: OmRPA(:)
   double precision,allocatable  :: XpY_RPA(:,:)
@@ -46,14 +52,10 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s
   double precision,allocatable  :: rho_RPA(:,:,:,:)
 
   integer                       :: nS_aa,nS_bb,nS_sc
-  double precision,allocatable  :: Om_sc(:)
-  double precision,allocatable  :: XpY_sc(:,:)
-  double precision,allocatable  :: XmY_sc(:,:)
-
   integer                       :: nS_ab,nS_ba,nS_sf
-  double precision,allocatable  :: Om_sf(:)
+  double precision,allocatable  :: Om(:)
-  double precision,allocatable  :: XpY_sf(:,:)
+  double precision,allocatable  :: XpY(:,:)
-  double precision,allocatable  :: XmY_sf(:,:)
+  double precision,allocatable  :: XmY(:,:)
 
 ! Output variables
 
@@ -63,7 +65,6 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s
 
   allocate(Ec(nAC,nspin))
 
-
 ! Hello World
 
   write(*,*) '-----------------------------------------------------------'
@@ -109,18 +110,26 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s
   nS_sf = nS_ab + nS_ba
 
   allocate(OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin))
+  allocate(Aph(nS_sc,nS_sc),Bph(nS_sc,nS_sc),KA(nS_sc,nS_sc),KB(nS_sc,nS_sc))
 
-  call phULR(isp_W,.true.,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,eW, &
+               call phULR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
-             ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
+  if(.not.TDA) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
+
+  call phULR(TDA_W,nS_aa,nS_bb,nS_sc,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
   call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
 
+  call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & 
+                                 ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,KA)
+  call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & 
+                                 ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,KB)
+
 ! Spin-conserved manifold
 
   if(spin_conserved) then
 
     ispin = 1
 
-    allocate(Om_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc))
+    allocate(Om(nS_sc),XpY(nS_sc,nS_sc),XmY(nS_sc,nS_sc))
 
     write(*,*) '------------------------'
     write(*,*) 'Spin-conserved manifold '
@@ -137,17 +146,29 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s
 
       if(doXBS) then
 
-        call phULR(isp_W,.true.,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, &
+                     call phULR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
-                   ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
+        if(.not.TDA) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
+   
+        call phULR(TDA_W,nS_aa,nS_bb,nS_sc,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
         call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
 
+        call UGW_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda, & 
+                                       ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,KA)
+        call UGW_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda, & 
+                                       ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,KB)
+
       end if
 
-      call phULR(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,e, &
+                   call phULR_A(isp_W,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
-                 ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcAC(ispin),Om_sc,XpY_sc,XmY_sc)
+      if(.not.TDA) call phULR_B(isp_W,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
+
+                   Aph(:,:) = Aph(:,:) + KA(:,:)
+      if(.not.TDA) Bph(:,:) = Bph(:,:) + KB(:,:)
+
+      call phULR(TDA_W,nS_aa,nS_bb,nS_sc,Aph,Bph,EcAC(ispin),OmRPA,XpY_RPA,XmY_RPA)
 
       call phUACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc, &
-                                       ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_sc,XmY_sc,Ec(iAC,ispin))
+                                       ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,XmY,Ec(iAC,ispin))
 
       write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)
 
@@ -162,7 +183,7 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s
     write(*,*) '-----------------------------------------------------------------------------------'
     write(*,*)
 
-    deallocate(Om_sc,XpY_sc,XmY_sc)
+    deallocate(Om,XpY,XmY)
 
   end if
  
@@ -174,7 +195,7 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s
 
     ! Memory allocation
 
-    allocate(Om_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf))
+    allocate(Om(nS_sf),XpY(nS_sf,nS_sf),XmY(nS_sf,nS_sf))
 
     write(*,*) '--------------------'
     write(*,*) ' Spin-flip manifold '
@@ -191,17 +212,24 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s
 
       if(doXBS) then
 
-        call phULR(isp_W,.true.,TDA_W,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eW, &
+                     call phULR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
-                   ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
+        if(.not.TDA) call phULR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
+
+        call phULR(TDA_W,nS_aa,nS_bb,nS_sc,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
         call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
 
       end if
 
-      call phULR(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sc,lambda,e, &
+                   call phULR_A(isp_W,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,lambda,eW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
-                 ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcAC(ispin),Om_sf,XpY_sf,XmY_sf)
+      if(.not.TDA) call phULR_B(isp_W,dRPA,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
+
+                   Aph(:,:) = Aph(:,:) + KA(:,:)
+      if(.not.TDA) Bph(:,:) = Bph(:,:) + KB(:,:)
+
+      call phULR(TDA_W,nS_aa,nS_bb,nS_sf,Aph,Bph,EcAC(ispin),OmRPA,XpY_RPA,XmY_RPA)
 
       call phUACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf, &
-                                       ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_sf,XmY_sf,Ec(iAC,ispin))
+                                       ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,XmY,Ec(iAC,ispin))
 
       write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)
 
@@ -216,7 +244,7 @@ subroutine UGW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,spin_conserved,s
     write(*,*) '-----------------------------------------------------------------------------------'
     write(*,*)
 
-    deallocate(Om_sf,XpY_sf,XmY_sf)
+    deallocate(Om,XpY,XmY)
 
   end if
 

src/GW/UGW_phBSE.f90

@@ -38,21 +38,17 @@ subroutine UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_fli
   integer                       :: ispin
   integer                       :: isp_W
 
-  integer                       :: nS_aa,nS_bb,nS_sc
   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_sc(:)
+  integer                       :: nS_aa,nS_bb,nS_sc
-  double precision,allocatable  :: XpY_BSE_sc(:,:)
-  double precision,allocatable  :: XmY_BSE_sc(:,:)
-
   integer                       :: nS_ab,nS_ba,nS_sf
-  double precision,allocatable  :: OmBSE_sf(:)
+  double precision,allocatable  :: OmBSE(:)
-  double precision,allocatable  :: XpY_BSE_sf(:,:)
+  double precision,allocatable  :: XpY_BSE(:,:)
-  double precision,allocatable  :: XmY_BSE_sf(:,:)
+  double precision,allocatable  :: XmY_BSE(:,:)
 
 ! Output variables
 
@@ -93,7 +89,6 @@ subroutine UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_fli
 
   call UGW_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA)
 
-
 !----------------------------!
 ! Spin-conserved excitations !
 !----------------------------!
@@ -103,15 +98,15 @@ subroutine UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_fli
     ispin = 1
     EcBSE(ispin) = 0d0
 
-    allocate(OmBSE_sc(nS_sc),XpY_BSE_sc(nS_sc,nS_sc),XmY_BSE_sc(nS_sc,nS_sc))
+    allocate(OmBSE(nS_sc),XpY_BSE(nS_sc,nS_sc),XmY_BSE(nS_sc,nS_sc))
 
     ! Compute spin-conserved BSE excitation energies
 
     call phULR(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, &
-               eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcBSE(ispin),OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc)
+               eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
-    call print_excitation_energies('phBSE@GW@UHF','spin-conserved',nS_sc,OmBSE_sc)
+    call print_excitation_energies('phBSE@GW@UHF','spin-conserved',nS_sc,OmBSE)
     call phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, &
-                                  cW,S,OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc)
+                                  cW,S,OmBSE,XpY_BSE,XmY_BSE)
 
     !-------------------------------------------------
     ! Compute the dynamical screening at the BSE level
@@ -120,9 +115,9 @@ subroutine UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_fli
     if(dBSE) &
       call UGW_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,nS_sc,    &
                                           eW,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, &
-                                          OmRPA,rho_RPA,OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc)
+                                          OmRPA,rho_RPA,OmBSE,XpY_BSE,XmY_BSE)
 
-    deallocate(OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc)
+    deallocate(OmBSE,XpY_BSE,XmY_BSE)
 
   end if
 
@@ -137,17 +132,17 @@ subroutine UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_fli
 
     ! Memory allocation
   
-    allocate(OmBSE_sf(nS_sf),XpY_BSE_sf(nS_sf,nS_sf),XmY_BSE_sf(nS_sf,nS_sf))
+    allocate(OmBSE(nS_sf),XpY_BSE(nS_sf,nS_sf),XmY_BSE(nS_sf,nS_sf))
 
     ! Compute spin-flip BSE excitation energies
 
     call phULR(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sc,1d0, &
                eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcBSE(ispin),       & 
-               OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf)
+               OmBSE,XpY_BSE,XmY_BSE)
 
-    call print_excitation_energies('phBSE@GW@UHF','spin-flip',nS_sf,OmBSE_sf)
+    call print_excitation_energies('phBSE@GW@UHF','spin-flip',nS_sf,OmBSE)
     call phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, &
-                                  cW,S,OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf)
+                                  cW,S,OmBSE,XpY_BSE,XmY_BSE)
 
     !-------------------------------------------------
     ! Compute the dynamical screening at the BSE level
@@ -156,9 +151,9 @@ subroutine UGW_phBSE(exchange_kernel,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_fli
     if(dBSE) &
       call UGW_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,nS_sc,    &
                                           eW,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, &
-                                          OmRPA,rho_RPA,OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf)
+                                          OmRPA,rho_RPA,OmBSE,XpY_BSE,XmY_BSE)
 
-    deallocate(OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf)
+    deallocate(OmBSE,XpY_BSE,XmY_BSE)
  
   end if