Merge branch 'master' of https://github.com/pfloos/QuAcK

mol/GW20/Ar.xyz

@@ -1,3 +0,0 @@
-1
-Argon; atom; s
-Ar 0.0 0.0 0.0

mol/GW20/BF.xyz

@@ -1,4 +0,0 @@
-2
-Fluoroborane; experimental structure from HCP92; s
-B 0.0000 0.0000 0.0000
-F 0.0000 0.0000 1.2626

mol/GW20/BH3.xyz

@@ -1,6 +0,0 @@
-4
-Borane; experimental structure from HCP92; s
-B  0.0000  0.0000  0.0000
-H   0.0000  0.0000  1.19
-H   0.0000  1.0306 -0.595
-H   0.0000 -1.0306 -0.595

mol/GW20/BN.xyz

@@ -1,4 +0,0 @@
-2
-Boron nitride; experimental structure from HCP92; s
-B 0.0000 0.0000 0.0000
-N 0.0000 0.0000 1.281

mol/GW20/BeO.xyz

@@ -1,4 +0,0 @@
-2
-Beryllium monoxide; experimental structure from HCP92; s
-Be 0.0000 0.0000 0.0000
-O  0.0000 0.0000 1.3308

mol/GW20/CH4.xyz

@@ -1,7 +0,0 @@
-5
-Methane; experimental structure from HCP92; m
-C  0.0000  0.0000  0.0000
-H  0.6276 -0.6275  0.6276
-H -0.6276  0.6276  0.6276
-H -0.6276 -0.6276 -0.6276
-H  0.6276  0.6276 -0.6276

mol/GW20/CO.xyz

@@ -1,4 +0,0 @@
-2
-Carbon monoxide; experimental structure from HCP92; s
-C 0.0000 0.0000 0.0000
-O 0.0000 0.0000 1.283

mol/GW20/Cu2.xyz

@@ -1,4 +0,0 @@
-2
-Copper dimer; experimental structure form HCP92; s
-Cu 0.0 0.0 0.0
-Cu 0.0 0.0 2.2197

mol/GW20/F2.xyz

@@ -1,4 +0,0 @@
-2
-Fluorine; experimental structure from HCP92; s
-F 0.0000 0.0000 0.0000
-F 0.0000 0.0000 1.4119

mol/GW20/H2.xyz

@@ -1,4 +0,0 @@
-2
-Hydrogen; experimental structure from HCP92; s
-H 0.0000 0.0000 0.0000
-H 0.0000 0.0000 0.74144

mol/GW20/H2O.xyz

@@ -1,5 +0,0 @@
-3
-Water; experimental structure from HCP92; s
-O  0.0000 0.0000 0.0000
-H  0.7571 0.0000 0.5861
-H -0.7571 0.0000 0.5861

mol/GW20/HCl.xyz

@@ -1,4 +0,0 @@
-2
-Hydrogen chloride; experimental structure from HCP92; s
-H  0.0000 0.0000 0.0000
-Cl 0.0000 0.0000 1.2746

mol/GW20/HF.xyz

@@ -1,4 +0,0 @@
-2
-Hydrogen fluoride; experimental structure from HCP92; s
-H 0.0000 0.0000 0.0000
-F 0.0000 0.0000 0.9169

mol/GW20/He.xyz

@@ -1,3 +0,0 @@
-1
-Helium; atom; s
-He 0.0 0.0 0.0

mol/GW20/Li2.xyz

@@ -1,4 +0,0 @@
-2
-Lithium dimer; experimental structure from HCP92; s
-Li 0.0000 0.0000 0.0000
-Li 0.0000 0.0000 2.6729

mol/GW20/LiF.xyz

@@ -1,4 +0,0 @@
-2
-Lithium fluoride; experimental structure from HCP92; s
-Li 0.0000 0.0000 0.0000
-F  0.0000 0.0000 1.5639

mol/GW20/LiH.xyz

@@ -1,4 +0,0 @@
-2
-Lithium hydride; experimental structure from HCP92; s
-Li 0.0000 0.0000 0.0000
-H  0.0000 0.0000 1.5949

mol/GW20/N2.xyz

@@ -1,4 +0,0 @@
-2
-Nitrogen; experimental structure from HCP92; s
-N 0.0000 0.0000 0.0000
-N 0.0000 0.0000 1.0977

mol/GW20/NH3.xyz

@@ -1,6 +0,0 @@
-4
-Amonia; experimental structure from HCP92; s
-N  0.0000  0.0000  0.0000
-H  0.0000 -0.9377 -0.3816
-H  0.8121  0.4689 -0.3816
-H -0.8121  0.4689 -0.3816

mol/GW20/Ne.xyz

@@ -1,3 +0,0 @@
-1
-Neon; atom; s
-Ne 0.0 0.0 0.0

mol/GW20/O3.xyz

@@ -1,5 +0,0 @@
-3
-Ozon; experimental structure from HCP92; s
-O  0.0000 0.0000 0.0000
-O  1.0869 0.0000 0.6600
-O -1.0869 0.0000 0.6600

mol/GW20/SH2.xyz

@@ -1,5 +0,0 @@
-3
-Hydrogen sulfide; experimental structure from HCP92; s
-S  0.0000 0.0000 0.0000
-H  0.9617 0.0000 0.9268
-H -0.9617 0.0000 0.9268

mol/LiF.xyz

@@ -1,4 +1,4 @@
 2
 
 Li 0.0000  0.0000  0.0000
-F  0.0000  0.0000  1.5783
+F  0.0000  0.0000  1.58753

src/GF/GG0F2.f90

@@ -1,5 +1,5 @@
 subroutine GG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize, &
-                 nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
+                 nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
 
 ! Perform a one-shot second-order Green function calculation
 
@@ -25,7 +25,7 @@ subroutine GG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize,
   integer,intent(in)            :: nR
   integer,intent(in)            :: nS
   double precision,intent(in)   :: ENuc
-  double precision,intent(in)   :: ERHF
+  double precision,intent(in)   :: EGHF
   double precision,intent(in)   :: eHF(nBas)
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: dipole_int(nBas,nBas,ncart)
@@ -84,13 +84,13 @@ subroutine GG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize,
   ! Print results
 
   call GMP2(.false.,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,eGF,Ec)
-  call print_RG0F2(nBas,nO,eHF,SigC,eGF,Z,ENuc,ERHF,Ec)
+  call print_RG0F2(nBas,nO,eHF,SigC,eGF,Z,ENuc,EGHF,Ec)
 
-! Perform BSE2 calculation
+! Perform BSE@GF2 calculation
 
   if(dophBSE) then 
   
-    call GGF2_phBSE2(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE)
+    call GGF2_phBSE(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE)
 
     write(*,*)
     write(*,*)'-------------------------------------------------------------------------------'
@@ -101,20 +101,20 @@ subroutine GG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize,
 
   end if
 
-! Perform ppBSE2 calculation
+! Perform ppBSE@GF2 calculation
 
-! if(doppBSE) then 
+  if(doppBSE) then 
-!  
-!   call GGF2_ppBSE2(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE)
    
-!   write(*,*)
+    call GGF2_ppBSE(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE)
-!   write(*,*)'-------------------------------------------------------------------------------'
-!   write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@GG0F2 correlation energy          =',EcBSE,' au'
-!   write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@GG0F2 total energy                =',ENuc + ERHF + EcBSE,' au'
-!   write(*,*)'-------------------------------------------------------------------------------'
-!   write(*,*)
 
-! end if
+    write(*,*)
+    write(*,*)'-------------------------------------------------------------------------------'
+    write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@GG0F2 correlation energy          =',EcBSE,' au'
+    write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@GG0F2 total energy                =',ENuc + EGHF + EcBSE,' au'
+    write(*,*)'-------------------------------------------------------------------------------'
+    write(*,*)
+
+  end if
 
 ! Testing zone
 

src/GF/GGF2_phBSE.f90

@@ -1,4 +1,4 @@
-subroutine GGF2_phBSE2(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE)
+subroutine GGF2_phBSE(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE)
 
 ! Compute the second-order Bethe-Salpeter excitation energies
 
@@ -25,7 +25,6 @@ subroutine GGF2_phBSE2(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,
 ! Local variables
 
   logical                       :: dRPA = .false.
-  integer                       :: ispin
   double precision,allocatable  :: OmBSE(:)
   double precision,allocatable  :: XpY(:,:)
   double precision,allocatable  :: XmY(:,:)
@@ -43,29 +42,28 @@ subroutine GGF2_phBSE2(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,
   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))
 
-  ispin = 3
   EcBSE = 0d0
 
-               call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI,A_sta)
+               call phGLR_A(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)
+  if(.not.TDA) call phGLR_B(dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,B_sta)
 
   ! Compute static kernel
 
-               call RGF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KA_sta)
+               call GGF2_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KA_sta)
-  if(.not.TDA) call RGF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KB_sta)
+  if(.not.TDA) call GGF2_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KB_sta)
 
                A_sta(:,:) = A_sta(:,:) + KA_sta(:,:)
   if(.not.TDA) B_sta(:,:) = B_sta(:,:) + KB_sta(:,:)
 
-  ! Compute phBSE2@GF2 excitation energies
+  ! Compute phBSE@GF2 excitation energies
 
   call phLR(TDA,nS,A_sta,B_sta,EcBSE,OmBSE,XpY,XmY)
-  call print_excitation_energies('phBSE2@GGF2','spinorbital',nS,OmBSE)
+  call print_excitation_energies('phBSE@GGF2','generalized',nS,OmBSE)
   call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY)
 
   ! Compute dynamic correction for BSE via perturbation theory
 
 ! if(dBSE) &
-!   call GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY)
+!   call GGF2_phBSE_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY)
 
 end subroutine 

src/GF/GGF2_phBSE_static_kernel_A.f90

@@ -0,0 +1,94 @@
+subroutine GGF2_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,KA_sta)
+
+! Compute the resonant part of the static BSE@GF2 matrix
+
+  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)   :: eGF(nBas)
+  
+! Local variables
+
+  double precision              :: dem,num
+  integer                       :: i,j,k,l
+  integer                       :: a,b,c,d
+  integer                       :: ia,jb
+
+! Output variables
+
+  double precision,intent(out)  :: KA_sta(nS,nS)
+
+! Initialization
+
+  KA_sta(:,:) = 0d0
+
+! Second-order correlation kernel for the block A of the spinorbital manifold
+
+    jb = 0
+
+!$omp parallel do default(private) shared(KA_sta,ERI,num,dem,eGF,nO,nBas,eta,nC,nR)
+  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)
+
+          do k=nC+1,nO
+            do c=nO+1,nBas-nR
+   
+              dem = - (eGF(c) - eGF(k))
+              num = ERI(j,k,i,c)*ERI(a,c,b,k) - ERI(j,k,i,c)*ERI(a,c,k,b) & 
+                  - ERI(j,k,c,i)*ERI(a,c,b,k) + ERI(j,k,c,i)*ERI(a,c,k,b)
+
+              KA_sta(ia,jb) = KA_sta(ia,jb) - num*dem/(dem**2 + eta**2)
+          
+              dem = + (eGF(c) - eGF(k))
+              num = ERI(j,c,i,k)*ERI(a,k,b,c) - ERI(j,c,i,k)*ERI(a,k,c,b) & 
+                  - ERI(j,c,k,i)*ERI(a,k,b,c) + ERI(j,c,k,i)*ERI(a,k,c,b)
+
+              KA_sta(ia,jb) = KA_sta(ia,jb) + num*dem/(dem**2 + eta**2)
+          
+            end do
+          end do
+
+          do c=nO+1,nBas-nR
+            do d=nO+1,nBas-nR
+   
+              dem = - (eGF(c) + eGF(d))
+              num = ERI(a,j,c,d)*ERI(c,d,i,b) - ERI(a,j,c,d)*ERI(c,d,b,i) & 
+                  - ERI(a,j,d,c)*ERI(c,d,i,b) + ERI(a,j,d,c)*ERI(c,d,b,i)
+
+              KA_sta(ia,jb) = KA_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2)
+          
+            end do
+          end do
+
+          do k=nC+1,nO
+            do l=nC+1,nO
+
+              dem = - (eGF(k) + eGF(l))
+              num = ERI(a,j,k,l)*ERI(k,l,i,b) - ERI(a,j,k,l)*ERI(k,l,b,i) & 
+                  - ERI(a,j,l,k)*ERI(k,l,i,b) + ERI(a,j,l,k)*ERI(k,l,b,i)
+
+              KA_sta(ia,jb) = KA_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2)
+          
+            end do
+          end do
+
+        end do
+      end do
+
+    end do
+  end do
+!$omp end parallel do
+
+end subroutine 

src/GF/GGF2_phBSE_static_kernel_B.f90

@@ -0,0 +1,94 @@
+subroutine GGF2_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,KB_sta)
+
+! Compute the anti-resonant part of the static BSE@GF2 matrix
+
+  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)   :: eGF(nBas)
+  
+! Local variables
+
+  double precision              :: dem,num
+  integer                       :: i,j,k,l
+  integer                       :: a,b,c,d
+  integer                       :: ia,jb
+
+! Output variables
+
+  double precision,intent(out)  :: KB_sta(nS,nS)
+
+! Initialization
+
+  KB_sta(:,:) = 0d0
+
+! Second-order correlation kernel for the block A of the spinorbital manifold
+
+    jb = 0
+
+!$omp parallel do default(private) shared(KB_sta,ERI,num,dem,eGF,nO,nBas,eta,nC,nR)
+  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) 
+
+          do k=nC+1,nO
+            do c=nO+1,nBas-nR
+   
+              dem = + eGF(k) - eGF(c)
+              num = ERI(b,k,i,c)*ERI(a,c,j,k) - ERI(b,k,i,c)*ERI(a,c,k,j) & 
+                  - ERI(b,k,c,i)*ERI(a,c,j,k) + ERI(b,k,c,i)*ERI(a,c,k,j)
+
+              KB_sta(ia,jb) = KB_sta(ia,jb) - num*dem/(dem**2 + eta**2)
+          
+              dem = - eGF(c) + eGF(k)
+              num = ERI(b,c,i,k)*ERI(a,k,j,c) - ERI(b,c,i,k)*ERI(a,k,c,j) & 
+                  - ERI(b,c,k,i)*ERI(a,k,j,c) + ERI(b,c,k,i)*ERI(a,k,c,j)
+
+              KB_sta(ia,jb) = KB_sta(ia,jb) - num*dem/(dem**2 + eta**2)
+          
+            end do
+          end do
+
+          do c=nO+1,nBas-nR
+            do d=nO+1,nBas-nR
+   
+              dem = - eGF(c) - eGF(d)
+              num = ERI(a,b,c,d)*ERI(c,d,i,j) - ERI(a,b,c,d)*ERI(c,d,j,i) & 
+                  - ERI(a,b,d,c)*ERI(c,d,i,j) + ERI(a,b,d,c)*ERI(c,d,j,i)
+
+              KB_sta(ia,jb) = KB_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2)
+          
+            end do
+          end do
+
+          do k=nC+1,nO
+            do l=nC+1,nO
+
+              dem = + eGF(k) + eGF(l)
+              num = ERI(a,b,k,l)*ERI(k,l,i,j) - ERI(a,b,k,l)*ERI(k,l,j,i) & 
+                  - ERI(a,b,l,k)*ERI(k,l,i,j) + ERI(a,b,l,k)*ERI(k,l,j,i)
+
+              KB_sta(ia,jb) = KB_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2)
+          
+            end do
+          end do
+
+        end do
+      end do
+
+    end do
+  end do
+!$omp end parallel do
+
+end subroutine 

src/GF/GGF2_ppBSE.f90

@@ -0,0 +1,89 @@
+subroutine GGF2_ppBSE(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE)
+
+! Compute the Bethe-Salpeter excitation energies at the pp level
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  logical,intent(in)            :: TDA
+  logical,intent(in)            :: dBSE
+  logical,intent(in)            :: dTDA
+
+  double precision,intent(in)   :: eta
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  double precision,intent(in)   :: eGF(nBas)
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: dipole_int(nBas,nBas,ncart)
+
+! Local variables
+
+  integer                       :: nOO
+  integer                       :: nVV
+
+  double precision,allocatable  :: Bpp(:,:)
+  double precision,allocatable  :: Cpp(:,:)
+  double precision,allocatable  :: Dpp(:,:)
+
+  double precision,allocatable  :: Om1(:)
+  double precision,allocatable  :: X1(:,:)
+  double precision,allocatable  :: Y1(:,:)
+
+  double precision,allocatable  :: Om2(:)
+  double precision,allocatable  :: X2(:,:)
+  double precision,allocatable  :: Y2(:,:)
+
+  double precision,allocatable  :: KB_sta(:,:)
+  double precision,allocatable  :: KC_sta(:,:)
+  double precision,allocatable  :: KD_sta(:,:)
+
+! Output variables
+
+  double precision,intent(out)  :: EcBSE
+
+! Initialization
+
+  EcBSE = 0d0
+
+  nOO = nO*(nO-1)/2
+  nVV = nV*(nV-1)/2
+
+  allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),       &
+           Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO),       &
+           Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO), &
+           KB_sta(nVV,nOO),KC_sta(nVV,nVV),KD_sta(nOO,nOO))
+
+  ! Compute BSE excitation energies
+
+  if(.not.TDA) call GGF2_ppBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,KB_sta)
+               call GGF2_ppBSE_static_kernel_C(eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,KC_sta)
+               call GGF2_ppBSE_static_kernel_D(eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,KD_sta)
+
+  if(.not.TDA) call ppGLR_B(nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
+               call ppGLR_C(nBas,nC,nO,nV,nR,nVV,1d0,eGF,ERI,Cpp)
+               call ppGLR_D(nBas,nC,nO,nV,nR,nOO,1d0,eGF,ERI,Dpp)
+
+  Bpp(:,:) = Bpp(:,:) + KB_sta(:,:)
+  Cpp(:,:) = Cpp(:,:) + KC_sta(:,:)
+  Dpp(:,:) = Dpp(:,:) + KD_sta(:,:)
+
+  call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE)
+
+  call ppLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)
+
+  !----------------------------------------------------!
+  ! Compute the dynamical screening at the ppBSE level !
+  !----------------------------------------------------!
+
+! if(dBSE) &
+!     call GGF2_ppBSE_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nOO,nVV,eGF,ERI,dipole_int, &
+!                                          Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta)
+
+  deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta)
+
+end subroutine 

src/GF/GGF2_ppBSE_static_kernel_B.f90

@@ -0,0 +1,68 @@
+subroutine GGF2_ppBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,eGF,KB_sta)
+
+! Compute the resonant part of the dynamic BSE@GF2 matrix
+
+  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)            :: 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)   :: eGF(nBas)
+  
+! Local variables
+
+  double precision,external     :: Kronecker_delta
+  double precision              :: dem,num
+  integer                       :: i,j,a,b,m,e
+  integer                       :: ab,ij
+
+! Output variables
+
+  double precision,intent(out)  :: KB_sta(nVV,nOO)
+
+! Initialization
+
+  KB_sta(:,:) = 0d0
+
+! Second-order correlation kernel for the block B of the spinorbital manifold
+
+  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
+
+          do m=nC+1,nO
+            do e=nO+1,nBas-nR
+   
+              dem = eGF(m) - eGF(e)
+              num =       (ERI(a,m,i,e) - ERI(a,m,e,i)) * (ERI(e,b,m,j) - ERI(e,b,j,m))
+              num = num + (ERI(a,e,i,m) - ERI(a,e,m,i)) * (ERI(m,b,e,j) - ERI(m,b,j,e))
+              num = num - (ERI(b,m,i,e) - ERI(b,m,e,i)) * (ERI(e,a,m,j) - ERI(e,a,j,m))
+              num = num - (ERI(b,e,i,m) - ERI(b,e,m,i)) * (ERI(m,a,e,j) - ERI(m,a,j,e))                                                                                
+              KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2)
+              
+            end do
+          end do
+
+        end do
+      end do
+
+    end do
+  end do
+
+end subroutine 

src/GF/GGF2_ppBSE_static_kernel_C.f90

@@ -0,0 +1,69 @@
+subroutine GGF2_ppBSE_static_kernel_C(eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI,eGF,KC_sta)
+
+! Compute the resonant part of the static BSE@GF2 matrix
+
+  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)            :: nVV
+  double precision,intent(in)   :: eta
+  double precision,intent(in)   :: lambda
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: eGF(nBas)
+  
+! Local variables
+
+  double precision,external     :: Kronecker_delta
+  double precision              :: dem,num
+  integer                       :: m
+  integer                       :: a,b,c,d,e
+  integer                       :: a0,aa,ab,cd
+  
+! Output variables
+
+  double precision,intent(out)  :: KC_sta(nVV,nVV)
+
+! Initialization
+
+  KC_sta(:,:) = 0d0
+  
+! Second-order correlation kernel for the block C of the singlet manifold
+
+  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
+ 
+          do m=nC+1,nO
+            do e=nO+1,nBas-nR
+  
+              dem = eGF(m) - eGF(e)
+              num =       (ERI(a,m,c,e) - ERI(a,m,e,c)) * (ERI(e,b,m,d) - ERI(e,b,d,m))
+              num = num + (ERI(a,e,c,m) - ERI(a,e,m,c)) * (ERI(m,b,e,d) - ERI(m,b,d,e))
+              num = num - (ERI(b,m,c,e) - ERI(b,m,e,c)) * (ERI(e,a,m,d) - ERI(e,a,d,m))
+              num = num - (ERI(b,e,c,m) - ERI(b,e,m,c)) * (ERI(m,a,e,d) - ERI(m,a,d,e))
+                                                                            
+              KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2)
+        
+            end do
+          end do
+ 
+        end do
+      end do
+ 
+    end do
+  end do
+
+end subroutine 

src/GF/GGF2_ppBSE_static_kernel_D.f90

@@ -0,0 +1,68 @@
+subroutine GGF2_ppBSE_static_kernel_D(eta,nBas,nC,nO,nV,nR,nOO,lambda,ERI,eGF,KD_sta)
+
+! Compute the resonant part of the static BSE@GF2 matrix
+
+  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)            :: nOO
+  double precision,intent(in)   :: eta
+  double precision,intent(in)   :: lambda
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: eGF(nBas)
+  
+! Local variables
+
+  double precision,external     :: Kronecker_delta
+  double precision              :: dem,num
+  integer                       :: i,j,k,l,m
+  integer                       :: e
+  integer                       :: ij,kl
+
+! Output variables
+
+  double precision,intent(out)  :: KD_sta(nOO,nOO)
+
+! Initialization
+
+  KD_sta(:,:) = 0d0
+
+! Second-order correlation kernel for the block D of the spinorbital manifold
+
+  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
+
+          do m=nC+1,nO
+            do e=nO+1,nBas-nR
+   
+              dem = eGF(m) - eGF(e)
+              num =       (ERI(i,m,k,e) - ERI(i,m,e,k)) * (ERI(e,j,m,l) - ERI(e,j,l,m))
+              num = num + (ERI(i,e,k,m) - ERI(i,e,m,k)) * (ERI(m,j,e,l) - ERI(m,j,l,e))
+              num = num - (ERI(j,m,k,e) - ERI(j,m,e,k)) * (ERI(e,i,m,l) - ERI(e,i,l,m))
+              num = num - (ERI(j,e,k,m) - ERI(j,e,m,k)) * (ERI(m,i,e,l) - ERI(m,i,l,e))
+                                                                               
+              KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2)
+           end do
+          end do
+
+        end do
+      end do
+
+    end do
+  end do
+
+end subroutine 

src/GF/RG0F2.f90

@@ -87,11 +87,11 @@ subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,
   call RMP2(.false.,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec)
   call print_RG0F2(nBas,nO,eHF,SigC,eGF,Z,ENuc,ERHF,Ec)
 
-! Perform BSE2 calculation
+! Perform BSE@GF2 calculation
 
   if(dophBSE) then 
   
-    call RGF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE)
+    call RGF2_phBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE)
 
     write(*,*)
     write(*,*)'-------------------------------------------------------------------------------'
@@ -104,11 +104,11 @@ subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,
 
   end if
 
-! Perform ppBSE2 calculation
+! Perform ppBSE@GF2 calculation
 
   if(doppBSE) then 
    
-    call RGF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE)
+    call RGF2_ppBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE)
 
     EcBSE(2) = 3d0*EcBSE(2)
 

src/GF/RGF2_phBSE.f90

@@ -1,4 +1,4 @@
-subroutine RGF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE)
+subroutine RGF2_phBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE)
 
 ! Compute the second-order Bethe-Salpeter excitation energies
 
@@ -59,22 +59,22 @@ subroutine RGF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI
 
     ! Compute static kernel
 
-                 call RGF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KA_sta)
+                 call RGF2_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KA_sta)
-    if(.not.TDA) call RGF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KB_sta)
+    if(.not.TDA) call RGF2_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KB_sta)
 
                  A_sta(:,:) = A_sta(:,:) + KA_sta(:,:)
     if(.not.TDA) B_sta(:,:) = B_sta(:,:) + KB_sta(:,:)
 
-    ! Compute phBSE2@GF2 excitation energies
+    ! Compute phBSE@GF2 excitation energies
 
     call phLR(TDA,nS,A_sta,B_sta,EcBSE(ispin),OmBSE,XpY,XmY)
-    call print_excitation_energies('phBSE2@GF2','singlet',nS,OmBSE)
+    call print_excitation_energies('phBSE@GF2','singlet',nS,OmBSE)
     call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY)
 
     ! Compute dynamic correction for BSE via perturbation theory
 
     if(dBSE) &
-      call RGF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY)
+      call RGF2_phBSE_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY)
 
   end if
 
@@ -92,22 +92,22 @@ subroutine RGF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI
 
     ! Compute static kernel
 
-                 call RGF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KA_sta)
+                 call RGF2_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KA_sta)
-    if(.not.TDA) call RGF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KB_sta)
+    if(.not.TDA) call RGF2_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KB_sta)
 
                  A_sta(:,:) = A_sta(:,:) + KA_sta(:,:)
     if(.not.TDA) B_sta(:,:) = B_sta(:,:) + KB_sta(:,:)
 
-    ! Compute phBSE2@GF2 excitation energies
+    ! Compute phBSE@GF2 excitation energies
 
     call phLR(TDA,nS,A_sta,B_sta,EcBSE(ispin),OmBSE,XpY,XmY)
-    call print_excitation_energies('phBSE2@GF2','triplet',nS,OmBSE)
+    call print_excitation_energies('phBSE@GF2','triplet',nS,OmBSE)
     call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY)
 
     ! Compute dynamic correction for BSE via perturbation theory
 
     if(dBSE) &
-      call RGF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY)
+      call RGF2_phBSE_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY)
 
   end if
 

src/GF/RGF2_phBSE_dynamic_kernel_A.f90

@@ -1,6 +1,6 @@
-subroutine RGF2_phBSE2_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,OmBSE,KA_dyn,ZA_dyn)
+subroutine RGF2_phBSE_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,OmBSE,KA_dyn,ZA_dyn)
 
-! Compute the resonant part of the dynamic BSE2 matrix
+! Compute the resonant part of the dynamic BSE@GF2 matrix
 
   implicit none
   include 'parameters.h'

src/GF/RGF2_phBSE_dynamic_kernel_B.f90

@@ -1,6 +1,6 @@
-subroutine RGF2_phBSE2_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,KB_dyn)
+subroutine RGF2_phBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,KB_dyn)
 
-! Compute the anti-resonant part of the dynamic BSE2 matrix
+! Compute the anti-resonant part of the dynamic BSE@GF2 matrix
 
   implicit none
   include 'parameters.h'

src/GF/RGF2_phBSE_dynamic_perturbation.f90

@@ -1,4 +1,4 @@
-subroutine RGF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY)
+subroutine RGF2_phBSE_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY)
 
 ! Compute dynamical effects via perturbation theory for BSE
 
@@ -76,7 +76,7 @@ subroutine RGF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,E
 
     ! Resonant part of the BSE correction for dynamical TDA
 
-    call RGF2_phBSE2_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,+OmBSE(ia),KAp_dyn,ZAp_dyn)
+    call RGF2_phBSE_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,+OmBSE(ia),KAp_dyn,ZAp_dyn)
 
     if(dTDA) then 
 
@@ -87,9 +87,9 @@ subroutine RGF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,E
 
       ! Second part of the resonant and anti-resonant part of the BSE correction (frequency independent)
 
-      call RGF2_phBSE2_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,-OmBSE(ia),KAm_dyn,ZAm_dyn)
+      call RGF2_phBSE_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,-OmBSE(ia),KAm_dyn,ZAm_dyn)
 
-      call RGF2_phBSE2_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KB_dyn)
+      call RGF2_phBSE_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))  

src/GF/RGF2_phBSE_static_kernel_A.f90

@@ -1,6 +1,6 @@
-subroutine RGF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,KA_sta)
+subroutine RGF2_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,KA_sta)
 
-! Compute the resonant part of the static BSE2 matrix
+! Compute the resonant part of the static BSE@GF2 matrix
 
   implicit none
   include 'parameters.h'

src/GF/RGF2_phBSE_static_kernel_B.f90

@@ -1,6 +1,6 @@
-subroutine RGF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,KB_sta)
+subroutine RGF2_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,KB_sta)
 
-! Compute the anti-resonant part of the static BSE2 matrix
+! Compute the anti-resonant part of the static BSE@GF2 matrix
 
   implicit none
   include 'parameters.h'

src/GF/RGF2_ppBSE.f90

@@ -1,4 +1,4 @@
-subroutine RGF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE)
+subroutine RGF2_ppBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE)
 
 ! Compute the Bethe-Salpeter excitation energies at the pp level
 
@@ -74,9 +74,9 @@ subroutine RGF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,di
 
     ! Compute BSE excitation energies
 
-    if(.not.TDA) call RGF2_ppBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,KB_sta)
+    if(.not.TDA) call RGF2_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,KB_sta)
-                 call RGF2_ppBSE2_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,KC_sta)
+                 call RGF2_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,KC_sta)
-                 call RGF2_ppBSE2_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,KD_sta)
+                 call RGF2_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,KD_sta)
 
     if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
                  call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGF,ERI,Cpp)
@@ -95,7 +95,7 @@ subroutine RGF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,di
     !----------------------------------------------------!
 
     if(dBSE) &
-        call RGF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nOO,nVV,eGF,ERI,dipole_int, &
+        call RGF2_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nOO,nVV,eGF,ERI,dipole_int, &
                                              Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta)
 
     deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta)
@@ -126,9 +126,9 @@ subroutine RGF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,di
 
     ! Compute BSE excitation energies
 
-    if(.not.TDA) call RGF2_ppBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,KB_sta)
+    if(.not.TDA) call RGF2_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,KB_sta)
-                 call RGF2_ppBSE2_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,KC_sta)
+                 call RGF2_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,KC_sta)
-                 call RGF2_ppBSE2_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,KD_sta)
+                 call RGF2_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,KD_sta)
 
 
     if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
@@ -148,7 +148,7 @@ subroutine RGF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,di
     !----------------------------------------------------!
 
     if(dBSE) &
-        call RGF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nOO,nVV,eGF,ERI,dipole_int, &
+        call RGF2_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nOO,nVV,eGF,ERI,dipole_int, &
                                              Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta)
 
     deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta)

src/GF/RGF2_ppBSE2_static_kernel_B.f90

@@ -1,119 +0,0 @@
-subroutine RGF2_ppBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,eGF,KB_sta)
-
-! Compute the resonant part of the dynamic BSE2 matrix
-
-  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)            :: 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)   :: eGF(nBas)
-  
-! Local variables
-
-  double precision,external     :: Kronecker_delta
-  double precision              :: dem,num
-  integer                       :: i,j,k,a,b,c
-  integer                       :: ab,ij
-
-! Output variables
-
-  double precision,intent(out)  :: KB_sta(nVV,nOO)
-
-! Initialization
-
-  KB_sta(:,:) = 0d0
-
-! Second-order correlation kernel for the block B of the singlet manifold
-
-  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
-  
-            do k=nC+1,nO
-              do c=nO+1,nBas-nR
-     
-               dem = eGF(k) - eGF(c)
-                num = 2d0*ERI(a,k,i,c)*ERI(b,c,j,k) -     ERI(a,k,i,c)*ERI(b,c,k,j) & 
-                    -     ERI(a,k,c,i)*ERI(b,c,j,k) -     ERI(a,k,c,i)*ERI(b,c,k,j)
-
-                KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2)
-            
-                dem = eGF(k) - eGF(c)
-                num = 2d0*ERI(b,k,i,c)*ERI(a,c,j,k) -     ERI(b,k,i,c)*ERI(a,c,k,j) & 
-                    -     ERI(b,k,c,i)*ERI(a,c,j,k) -     ERI(b,k,c,i)*ERI(a,c,k,j)
-
-                KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2)
-            
-              end do
-            end do
-
-            KB_sta(ab,ij) = 2d0*lambda*KB_sta(ab,ij)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j)))
-          end do
-        end do
-
-      end do
-    end do
-
-  end if
-
-! Second-order correlation kernel for the block B of the triplet manifold
-
-  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
-  
-            do k=nC+1,nO
-              do c=nO+1,nBas-nR
-     
-                dem = eGF(k) - eGF(c)
-                num = 2d0*ERI(a,k,i,c)*ERI(b,c,j,k) -     ERI(a,k,i,c)*ERI(b,c,k,j) & 
-                    -     ERI(a,k,c,i)*ERI(b,c,j,k) +     ERI(a,k,c,i)*ERI(b,c,k,j)
-
-                KB_sta(ab,ij) = KB_sta(ab,ij) + 2d0*num*dem/(dem**2 + eta**2)
-            
-                dem = eGF(k) - eGF(c)
-                num = 2d0*ERI(b,k,i,c)*ERI(a,c,j,k) -     ERI(b,k,i,c)*ERI(a,c,k,j) & 
-                    -     ERI(b,k,c,i)*ERI(a,c,j,k) +     ERI(b,k,c,i)*ERI(a,c,k,j)
-
-                KB_sta(ab,ij) = KB_sta(ab,ij) - 2d0*num*dem/(dem**2 + eta**2)
-            
-              end do
-            end do
-
-          end do
-        end do
-
-      end do
-    end do
-
-  end if
-
-end subroutine 

src/GF/RGF2_ppBSE2_static_kernel_C.f90

@@ -1,120 +0,0 @@
-subroutine RGF2_ppBSE2_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI,eGF,KC_sta)
-
-! Compute the resonant part of the static BSE2 matrix
-
-  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)            :: nVV
-  double precision,intent(in)   :: eta
-  double precision,intent(in)   :: lambda
-  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
-  double precision,intent(in)   :: eGF(nBas)
-  
-! Local variables
-
-  double precision,external     :: Kronecker_delta
-  double precision              :: dem,num
-  integer                       :: m
-  integer                       :: a,b,c,d,e
-  integer                       :: ab,cd
-
-! Output variables
-
-  double precision,intent(out)  :: KC_sta(nVV,nVV)
-
-! Initialization
-
-  KC_sta(:,:) = 0d0
-
-! Second-order correlation kernel for the block C of the singlet manifold
-
-  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
-  
-            do m=nC+1,nO
-              do e=nO+1,nBas-nR
-     
-                dem = eGF(m) - eGF(e)
-                num = 2d0*ERI(a,m,c,e)*ERI(b,e,d,m) -     ERI(a,m,c,e)*ERI(b,e,m,d) & 
-                    -     ERI(a,m,e,c)*ERI(b,e,d,m) -     ERI(a,m,e,c)*ERI(b,e,m,d)
-
-                KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2)
-            
-                dem = eGF(m) - eGF(e)
-                num = 2d0*ERI(b,m,c,e)*ERI(a,e,d,m) -     ERI(b,m,c,e)*ERI(a,e,m,d) & 
-                    -     ERI(b,m,e,c)*ERI(a,e,d,m) -     ERI(b,m,e,c)*ERI(a,e,m,d)
-
-                KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2)
-            
-              end do
-           end do
-
-           KC_sta(ab,cd) = 2d0*lambda*KC_sta(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d)))
-            
-          end do
-        end do
-
-      end do
-    end do
-
-  end if
-
-! Second-order correlation kernel for the block C of the triplet manifold
-
-  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
-  
-            do m=nC+1,nO
-              do e=nO+1,nBas-nR
-     
-                dem = eGF(m) - eGF(e)
-                num = 2d0*ERI(a,m,c,e)*ERI(b,e,d,m) -     ERI(a,m,c,e)*ERI(b,e,m,d) &
-                    -     ERI(a,m,e,c)*ERI(b,e,d,m) +     ERI(a,m,e,c)*ERI(b,e,m,d)
-
-                KC_sta(ab,cd) = KC_sta(ab,cd) + 2d0*num*dem/(dem**2 + eta**2)
-            
-                dem = eGF(m) - eGF(e)
-                num = 2d0*ERI(b,m,c,e)*ERI(a,e,d,m) -     ERI(b,m,c,e)*ERI(a,e,m,d) &
-                    -     ERI(b,m,e,c)*ERI(a,e,d,m) +     ERI(b,m,e,c)*ERI(a,e,m,d)
-
-                KC_sta(ab,cd) = KC_sta(ab,cd) - 2d0*num*dem/(dem**2 + eta**2)
-            
-              end do
-            end do
-
-          end do
-        end do
-
-      end do
-    end do
-
-  end if
-
-end subroutine 

src/GF/RGF2_ppBSE2_static_kernel_D.f90

@@ -1,119 +0,0 @@
-subroutine RGF2_ppBSE2_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,lambda,ERI,eGF,KD_sta)
-
-! Compute the resonant part of the static BSE2 matrix
-
-  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)            :: nOO
-  double precision,intent(in)   :: eta
-  double precision,intent(in)   :: lambda
-  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
-  double precision,intent(in)   :: eGF(nBas)
-  
-! Local variables
-
-  double precision,external     :: Kronecker_delta
-  double precision              :: dem,num
-  integer                       :: i,j,k,l,m
-  integer                       :: e
-  integer                       :: ij,kl
-
-! Output variables
-
-  double precision,intent(out)  :: KD_sta(nOO,nOO)
-
-! Initialization
-
-  KD_sta(:,:) = 0d0
-
-! Second-order correlation kernel for the block D of the singlet manifold
-
-  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
-  
-            do m=nC+1,nO
-              do e=nO+1,nBas-nR
-     
-                dem = - eGF(e) + eGF(m)
-                num = 2d0*ERI(i,e,k,m)*ERI(j,m,l,e) -     ERI(i,e,k,m)*ERI(j,m,e,l) & 
-                    -     ERI(i,e,m,k)*ERI(j,m,l,e) -     ERI(i,e,m,k)*ERI(j,m,e,l)
-
-                KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2)
-            
-                dem = - eGF(e) + eGF(m)
-                num = 2d0*ERI(j,e,k,m)*ERI(i,m,l,e) -     ERI(j,e,k,m)*ERI(i,m,e,l) & 
-                    -     ERI(j,e,m,k)*ERI(i,m,l,e) -     ERI(j,e,m,k)*ERI(i,m,e,l)
-
-                KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2)
-            
-              end do
-            end do
-
-            KD_sta(ij,kl) = 2d0*lambda*KD_sta(ij,kl)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l)))
-          end do
-        end do
-
-      end do
-    end do
-
-  end if
-
-! Second-order correlation kernel for the block D of the triplet manifold
-
-  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
-  
-            do m=nC+1,nO
-              do e=nO+1,nBas-nR
-     
-                dem = - eGF(e) + eGF(m)
-                num = 2d0*ERI(i,e,k,m)*ERI(j,m,l,e) -     ERI(i,e,k,m)*ERI(j,m,e,l) &
-                    -     ERI(i,e,m,k)*ERI(j,m,l,e) +     ERI(i,e,m,k)*ERI(j,m,e,l)
-
-                KD_sta(ij,kl) = KD_sta(ij,kl) + 2d0*num*dem/(dem**2 + eta**2)
-
-                dem = - eGF(e) + eGF(m)
-                num = 2d0*ERI(j,e,k,m)*ERI(i,m,l,e) -     ERI(j,e,k,m)*ERI(i,m,e,l) &
-                    -     ERI(j,e,m,k)*ERI(i,m,l,e) +     ERI(j,e,m,k)*ERI(i,m,e,l)
-
-                KD_sta(ij,kl) = KD_sta(ij,kl) - 2d0*num*dem/(dem**2 + eta**2)
-
-              end do
-            end do
-
-          end do
-        end do
-
-      end do
-    end do
-
-  end if
-
-end subroutine 

src/GF/RGF2_ppBSE_dynamic_kernel_B.f90

@@ -1,6 +1,6 @@
-subroutine RGF2_ppBSE2_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,eGF,KB_dyn)
+subroutine RGF2_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,eGF,KB_dyn)
 
-! Compute the resonant part of the dynamic BSE2 matrix
+! Compute the resonant part of the dynamic BSE@GF2 matrix
 
   implicit none
   include 'parameters.h'

src/GF/RGF2_ppBSE_dynamic_kernel_C.f90

@@ -1,6 +1,6 @@
-subroutine RGF2_ppBSE2_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI,eGF,OmBSE,KC_dyn,ZC_dyn)
+subroutine RGF2_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI,eGF,OmBSE,KC_dyn,ZC_dyn)
 
-! Compute the resonant part of the dynamic BSE2 matrix
+! Compute the resonant part of the dynamic BSE@GF2 matrix
 
   implicit none
   include 'parameters.h'

src/GF/RGF2_ppBSE_dynamic_kernel_D.f90

@@ -1,6 +1,6 @@
-subroutine RGF2_ppBSE2_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,lambda,ERI,eGF,OmBSE,KD_dyn,ZD_dyn)
+subroutine RGF2_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,lambda,ERI,eGF,OmBSE,KD_dyn,ZD_dyn)
 
-! Compute the resonant part of the dynamic BSE2 matrix
+! Compute the resonant part of the dynamic BSE@GF2 matrix
 
   implicit none
   include 'parameters.h'

src/GF/RGF2_ppBSE_dynamic_perturbation.f90

@@ -1,4 +1,4 @@
-subroutine RGF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nOO,nVV,eGF,ERI,dipole_int, & 
+subroutine RGF2_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nOO,nVV,eGF,ERI,dipole_int, & 
                                            Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta)
 
 ! Compute dynamical effects via perturbation theory for BSE
@@ -63,7 +63,7 @@ subroutine RGF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nOO,
   end if
 
   write(*,*) '---------------------------------------------------------------------------------------------------'
-  write(*,*) ' First-order dynamical correction to static ppBSE2 double electron attachment energies             '
+  write(*,*) ' First-order dynamical correction to static ppBSE@GF2 double electron attachment energies             '
   write(*,*) '---------------------------------------------------------------------------------------------------'
   write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)'
   write(*,*) '---------------------------------------------------------------------------------------------------'
@@ -72,16 +72,16 @@ subroutine RGF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nOO,
 
     if(dTDA) then
 
-      call RGF2_ppBSE2_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,Om1(ab),KC_dyn,ZC_dyn)
+      call RGF2_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,Om1(ab),KC_dyn,ZC_dyn)
   
       Z1_dyn(ab)  = dot_product(X1(:,ab),matmul(ZC_dyn,X1(:,ab)))
       Om1_dyn(ab) = dot_product(X1(:,ab),matmul(KC_dyn - KC_sta,X1(:,ab)))
 
     else
 
-      call RGF2_ppBSE2_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,KB_dyn)
+      call RGF2_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,KB_dyn)
-      call RGF2_ppBSE2_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,Om1(ab),KC_dyn,ZC_dyn)
+      call RGF2_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,Om1(ab),KC_dyn,ZC_dyn)
-      call RGF2_ppBSE2_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,Om1(ab),KD_dyn,ZD_dyn)
+      call RGF2_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,Om1(ab),KD_dyn,ZD_dyn)
 
       Z1_dyn(ab)  = dot_product(X1(:,ab),matmul(ZC_dyn,X1(:,ab))) &
                   + dot_product(Y1(:,ab),matmul(ZD_dyn,Y1(:,ab)))
@@ -104,7 +104,7 @@ subroutine RGF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nOO,
   write(*,*) 
 
   write(*,*) '---------------------------------------------------------------------------------------------------'
-  write(*,*) ' First-order dynamical correction to static ppBSE2 double electron detachment energies             '
+  write(*,*) ' First-order dynamical correction to static ppBSE@GF2 double electron detachment energies             '
   write(*,*) '---------------------------------------------------------------------------------------------------'
   write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)'
   write(*,*) '---------------------------------------------------------------------------------------------------'
@@ -115,16 +115,16 @@ subroutine RGF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nOO,
 
     if(dTDA) then
 
-      call RGF2_ppBSE2_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,-Om2(ij),KD_dyn,ZD_dyn)
+      call RGF2_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,-Om2(ij),KD_dyn,ZD_dyn)
 
       Z2_dyn(kl)  = dot_product(Y2(:,ij),matmul(ZD_dyn,Y2(:,ij)))
       Om2_dyn(kl) = dot_product(Y2(:,ij),matmul(KD_dyn - KD_sta,Y2(:,ij)))
 
     else
 
-      call RGF2_ppBSE2_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,KB_dyn)
+      call RGF2_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,KB_dyn)
-      call RGF2_ppBSE2_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,-Om2(ij),KC_dyn,ZC_dyn)
+      call RGF2_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,-Om2(ij),KC_dyn,ZC_dyn)
-      call RGF2_ppBSE2_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,-Om2(ij),KD_dyn,ZD_dyn)
+      call RGF2_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,-Om2(ij),KD_dyn,ZD_dyn)
 
       Z2_dyn(kl)  = dot_product(X2(:,ij),matmul(ZC_dyn,X2(:,ij))) &
                   + dot_product(Y2(:,ij),matmul(ZD_dyn,Y2(:,ij)))

src/GF/RGF2_ppBSE_static_kernel_B.f90

@@ -0,0 +1,173 @@
+subroutine RGF2_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,eGF,KB_sta)
+
+! Compute the resonant part of the dynamic BSE@GF2 matrix
+
+  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)            :: 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)   :: eGF(nBas)
+  
+! Local variables
+
+  double precision,external     :: Kronecker_delta
+  double precision              :: dem,num
+  integer                       :: i,j,a,b,m,e
+  integer                       :: ab,ij
+
+! Output variables
+
+  double precision,intent(out)  :: KB_sta(nVV,nOO)
+
+! Initialization
+
+  KB_sta(:,:) = 0d0
+
+! Second-order correlation kernel for the block B of the singlet manifold
+
+  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
+  
+            do m=nC+1,nO
+              do e=nO+1,nBas-nR
+     
+                dem = eGF(m) - eGF(e)
+                num = 2d0*ERI(a,m,i,e)*ERI(e,b,m,j) - ERI(a,m,i,e)*ERI(e,b,j,m)  &
+                     -    ERI(a,m,e,i)*ERI(e,b,m,j) - ERI(a,m,e,i)*ERI(e,b,j,m)
+                                                                                 
+                KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2)
+
+                num = 2d0*ERI(a,e,i,m)*ERI(m,b,e,j) - ERI(a,e,i,m)*ERI(m,b,j,e)  &
+                     -    ERI(a,e,m,i)*ERI(m,b,e,j) - ERI(a,e,m,i)*ERI(m,b,j,e)
+                                                                                 
+                KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2) 
+                                                                                 
+                num = 2d0*ERI(b,m,i,e)*ERI(e,a,m,j) - ERI(b,m,i,e)*ERI(e,a,j,m)  &
+                     -    ERI(b,m,e,i)*ERI(e,a,m,j) - ERI(b,m,e,i)*ERI(e,a,j,m)
+                                                                                 
+                KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2)
+
+                num = 2d0*ERI(b,e,i,m)*ERI(m,a,e,j) - ERI(b,e,i,m)*ERI(m,a,j,e)  &
+                     -    ERI(b,e,m,i)*ERI(m,a,e,j) - ERI(b,e,m,i)*ERI(m,a,j,e)
+                                                                                 
+                KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2) 
+            
+              end do
+            end do
+
+            KB_sta(ab,ij) = lambda*KB_sta(ab,ij)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j)))
+
+          end do
+        end do
+
+      end do
+    end do
+
+  end if
+
+! Second-order correlation kernel for the block B of the triplet manifold
+
+  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
+  
+            do m=nC+1,nO
+              do e=nO+1,nBas-nR
+     
+                dem = eGF(m) - eGF(e)
+                num = 2d0*ERI(a,m,i,e)*ERI(e,b,m,j) - ERI(a,m,i,e)*ERI(e,b,j,m)  &
+                     -    ERI(a,m,e,i)*ERI(e,b,m,j) + ERI(a,m,e,i)*ERI(e,b,j,m)
+                                                                                 
+                KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2)
+
+                num = 2d0*ERI(a,e,i,m)*ERI(m,b,e,j) - ERI(a,e,i,m)*ERI(m,b,j,e)  &
+                     -    ERI(a,e,m,i)*ERI(m,b,e,j) + ERI(a,e,m,i)*ERI(m,b,j,e)
+                                                                                 
+                KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2) 
+                                                                                 
+                num = 2d0*ERI(b,m,i,e)*ERI(e,a,m,j) - ERI(b,m,i,e)*ERI(e,a,j,m)  &
+                     -    ERI(b,m,e,i)*ERI(e,a,m,j) + ERI(b,m,e,i)*ERI(e,a,j,m)
+                                                                                 
+                KB_sta(ab,ij) = KB_sta(ab,ij) - num*dem/(dem**2 + eta**2)
+
+                num = 2d0*ERI(b,e,i,m)*ERI(m,a,e,j) - ERI(b,e,i,m)*ERI(m,a,j,e)  &
+                     -    ERI(b,e,m,i)*ERI(m,a,e,j) + ERI(b,e,m,i)*ERI(m,a,j,e)
+                                                                                 
+                KB_sta(ab,ij) = KB_sta(ab,ij) - num*dem/(dem**2 + eta**2) 
+            
+              end do
+            end do
+
+          end do
+        end do
+
+      end do
+    end do
+
+  end if
+
+! Second-order correlation kernel for the block B of the spinorbital manifold
+
+  if(ispin == 4) 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
+  
+            do m=nC+1,nO
+              do e=nO+1,nBas-nR
+     
+                dem = eGF(m) - eGF(e)
+                num =       (ERI(a,m,i,e) - ERI(a,m,e,i)) * (ERI(e,b,m,j) - ERI(e,b,j,m))
+                num = num + (ERI(a,e,i,m) - ERI(a,e,m,i)) * (ERI(m,b,e,j) - ERI(m,b,j,e))
+                num = num - (ERI(b,m,i,e) - ERI(b,m,e,i)) * (ERI(e,a,m,j) - ERI(e,a,j,m))
+                num = num - (ERI(b,e,i,m) - ERI(b,e,m,i)) * (ERI(m,a,e,j) - ERI(m,a,j,e))                                                                                
+                KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2)
+                
+              end do
+            end do
+
+          end do
+        end do
+
+      end do
+    end do
+
+  end if
+
+end subroutine 

src/GF/RGF2_ppBSE_static_kernel_C.f90

@@ -0,0 +1,304 @@
+subroutine RGF2_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI,eGF,KC_sta)
+
+! Compute the resonant part of the static BSE@GF2 matrix
+
+  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)            :: nVV
+  double precision,intent(in)   :: eta
+  double precision,intent(in)   :: lambda
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: eGF(nBas)
+  
+! Local variables
+
+  double precision,external     :: Kronecker_delta
+  double precision              :: dem,num
+  integer                       :: m
+  integer                       :: a,b,c,d,e
+  integer                       :: a0,aa,ab,cd
+  
+! Output variables
+
+  double precision,intent(out)  :: KC_sta(nVV,nVV)
+
+! Initialization
+
+  KC_sta(:,:) = 0d0
+  
+! Second-order correlation kernel for the block C of the singlet manifold
+
+! --- --- ---
+! OpenMP implementation
+! --- --- ---
+
+if(ispin == 1) then
+
+  a0 = nBas - nR - nO
+  !$OMP PARALLEL DEFAULT(NONE)                              &
+  !$OMP          PRIVATE(a, b, aa, ab, c, d, cd, m, e, num, dem) &
+  !$OMP          SHARED(nO, nBas, nR, nC, a0, ERI, KC_sta, lambda, eGF, eta)
+  !$OMP DO
+  do a=nO+1,nBas-nR
+    aa = a0 * (a - nO - 1) - (a - nO - 1) * (a - nO) / 2 - nO
+    do b=a,nBas-nR
+      ab = aa + b
+
+      cd = 0
+      do c=nO+1,nBas-nR
+        do d=c,nBas-nR
+          cd = cd + 1
+
+          do m=nC+1,nO
+            do e=nO+1,nBas-nR
+   
+                dem = eGF(m) - eGF(e)
+                num = 2d0*ERI(a,m,c,e)*ERI(e,b,m,d) - ERI(a,m,c,e)*ERI(e,b,d,m)  &
+                     -    ERI(a,m,e,c)*ERI(e,b,m,d) - ERI(a,m,e,c)*ERI(e,b,d,m)
+                                                                                 
+                KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2)
+
+                num = 2d0*ERI(a,e,c,m)*ERI(m,b,e,d) - ERI(a,e,c,m)*ERI(m,b,d,e)  &
+                     -    ERI(a,e,m,c)*ERI(m,b,e,d) - ERI(a,e,m,c)*ERI(m,b,d,e)
+                                                                                 
+                KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2) 
+                                                                                 
+                num = 2d0*ERI(b,m,c,e)*ERI(e,a,m,d) - ERI(b,m,c,e)*ERI(e,a,d,m)  &
+                     -    ERI(b,m,e,c)*ERI(e,a,m,d) - ERI(b,m,e,c)*ERI(e,a,d,m)
+                                                                                 
+                KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2)
+
+                num = 2d0*ERI(b,e,c,m)*ERI(m,a,e,d) - ERI(b,e,c,m)*ERI(m,a,d,e)  &
+                     -    ERI(b,e,m,c)*ERI(m,a,e,d) - ERI(b,e,m,c)*ERI(m,a,d,e)
+                                                                                 
+                KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2) 
+          
+            end do
+          end do
+
+          KC_sta(ab,cd) = lambda*KC_sta(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d)))
+          
+         end do
+       end do
+
+    end do
+  end do
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+end if
+
+!  --- --- ---
+!  Naive implementation
+!  --- --- ---
+
+!   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
+   
+!             do m=nC+1,nO
+!               do e=nO+1,nBas-nR
+   
+!                 dem = eGF(m) - eGF(e)
+!                 num = 2d0*ERI(a,m,c,e)*ERI(e,b,m,d) - ERI(a,m,c,e)*ERI(e,b,d,m)  &
+!                      -    ERI(a,m,e,c)*ERI(e,b,m,d) - ERI(a,m,e,c)*ERI(e,b,d,m)
+                                                                                 
+!                 KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2)
+   
+!                 num = 2d0*ERI(a,e,c,m)*ERI(m,b,e,d) - ERI(a,e,c,m)*ERI(m,b,d,e)  &
+!                      -    ERI(a,e,m,c)*ERI(m,b,e,d) - ERI(a,e,m,c)*ERI(m,b,d,e)
+                                                                                 
+!                 KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2) 
+                                                                                 
+!                 num = 2d0*ERI(b,m,c,e)*ERI(e,a,m,d) - ERI(b,m,c,e)*ERI(e,a,d,m)  &
+!                      -    ERI(b,m,e,c)*ERI(e,a,m,d) - ERI(b,m,e,c)*ERI(e,a,d,m)
+                                                                                 
+!                 KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2)
+   
+!                 num = 2d0*ERI(b,e,c,m)*ERI(m,a,e,d) - ERI(b,e,c,m)*ERI(m,a,d,e)  &
+!                      -    ERI(b,e,m,c)*ERI(m,a,e,d) - ERI(b,e,m,c)*ERI(m,a,d,e)
+                                                                                 
+!                 KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2) 
+          
+!               end do
+!            end do
+   
+!            KC_sta(ab,cd) = lambda*KC_sta(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d)))
+          
+!           end do
+!         end do
+   
+!       end do
+!     end do
+   
+!   end if
+
+! Second-order correlation kernel for the block C of the triplet manifold
+
+!  --- --- ---
+!  OpenMP implementation
+!  --- --- ---
+  
+if(ispin == 2) then
+
+  a0 = nBas - nR - nO - 1
+  !$OMP PARALLEL DEFAULT(NONE)                              &
+  !$OMP          PRIVATE(a, b, aa, ab, c, d, cd, m, e, num, dem) &
+  !$OMP          SHARED(nO, nBas, nR, nC, a0, ERI, KC_sta, eGF, eta)
+  !$OMP DO
+  do a = nO+1, nBas-nR
+    aa = a0 * (a - nO - 1) - (a - nO - 1) * (a - nO) / 2 - nO - 1
+    do b = a+1, nBas-nR
+      ab = aa + b
+
+      cd = 0
+      do c=nO+1,nBas-nR
+        do d=c+1,nBas-nR
+          cd = cd + 1
+
+          do m=nC+1,nO
+            do e=nO+1,nBas-nR
+   
+                dem = eGF(m) - eGF(e)
+                num = 2d0*ERI(a,m,c,e)*ERI(e,b,m,d) - ERI(a,m,c,e)*ERI(e,b,d,m)  &
+                     -    ERI(a,m,e,c)*ERI(e,b,m,d) + ERI(a,m,e,c)*ERI(e,b,d,m)
+                                                                                 
+                KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2)
+
+                num = 2d0*ERI(a,e,c,m)*ERI(m,b,e,d) - ERI(a,e,c,m)*ERI(m,b,d,e)  &
+                     -    ERI(a,e,m,c)*ERI(m,b,e,d) + ERI(a,e,m,c)*ERI(m,b,d,e)
+                                                                                 
+                KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2) 
+                                                                                 
+                num = 2d0*ERI(b,m,c,e)*ERI(e,a,m,d) - ERI(b,m,c,e)*ERI(e,a,d,m)  &
+                     -    ERI(b,m,e,c)*ERI(e,a,m,d) + ERI(b,m,e,c)*ERI(e,a,d,m)
+                                                                                 
+                KC_sta(ab,cd) = KC_sta(ab,cd) - num*dem/(dem**2 + eta**2)
+
+                num = 2d0*ERI(b,e,c,m)*ERI(m,a,e,d) - ERI(b,e,c,m)*ERI(m,a,d,e)  &
+                     -    ERI(b,e,m,c)*ERI(m,a,e,d) + ERI(b,e,m,c)*ERI(m,a,d,e)
+                                                                                 
+                KC_sta(ab,cd) = KC_sta(ab,cd) - num*dem/(dem**2 + eta**2) 
+          
+            end do
+          end do
+
+        end do
+      end do
+
+    end do
+  end do
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+end if
+
+!  --- --- ---
+!  Naive implementation
+!  --- --- ---
+
+! 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
+  
+!            do m=nC+1,nO
+!              do e=nO+1,nBas-nR
+  
+!                 dem = eGF(m) - eGF(e)
+!                 num = 2d0*ERI(a,m,c,e)*ERI(e,b,m,d) - ERI(a,m,c,e)*ERI(e,b,d,m)  &
+!                      -    ERI(a,m,e,c)*ERI(e,b,m,d) + ERI(a,m,e,c)*ERI(e,b,d,m)
+                                                                                 
+!                 KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2)
+  
+!                 num = 2d0*ERI(a,e,c,m)*ERI(m,b,e,d) - ERI(a,e,c,m)*ERI(m,b,d,e)  &
+!                      -    ERI(a,e,m,c)*ERI(m,b,e,d) + ERI(a,e,m,c)*ERI(m,b,d,e)
+                                                                                 
+!                 KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2) 
+                                                                                 
+!                 num = 2d0*ERI(b,m,c,e)*ERI(e,a,m,d) - ERI(b,m,c,e)*ERI(e,a,d,m)  &
+!                      -    ERI(b,m,e,c)*ERI(e,a,m,d) + ERI(b,m,e,c)*ERI(e,a,d,m)
+                                                                                 
+!                 KC_sta(ab,cd) = KC_sta(ab,cd) - num*dem/(dem**2 + eta**2)
+  
+!                 num = 2d0*ERI(b,e,c,m)*ERI(m,a,e,d) - ERI(b,e,c,m)*ERI(m,a,d,e)  &
+!                      -    ERI(b,e,m,c)*ERI(m,a,e,d) + ERI(b,e,m,c)*ERI(m,a,d,e)
+                                                                                 
+!                 KC_sta(ab,cd) = KC_sta(ab,cd) - num*dem/(dem**2 + eta**2) 
+         
+!              end do
+!            end do
+  
+!          end do
+!        end do
+  
+!      end do
+!    end do
+  
+!   end if
+  
+  if(ispin == 4) 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
+ 
+            do m=nC+1,nO
+              do e=nO+1,nBas-nR
+   
+                dem = eGF(m) - eGF(e)
+                num =       (ERI(a,m,c,e) - ERI(a,m,e,c)) * (ERI(e,b,m,d) - ERI(e,b,d,m))
+                num = num + (ERI(a,e,c,m) - ERI(a,e,m,c)) * (ERI(m,b,e,d) - ERI(m,b,d,e))
+                num = num - (ERI(b,m,c,e) - ERI(b,m,e,c)) * (ERI(e,a,m,d) - ERI(e,a,d,m))
+                num = num - (ERI(b,e,c,m) - ERI(b,e,m,c)) * (ERI(m,a,e,d) - ERI(m,a,d,e))
+                                                                              
+                KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2)
+          
+              end do
+            end do
+ 
+          end do
+        end do
+ 
+      end do
+    end do
+
+  end if
+  
+! Second-order correlation kernel for the block C of the spinorbital manifold
+
+  
+
+! deallocate(Om_tmp)
+
+end subroutine 

src/GF/RGF2_ppBSE_static_kernel_D.f90

@@ -0,0 +1,173 @@
+subroutine RGF2_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,lambda,ERI,eGF,KD_sta)
+
+! Compute the resonant part of the static BSE@GF2 matrix
+
+  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)            :: nOO
+  double precision,intent(in)   :: eta
+  double precision,intent(in)   :: lambda
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: eGF(nBas)
+  
+! Local variables
+
+  double precision,external     :: Kronecker_delta
+  double precision              :: dem,num
+  integer                       :: i,j,k,l,m
+  integer                       :: e
+  integer                       :: ij,kl
+
+! Output variables
+
+  double precision,intent(out)  :: KD_sta(nOO,nOO)
+
+! Initialization
+
+  KD_sta(:,:) = 0d0
+
+! Second-order correlation kernel for the block D of the singlet manifold
+
+  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
+  
+            do m=nC+1,nO
+              do e=nO+1,nBas-nR
+     
+                dem = eGF(m) - eGF(e)
+                num = 2d0*ERI(i,m,k,e)*ERI(e,j,m,l) - ERI(i,m,k,e)*ERI(e,j,l,m)  &
+                     -    ERI(i,m,e,k)*ERI(e,j,m,l) - ERI(i,m,e,k)*ERI(e,j,l,m)
+                                                                                 
+                KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2)
+
+                num = 2d0*ERI(i,e,k,m)*ERI(m,j,e,l) - ERI(i,e,k,m)*ERI(m,j,l,e)  &
+                     -    ERI(i,e,m,k)*ERI(m,j,e,l) - ERI(i,e,m,k)*ERI(m,j,l,e)
+                                                                                 
+                KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2) 
+                                                                                 
+                num = 2d0*ERI(j,m,k,e)*ERI(e,i,m,l) - ERI(j,m,k,e)*ERI(e,i,l,m)  &
+                     -    ERI(j,m,e,k)*ERI(e,i,m,l) - ERI(j,m,e,k)*ERI(e,i,l,m)
+                                                                                 
+                KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2)
+
+                num = 2d0*ERI(j,e,k,m)*ERI(m,i,e,l) - ERI(j,e,k,m)*ERI(m,i,l,e)  &
+                     -    ERI(j,e,m,k)*ERI(m,i,e,l) - ERI(j,e,m,k)*ERI(m,i,l,e)
+                                                                                 
+                KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2) 
+            
+              end do
+            end do
+
+            KD_sta(ij,kl) = lambda*KD_sta(ij,kl)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l)))
+
+          end do
+        end do
+
+      end do
+    end do
+
+  end if
+
+! Second-order correlation kernel for the block D of the triplet manifold
+
+  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
+  
+            do m=nC+1,nO
+              do e=nO+1,nBas-nR
+     
+                dem = eGF(m) - eGF(e)
+                num = 2d0*ERI(i,m,k,e)*ERI(e,j,m,l) - ERI(i,m,k,e)*ERI(e,j,l,m)  &
+                     -    ERI(i,m,e,k)*ERI(e,j,m,l) + ERI(i,m,e,k)*ERI(e,j,l,m)
+                                                                                 
+                KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2)
+
+                num = 2d0*ERI(i,e,k,m)*ERI(m,j,e,l) - ERI(i,e,k,m)*ERI(m,j,l,e)  &
+                     -    ERI(i,e,m,k)*ERI(m,j,e,l) + ERI(i,e,m,k)*ERI(m,j,l,e)
+                                                                                 
+                KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2) 
+                                                                                 
+                num = 2d0*ERI(j,m,k,e)*ERI(e,i,m,l) - ERI(j,m,k,e)*ERI(e,i,l,m)  &
+                     -    ERI(j,m,e,k)*ERI(e,i,m,l) + ERI(j,m,e,k)*ERI(e,i,l,m)
+                                                                                 
+                KD_sta(ij,kl) = KD_sta(ij,kl) - num*dem/(dem**2 + eta**2)
+
+                num = 2d0*ERI(j,e,k,m)*ERI(m,i,e,l) - ERI(j,e,k,m)*ERI(m,i,l,e)  &
+                     -    ERI(j,e,m,k)*ERI(m,i,e,l) + ERI(j,e,m,k)*ERI(m,i,l,e)
+                                                                                 
+                KD_sta(ij,kl) = KD_sta(ij,kl) - num*dem/(dem**2 + eta**2) 
+
+              end do
+            end do
+
+          end do
+        end do
+
+      end do
+    end do
+
+  end if
+
+! Second-order correlation kernel for the block D of the spinorbital manifold
+
+  if(ispin == 4) 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
+  
+            do m=nC+1,nO
+              do e=nO+1,nBas-nR
+     
+                dem = eGF(m) - eGF(e)
+                num =       (ERI(i,m,k,e) - ERI(i,m,e,k)) * (ERI(e,j,m,l) - ERI(e,j,l,m))
+                num = num + (ERI(i,e,k,m) - ERI(i,e,m,k)) * (ERI(m,j,e,l) - ERI(m,j,l,e))
+                num = num - (ERI(j,m,k,e) - ERI(j,m,e,k)) * (ERI(e,i,m,l) - ERI(e,i,l,m))
+                num = num - (ERI(j,e,k,m) - ERI(j,e,m,k)) * (ERI(m,i,e,l) - ERI(m,i,l,e))
+                                                                                 
+                KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2)
+             end do
+            end do
+
+          end do
+        end do
+
+      end do
+    end do
+
+  end if
+
+end subroutine 

src/GF/UG0F2.f90

@@ -124,7 +124,7 @@ subroutine UG0F2(dotest,BSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,linearize,eta
 
   if(BSE) then
 
-    print*,'!!! BSE2 NYI for UG0F2 !!!'
+    print*,'!!! BSE@GF2 NYI for UG0F2 !!!'
 
   end if
 

src/GF/evGGF2.f90

@@ -142,11 +142,11 @@ subroutine evGGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, &
 
   end if
 
-! Perform BSE2 calculation
+! Perform BSE@GF2 calculation
 
   if(dophBSE) then 
   
-    call GGF2_phBSE2(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE)
+    call GGF2_phBSE(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE)
 
     write(*,*)
     write(*,*)'-------------------------------------------------------------------------------'
@@ -157,11 +157,11 @@ subroutine evGGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, &
 
   end if
 
-! Perform ppBSE2 calculation
+! Perform ppBSE@GF2 calculation
 
 ! if(doppBSE) then
 
-!   call GGF2_ppBSE2(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE)
+!   call GGF2_ppBSE(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE)
 
 !   write(*,*)
 !   write(*,*)'-------------------------------------------------------------------------------'

src/GF/evRGF2.f90

@@ -145,11 +145,11 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si
 
   end if
 
-! Perform BSE2 calculation
+! Perform BSE@GF2 calculation
 
   if(dophBSE) then 
   
-    call RGF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE)
+    call RGF2_phBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE)
 
     write(*,*)
     write(*,*)'-------------------------------------------------------------------------------'
@@ -162,11 +162,11 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si
 
   end if
 
-! Perform ppBSE2 calculation
+! Perform ppBSE@GF2 calculation
 
   if(doppBSE) then
 
-    call RGF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE)
+    call RGF2_ppBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE)
 
     write(*,*)
     write(*,*)'-------------------------------------------------------------------------------'

src/GF/evUGF2.f90

@@ -197,7 +197,7 @@ subroutine evUGF2(dotest,maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved
 
   if(BSE) then
 
-    print*,'!!! BSE2 NYI for evUGF2 !!!'
+    print*,'!!! BSE@GF2 NYI for evUGF2 !!!'
 
   end if
 

src/GF/qsRGF2.f90

@@ -291,11 +291,11 @@ subroutine qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, TDA,   &
 
   deallocate(c, cp, P, F, Fp, J, K, SigC, SigCp, Z, error, error_diis, F_diis)
 
-! Perform BSE calculation
+! Perform phBSE@GF2 calculation
 
   if(dophBSE) then
 
-    call RGF2_phBSE2(TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, nC, nO, &
+    call RGF2_phBSE(TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, nC, nO, &
                     nV, nR, nS, ERI_MO, dipole_int_MO, eGF, EcBSE)
 
     write(*,*)
@@ -310,11 +310,11 @@ subroutine qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, TDA,   &
   end if
 
 
-! Perform ppBSE2 calculation
+! Perform ppBSE@GF2 calculation
 
   if(doppBSE) then
 
-    call RGF2_ppBSE2(TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, &
+    call RGF2_ppBSE(TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, &
                     nC, nO, nV, nR, ERI_MO, dipole_int_MO, eGF, EcBSE)
 
     write(*,*)

src/GF/qsUGF2.f90

@@ -333,7 +333,7 @@ subroutine qsUGF2(dotest,maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved
 
   if(BSE) then
 
-    print*,'!!! BSE2 NYI for qsUGF2 !!!'
+    print*,'!!! BSE@GF(2) NYI for qsUGF2 !!!'
 
   end if
 

src/GF/ufRG0F02.f90

@@ -71,7 +71,7 @@ subroutine ufRG0F02(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
 ! Main loop over orbitals !
 !-------------------------!
 
-  do p=nO-1,nO
+  do p=nO-3,nO
 
      H(:,:) = 0d0
      Reigv(:,:) = 0d0

src/GW/GG0W0.f90

@@ -1,5 +1,5 @@
 subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & 
-                 linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
+                 linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
 ! Perform G0W0 calculation
 
   implicit none
@@ -31,7 +31,7 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
   integer,intent(in)            :: nR
   integer,intent(in)            :: nS
   double precision,intent(in)   :: ENuc
-  double precision,intent(in)   :: ERHF
+  double precision,intent(in)   :: EGHF
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: dipole_int(nBas,nBas,ncart)
   double precision,intent(in)   :: eHF(nBas)
@@ -40,7 +40,6 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
 
   logical                       :: print_W = .true.
   logical                       :: dRPA
-  integer                       :: ispin
   double precision              :: EcRPA
   double precision              :: EcBSE
   double precision              :: EcGM
@@ -85,10 +84,6 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
     write(*,*)
   end if
 
-! Spin manifold 
-
-  ispin = 3
-
 ! Memory allocation
 
   allocate(Aph(nS,nS),Bph(nS,nS),SigC(nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS), & 
@@ -98,12 +93,12 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
 ! Compute screening !
 !-------------------!
 
-                 call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
+                 call phGLR_A(dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
-  if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+  if(.not.TDA_W) call phGLR_B(dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
-  call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
+  call phGLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
 
-  if(print_W) call print_excitation_energies('phRPA@GHF','spinorbital',nS,Om)
+  if(print_W) call print_excitation_energies('phRPA@GHF','generalized',nS,Om)
 
 !--------------------------!
 ! Compute spectral weights !
@@ -147,16 +142,16 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
 
 ! Compute the RPA correlation energy
 
-  call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
+                 call phGLR_A(dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
-  if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+  if(.not.TDA_W) call phGLR_B(dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
-  call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
+  call phGLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
 
 !--------------!
 ! Dump results !
 !--------------!
 
-  call print_GG0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM)
+  call print_GG0W0(nBas,nO,eHF,ENuc,EGHF,SigC,Z,eGW,EcRPA,EcGM)
 
 ! Deallocate memory
 
@@ -171,7 +166,7 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
     write(*,*)
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@GHF correlation energy = ',EcBSE,' au'
-    write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@GHF total energy       = ',ENuc + ERHF + EcBSE,' au'
+    write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@GHF total energy       = ',ENuc + EGHF + EcBSE,' au'
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,*)
 
@@ -198,7 +193,7 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
 !     write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy (singlet) =',EcBSE(1),' au'
 !     write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy (triplet) =',EcBSE(2),' au'
 !     write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy           =',EcBSE(1) + EcBSE(2),' au'
-!     write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 total energy                 =',ENuc + ERHF + EcBSE(1) + EcBSE(2),' au'
+!     write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 total energy                 =',ENuc + EGHF + EcBSE(1) + EcBSE(2),' au'
 !     write(*,*)'-------------------------------------------------------------------------------'
 !     write(*,*)
 
@@ -206,26 +201,25 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
 
   end if
 
-! if(doppBSE) then
+  if(doppBSE) then
 
-!   call GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE)
+    call GGW_ppBSE(TDA_W,TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE)
 
-!   write(*,*)
+    write(*,*)
-!   write(*,*)'-------------------------------------------------------------------------------'
+    write(*,*)'-------------------------------------------------------------------------------'
-!   write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy (singlet) =',EcBSE(1),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@GG0W0 correlation energy          =',EcBSE,' au'
-!   write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy (triplet) =',3d0*EcBSE(2),' au'
+    write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@GG0W0 total energy                =',ENuc + EGHF + EcBSE,' au'
-!   write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy           =',EcBSE(1) + 3d0*EcBSE(2),' au'
+    write(*,*)'-------------------------------------------------------------------------------'
-!   write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 total energy                 =',ENuc + ERHF + EcBSE(1) + 3d0*EcBSE(2),' au'
+    write(*,*)
-!   write(*,*)'-------------------------------------------------------------------------------'
-!   write(*,*)
 
-! end if
+  end if
 
 ! Testing zone
 
   if(dotest) then
 
-    call dump_test_value('G','G0W0 correlation energy',EcRPA)
+    call dump_test_value('G','RPA@G0W0 correlation energy',EcRPA)
+    call dump_test_value('G','Tr@ppBSE@G0W0 correlation energy',EcBSE)
     call dump_test_value('G','G0W0 HOMO energy',eGW(nO))
     call dump_test_value('G','G0W0 LUMO energy',eGW(nO+1))
 

src/GW/GGW_ppBSE.f90

@@ -0,0 +1,125 @@
+subroutine GGW_ppBSE(TDA_W,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE)
+
+! Compute the Bethe-Salpeter excitation energies at the pp level based on a GHF reference
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  logical,intent(in)            :: TDA_W
+  logical,intent(in)            :: TDA
+  logical,intent(in)            :: dBSE
+  logical,intent(in)            :: dTDA
+
+  double precision,intent(in)   :: eta
+  integer,intent(in)            :: 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)   :: eW(nOrb)
+  double precision,intent(in)   :: eGW(nOrb)
+  double precision,intent(in)   :: ERI(nOrb,nOrb,nOrb,nOrb)
+  double precision,intent(in)   :: dipole_int(nOrb,nOrb,ncart)
+
+! Local variables
+
+  integer                       :: isp_W
+
+  logical                       :: dRPA   = .false.
+  logical                       :: dRPA_W = .true.
+
+  integer                       :: nOO
+  integer                       :: nVV
+
+  double precision,allocatable  :: Aph(:,:)
+  double precision,allocatable  :: Bph(:,:)
+
+  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  :: Bpp(:,:)
+  double precision,allocatable  :: Cpp(:,:)
+  double precision,allocatable  :: Dpp(:,:)
+
+  double precision,allocatable  :: Om1(:)
+  double precision,allocatable  :: X1(:,:)
+  double precision,allocatable  :: Y1(:,:)
+
+  double precision,allocatable  :: Om2(:)
+  double precision,allocatable  :: X2(:,:)
+  double precision,allocatable  :: Y2(:,:)
+
+  double precision,allocatable  :: KB_sta(:,:)
+  double precision,allocatable  :: KC_sta(:,:)
+  double precision,allocatable  :: KD_sta(:,:)
+  
+! Output variables
+
+  double precision,intent(out)  :: EcBSE
+
+!-----------------------!
+! Compute RPA screening !
+!-----------------------!
+
+  isp_W = 3
+  EcRPA = 0d0
+
+  allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nOrb,nOrb,nS), &
+           Aph(nS,nS),Bph(nS,nS))
+ 
+  call phLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph)
+
+  if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+
+  call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
+
+  call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
+
+  deallocate(XpY_RPA,XmY_RPA,Aph,Bph)
+
+  EcBSE = 0d0
+
+  nOO = nO*(nO-1)/2
+  nVV = nV*(nV-1)/2
+
+  allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),       &
+           Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO),       &
+           Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO), &
+           KB_sta(nVV,nOO),KC_sta(nVV,nVV),KD_sta(nOO,nOO))
+
+  ! Compute BSE excitation energies
+
+               call GGW_ppBSE_static_kernel_C(eta,nOrb,nC,nO,nV,nR,nS,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta)
+               call GGW_ppBSE_static_kernel_D(eta,nOrb,nC,nO,nV,nR,nS,nOO,1d0,ERI,OmRPA,rho_RPA,KD_sta)
+  if(.not.TDA) call GGW_ppBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta)
+  
+               call ppGLR_C(nOrb,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp)
+               call ppGLR_D(nOrb,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp)
+  if(.not.TDA) call ppGLR_B(nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
+
+  Bpp(:,:) = Bpp(:,:) + KB_sta(:,:)
+  Cpp(:,:) = Cpp(:,:) + KC_sta(:,:)
+  Dpp(:,:) = Dpp(:,:) + KD_sta(:,:)
+
+  call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE)
+
+  call ppLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)
+
+  !----------------------------------------------------!
+  ! Compute the dynamical screening at the ppBSE level !
+  !----------------------------------------------------!
+
+! if(dBSE) &
+!     call GGW_ppBSE_dynamic_perturbation(dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, &
+!                                        Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta)
+
+
+  deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta)
+
+end subroutine 

src/GW/GGW_ppBSE_static_kernel_B.f90

@@ -0,0 +1,66 @@
+subroutine GGW_ppBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Om,rho,KB)
+
+! Compute the VVOO block of the static screening W for the pp-BSE
+
+  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
+  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)   :: Om(nS)
+  double precision,intent(in)   :: rho(nBas,nBas,nS)
+
+! Local variables
+
+  double precision,external     :: Kronecker_delta
+  double precision              :: chi
+  double precision              :: eps
+  integer                       :: a,b,i,j,ab,ij,m
+
+! Output variables
+
+  double precision,intent(out)  :: KB(nVV,nOO)
+
+! Initialization
+
+  KB(:,:) = 0d0
+
+!---------------!
+! SpinOrb block !
+!---------------!
+
+  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 = Om(m)**2 + eta**2
+            chi = chi - rho(i,a,m)*rho(j,b,m)*Om(m)/eps &
+                      + rho(i,b,m)*rho(a,j,m)*Om(m)/eps
+          end do
+
+          KB(ab,ij) = 2d0*lambda*chi
+
+        end do
+      end do
+    end do
+  end do
+
+end subroutine 

src/GW/GGW_ppBSE_static_kernel_C.f90

@@ -0,0 +1,67 @@
+subroutine GGW_ppBSE_static_kernel_C(eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ERI,Om,rho,KC)
+
+! Compute the VVVV block of the static screening W for the pp-BSE
+
+  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
+  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)   :: Om(nS)
+  double precision,intent(in)   :: rho(nBas,nBas,nS)
+
+! Local variables
+
+  double precision,external     :: Kronecker_delta
+  double precision              :: chi
+  double precision              :: eps
+  double precision              :: tmp_ab, lambda4, eta2
+  integer                       :: a,b,c,d,ab,cd,m
+  integer                       :: a0, aa
+
+  double precision, allocatable :: Om_tmp(:)
+  double precision, allocatable :: tmp_m(:,:,:)
+  double precision, allocatable :: tmp(:,:,:,:)
+
+! Output variables
+
+  double precision,intent(out)  :: KC(nVV,nVV)
+
+!---------------!
+! SpinOrb block !
+!---------------!
+
+  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 = Om(m)**2 + eta**2
+            chi = chi - rho(a,c,m)*rho(b,d,m)*Om(m)/eps &
+                      + rho(a,d,m)*rho(b,c,m)*Om(m)/eps
+          end do
+         
+          KC(ab,cd) = 2d0*lambda*chi
+
+        end do
+      end do
+    end do
+  end do
+
+end subroutine 

src/GW/GGW_ppBSE_static_kernel_D.f90

@@ -0,0 +1,65 @@
+subroutine GGW_ppBSE_static_kernel_D(eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,ERI,Om,rho,KD)
+
+! Compute the OOOO block of the static screening W for the pp-BSE
+
+  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
+  integer,intent(in)            :: nOO
+  double precision,intent(in)   :: eta
+  double precision,intent(in)   :: lambda
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: Om(nS)
+  double precision,intent(in)   :: rho(nBas,nBas,nS)
+
+! Local variables
+
+  double precision,external     :: Kronecker_delta
+  double precision              :: chi
+  double precision              :: eps
+  integer                       :: i,j,k,l,ij,kl,m
+
+! Output variables
+
+  double precision,intent(out)  :: KD(nOO,nOO)
+
+! Initialization
+
+  KD(:,:) = 0d0
+
+!---------------!
+! SpinOrb block !
+!---------------!
+
+  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 = Om(m)**2 + eta**2
+            chi = chi - rho(i,k,m)*rho(j,l,m)*Om(m)/eps &
+                      + rho(i,l,m)*rho(j,k,m)*Om(m)/eps
+          end do
+
+          KD(ij,kl) = 2d0*lambda*chi
+
+        end do
+      end do
+    end do
+  end do
+
+end subroutine 

src/GW/evGGW.f90

@@ -43,7 +43,6 @@ subroutine evGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
 
   logical                       :: linear_mixing
   logical                       :: dRPA = .true.
-  integer                       :: ispin
   integer                       :: nSCF
   integer                       :: n_diis
   double precision              :: rcond
@@ -100,7 +99,6 @@ subroutine evGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
 ! Initialization
 
   nSCF            = 0
-  ispin           = 3
   n_diis          = 0
   Conv            = 1d0
   e_diis(:,:)     = 0d0
@@ -118,10 +116,10 @@ subroutine evGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
 
    ! Compute screening
 
-    call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
+                   call phGLR_A(dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
-    if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+    if(.not.TDA_W) call phGLR_B(dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
-    call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
+    call phGLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
 
     ! Compute spectral weights
 

src/GW/qsGGW.f90

@@ -58,7 +58,6 @@ subroutine qsGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
   integer                       :: nSCF
   integer                       :: nBasSq
   integer                       :: nBas2Sq
-  integer                       :: ispin
   integer                       :: ixyz
   integer                       :: n_diis
   double precision              :: ET,ETaa,ETbb
@@ -154,7 +153,6 @@ subroutine qsGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
   
   nSCF          = -1
   n_diis        = 0
-  ispin         = 3
   Conv          = 1d0
   P(:,:)        = PHF(:,:)
   eGW(:)        = eHF(:)
@@ -253,11 +251,11 @@ subroutine qsGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
 
 !   Compute linear response
 
-    call phLR_A(ispin,dRPA,nBas2,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph)
+                   call phGLR_A(dRPA,nBas2,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph)
-    if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas2,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
+    if(.not.TDA_W) call phGLR_B(dRPA,nBas2,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
 
-    call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
+    call phGLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
-    if(print_W) call print_excitation_energies('phRPA@GW@GHF','spinorbital',nS,Om)
+    if(print_W) call print_excitation_energies('phRPA@GW@GHF','generalized',nS,Om)
 
 !   Compute correlation part of the self-energy 
 

src/LR/phGLR.f90

@@ -0,0 +1,82 @@
+subroutine phGLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
+
+! Compute linear response
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  logical,intent(in)            :: TDA
+  integer,intent(in)            :: nS
+  double precision,intent(in)   :: Aph(nS,nS)
+  double precision,intent(in)   :: Bph(nS,nS)
+
+  ! Local variables
+
+  double precision              :: trace_matrix
+  double precision,allocatable  :: ApB(:,:)
+  double precision,allocatable  :: AmB(:,:)
+  double precision,allocatable  :: AmBSq(:,:)
+  double precision,allocatable  :: AmBIv(:,:)
+  double precision,allocatable  :: Z(:,:)
+  double precision,allocatable  :: tmp(:,:)
+
+! Output variables
+
+  double precision,intent(out)  :: EcRPA
+  double precision,intent(out)  :: Om(nS)
+  double precision,intent(out)  :: XpY(nS,nS)
+  double precision,intent(out)  :: XmY(nS,nS)
+
+! Memory allocation
+
+  allocate(ApB(nS,nS),AmB(nS,nS),AmBSq(nS,nS),AmBIv(nS,nS),Z(nS,nS),tmp(nS,nS))
+
+! Tamm-Dancoff approximation
+
+  if(TDA) then
+ 
+    XpY(:,:) = Aph(:,:)
+    call diagonalize_matrix(nS,XpY,Om)
+    XpY(:,:) = transpose(XpY(:,:))
+    XmY(:,:) = XpY(:,:)
+
+  else
+
+    ApB(:,:) = Aph(:,:) + Bph(:,:)
+    AmB(:,:) = Aph(:,:) - Bph(:,:)
+
+  ! Diagonalize linear response matrix
+
+    call diagonalize_matrix(nS,AmB,Om)
+
+    if(minval(Om) < 0d0) &
+      call print_warning('You may have instabilities in linear response: A-B is not positive definite!!')
+
+    call ADAt(nS,AmB,1d0*dsqrt(Om),AmBSq)
+    call ADAt(nS,AmB,1d0/dsqrt(Om),AmBIv)
+
+    call dgemm('N','N',nS,nS,nS,1d0,ApB,size(ApB,1),AmBSq,size(AmBSq,1),0d0,tmp,size(tmp,1))
+    call dgemm('N','N',nS,nS,nS,1d0,AmBSq,size(AmBSq,1),tmp,size(tmp,1),0d0,Z,size(Z,1))
+
+    call diagonalize_matrix(nS,Z,Om)
+
+    if(minval(Om) < 0d0) & 
+      call print_warning('You may have instabilities in linear response: negative excitations!!')
+ 
+    Om = sqrt(Om)
+
+    call dgemm('T','N',nS,nS,nS,1d0,Z,size(Z,1),AmBSq,size(AmBSq,1),0d0,XpY,size(XpY,1))
+    call DA(nS,1d0/dsqrt(Om),XpY)
+
+    call dgemm('T','N',nS,nS,nS,1d0,Z,size(Z,1),AmBIv,size(AmBIv,1),0d0,XmY,size(XmY,1))
+    call DA(nS,1d0*dsqrt(Om),XmY)
+ 
+  end if
+
+  ! Compute the RPA correlation energy
+
+  EcRPA = 0.5d0*(sum(Om) - trace_matrix(nS,Aph))
+
+end subroutine 

src/LR/phGLR_A.f90

@@ -0,0 +1,56 @@
+subroutine phGLR_A(dRPA,nOrb,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
+
+! Compute resonant block of the ph channel
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  logical,intent(in)            :: dRPA
+  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)   :: lambda
+  double precision,intent(in)   :: e(nOrb)
+  double precision,intent(in)   :: ERI(nOrb,nOrb,nOrb,nOrb) 
+  
+! Local variables
+
+  double precision              :: delta_dRPA
+  double precision,external     :: Kronecker_delta
+
+  integer                       :: i,j,a,b,ia,jb
+
+! Output variables
+
+  double precision,intent(out)  :: Aph(nS,nS)
+
+! Direct RPA
+
+  delta_dRPA = 0d0
+  if(dRPA) delta_dRPA = 1d0
+
+! Build A matrix for spin orbitals
+
+  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
+
+          Aph(ia,jb) = (e(a) - e(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
+                     + lambda*ERI(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a)
+
+        end  do
+      end  do
+    end  do
+  end  do
+
+end subroutine 

src/LR/phGLR_B.f90

@@ -0,0 +1,53 @@
+subroutine phGLR_B(dRPA,nOrb,nC,nO,nV,nR,nS,lambda,ERI,Bph)
+
+! Compute the coupling block of the ph channel
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  logical,intent(in)            :: dRPA
+  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)   :: lambda
+  double precision,intent(in)   :: ERI(nOrb,nOrb,nOrb,nOrb)
+  
+! Local variables
+
+  double precision              :: delta_dRPA
+
+  integer                       :: i,j,a,b,ia,jb
+
+! Output variables
+
+  double precision,intent(out)  :: Bph(nS,nS)
+
+! Direct RPA
+
+  delta_dRPA = 0d0
+  if(dRPA) delta_dRPA = 1d0
+
+! Build B matrix for spin orbitals 
+
+  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
+
+          Bph(ia,jb) = lambda*ERI(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a)
+
+        end do
+      end do
+    end do
+  end do
+
+end subroutine 

src/LR/phLR_oscillator_strength.f90

@@ -1,4 +1,4 @@
-subroutine phLR_oscillator_strength(nBas,nC,nO,nV,nR,nS,maxS,dipole_int,Om,XpY,XmY,os)
+subroutine phLR_oscillator_strength(nOrb,nC,nO,nV,nR,nS,maxS,dipole_int,Om,XpY,XmY,os)
 
 ! Compute linear response
 
@@ -7,14 +7,14 @@ subroutine phLR_oscillator_strength(nBas,nC,nO,nV,nR,nS,maxS,dipole_int,Om,XpY,X
 
 ! Input variables
 
-  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
   integer,intent(in)            :: maxS
-  double precision              :: dipole_int(nBas,nBas,ncart)
+  double precision              :: dipole_int(nOrb,nOrb,ncart)
   double precision,intent(in)   :: Om(nS)
   double precision,intent(in)   :: XpY(nS,nS)
   double precision,intent(in)   :: XmY(nS,nS)
@@ -44,7 +44,7 @@ subroutine phLR_oscillator_strength(nBas,nC,nO,nV,nR,nS,maxS,dipole_int,Om,XpY,X
     do ixyz=1,ncart
       jb = 0
       do j=nC+1,nO
-        do b=nO+1,nBas-nR
+        do b=nO+1,nOrb-nR
           jb = jb + 1
           f(m,ixyz) = f(m,ixyz) + dipole_int(j,b,ixyz)*XpY(m,jb)
         end do
@@ -68,8 +68,4 @@ subroutine phLR_oscillator_strength(nBas,nC,nO,nV,nR,nS,maxS,dipole_int,Om,XpY,X
   write(*,*) '---------------------------------------------------------------'
   write(*,*)
 
-! do m=1,maxS
-!   write(*,'(I3,3F12.6)') m,Om(m),os(m)
-! end do
-
 end subroutine 

src/LR/ppGLR.f90

@@ -0,0 +1,125 @@
+subroutine ppGLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
+
+  ! 
+  ! Solve the pp-RPA linear eigenvalue problem
+  !
+  ! right eigen-problem: H   R = R w
+  ! left  eigen-problem: H.T L = L w
+  !
+  !     where L.T R = 1
+  !
+  !
+  !        (+C    +B) 
+  !   H =  (        )  where C = C.T and D = D.T
+  !        (-B.T  -D)
+  !
+  !      (w1   0)        (X1   X2)          (+X1  +X2)
+  !  w = (      ),   R = (       )  and L = (        )
+  !      (0   w2)        (Y1   Y2)          (-Y1  -Y2)
+  ! 
+  !
+  !  the normalisation condition reduces to
+  !
+  !    X1.T X2 - Y1.T Y2 = 0
+  !    X1.T X1 - Y1.T Y1 = 1
+  !    X2.T X2 - Y2.T Y2 = 1
+  !
+
+  implicit none
+  include 'parameters.h'
+
+  logical,          intent(in)  :: TDA
+  integer,          intent(in)  :: nOO, nVV
+  double precision, intent(in)  :: Bpp(nVV,nOO), Cpp(nVV,nVV), Dpp(nOO,nOO)
+  double precision, intent(out) :: Om1(nVV), X1(nVV,nVV), Y1(nOO,nVV)
+  double precision, intent(out) :: Om2(nOO), X2(nVV,nOO), Y2(nOO,nOO)
+  double precision, intent(out) :: EcRPA
+  
+  logical                       :: imp_bio, verbose
+  integer                       :: i, j, N
+  double precision              :: EcRPA1, EcRPA2
+  double precision              :: thr_d, thr_nd, thr_deg
+  double precision,allocatable  :: M(:,:), Z(:,:), Om(:)
+
+  double precision, external    :: trace_matrix
+
+
+
+  N = nOO + nVV
+
+  allocate(M(N,N), Z(N,N), Om(N))
+
+  if(TDA) then
+
+    X1(:,:) = +Cpp(:,:)
+    Y1(:,:) = 0d0
+    if(nVV > 0) call diagonalize_matrix(nVV, X1, Om1)
+
+    X2(:,:) = 0d0
+    Y2(:,:) = -Dpp(:,:)
+    if(nOO > 0) call diagonalize_matrix(nOO, Y2, Om2)
+
+  else
+
+    ! Diagonal blocks 
+    M(    1:nVV    ,    1:nVV)     = + Cpp(1:nVV,1:nVV)
+    M(nVV+1:nVV+nOO,nVV+1:nVV+nOO) = - Dpp(1:nOO,1:nOO)
+
+    ! Off-diagonal blocks
+    M(    1:nVV    ,nVV+1:nOO+nVV) = -           Bpp(1:nVV,1:nOO)
+    M(nVV+1:nOO+nVV,    1:nVV)     = + transpose(Bpp(1:nVV,1:nOO))
+
+    if((nOO .eq. 0) .or. (nVV .eq. 0)) then
+
+      ! Diagonalize the p-p matrix
+      if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV, M, Om, Z)
+      ! Split the various quantities in p-p and h-h parts
+      call sort_ppRPA(nOO, nVV, Om, Z, Om1, X1, Y1, Om2, X2, Y2)
+
+    else
+
+        thr_d   = 1d-6   ! to determine if     diagonal elements of L.T x R are close enouph to 1
+        thr_nd  = 1d-6   ! to determine if non-diagonal elements of L.T x R are close enouph to 1
+        thr_deg = 1d-8   ! to determine if two eigenvectors are degenerate or not
+        imp_bio = .True. ! impose bi-orthogonality
+        verbose = .False.
+        call diagonalize_nonsym_matrix(N, M, Z, Om, thr_d, thr_nd, thr_deg, imp_bio, verbose)
+    
+        do i = 1, nOO
+          Om2(i) = Om(i)
+          do j = 1, nVV
+            X2(j,i) = Z(j,i)
+          enddo
+          do j = 1, nOO
+            Y2(j,i) = Z(nVV+j,i)
+          enddo
+        enddo
+    
+        do i = 1, nVV
+          Om1(i) = Om(nOO+i)
+          do j = 1, nVV
+            X1(j,i) = M(j,nOO+i)
+          enddo
+          do j = 1, nOO
+            Y1(j,i) = M(nVV+j,nOO+i)
+          enddo
+        enddo
+
+    endif
+
+  end if
+
+  ! Compute the RPA correlation energy
+  EcRPA = 0.5d0 * (sum(Om1) - sum(Om2) - trace_matrix(nVV, Cpp) - trace_matrix(nOO, Dpp))
+  EcRPA1 = +sum(Om1) - trace_matrix(nVV, Cpp)
+  EcRPA2 = -sum(Om2) - trace_matrix(nOO, Dpp)
+
+  if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) then
+    print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!'
+  endif
+
+  deallocate(M, Z, Om)
+
+end subroutine 
+
+

src/LR/ppGLR_B.f90

@@ -0,0 +1,47 @@
+subroutine ppGLR_B(nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp)
+
+! Compute the B matrix of the pp channel
+
+  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)            :: nOO
+  integer,intent(in)            :: nVV
+  double precision,intent(in)   :: lambda
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas) 
+  
+! Local variables
+
+  double precision,external     :: Kronecker_delta
+  integer                       :: a,b,i,j,ab,ij
+
+! Output variables
+
+  double precision,intent(out)  :: Bpp(nVV,nOO)
+
+! Build the B matrix for the spin-orbital basis
+
+  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
+ 
+          Bpp(ab,ij) = lambda*(ERI(a,b,i,j) - ERI(a,b,j,i))
+ 
+        end do
+      end do
+    end do
+  end do
+
+end subroutine 

src/LR/ppGLR_C.f90

@@ -0,0 +1,61 @@
+subroutine ppGLR_C(nBas,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp)
+
+! Compute the C matrix of the pp channel
+
+  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)            :: nVV
+  double precision,intent(in)   :: lambda
+  double precision,intent(in)   :: e(nBas),ERI(nBas,nBas,nBas,nBas) 
+  
+! Local variables
+
+  double precision              :: eF
+  double precision,external     :: Kronecker_delta
+
+  integer                       :: a,b,c,d,ab,cd
+  integer                       :: a0, aa
+  double precision              :: e_ab, tmp_ab, delta_ac, tmp_cd
+
+! Output variables
+
+  double precision,intent(out)  :: Cpp(nVV,nVV)
+
+! Define the chemical potential
+ 
+! eF = e(nO) + e(nO+1)
+  eF = 0d0
+
+! Build C matrix for the singlet manifold
+
+  !$OMP PARALLEL &
+  !$OMP SHARED(Cpp,lambda,ERI,e,eF,nC,nO,nBas,nR) &
+  !$OMP PRIVATE(c,d,a,b,ab,cd) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO
+  do c=nO+1,nBas-nR
+    do d=c+1,nBas-nR
+      cd = (c-(nO+1))*(nBas-nR-(nO+1)) - (c-1-(nO+1))*(c-(nO+1))/2 + d - c
+        do a=nO+1,nBas-nR
+          do b=a+1,nBas-nR
+            ab = (a-(nO+1))*(nBas-nR-(nO+1)) - (a-1-(nO+1))*(a-(nO+1))/2 + b - a
+
+            Cpp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & 
+                 + lambda*(ERI(a,b,c,d) - ERI(a,b,d,c))
+
+        end do
+      end do
+    end do
+ end do
+ !$OMP END DO
+ !$OMP END PARALLEL
+
+end subroutine 

src/LR/ppGLR_D.f90

@@ -0,0 +1,54 @@
+subroutine ppGLR_D(nBas,nC,nO,nV,nR,nOO,lambda,e,ERI,Dpp)
+
+! Compute the D matrix of the pp channel
+
+  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)            :: nOO
+  double precision,intent(in)   :: lambda
+  double precision,intent(in)   :: e(nBas),ERI(nBas,nBas,nBas,nBas) 
+  
+! Local variables
+
+  double precision              :: eF
+  double precision,external     :: Kronecker_delta
+
+  integer                       :: i,j,k,l,ij,kl
+
+! Output variables
+
+  double precision,intent(out)  :: Dpp(nOO,nOO)
+
+! Define the chemical potential
+
+! eF = e(nO) + e(nO+1)
+  eF = 0d0
+ 
+! Build the D matrix for the spin-orbital basis 
+
+  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
+
+          Dpp(ij,kl) = - (e(i) + e(j) - eF)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & 
+                       + lambda*(ERI(i,j,k,l) - ERI(i,j,l,k))
+
+        end do
+      end do
+    end do
+  end do
+
+end subroutine 

src/LR/ppLR.f90

@@ -1,7 +1,4 @@
-
+subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
-! ---
-
-subroutine ppLR(TDA, nOO, nVV, Bpp, Cpp, Dpp, Om1, X1, Y1, Om2, X2, Y2, EcRPA)
 
   ! 
   ! Solve the pp-RPA linear eigenvalue problem

src/LR/ppLR_transition_vectors.f90

@@ -117,8 +117,8 @@ subroutine ppLR_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nOO,nVV,dipole_
 
 ! Thomas-Reiche-Kuhn sum rule
 
-  if(nVV > 0) write(*,'(A50,F10.6)') 'Thomas-Reiche-Kuhn sum rule for p-p sector = ',sum(os1(:))
+!  if(nVV > 0) write(*,'(A50,F10.6)') 'Thomas-Reiche-Kuhn sum rule for p-p sector = ',sum(os1(:))
-  write(*,*)
+!  write(*,*)
 
 !-----------------------------------------------!
 ! Print details about excitations for hh sector !
@@ -188,7 +188,7 @@ subroutine ppLR_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nOO,nVV,dipole_
 
 ! Thomas-Reiche-Kuhn sum rule
 
-  if(nOO > 0) write(*,'(A50,F10.6)') 'Thomas-Reiche-Kuhn sum rule for h-h sector = ',sum(os2(:))
+!  if(nOO > 0) write(*,'(A50,F10.6)') 'Thomas-Reiche-Kuhn sum rule for h-h sector = ',sum(os2(:))
-  write(*,*)
+!  write(*,*)
 
 end subroutine 

src/RPA/crGRPA.f90

@@ -1,4 +1,4 @@
-subroutine crGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
+subroutine crGRPA(dotest,TDA,nOrb,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
 
 ! Crossed-ring channel of the random phase approximation 
 
@@ -11,7 +11,7 @@ subroutine crGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
   logical,intent(in)            :: dotest
 
   logical,intent(in)            :: TDA
-  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nOrb
   integer,intent(in)            :: nC
   integer,intent(in)            :: nO
   integer,intent(in)            :: nV
@@ -19,13 +19,12 @@ subroutine crGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
   integer,intent(in)            :: nS
   double precision,intent(in)   :: ENuc
   double precision,intent(in)   :: EGHF
-  double precision,intent(in)   :: eHF(nBas)
+  double precision,intent(in)   :: eHF(nOrb)
-  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI(nOrb,nOrb,nOrb,nOrb)
-  double precision,intent(in)   :: dipole_int(nBas,nBas,ncart)
+  double precision,intent(in)   :: dipole_int(nOrb,nOrb,ncart)
 
 ! Local variables
 
-  integer                       :: ispin
   logical                       :: dRPA
   double precision,allocatable  :: Aph(:,:)
   double precision,allocatable  :: Bph(:,:)
@@ -59,14 +58,12 @@ subroutine crGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
 
   allocate(Om(nS),XpY(nS,nS),XmY(nS,nS),Aph(nS,nS),Bph(nS,nS))
 
-  ispin = 3
+               call phLR_A(dRPA,nOrb,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Aph)
+  if(.not.TDA) call phLR_B(dRPA,nOrb,nC,nO,nV,nR,nS,-1d0,ERI,Bph)
 
-  call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Aph)
+  call phGLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
-  if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,ERI,Bph)
-
-  call phLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
   call print_excitation_energies('crRPA@GHF','spinorbital',nS,Om)
-  call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
+  call phLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
 
   write(*,*)
   write(*,*)'-------------------------------------------------------------------------------'

src/RPA/phGRPA.f90

@@ -1,4 +1,4 @@
-subroutine phGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
+subroutine phGRPA(dotest,TDA,nOrb,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
 
 ! Perform a direct random phase approximation calculation
 
@@ -11,7 +11,7 @@ subroutine phGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
   logical,intent(in)            :: dotest
 
   logical,intent(in)            :: TDA
-  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nOrb
   integer,intent(in)            :: nC
   integer,intent(in)            :: nO
   integer,intent(in)            :: nV
@@ -19,13 +19,12 @@ subroutine phGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
   integer,intent(in)            :: nS
   double precision,intent(in)   :: ENuc
   double precision,intent(in)   :: EGHF
-  double precision,intent(in)   :: eHF(nBas)
+  double precision,intent(in)   :: eHF(nOrb)
-  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI(nOrb,nOrb,nOrb,nOrb)
-  double precision,intent(in)   :: dipole_int(nBas,nBas,ncart)
+  double precision,intent(in)   :: dipole_int(nOrb,nOrb,ncart)
 
 ! Local variables
 
-  integer                       :: ispin
   logical                       :: dRPA
   double precision              :: lambda
 
@@ -62,14 +61,12 @@ subroutine phGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
 
   allocate(Om(nS),XpY(nS,nS),XmY(nS,nS),Aph(nS,nS),Bph(nS,nS))
 
-  ispin = 3
+               call phGLR_A(dRPA,nOrb,nC,nO,nV,nR,nS,lambda,eHF,ERI,Aph)
+  if(.not.TDA) call phGLR_B(dRPA,nOrb,nC,nO,nV,nR,nS,lambda,ERI,Bph)
 
-  call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,Aph)
+  call phGLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
-  if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
-
-  call phLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
   call print_excitation_energies('phRPA@GHF','spinorbital',nS,Om)
-  call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
+  call phLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
 
   write(*,*)
   write(*,*)'-------------------------------------------------------------------------------'

src/RPA/phGRPAx.f90

@@ -1,4 +1,4 @@
-subroutine phGRPAx(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
+subroutine phGRPAx(dotest,TDA,nOrb,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
 
 ! Perform random phase approximation calculation with exchange (aka TDHF)
 
@@ -11,7 +11,7 @@ subroutine phGRPAx(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
   logical,intent(in)            :: dotest
 
   logical,intent(in)            :: TDA
-  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nOrb
   integer,intent(in)            :: nC
   integer,intent(in)            :: nO
   integer,intent(in)            :: nV
@@ -19,13 +19,12 @@ subroutine phGRPAx(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
   integer,intent(in)            :: nS
   double precision,intent(in)   :: ENuc
   double precision,intent(in)   :: EGHF
-  double precision,intent(in)   :: eHF(nBas)
+  double precision,intent(in)   :: eHF(nOrb)
-  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI(nOrb,nOrb,nOrb,nOrb)
-  double precision,intent(in)   :: dipole_int(nBas,nBas,ncart)
+  double precision,intent(in)   :: dipole_int(nOrb,nOrb,ncart)
 
 ! Local variables
 
-  integer                       :: ispin
   logical                       :: dRPA
   double precision,allocatable  :: Aph(:,:)
   double precision,allocatable  :: Bph(:,:)
@@ -59,14 +58,12 @@ subroutine phGRPAx(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
 
   allocate(Om(nS),XpY(nS,nS),XmY(nS,nS),Aph(nS,nS),Bph(nS,nS))
 
-  ispin = 3
+               call phGLR_A(dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
+  if(.not.TDA) call phLR_B(dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
-  call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
+  call phGLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
-  if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
-
-  call phLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
   call print_excitation_energies('phRPAx@GHF','spinorbital',nS,Om)
-  call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
+  call phLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
 
   write(*,*)
   write(*,*)'-------------------------------------------------------------------------------'

src/RPA/ppGRPA.f90

@@ -23,7 +23,6 @@ subroutine ppGRPA(dotest,TDA,nBas,nC,nO,nV,nR,ENuc,EGHF,ERI,dipole_int,eHF)
 
 ! Local variables
 
-  integer                       :: ispin
   integer                       :: nOO
   integer                       :: nVV
   double precision,allocatable  :: Bpp(:,:)
@@ -50,19 +49,17 @@ subroutine ppGRPA(dotest,TDA,nBas,nC,nO,nV,nR,ENuc,EGHF,ERI,dipole_int,eHF)
 
   EcRPA = 0d0
 
-  ispin = 4
-
   nOO = nO*(nO-1)/2
   nVV = nV*(nV-1)/2
 
   allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO), &
            Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO))
 
-  if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
+  if(.not.TDA) call ppGLR_B(nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
-               call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eHF,ERI,Cpp)
+               call ppGLR_C(nBas,nC,nO,nV,nR,nVV,1d0,eHF,ERI,Cpp)
-               call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eHF,ERI,Dpp)
+               call ppGLR_D(nBas,nC,nO,nV,nR,nOO,1d0,eHF,ERI,Dpp)
 
-  call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
+  call ppGLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
 
 !   call print_transition_vectors_pp(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)
 

src/make_ninja.py

@@ -89,8 +89,8 @@ FIX_ORDER_OF_LIBS=-Wl,--start-group
 
 if sys.platform in ["linux", "linux2"]:
 #   compiler = compile_gfortran_linux
-   compiler = compile_ifort_linux
+#   compiler = compile_ifort_linux 
-#   compiler = compile_olympe
+    compiler = compile_olympe
 elif sys.platform == "darwin":
   compiler = compile_gfortran_mac
 else: