clean up LR in GW mess

src/CC/CCD.f90

@@ -1,4 +1,4 @@
-subroutine CCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF)
+subroutine CCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF)
 
 ! CCD module
 
@@ -6,7 +6,6 @@ subroutine CCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ER
 
 ! Input variables
 
-  logical,intent(in)            :: BSE
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
   double precision,intent(in)   :: thresh
@@ -93,15 +92,7 @@ subroutine CCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ER
 
   allocate(dbERI(nBas,nBas,nBas,nBas))
 
-  if(BSE) then 
+  call antisymmetrize_ERI(2,nBas,sERI,dbERI)
-    
-    call static_screening(nBas,nC,nO,nV,nR,seHF,sERI,dbERI)
-    
-  else
-
-    call antisymmetrize_ERI(2,nBas,sERI,dbERI)
-
-  end if
 
   deallocate(sERI)
 

src/CC/CCSD.f90

@@ -1,4 +1,4 @@
-subroutine CCSD(BSE,maxSCF,thresh,max_diis,doCCSDT,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF)
+subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF)
 
 ! CCSD module
 
@@ -6,7 +6,6 @@ subroutine CCSD(BSE,maxSCF,thresh,max_diis,doCCSDT,nBasin,nCin,nOin,nVin,nRin,ER
 
 ! Input variables
 
-  logical,intent(in)            :: BSE
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
   double precision,intent(in)   :: thresh
@@ -105,15 +104,7 @@ subroutine CCSD(BSE,maxSCF,thresh,max_diis,doCCSDT,nBasin,nCin,nOin,nVin,nRin,ER
 
   allocate(dbERI(nBas,nBas,nBas,nBas))
 
-  if(BSE) then
+  call antisymmetrize_ERI(2,nBas,sERI,dbERI)
-
-    call static_screening(nBas,nC,nO,nV,nR,seHF,sERI,dbERI)
-
-  else
-
-    call antisymmetrize_ERI(2,nBas,sERI,dbERI)
-
-  end if
 
   deallocate(sERI)
 

src/CC/rCCD.f90

@@ -1,4 +1,4 @@
-subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF,eGW)
+subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF,eGW)
 
 ! Ring CCD module
 
@@ -6,7 +6,6 @@ subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,E
 
 ! Input variables
 
-  logical,intent(in)            :: BSE
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
   double precision,intent(in)   :: thresh
@@ -77,23 +76,11 @@ subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,E
   call spatial_to_spin_MO_energy(nBasin,eGW,nBas,seGW)
   call spatial_to_spin_ERI(nBasin,ERI,nBas,sERI)
 
-! call soRPAx(.false.,nBas,nC,nO,nV,nR,nO*nV,ENuc,ERHF,sERI,seHF)
-
-  call soBSE(.false.,.false.,0.0,nBas,nC,nO,nV,nR,nO*nV,sERI,seHF,seGW)
-
 ! Antysymmetrize ERIs
 
   allocate(dbERI(nBas,nBas,nBas,nBas))
 
-  if(BSE) then
+  call antisymmetrize_ERI(2,nBas,sERI,dbERI)
-
-    call static_screening(nBas,nC,nO,nV,nR,seHF,sERI,dbERI)
-
-  else
-
-    call antisymmetrize_ERI(2,nBas,sERI,dbERI)
-
-  end if
 
   deallocate(sERI)
 

src/GF/GF2_phBSE2.f90

@@ -28,13 +28,15 @@ subroutine GF2_phBSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,n
 
 ! Local variables
 
+  logical                       :: dRPA = .false.
   integer                       :: ispin
   double precision,allocatable  :: OmBSE(:)
   double precision,allocatable  :: XpY(:,:)
   double precision,allocatable  :: XmY(:,:)
-  double precision              :: rho
   double precision,allocatable  :: A_sta(:,:)
   double precision,allocatable  :: B_sta(:,:)
+  double precision,allocatable  :: KA_sta(:,:)
+  double precision,allocatable  :: KB_sta(:,:)
 
 ! Output variables
 
@@ -42,7 +44,8 @@ subroutine GF2_phBSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,n
 
 ! Memory allocation
 
-  allocate(OmBSE(nS),XpY(nS,nS),XmY(nS,nS),A_sta(nS,nS),B_sta(nS,nS))
+  allocate(OmBSE(nS),XpY(nS,nS),XmY(nS,nS),A_sta(nS,nS),KA_sta(nS,nS))
+  allocate(B_sta(nS,nS),KB_sta(nS,nS))
 
 !-------------------
 ! Singlet manifold
@@ -53,18 +56,23 @@ subroutine GF2_phBSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,n
     ispin = 1
     EcBSE(ispin) = 0d0
 
+                 call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI,A_sta)
+    if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,B_sta)
+
     ! Compute static kernel
 
-                 call GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,A_sta)
+                 call GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KA_sta)
-    if(.not.TDA) call GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,B_sta)
+    if(.not.TDA) call GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KB_sta)
 
-    ! Compute BSE2 excitation energies
+                 A_sta(:,:) = A_sta(:,:) + KA_sta(:,:)
+    if(.not.TDA) B_sta(:,:) = B_sta(:,:) + KB_sta(:,:)
 
-    call linear_response_BSE(ispin,.false.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI,-A_sta,-B_sta,EcBSE(ispin),OmBSE,XpY,XmY)
+    ! Compute phBSE2@GF2 excitation energies
-    call print_excitation('BSE2        ',ispin,nS,OmBSE)
+
+    call phLR(TDA,nS,A_sta,B_sta,EcBSE(ispin),OmBSE,XpY,XmY)
+    call print_excitation('phBSE2@GF2  ',ispin,nS,OmBSE)
     call print_transition_vectors_ph(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY)
 
-
     ! Compute dynamic correction for BSE via perturbation theory
 
     if(dBSE) then
@@ -72,10 +80,10 @@ subroutine GF2_phBSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,n
      if(evDyn) then
 
       call GF2_phBSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & 
-                                                     A_sta,B_sta,OmBSE,XpY,XmY)
+                                                     KA_sta,KB_sta,OmBSE,XpY,XmY)
      else
  
-      call GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,A_sta,B_sta,OmBSE,XpY,XmY)
+      call GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY)
  
       end if
 
@@ -92,15 +100,21 @@ subroutine GF2_phBSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,n
     ispin = 2
     EcBSE(ispin) = 0d0
 
+                 call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI,A_sta)
+    if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,B_sta)
+
     ! Compute static kernel
 
-                 call GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,A_sta)
+                 call GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KA_sta)
-    if(.not.TDA) call GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,B_sta)
+    if(.not.TDA) call GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KB_sta)
 
-    ! Compute BSE2 excitation energies
+                 A_sta(:,:) = A_sta(:,:) + KA_sta(:,:)
+    if(.not.TDA) B_sta(:,:) = B_sta(:,:) + KB_sta(:,:)
 
-    call linear_response_BSE(ispin,.false.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI,-A_sta,-B_sta,EcBSE(ispin),OmBSE,XpY,XmY)
+    ! Compute phBSE2@GF2 excitation energies
-    call print_excitation('BSE2        ',ispin,nS,OmBSE)
+
+    call phLR(TDA,nS,A_sta,B_sta,EcBSE(ispin),OmBSE,XpY,XmY)
+    call print_excitation('phBSE2@GF2  ',ispin,nS,OmBSE)
     call print_transition_vectors_ph(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY)
 
     ! Compute dynamic correction for BSE via perturbation theory
@@ -110,10 +124,10 @@ subroutine GF2_phBSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,n
      if(evDyn) then
 
       call GF2_phBSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & 
-                                                     A_sta,B_sta,OmBSE,XpY,XmY)
+                                                     KA_sta,KB_sta,OmBSE,XpY,XmY)
      else
  
-      call GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,A_sta,B_sta,OmBSE,XpY,XmY)
+      call GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY)
  
       end if
 

src/GF/GF2_phBSE2_dynamic_perturbation.f90

@@ -1,4 +1,4 @@
-subroutine GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,A_sta,B_sta,OmBSE,XpY,XmY)
+subroutine GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY)
 
 ! Compute dynamical effects via perturbation theory for BSE
 
@@ -21,8 +21,8 @@ subroutine GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ER
   double precision,intent(in)   :: dipole_int(nBas,nBas,ncart)
   double precision,intent(in)   :: eHF(nBas)
   double precision,intent(in)   :: eGF(nBas)
-  double precision,intent(in)   :: A_sta(nS,nS)
+  double precision,intent(in)   :: KA_sta(nS,nS)
-  double precision,intent(in)   :: B_sta(nS,nS)
+  double precision,intent(in)   :: KB_sta(nS,nS)
   double precision,intent(in)   :: OmBSE(nS)
   double precision,intent(in)   :: XpY(nS,nS)
   double precision,intent(in)   :: XmY(nS,nS)
@@ -38,18 +38,17 @@ subroutine GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ER
   double precision,allocatable  :: X(:)
   double precision,allocatable  :: Y(:)
 
-  double precision,allocatable  ::  Ap_dyn(:,:)
+  double precision,allocatable  :: KAp_dyn(:,:)
-  double precision,allocatable  ::  Am_dyn(:,:)
+  double precision,allocatable  :: KAm_dyn(:,:)
   double precision,allocatable  :: ZAp_dyn(:,:)
   double precision,allocatable  :: ZAm_dyn(:,:)
 
-  double precision,allocatable  ::  B_dyn(:,:)
+  double precision,allocatable  ::  KB_dyn(:,:)
 
 ! Memory allocation
 
-  allocate(OmDyn(nS),ZDyn(nS),X(nS),Y(nS),Ap_dyn(nS,nS),ZAp_dyn(nS,nS))
+  allocate(OmDyn(nS),ZDyn(nS),X(nS),Y(nS),KAp_dyn(nS,nS),ZAp_dyn(nS,nS))
-
+  allocate(KAm_dyn(nS,nS),ZAm_dyn(nS,nS),KB_dyn(nS,nS))
-  if(.not.dTDA) allocate(Am_dyn(nS,nS),ZAm_dyn(nS,nS),B_dyn(nS,nS))
 
   if(dTDA) then
     write(*,*)
@@ -78,28 +77,28 @@ subroutine GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ER
 
     ! Resonant part of the BSE correction for dynamical TDA
 
-    call GF2_phBSE2_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,+OmBSE(ia),Ap_dyn,ZAp_dyn)
+    call GF2_phBSE2_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,+OmBSE(ia),KAp_dyn,ZAp_dyn)
 
     if(dTDA) then 
 
       ZDyn(ia)  = dot_product(X,matmul(ZAp_dyn,X))
-      OmDyn(ia) = dot_product(X,matmul(Ap_dyn - A_sta,X))
+      OmDyn(ia) = dot_product(X,matmul(KAp_dyn - KA_sta,X))
 
     else
 
       ! Second part of the resonant and anti-resonant part of the BSE correction (frequency independent)
 
-      call GF2_phBSE2_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,-OmBSE(ia),Am_dyn,ZAm_dyn)
+      call GF2_phBSE2_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,-OmBSE(ia),KAm_dyn,ZAm_dyn)
 
-      call GF2_phBSE2_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,B_dyn)
+      call GF2_phBSE2_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KB_dyn)
 
       ZDyn(ia)  = dot_product(X,matmul(ZAp_dyn,X)) &
                 + dot_product(Y,matmul(ZAm_dyn,Y))  
 
-      OmDyn(ia) = dot_product(X,matmul(Ap_dyn - A_sta,X)) &
+      OmDyn(ia) = dot_product(X,matmul(KAp_dyn - KA_sta,X)) &
-                - dot_product(Y,matmul(Am_dyn - A_sta,Y)) &
+                - dot_product(Y,matmul(KAm_dyn - KA_sta,Y)) &
-                + dot_product(X,matmul(B_dyn  - B_sta,Y)) & 
+                + dot_product(X,matmul(KB_dyn  - KB_sta,Y)) & 
-                - dot_product(Y,matmul(B_dyn  - B_sta,X))  
+                - dot_product(Y,matmul(KB_dyn  - KB_sta,X))  
 
     end if
 

src/GW/Bethe_Salpeter_A_matrix.f90

@@ -1,53 +0,0 @@
-subroutine Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,A_lr)
-
-! Compute the extra term for Bethe-Salpeter equation for linear response 
-
-  implicit none
-  include 'parameters.h'
-
-! Input variables
-
-  integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
-  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,a,b,ia,jb,kc
-
-! Output variables
-
-  double precision,intent(out)  :: A_lr(nS,nS)
-
-  jb = 0
-!$omp parallel do default(private) shared(A_lr,ERI,Omega,rho,nO,nBas,nS,chi,eps,eta,nC,nR,lambda)
-  do j=nC+1,nO
-    do b=nO+1,nBas-nR
-      jb = (b-nO) + (j-1)*(nBas-nO)
-
-      ia = 0
-      do i=nC+1,nO
-        do a=nO+1,nBas-nR
-          ia = (a-nO) + (i-1)*(nBas-nO) 
- 
-          chi = 0d0
-          do kc=1,nS
-            eps = Omega(kc)**2 + eta**2
-!           chi = chi + lambda*rho(i,j,kc)*rho(a,b,kc)*Omega(kc)/eps
-            chi = chi + rho(i,j,kc)*rho(a,b,kc)*Omega(kc)/eps
-          enddo
-
-          A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI(i,b,j,a) + 4d0*lambda*chi
-
-        enddo
-      enddo
-    enddo
-  enddo
-!$omp end parallel do
-
-end subroutine Bethe_Salpeter_A_matrix

src/GW/Bethe_Salpeter_B_matrix.f90

@@ -1,50 +0,0 @@
-subroutine Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,B_lr)
-
-! Compute the extra term for Bethe-Salpeter equation for linear response 
-
-  implicit none
-  include 'parameters.h'
-
-! Input variables
-
-  integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
-  double precision,intent(in)   :: eta
-  double precision,intent(in)   :: lambda
-  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
-  double precision,intent(in)   :: OmRPA(nS)
-  double precision,intent(in)   :: rho_RPA(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)  :: B_lr(nS,nS)
-
-  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 = OmRPA(kc)**2 + eta**2
-!           chi = chi + lambda*rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*OmRPA(kc)/eps
-            chi = chi + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*OmRPA(kc)/eps
-          enddo
-
-          B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI(i,j,b,a) + 4d0*lambda*chi
-
-        enddo
-      enddo
-    enddo
-  enddo
-
-end subroutine Bethe_Salpeter_B_matrix

src/GW/G0W0.f90

@@ -107,7 +107,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT
 ! Compute screening !
 !-------------------!
 
-  call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI_MO,Aph)
+                 call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI_MO,Aph)
   if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
 
   call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)

src/GW/GW_phACFDT.f90

@@ -30,6 +30,7 @@ subroutine GW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,e
 
 ! Local variables
 
+  logical                       :: dRPA_W = .true.
   integer                       :: ispin
   integer                       :: isp_W
   integer                       :: iAC
@@ -57,8 +58,9 @@ subroutine GW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,e
 ! Memory allocation
 
   allocate(Ec(nAC,nspin))
-  allocate(KA(nS,nS),KB(nS,nS),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS))
+  allocate(Aph(nS,nS),KA(nS,nS),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), & 
-  allocate(Om(nS),XpY(nS,nS),XmY(nS,nS))
+           rho_RPA(nBas,nBas,nS),Om(nS),XpY(nS,nS),XmY(nS,nS))
+  if(.not.TDA) allocate(Aph(nS,nS),KB(nS,nS))
 
 ! Antisymmetrized kernel version
 
@@ -78,11 +80,14 @@ subroutine GW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,e
   isp_W = 1
   EcRPA = 0d0
 
-  call phLR(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
+                 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 BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA)
+  call GW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA)
-  call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB)
+  call GW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB)
 
 ! Singlet manifold
 
@@ -105,15 +110,21 @@ subroutine GW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,e
 
       if(doXBS) then
 
-        call phLR(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
+                       call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,Aph)
+        if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,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 BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KA)
+        call GW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KA)
-        call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KB)
+        call GW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KB)
 
       end if
 
-      call linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,KA,KB,EcAC(ispin),Om,XpY,XmY)
+                   Aph(:,:) = Aph(:,:) + KA(:,:)
+      if(.not.TDA) Bph(:,:) = Bph(:,:) + KB(:,:)
+
+      call phLR(TDA,nS,Aph,Bph,EcAC(ispin),Om,XpY,XmY)
 
       call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin))
 
@@ -153,15 +164,21 @@ subroutine GW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,e
 
       if(doXBS) then
 
-        call phLR(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
+                       call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,Aph)
+        if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,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 BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KA)
+        call GW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KA)
-        call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KB)
+        call GW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KB)
 
       end if  
 
-      call linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,KA,KB,EcAC(ispin),Om,XpY,XmY)
+                   Aph(:,:) = Aph(:,:) + KA(:,:)
+      if(.not.TDA) Bph(:,:) = Bph(:,:) + KB(:,:)
+
+      call phLR(TDA,nS,Aph,Bph,EcAC(ispin),Om,XpY,XmY)
 
       call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin))
 

src/GW/GW_phBSE.f90

@@ -30,6 +30,9 @@ subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,n
 
 ! Local variables
 
+  logical                       :: dRPA   = .false.
+  logical                       :: dRPA_W = .true.
+
   integer                       :: ispin
   integer                       :: isp_W
 
@@ -43,6 +46,9 @@ subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,n
   double precision,allocatable  :: XpY_BSE(:,:)
   double precision,allocatable  :: XmY_BSE(:,:)
 
+  double precision,allocatable  :: A_sta(:,:)
+  double precision,allocatable  :: B_sta(:,:)
+
   double precision,allocatable  :: KA_sta(:,:)
   double precision,allocatable  :: KB_sta(:,:)
 
@@ -57,7 +63,8 @@ subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,n
 ! Memory allocation
 
   allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), &
-           KA_sta(nS,nS),KB_sta(nS,nS),OmBSE(nS),XpY_BSE(nS,nS),XmY_BSE(nS,nS))
+           A_sta(nS,nS),KA_sta(nS,nS),KA2_sta(nS,nS),OmBSE(nS),XpY_BSE(nS,nS),XmY_BSE(nS,nS))
+  allocate(B_sta(nS,nS),KB_sta(nS,nS),KB2_sta(nS,nS))
 
 !---------------------------------
 ! Compute (singlet) RPA screening 
@@ -66,11 +73,14 @@ subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,n
   isp_W = 1
   EcRPA = 0d0
 
-  call phLR(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
+                 call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,A_sta)
+  if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,ERI,B_sta)
+
+  call phLR(TDA_W,nS,A_sta,B_sta,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
   call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
 
-  call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA_sta)
+  call GW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA_sta)
-  call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB_sta)
+  call GW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB_sta)
 
 !-------------------
 ! Singlet manifold
@@ -81,27 +91,39 @@ subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,n
     ispin = 1
     EcBSE(ispin) = 0d0
 
+    ! Compute BSE excitation energies
+
+                 call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,A_sta)
+    if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,B_sta)
+
+                 A_sta(:,:) = A_sta(:,:) + KA_sta(:,:)
+    if(.not.TDA) B_sta(:,:) = B_sta(:,:) + KB_sta(:,:)
+
     ! Second-order BSE static kernel
   
-    allocate(W(nBas,nBas,nBas,nBas),KA2_sta(nS,nS),KB2_sta(nS,nS))
-    KA2_sta(:,:) = 0d0
-    KB2_sta(:,:) = 0d0
-
     if(BSE2) then 
 
-      write(*,*) '*** Second-order BSE static kernel activated! ***'
+      allocate(W(nBas,nBas,nBas,nBas))
-      call static_kernel_W(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,W)
-      call BSE2_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,1d0,eW,W,KA2_sta)
 
-      if(.not.TDA) call BSE2_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,1d0,eW,W,KB2_sta)
+      write(*,*) 
+      write(*,*) '*** Second-order BSE static kernel activated! ***'
+      write(*,*) 
+
+      call static_kernel_W(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,KA2_sta)
+
+      if(.not.TDA) call GW_phBSE2_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,eW,W,KB2_sta)
+
+      deallocate(W)
+
+                   A_sta(:,:) = A_sta(:,:) + KA2_sta(:,:)
+      if(.not.TDA) B_sta(:,:) = B_sta(:,:) + KB2_sta(:,:)
 
     end if
 
-    ! Compute BSE excitation energies
+    call phLR(TDA,nS,A_sta,B_sta,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
 
-    call linear_response_BSE(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,KA_sta-KA2_sta,KB_sta-KB2_sta, &
+    call print_excitation('phBSE@GW    ',ispin,nS,OmBSE)
-                             EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
-    call print_excitation('BSE@GW      ',ispin,nS,OmBSE)
     call print_transition_vectors_ph(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)
 
     !-------------------------------------------------
@@ -114,12 +136,12 @@ subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,n
  
       if(evDyn) then
  
-        call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int,OmRPA,rho_RPA, &
+        call GW_phBSE_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int,OmRPA,rho_RPA, &
-                                                           OmBSE,XpY_BSE,XmY_BSE)
+                                                     OmBSE,XpY_BSE,XmY_BSE)
       else
 
-        call Bethe_Salpeter_dynamic_perturbation(BSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,dipole_int,OmRPA,rho_RPA, &
+        call GW_phBSE_dynamic_perturbation(BSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,dipole_int,OmRPA,rho_RPA, &
-                            OmBSE,XpY_BSE,XmY_BSE,W,KA2_sta)
+                                           OmBSE,XpY_BSE,XmY_BSE,W,KA2_sta)
       end if
 
     end if
@@ -137,8 +159,15 @@ subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,n
 
     ! Compute BSE excitation energies
 
-    call linear_response_BSE(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,KA_sta,KB_sta,EcBSE,OmBSE,XpY_BSE,XmY_BSE)
+                 call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,A_sta)
-    call print_excitation('BSE@GW      ',ispin,nS,OmBSE)
+    if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,B_sta)
+
+                 A_sta(:,:) = A_sta(:,:) + KA_sta(:,:)
+    if(.not.TDA) B_sta(:,:) = B_sta(:,:) + KB_sta(:,:)
+
+    call phLR(TDA,nS,A_sta,B_sta,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
+
+    call print_excitation('phBSE@GW    ',ispin,nS,OmBSE)
     call print_transition_vectors_ph(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)
 
     !-------------------------------------------------
@@ -151,12 +180,12 @@ subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,n
 
       if(evDyn) then
      
-        call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int,OmRPA,rho_RPA, &
+        call GW_phBSE_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int,OmRPA,rho_RPA, &
-                                                           OmBSE,XpY_BSE,XmY_BSE)
+                                                     OmBSE,XpY_BSE,XmY_BSE)
       else
      
-        call Bethe_Salpeter_dynamic_perturbation(BSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,dipole_int,OmRPA,rho_RPA, &
+        call GW_phBSE_dynamic_perturbation(BSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,dipole_int,OmRPA,rho_RPA, &
-                             OmBSE,XpY_BSE,XmY_BSE,W,KA2_sta)
+                                           OmBSE,XpY_BSE,XmY_BSE,W,KA2_sta)
       end if
 
     end if

src/GW/GW_phBSE2_dynamic_kernel_A.f90

@@ -1,4 +1,4 @@
-subroutine BSE2_GW_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE,A_dyn,ZA_dyn)
+subroutine GW_phBSE2_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE,A_dyn,ZA_dyn)
 
 ! Compute the dynamic part of the Bethe-Salpeter equation matrices
 
@@ -92,4 +92,4 @@ subroutine BSE2_GW_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE,A_dyn,ZA
   enddo
 !$omp end parallel do
 
-end subroutine BSE2_GW_A_matrix_dynamic
+end subroutine 

src/GW/GW_phBSE2_static_kernel_A.f90

@@ -1,4 +1,4 @@
-subroutine BSE2_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,eW,W,KA2_sta)
+subroutine GW_phBSE2_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,eW,W,KA2_sta)
 
 ! Compute the second-order static BSE kernel for the resonant block (only for singlets!)
 
@@ -91,4 +91,4 @@ subroutine BSE2_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,eW,W,KA2_sta)
 
 !$omp end parallel do
 
-end subroutine BSE2_static_kernel_KA
+end subroutine 

src/GW/GW_phBSE2_static_kernel_B.f90

@@ -1,4 +1,4 @@
-subroutine BSE2_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,eW,W,KB2_sta)
+subroutine GW_phBSE2_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,lambda,eW,W,KB2_sta)
 
 ! Compute the second-order static BSE kernel for the coupling block (only for singlets!)
 
@@ -91,4 +91,4 @@ subroutine BSE2_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,eW,W,KB2_sta)
 
 !$omp end parallel do
 
-end subroutine BSE2_static_kernel_KB
+end subroutine 

src/GW/GW_phBSE_dynamic_kernel.f90

@@ -1,4 +1,4 @@
-subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rhO_RPA,OmBSE,Ap,Am,Bp,Bm)
+subroutine GW_phBSE_dynamic_kernel(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rhO_RPA,OmBSE,Ap,Am,Bp,Bm)
 
 ! Compute the dynamic part of the Bethe-Salpeter equation matrices
 
@@ -115,4 +115,4 @@ subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O
     enddo
   enddo
 
-end subroutine Bethe_Salpeter_AB_matrix_dynamic
+end subroutine 

src/GW/GW_phBSE_dynamic_kernel_A.f90

@@ -1,4 +1,4 @@
-subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE,A_dyn,ZA_dyn)
+subroutine GW_phBSE_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE,A_dyn,ZA_dyn)
 
 ! Compute the dynamic part of the Bethe-Salpeter equation matrices
 
@@ -92,4 +92,4 @@ subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,Om
 
 !$omp end parallel do
 
-end subroutine Bethe_Salpeter_A_matrix_dynamic
+end subroutine 

src/GW/GW_phBSE_dynamic_kernel_Z.f90

@@ -1,5 +1,4 @@
-subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE, & 
+subroutine GW_phBSE_dynamic_kernel_Z(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE,ZAp,ZAm,ZBp,ZBm)
-                                             ZAp,ZAm,ZBp,ZBm)
 
 ! Compute the dynamic part of the renormalization for the Bethe-Salpeter equation matrices
 
@@ -99,4 +98,4 @@ subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,
     enddo
   enddo
 
-end subroutine Bethe_Salpeter_ZAB_matrix_dynamic
+end subroutine 

src/GW/GW_phBSE_dynamic_kernel_ZA.f90

@@ -1,4 +1,4 @@
-subroutine Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE,ZA_dyn)
+subroutine GW_phBSE_dynamic_kernel_ZA(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE,ZA_dyn)
 
 ! Compute the dynamic part of the Bethe-Salpeter equation matrices
 
@@ -63,4 +63,4 @@ subroutine Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O
     enddo
   enddo
 
-end subroutine Bethe_Salpeter_ZA_matrix_dynamic
+end subroutine 

src/GW/GW_phBSE_dynamic_perturbation.f90

@@ -1,5 +1,5 @@
-subroutine Bethe_Salpeter_dynamic_perturbation(BSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW, &
+subroutine GW_phBSE_dynamic_perturbation(BSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,dipole_int, & 
-dipole_int,OmRPA,rho_RPA,OmBSE,XpY,XmY,W,A_stat)
+                                         OmRPA,rho_RPA,OmBSE,XpY,XmY,W,A_stat)
 
 ! Compute dynamical effects via perturbation theory for BSE
 
@@ -57,8 +57,7 @@ dipole_int,OmRPA,rho_RPA,OmBSE,XpY,XmY,W,A_stat)
 
   maxS = min(nS,maxS)
   allocate(OmDyn(maxS),ZDyn(maxS),X(nS),Y(nS),Ap_dyn(nS,nS),ZAp_dyn(nS,nS))
-
+  allocate(Am_dyn(nS,nS),ZAm_dyn(nS,nS),Bp_dyn(nS,nS),ZBp_dyn(nS,nS),Bm_dyn(nS,nS),ZBm_dyn(nS,nS))
-  if(.not.dTDA) allocate(Am_dyn(nS,nS),ZAm_dyn(nS,nS),Bp_dyn(nS,nS),ZBp_dyn(nS,nS),Bm_dyn(nS,nS),ZBm_dyn(nS,nS))
 
   if(dTDA) then 
     write(*,*)
@@ -87,9 +86,9 @@ dipole_int,OmRPA,rho_RPA,OmBSE,XpY,XmY,W,A_stat)
 
       ! Resonant part of the BSE correction for dynamical TDA
 
-      call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia),Ap_dyn,ZAp_dyn)
+      call GW_phBSE_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia),Ap_dyn,ZAp_dyn)
 
-if(BSE2) call BSE2_GW_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE(ia),Ap_dyn,ZAp_dyn,W)
+      if(BSE2) call GW_phBSE2_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE(ia),Ap_dyn,ZAp_dyn,W)
 
       ZDyn(ia)  = dot_product(X,matmul(ZAp_dyn,X))
       OmDyn(ia) = dot_product(X,matmul( Ap_dyn - A_stat,X))
@@ -98,13 +97,11 @@ if(BSE2) call BSE2_GW_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE(ia),A
 
       ! Resonant and anti-resonant part of the BSE correction
 
-      call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia), &
+      call GW_phBSE_dynamic_kernel(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia),Ap_dyn,Am_dyn,Bp_dyn,Bm_dyn)
-                                            Ap_dyn,Am_dyn,Bp_dyn,Bm_dyn)
 
       ! Renormalization factor of the resonant and anti-resonant parts
 
-      call Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia), &
+      call GW_phBSE_dynamic_kernel_Z(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia),ZAp_dyn,ZAm_dyn,ZBp_dyn,ZBm_dyn)
-                                             ZAp_dyn,ZAm_dyn,ZBp_dyn,ZBm_dyn)
 
       ZDyn(ia)  = dot_product(X,matmul(ZAp_dyn,X)) &
                 - dot_product(Y,matmul(ZAm_dyn,Y)) &
@@ -128,4 +125,4 @@ if(BSE2) call BSE2_GW_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE(ia),A
   write(*,*) '---------------------------------------------------------------------------------------------------'
   write(*,*) 
 
-end subroutine Bethe_Salpeter_dynamic_perturbation
+end subroutine 

src/GW/GW_phBSE_dynamic_perturbation_iterative.f90

@@ -1,5 +1,4 @@
-subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int, & 
+subroutine GW_phBSE_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int,OmRPA,rho_RPA,OmBSE,XpY,XmY)
-                                                         OmRPA,rho_RPA,OmBSE,XpY,XmY)
 
 ! Compute self-consistently the dynamical effects via perturbation theory for BSE
 
@@ -57,9 +56,7 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,
 
   maxS = min(nS,maxS)
   allocate(OmDyn(maxS),OmOld(maxS),X(nS),Y(nS),Ap_dyn(nS,nS),ZAp_dyn(nS,nS))
-
+  allocate(Am_dyn(nS,nS),ZAm_dyn(nS,nS),Bp_dyn(nS,nS),Bm_dyn(nS,nS))
-  if(.not.dTDA) &
-    allocate(Am_dyn(nS,nS),ZAm_dyn(nS,nS),Bp_dyn(nS,nS),Bm_dyn(nS,nS))
 
   if(dTDA) then
     write(*,*)
@@ -101,7 +98,7 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,
 
        ! Resonant part of the BSE correction
  
-        call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmOld(ia),Ap_dyn,ZAp_dyn)
+        call GW_phBSE_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmOld(ia),Ap_dyn,ZAp_dyn)
  
         OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:)))
  
@@ -109,8 +106,7 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,
  
         ! Anti-resonant part of the BSE correction
  
-        call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmOld(ia), &
+        call GW_phBSE_dynamic_kernel(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmOld(ia),Ap_dyn,Am_dyn,Bp_dyn,Bm_dyn)
-                                              Ap_dyn,Am_dyn,Bp_dyn,Bm_dyn)
  
         OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) &
                   - dot_product(Y(:),matmul(Am_dyn(:,:),Y(:))) &
@@ -146,4 +142,4 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,
 
   endif
 
-end subroutine Bethe_Salpeter_dynamic_perturbation_iterative
+end subroutine 

src/GW/GW_phBSE_static_kernel_A.f90

@@ -1,4 +1,4 @@
-subroutine BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,WA)
+subroutine GW_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
 
@@ -11,7 +11,7 @@ subroutine BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,WA)
   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)   :: Om(nS)
   double precision,intent(in)   :: rho(nBas,nBas,nS)
 
 ! Local variables
@@ -22,11 +22,9 @@ subroutine BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,WA)
 
 ! Output variables
 
-  double precision,intent(out)  :: WA(nS,nS)
+  double precision,intent(out)  :: KA(nS,nS)
 
-! Initialize 
+! Compute static kernel
-
-  WA(:,:) = 0d0
 
   ia = 0
   do i=nC+1,nO
@@ -39,15 +37,15 @@ subroutine BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,WA)
 
           chi = 0d0
           do kc=1,nS
-            eps = Omega(kc)**2 + eta**2
+            eps = Om(kc)**2 + eta**2
-            chi = chi + rho(i,j,kc)*rho(a,b,kc)*Omega(kc)/eps
+            chi = chi + rho(i,j,kc)*rho(a,b,kc)*Om(kc)/eps
           enddo
 
-          WA(ia,jb) = WA(ia,jb) + lambda*ERI(i,b,j,a) - 4d0*lambda*chi
+          KA(ia,jb) = 4d0*lambda*chi
 
         enddo
       enddo
     enddo
   enddo
 
-end subroutine BSE_static_kernel_KA
+end subroutine 

src/GW/GW_phBSE_static_kernel_B.f90

@@ -1,4 +1,4 @@
-subroutine BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,KB)
+subroutine GW_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
 
@@ -16,7 +16,7 @@ subroutine BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,KB)
   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)   :: Om(nS)
   double precision,intent(in)   :: rho(nBas,nBas,nS)
 
 ! Local variables
@@ -29,9 +29,7 @@ subroutine BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,KB)
 
   double precision,intent(out)  :: KB(nS,nS)
 
-! Initialize 
+! Compute static kernel
-
-  KB(:,:) = 0d0
 
   ia = 0
   do i=nC+1,nO
@@ -44,15 +42,15 @@ subroutine BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,KB)
 
           chi = 0d0
           do kc=1,nS
-            eps = Omega(kc)**2 + eta**2
+            eps = Om(kc)**2 + eta**2
-            chi = chi + rho(i,b,kc)*rho(a,j,kc)*Omega(kc)/eps
+            chi = chi + rho(i,b,kc)*rho(a,j,kc)*Om(kc)/eps
           enddo
 
-          KB(ia,jb) = KB(ia,jb) + lambda*ERI(i,j,b,a) - 4d0*lambda*chi
+          KB(ia,jb) = 4d0*lambda*chi
 
         enddo
       enddo
     enddo
   enddo
 
-end subroutine BSE_static_kernel_KB
+end subroutine 

src/GW/soBSE.f90

@@ -1,121 +0,0 @@
-subroutine soBSE(TDA_W,TDA,eta,nBas,nC,nO,nV,nR,nS,ERI,eW,eGW)
-
-! Compute the Bethe-Salpeter excitation energies
-
-  implicit none
-  include 'parameters.h'
-
-! Input variables
-
-  logical,intent(in)            :: TDA_W
-  logical,intent(in)            :: TDA
-
-  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)
-
-! Local variables
-
-  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  :: WA_sta(:,:)
-  double precision,allocatable  :: WB_sta(:,:)
-
-  double precision,allocatable  :: X(:,:)
-  double precision,allocatable  :: Y(:,:)
-  double precision,allocatable  :: Xinv(:,:)
-  double precision,allocatable  :: t(:,:,:,:)
-
-  integer                       ::i,a,j,b,k,c,ia,jb,kc
-
-  double precision              :: EcBSE
-
-! Memory allocation
-
-  allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), &
-           WA_sta(nS,nS),WB_sta(nS,nS),OmBSE(nS),XpY_BSE(nS,nS),XmY_BSE(nS,nS))
-
-!---------------------------------
-! Compute (singlet) RPA screening 
-!---------------------------------
-
-  isp_W = 3
-  EcRPA = 0d0
-
-  call phLR(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
-  call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
-
-  call static_screening_WA_so(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WA_sta)
-  call static_screening_WB_so(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WB_sta)
-
-
-  ispin = 3
-  EcBSE = 0d0
-
-  ! Compute BSE excitation energies
-
-  call linear_response_BSE(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,WA_sta,WB_sta, &
-                           EcBSE,OmBSE,XpY_BSE,XmY_BSE)
-  call print_excitation('soBSE@GW    ',ispin,nS,OmBSE)
-
-    write(*,*)
-    write(*,*)'-------------------------------------------------------------------------------'
-    write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0 correlation energy =',0.5d0*EcBSE,' au'
-    write(*,*)'-------------------------------------------------------------------------------'
-    write(*,*)
-
-! allocate(X(nS,nS),Y(nS,nS),Xinv(nS,nS),t(nO,nO,nV,nV))
-
-! X(:,:) = transpose(0.5d0*(XpY_BSE(:,:) + XmY_BSE(:,:)))
-! Y(:,:) = transpose(0.5d0*(XpY_BSE(:,:) - XmY_BSE(:,:)))
-
-! call matout(nS,nS,matmul(X,transpose(X))-matmul(Y,transpose(Y)))
-
-! call inverse_matrix(nS,X,Xinv)
-
-! t = 0d0
-! ia = 0
-! do i=1,nO
-!   do a=1,nV
-!     ia = ia + 1
-
-!     jb = 0
-!     do j=1,nO
-!       do b=1,nV
-!       jb = jb + 1
-
-!       kc = 0
-!       do k=1,nO
-!         do c=1,nV
-!         kc = kc + 1
-
-!           t(i,j,a,b) = t(i,j,a,b) + Y(ia,kc)*Xinv(kc,jb)
-
-!           end do
-!         end do
-!       end do
-!     end do
-!   end do
-! end do
-
-! call matout(nO*nO,nV*nV,t)
-
-end subroutine 

src/GW/static_screening_WA_so.f90

@@ -1,53 +0,0 @@
-subroutine static_screening_WA_so(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,WA)
-
-! Compute the OOVV block of the static screening W for the resonant block
-
-  implicit none
-  include 'parameters.h'
-
-! Input variables
-
-  integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
-  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,a,b,ia,jb,kc
-
-! Output variables
-
-  double precision,intent(out)  :: WA(nS,nS)
-
-! Initialize 
-
-  WA(:,:) = 0d0
-
-  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 = Omega(kc)**2 + eta**2
-            chi = chi + rho(i,j,kc)*rho(a,b,kc)*Omega(kc)/eps
-          enddo
-
-          WA(ia,jb) = WA(ia,jb) + lambda*ERI(i,b,j,a) - 2d0*lambda*chi
-
-        enddo
-      enddo
-    enddo
-  enddo
-
-end subroutine static_screening_WA_so

src/GW/static_screening_WB_so.f90

@@ -1,58 +0,0 @@
-subroutine static_screening_WB_so(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,WB)
-
-! Compute the static screening W 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)   :: Omega(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)  :: WB(nS,nS)
-
-! Initialize 
-
-  WB(:,:) = 0d0
-
-  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 = Omega(kc)**2 + eta**2
-            chi = chi + rho(i,b,kc)*rho(a,j,kc)*Omega(kc)/eps
-          enddo
-
-          WB(ia,jb) = WB(ia,jb) + lambda*ERI(i,j,b,a) - 2d0*lambda*chi
-
-        enddo
-      enddo
-    enddo
-  enddo
-
-end subroutine static_screening_WB_so

src/QuAcK/QuAcK.f90

@@ -384,7 +384,7 @@ program QuAcK
   if(doCCD) then
 
     call cpu_time(start_CC)
-    call CCD(.false.,maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF)
+    call CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF)
     call cpu_time(end_CC)
 
     t_CC = end_CC - start_CC
@@ -419,7 +419,7 @@ program QuAcK
   if(doCCSD) then
 
     call cpu_time(start_CC)
-    call CCSD(.false.,maxSCF_CC,thresh_CC,n_diis_CC,doCCSDT,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF)
+    call CCSD(maxSCF_CC,thresh_CC,n_diis_CC,doCCSDT,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF)
     call cpu_time(end_CC)
 
     t_CC = end_CC - start_CC
@@ -451,7 +451,7 @@ program QuAcK
   if(do_rCCD) then
 
     call cpu_time(start_CC)
-    call rCCD(.false.,maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF,epsHF)
+    call rCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF,epsHF)
     call cpu_time(end_CC)
 
     t_CC = end_CC - start_CC