diff --git a/src/GF/GF2_phBSE2.f90 b/src/GF/GF2_phBSE2.f90
index 3acba64..0ea1fd8 100644
--- a/src/GF/GF2_phBSE2.f90
+++ b/src/GF/GF2_phBSE2.f90
@@ -56,8 +56,8 @@ subroutine GF2_phBSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,n
 
     ! Compute static kernel
 
-                 call BSE2_A_matrix_static(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,A_sta(:,:,ispin))
-    if(.not.TDA) call BSE2_B_matrix_static(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,B_sta(:,:,ispin))
+                 call GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,A_sta(:,:,ispin))
+    if(.not.TDA) call GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,B_sta(:,:,ispin))
 
     ! Compute BSE2 excitation energies
 
@@ -73,12 +73,12 @@ subroutine GF2_phBSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,n
 
      if(evDyn) then
 
-      call BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & 
-                                               A_sta(:,:,ispin),B_sta(:,:,ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
+      call GF2_phBSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & 
+                                                     A_sta(:,:,ispin),B_sta(:,:,ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
      else
  
-      call BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & 
-                                     A_sta(:,:,ispin),B_sta(:,:,ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
+      call GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & 
+                                           A_sta(:,:,ispin),B_sta(:,:,ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
  
       end if
 
@@ -97,8 +97,8 @@ subroutine GF2_phBSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,n
 
     ! Compute static kernel
 
-                 call BSE2_A_matrix_static(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,A_sta(:,:,ispin))
-    if(.not.TDA) call BSE2_B_matrix_static(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,B_sta(:,:,ispin))
+                 call GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,A_sta(:,:,ispin))
+    if(.not.TDA) call GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,B_sta(:,:,ispin))
 
     ! Compute BSE2 excitation energies
 
@@ -113,12 +113,12 @@ subroutine GF2_phBSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,n
 
      if(evDyn) then
 
-      call BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & 
-                                               A_sta(:,:,ispin),B_sta(:,:,ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
+      call GF2_phBSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & 
+                                                     A_sta(:,:,ispin),B_sta(:,:,ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
      else
  
-      call BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, &
-                                     A_sta(:,:,ispin),B_sta(:,:,ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
+      call GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, &
+                                           A_sta(:,:,ispin),B_sta(:,:,ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
  
       end if
 
diff --git a/src/GF/BSE2_A_matrix_dynamic.f90 b/src/GF/GF2_phBSE2_dynamic_kernel_A.f90
similarity index 97%
rename from src/GF/BSE2_A_matrix_dynamic.f90
rename to src/GF/GF2_phBSE2_dynamic_kernel_A.f90
index 18216b4..31a01b7 100644
--- a/src/GF/BSE2_A_matrix_dynamic.f90
+++ b/src/GF/GF2_phBSE2_dynamic_kernel_A.f90
@@ -1,4 +1,4 @@
-subroutine BSE2_A_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,OmBSE,A_dyn,ZA_dyn)
+subroutine GF2_phBSE2_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,OmBSE,A_dyn,ZA_dyn)
 
 ! Compute the resonant part of the dynamic BSE2 matrix
 
@@ -165,4 +165,4 @@ subroutine BSE2_A_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,Om
   end if
 
 
-end subroutine BSE2_A_matrix_dynamic
+end subroutine 
diff --git a/src/GF/BSE2_B_matrix_dynamic.f90 b/src/GF/GF2_phBSE2_dynamic_kernel_B.f90
similarity index 97%
rename from src/GF/BSE2_B_matrix_dynamic.f90
rename to src/GF/GF2_phBSE2_dynamic_kernel_B.f90
index aefbea7..2ced77e 100644
--- a/src/GF/BSE2_B_matrix_dynamic.f90
+++ b/src/GF/GF2_phBSE2_dynamic_kernel_B.f90
@@ -1,4 +1,4 @@
-subroutine BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_dyn)
+subroutine GF2_phBSE2_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_dyn)
 
 ! Compute the anti-resonant part of the dynamic BSE2 matrix
 
@@ -154,4 +154,4 @@ subroutine BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_
   end if
 
 
-end subroutine BSE2_B_matrix_dynamic
+end subroutine 
diff --git a/src/GF/BSE2_dynamic_perturbation.f90 b/src/GF/GF2_phBSE2_dynamic_perturbation.f90
similarity index 89%
rename from src/GF/BSE2_dynamic_perturbation.f90
rename to src/GF/GF2_phBSE2_dynamic_perturbation.f90
index 50077f1..5705fdf 100644
--- a/src/GF/BSE2_dynamic_perturbation.f90
+++ b/src/GF/GF2_phBSE2_dynamic_perturbation.f90
@@ -1,4 +1,4 @@
-subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,A_sta,B_sta,OmBSE,XpY,XmY)
+subroutine GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,A_sta,B_sta,OmBSE,XpY,XmY)
 
 ! Compute dynamical effects via perturbation theory for BSE
 
@@ -78,7 +78,7 @@ subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipo
 
     ! Resonant part of the BSE correction for dynamical TDA
 
-    call BSE2_A_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,+OmBSE(ia),Ap_dyn,ZAp_dyn)
+    call GF2_phBSE2_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,+OmBSE(ia),Ap_dyn,ZAp_dyn)
 
     if(dTDA) then 
 
@@ -89,9 +89,9 @@ subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipo
 
       ! Second part of the resonant and anti-resonant part of the BSE correction (frequency independent)
 
-      call BSE2_A_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,-OmBSE(ia),Am_dyn,ZAm_dyn)
+      call GF2_phBSE2_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,-OmBSE(ia),Am_dyn,ZAm_dyn)
 
-      call BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,B_dyn)
+      call GF2_phBSE2_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,B_dyn)
 
       ZDyn(ia)  = dot_product(X,matmul(ZAp_dyn,X)) &
                 + dot_product(Y,matmul(ZAm_dyn,Y))  
@@ -113,4 +113,4 @@ subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipo
   write(*,*) '---------------------------------------------------------------------------------------------------'
   write(*,*) 
 
-end subroutine BSE2_dynamic_perturbation
+end subroutine 
diff --git a/src/GF/BSE2_dynamic_perturbation_iterative.f90 b/src/GF/GF2_phBSE2_dynamic_perturbation_iterative.f90
similarity index 90%
rename from src/GF/BSE2_dynamic_perturbation_iterative.f90
rename to src/GF/GF2_phBSE2_dynamic_perturbation_iterative.f90
index 5f72a5f..97114a2 100644
--- a/src/GF/BSE2_dynamic_perturbation_iterative.f90
+++ b/src/GF/GF2_phBSE2_dynamic_perturbation_iterative.f90
@@ -1,5 +1,5 @@
-subroutine BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int, & 
-                                               eHF,eGF,A_sta,B_sta,OmBSE,XpY,XmY)
+subroutine GF2_phBSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int, & 
+                                                     eHF,eGF,A_sta,B_sta,OmBSE,XpY,XmY)
 
 ! Compute self-consistently the dynamical effects via perturbation theory for BSE2
 
@@ -100,7 +100,7 @@ subroutine BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,n
       Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:))
  
      ! Resonant part of the BSE correction
-      call BSE2_A_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,+OmOld(ia),Ap_dyn,ZAp_dyn)
+      call GF2_phBSE2_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,+OmOld(ia),Ap_dyn,ZAp_dyn)
 
       if(dTDA) then 
  
@@ -111,9 +111,9 @@ subroutine BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,n
  
         ! Anti-resonant part of the BSE correction
  
-      call BSE2_A_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,-OmOld(ia),Am_dyn,ZAm_dyn)
+      call GF2_phBSE2_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,-OmOld(ia),Am_dyn,ZAm_dyn)
 
-      call BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,B_dyn)
+      call GF2_phBSE2_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,B_dyn)
 
       ZDyn(ia)  = dot_product(X,matmul(ZAp_dyn,X)) &
                 + dot_product(Y,matmul(ZAm_dyn,Y))  
@@ -152,4 +152,4 @@ subroutine BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,n
 
   endif
 
-end subroutine BSE2_dynamic_perturbation_iterative
+end subroutine 
diff --git a/src/GF/BSE2_A_matrix_static.f90 b/src/GF/GF2_phBSE2_static_kernel_A.f90
similarity index 97%
rename from src/GF/BSE2_A_matrix_static.f90
rename to src/GF/GF2_phBSE2_static_kernel_A.f90
index dee7202..ce8fede 100644
--- a/src/GF/BSE2_A_matrix_static.f90
+++ b/src/GF/GF2_phBSE2_static_kernel_A.f90
@@ -1,4 +1,4 @@
-subroutine BSE2_A_matrix_static(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,A_sta)
+subroutine GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,A_sta)
 
 ! Compute the resonant part of the static BSE2 matrix
 
@@ -154,4 +154,4 @@ subroutine BSE2_A_matrix_static(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,A_s
   end if
 
 
-end subroutine BSE2_A_matrix_static
+end subroutine 
diff --git a/src/GF/BSE2_B_matrix_static.f90 b/src/GF/GF2_phBSE2_static_kernel_B.f90
similarity index 97%
rename from src/GF/BSE2_B_matrix_static.f90
rename to src/GF/GF2_phBSE2_static_kernel_B.f90
index e93abba..011949c 100644
--- a/src/GF/BSE2_B_matrix_static.f90
+++ b/src/GF/GF2_phBSE2_static_kernel_B.f90
@@ -1,4 +1,4 @@
-subroutine BSE2_B_matrix_static(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_sta)
+subroutine GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_sta)
 
 ! Compute the anti-resonant part of the static BSE2 matrix
 
@@ -154,4 +154,4 @@ subroutine BSE2_B_matrix_static(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_s
   end if
 
 
-end subroutine BSE2_B_matrix_static
+end subroutine