diff --git a/plugins/FourIdx/four_index.irp.f b/plugins/FourIdx/four_index.irp.f
index e9e6e9c0..fcdad326 100644
--- a/plugins/FourIdx/four_index.irp.f
+++ b/plugins/FourIdx/four_index.irp.f
@@ -30,7 +30,10 @@ subroutine four_index_transform(map_a,map_c,matrix_B,LDB,            &
   integer(key_kind), allocatable :: key(:)
   real(integral_kind), allocatable :: value(:)
 
-
+  ASSERT (k_start == i_start)
+  ASSERT (l_start == j_start)
+  ASSERT (a_start == c_start)
+  ASSERT (b_start == d_start)
 
   i_min = min(i_start,a_start)
   i_max = max(i_end  ,a_end  )
@@ -50,42 +53,40 @@ subroutine four_index_transform(map_a,map_c,matrix_B,LDB,            &
   ASSERT (LDB >= k_max)
   ASSERT (LDB >= l_max)
 
+  !$OMP PARALLEL DEFAULT(PRIVATE) SHARED(                            &
+      !$OMP  a_start,a_end,b_start,b_end,c_start,c_end,d_start,d_end,&
+      !$OMP  i_start,i_end,j_start,j_end,k_start,k_end,l_start,l_end,&
+      !$OMP  i_min,i_max,j_min,j_max,k_min,k_max,l_min,l_max,        &
+      !$OMP  map_a,map_c,matrix_B)
   allocate( key(i_max*j_max*k_max), value(i_max*j_max*k_max) )
   allocate( U(a_start:a_end, b_start:b_end, c_start:c_end) )
 
+  !$OMP DO
   do d=d_start,d_end
     U = 0.d0
     print *,  d
-    do l=l_start,l_end
+    do l=1,l_end
       if (dabs(matrix_B(l,d)) < 1.d-10) then
         cycle
       endif
 
-      allocate( T(i_start:i_end, k_start:k_end, j_start:j_end) )
+      allocate( T(i_start:i_end, k_start:k_end, j_start:j_end), &
+                V(a_start:a_end, k_start:k_end, j_start:j_end) )
 
       do k=k_start,k_end
         do j=j_start,j_end
-          do i=i_start,i_end
+          do i=i_start,k
             call bielec_integrals_index(i,j,k,l,idx)
             call map_get(map_a,idx,tmp)
             T(i, k,j) = tmp
+            T(k, i,j) = tmp
           enddo
         enddo
       enddo
 
-      allocate( V(a_start:a_end, k_start:k_end, j_start:j_end) )
 
-!      V = 0.d0
-!      do a=a_start,a_end
-!        do k=k_start,k_end
-!          do j=j_start,j_end
-!            do i=i_start,i_end
-!              V(a, k,j) = V(a, k,j) + T(i, k,j)*matrix_B(i, a)
-!            enddo
-!          enddo
-!        enddo
-!      enddo
-      call DGEMM('T','N', (a_end-a_start+1), (k_end-k_start+1)*(j_end-j_start+1),&
+      call DGEMM('T','N', (a_end-a_start+1),                         &
+          (k_end-k_start+1)*(j_end-j_start+1),                       &
           (i_end-i_start+1), 1.d0,                                   &
           matrix_B(i_start,a_start), size(matrix_B,1),               &
           T(i_start,k_start,j_start), size(T,1),  0.d0,              &
@@ -94,50 +95,16 @@ subroutine four_index_transform(map_a,map_c,matrix_B,LDB,            &
       deallocate(T)
       allocate( T(a_start:a_end, k_start:k_end, b_start:b_end) )
 
-!      V = 0.d0
-!      do a=a_start,a_end
-!        do k=k_start,k_end
-!          do b=b_start,b_end
-!            do j=j_start,j_end
-!              V(a,k, b) = V(a,k, b) + T(a,k, j)*matrix_B(j, b)
-!            enddo
-!          enddo
-!        enddo
-!      enddo
-      call DGEMM('N','N', (a_end-a_start+1)*(k_end-k_start+1),(b_end-b_start+1),&
-          (j_end-j_start+1), 1.d0,                                   &
-          V(a_start,k_start,j_start), size(V,1)*size(V,2),           &
-          matrix_B(j_start,b_start), size(matrix_B,1),0.d0,          &
-          T(a_start,k_start,b_start), size(T,1)*size(T,2) )
+      call DGEMM('N','N', (a_end-a_start+1)*(k_end-k_start+1),       &
+              (b_end-b_start+1),                                     &
+              (j_end-j_start+1), 1.d0,                               &
+              V(a_start,k_start,j_start), size(V,1)*size(V,2),       &
+              matrix_B(j_start,b_start), size(matrix_B,1),0.d0,      &
+              T(a_start,k_start,b_start), size(T,1)*size(T,2) )
 
       deallocate(V)
-      allocate( V(a_start:a_end, k_start:k_end, b_start:b_end) )
-      V = T
-      deallocate(T)
-      allocate( T(a_start:a_end, k_start:k_end, b_start:b_end) )
 
-      do b=b_start,b_end
-        do a=a_start,a_end
-          do k=k_start,k_end
-              T(a, k,b) = V(a, k,b) 
-          enddo
-        enddo
-      enddo
-
-      deallocate(V)
       allocate( V(a_start:a_end, b_start:b_end, c_start:c_end) )
-
-!      V = 0.d0
-!      do b=b_start,b_end
-!        do c=c_start,c_end
-!          do a=a_start,a_end
-!            do k=k_start,k_end
-!              V(a,b,c) = V(a,b,c) + T(a,k ,b)*matrix_B(k, c)
-!            enddo
-!          enddo
-!        enddo
-!      enddo
-
       V = 0.d0
       do b=b_start,b_end
         call DGEMM('N','N', (a_end-a_start+1), (c_end-c_start+1),    &
@@ -147,22 +114,8 @@ subroutine four_index_transform(map_a,map_c,matrix_B,LDB,            &
             V(a_start,c_start,b), size(V,1) )
       enddo
 
-
-      deallocate(T)
       U = U + V*matrix_B(l, d)
-
-!      do a=a_start,a_end
-!        do b=b_start,b_end
-!          do c=c_start,c_end
-!    do c=c_start,c_end
-!      do b=b_start,d
-!        do a=a_start,min(b,c)
-!            U(a,b,c) = U(a,b,c) + V(a,b,c) * matrix_B(l, d)
-!          enddo
-!        enddo
-!      enddo
-
-      deallocate(V)
+      deallocate(T,V)
 
     enddo
 
@@ -182,11 +135,16 @@ subroutine four_index_transform(map_a,map_c,matrix_B,LDB,            &
         enddo
       enddo
     enddo
+    !$OMP CRITICAL
     call map_append(map_c, key, value, idx) 
     call map_sort(map_c)
     call map_unique(map_c)
+    !$OMP END CRITICAL
 
   enddo
+  !$OMP END DO
+
   deallocate(key,value)
+  !$OMP END PARALLEL
 
 end