diff --git a/src/QuAcK/Bethe_Salpeter.f90 b/src/QuAcK/Bethe_Salpeter.f90
index ba77a3a..1ffe1b0 100644
--- a/src/QuAcK/Bethe_Salpeter.f90
+++ b/src/QuAcK/Bethe_Salpeter.f90
@@ -57,8 +57,8 @@ subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, &
 
     ! Compute dynamic correction for BSE via perturbation theory
 
-    call BSE_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,OmRPA(:,ispin),OmBSE(:,ispin), & 
-                                  XpY(:,:,ispin),XmY(:,:,ispin),rho(:,:,:,ispin))
+    call Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,OmRPA(:,ispin),OmBSE(:,ispin), & 
+                                             XpY(:,:,ispin),XmY(:,:,ispin),rho(:,:,:,ispin))
 
   end if
 
@@ -82,8 +82,8 @@ subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, &
 
     ! Compute dynamic correction for BSE via perturbation theory
 
-    call BSE_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,OmRPA(:,ispin),OmBSE(:,ispin), & 
-                                  XpY(:,:,ispin),XmY(:,:,ispin),rho(:,:,:,ispin))
+    call Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,OmRPA(:,ispin),OmBSE(:,ispin), & 
+                                             XpY(:,:,ispin),XmY(:,:,ispin),rho(:,:,:,ispin))
 
 
   end if
diff --git a/src/QuAcK/Bethe_Salpeter_A_matrix_dynamic.f90 b/src/QuAcK/Bethe_Salpeter_A_matrix_dynamic.f90
index b456b66..e8d8310 100644
--- a/src/QuAcK/Bethe_Salpeter_A_matrix_dynamic.f90
+++ b/src/QuAcK/Bethe_Salpeter_A_matrix_dynamic.f90
@@ -11,7 +11,7 @@ subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,OmRPA,
   double precision,intent(in)   :: eta
   double precision,intent(in)   :: lambda
   double precision,intent(in)   :: OmRPA(nS)
-  double precision,intent(in)   :: OmBSE(nS)
+  double precision,intent(in)   :: OmBSE
   double precision,intent(in)   :: rho(nBas,nBas,nS)
   
 ! Local variables
@@ -48,11 +48,11 @@ subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,OmRPA,
           chi = 0d0
           do kc=1,nS
 
-            eps = (OmBSE(kc) - OmRPA(kc))**2 + eta**2
-            chi = chi + rho(i,j,kc)*rho(a,b,kc)*(OmBSE(kc) - OmRPA(kc))/eps
+            eps = (OmBSE - OmRPA(kc))**2 + eta**2
+            chi = chi + rho(i,j,kc)*rho(a,b,kc)*(OmBSE - OmRPA(kc))/eps
 
-            eps = (OmBSE(kc) + OmRPA(kc))**2 + eta**2
-            chi = chi - rho(i,j,kc)*rho(a,b,kc)*(OmBSE(kc) + OmRPA(kc))/eps
+            eps = (OmBSE + OmRPA(kc))**2 + eta**2
+            chi = chi - rho(i,j,kc)*rho(a,b,kc)*(OmBSE + OmRPA(kc))/eps
 
           enddo
 
diff --git a/src/QuAcK/Bethe_Salpeter_B_matrix_dynamic.f90 b/src/QuAcK/Bethe_Salpeter_B_matrix_dynamic.f90
index 8b46d84..d6950df 100644
--- a/src/QuAcK/Bethe_Salpeter_B_matrix_dynamic.f90
+++ b/src/QuAcK/Bethe_Salpeter_B_matrix_dynamic.f90
@@ -11,7 +11,7 @@ subroutine Bethe_Salpeter_B_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,OmRPA,
   double precision,intent(in)   :: eta
   double precision,intent(in)   :: lambda
   double precision,intent(in)   :: OmRPA(nS)
-  double precision,intent(in)   :: OmBSE(nS)
+  double precision,intent(in)   :: OmBSE
   double precision,intent(in)   :: rho(nBas,nBas,nS)
   
 ! Local variables
@@ -48,11 +48,11 @@ subroutine Bethe_Salpeter_B_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,OmRPA,
           chi = 0d0
           do kc=1,nS
 
-            eps = (OmBSE(kc) - OmRPA(kc))**2 + eta**2
-            chi = chi + rho(i,b,kc)*rho(a,j,kc)*(OmBSE(kc) - OmRPA(kc))/eps
+            eps = (OmBSE - OmRPA(kc))**2 + eta**2
+            chi = chi + rho(i,b,kc)*rho(a,j,kc)*(OmBSE - OmRPA(kc))/eps
 
-            eps = (OmBSE(kc) + OmRPA(kc))**2 + eta**2
-            chi = chi - rho(i,b,kc)*rho(a,j,kc)*(OmBSE(kc) + OmRPA(kc))/eps
+            eps = (OmBSE + OmRPA(kc))**2 + eta**2
+            chi = chi - rho(i,b,kc)*rho(a,j,kc)*(OmBSE + OmRPA(kc))/eps
 
 
           enddo
diff --git a/src/QuAcK/BSE_dynamic_perturbation.f90 b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90
similarity index 68%
rename from src/QuAcK/BSE_dynamic_perturbation.f90
rename to src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90
index 63dc325..befd12b 100644
--- a/src/QuAcK/BSE_dynamic_perturbation.f90
+++ b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90
@@ -1,4 +1,4 @@
-subroutine BSE_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,OmRPA,OmBSE,XpY,XmY,rho)
+subroutine Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,OmRPA,OmBSE,XpY,XmY,rho)
 
 ! Compute dynamical effects via perturbation theory for BSE
 
@@ -35,22 +35,27 @@ subroutine BSE_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,OmRPA,OmBSE,XpY,
   write(*,*) '-----------------------------------------------------------------------------------------------------------'
   write(*,'(2X,A5,1X,A30,1X,A30,1X,A30)') '#','Static excitation (eV)','Dynamic correction (eV)','Dynamic excitation (eV)'
   write(*,*) '-----------------------------------------------------------------------------------------------------------'
+  do ia=1,min(nS,maxS)
 
-  do ia=1,maxS
+    X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:))
+    Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:))
 
-    X(:) = 0.5d0*(XpY(:,ia) + XmY(:,ia))
-    Y(:) = 0.5d0*(XpY(:,ia) - XmY(:,ia))
+    call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,OmRPA(:),OmBSE(ia),rho(:,:,:),A_dyn(:,:))
 
-    call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,OmRPA(:),OmBSE(:),rho(:,:,:),A_dyn(:,:))
+  if(TDA) then 
 
-    if(TDA) then 
-      B_dyn(:,:) = 0d0
-    else
-      call Bethe_Salpeter_B_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,OmRPA(:),OmBSE(:),rho(:,:,:),B_dyn(:,:))
-    end if
+    OmDyn(ia) = dot_product(X(:),matmul(A_dyn(:,:),X(:)))
 
-    OmDyn(ia) = dot_product(X(:),matmul(A_dyn(:,:),X(:))) - dot_product(Y(:),matmul(A_dyn(:,:),Y(:))) &
-              + dot_product(X(:),matmul(B_dyn(:,:),Y(:))) - dot_product(Y(:),matmul(B_dyn(:,:),X(:)))  
+  else
+
+    call Bethe_Salpeter_B_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,OmRPA(:),OmBSE(ia),rho(:,:,:),B_dyn(:,:))
+
+    OmDyn(ia) = dot_product(X(:),matmul(A_dyn(:,:),X(:))) &
+              - dot_product(Y(:),matmul(A_dyn(:,:),Y(:))) &
+              + dot_product(X(:),matmul(B_dyn(:,:),Y(:))) & 
+              - dot_product(Y(:),matmul(B_dyn(:,:),X(:)))  
+
+  end if
 
     write(*,'(2X,I5,15X,F15.6,15X,F15.6,15X,F15.6)') ia,OmBSE(ia)*HaToeV,OmDyn(ia)*HaToeV,(OmBSE(ia)+OmDyn(ia))*HaToeV
 
@@ -58,4 +63,4 @@ subroutine BSE_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,OmRPA,OmBSE,XpY,
   write(*,*) '-----------------------------------------------------------------------------------------------------------'
   write(*,*) 
 
-end subroutine BSE_dynamic_perturbation
+end subroutine Bethe_Salpeter_dynamic_perturbation