diff --git a/input/methods b/input/methods
index 33c74d4..d9f3262 100644
--- a/input/methods
+++ b/input/methods
@@ -13,7 +13,7 @@
 # G0F2* evGF2* qsGF2* G0F3 evGF3
   F     F      F      F    F
 # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW
-  F     F     F     F        F      F
+  T     F     F     F        F      F
 # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh
-  F      F      F      F      F      F
+  T      F      F      F      F      F
 # * unrestricted version available
diff --git a/input/options b/input/options
index c978a04..e8d6778 100644
--- a/input/options
+++ b/input/options
@@ -15,4 +15,4 @@
 # ACFDT: AC Kx  XBS
          F  T   T
 # BSE: BSE dBSE dTDA evDyn ppBSE BSE2
-        F    T    T    F   F     F
+        T    T    T    F   F     F
diff --git a/src/GT/GTpp_excitation_density.f90 b/src/GT/GTpp_excitation_density.f90
index 14e7ed5..897af3e 100644
--- a/src/GT/GTpp_excitation_density.f90
+++ b/src/GT/GTpp_excitation_density.f90
@@ -44,8 +44,14 @@ subroutine GTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1
 
   if(ispin == 1) then
 
-    do p=nC+1,nBas-nR
-      do q=nC+1,nBas-nR
+    !$OMP PARALLEL &
+    !$OMP SHARED(nC,nBas,nR,nO,nVV,nOO,rho1,rho2,ERI,X1,Y1,X2,Y2) &
+    !$OMP PRIVATE(q,p,ab,cd,kl,ij) &
+    !$OMP DEFAULT(NONE)
+    !$OMP DO
+
+    do q=nC+1,nBas-nR
+      do p=nC+1,nBas-nR
 
         do ab=1,nVV
  
@@ -68,7 +74,9 @@ subroutine GTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1
           end do
  
         end do
- 
+    !$OMP END DO
+
+    !$OMP DO
         do ij=1,nOO
  
           cd = 0
@@ -93,7 +101,8 @@ subroutine GTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1
 
       end do
     end do
-
+  !$OMP END DO
+  !$OMP END PARALLEL
   end if
 
 !----------------------------------------------
@@ -102,8 +111,8 @@ subroutine GTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1
 
   if(ispin == 2 .or. ispin == 4) then
 
-    do p=nC+1,nBas-nR
-      do q=nC+1,nBas-nR
+    do q=nC+1,nBas-nR
+      do p=nC+1,nBas-nR
  
 !       do ab=1,nVV
         ab = 0
@@ -170,8 +179,14 @@ subroutine GTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1
 
   if(ispin == 3) then
 
-    do p=nC+1,nBas-nR
-      do q=nC+1,nBas-nR
+    !$OMP PARALLEL &
+    !$OMP SHARED(nC,nBas,nR,nO,nVV,nOO,rho1,rho2,ERI,X1,Y1,X2,Y2) &
+    !$OMP PRIVATE(q,p,ab,cd,kl,ij,c,d,k,l) &
+    !$OMP DEFAULT(NONE)
+    !$OMP DO
+
+    do q=nC+1,nBas-nR
+      do p=nC+1,nBas-nR
  
 !       do ab=1,nVV
         ab = 0
@@ -197,6 +212,9 @@ subroutine GTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1
  
         end do
         end do
+    !$OMP END DO 
+
+    !$OMP DO
  
 !       do ij=1,nOO
         ij = 0
@@ -225,6 +243,8 @@ subroutine GTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1
 
       end do
     end do
+  !$OMP END DO
+  !$OMP END PARALLEL
 
   end if
 
diff --git a/src/GT/dynamic_Tmatrix_A.f90 b/src/GT/dynamic_Tmatrix_A.f90
index 7476be2..65361d7 100644
--- a/src/GT/dynamic_Tmatrix_A.f90
+++ b/src/GT/dynamic_Tmatrix_A.f90
@@ -46,14 +46,16 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O
 
 ! Build dynamic A matrix
 
-  ia = 0
-  do i=nC+1,nO
-    do a=nO+1,nBas-nR
-      ia = ia + 1
-      jb = 0
-      do j=nC+1,nO
-        do b=nO+1,nBas-nR
-          jb = jb + 1
+  jb = 0
+!$omp parallel do default(private) shared(TA,ZA,Omega1,Omega2,OmBSE,eGT,rho1,rho2,nO,nBas,nVV,nOO,chi,eps,eta,nC,nR,lambda)
+  do j=nC+1,nO
+    do b=nO+1,nBas-nR
+      jb = (b-nO) + (j-1)*(nBas-nO) 
+
+      ia = 0
+      do i=nC+1,nO
+        do a=nO+1,nBas-nR
+          ia = (a-nO) + (i-1)*(nBas-nO) 
  
           chi = 0d0
 
@@ -88,4 +90,6 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O
     end do
   end do
 
+!$omp end parallel do
+
 end subroutine dynamic_Tmatrix_A
diff --git a/src/GT/static_Tmatrix_A.f90 b/src/GT/static_Tmatrix_A.f90
index 2d8917b..f4c42b2 100644
--- a/src/GT/static_Tmatrix_A.f90
+++ b/src/GT/static_Tmatrix_A.f90
@@ -34,14 +34,16 @@ subroutine static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,Omega1,rho1,O
 
   TA(:,:) = 0d0
 
-  ia = 0
-  do i=nC+1,nO
-    do a=nO+1,nBas-nR
-      ia = ia + 1
-      jb = 0
-      do j=nC+1,nO
-        do b=nO+1,nBas-nR
-          jb = jb + 1
+  jb = 0
+!$omp parallel do default(private) shared(TA,Omega1,Omega2,rho1,rho2,nO,nBas,nVV,nOO,chi,eps,eta,nC,nR,lambda)
+  do j=nC+1,nO
+    do b=nO+1,nBas-nR
+      jb = (b-nO) + (j-1)*(nBas-nO) 
+
+      ia = 0
+      do i=nC+1,nO
+        do a=nO+1,nBas-nR
+          ia = (a-nO) + (i-1)*(nBas-nO) 
 
           chi = 0d0
 
@@ -63,4 +65,6 @@ subroutine static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,Omega1,rho1,O
     enddo
   enddo
 
+!$omp end parallel do
+
 end subroutine static_Tmatrix_A