diff --git a/src/casscf/bielec_natorb.irp.f b/src/casscf/bielec_natorb.irp.f
index 07836591..cb09be3e 100644
--- a/src/casscf/bielec_natorb.irp.f
+++ b/src/casscf/bielec_natorb.irp.f
@@ -8,12 +8,18 @@
   double precision, allocatable  :: f(:,:,:), d(:,:,:)
 
 
-  bielec_PQxx_no(:,:,:,:) = bielec_PQxx(:,:,:,:)
+
+  !$OMP PARALLEL DEFAULT(NONE) &
+  !$OMP PRIVATE(j,k,l,p,pp,d,f) &
+  !$OMP SHARED(n_core_inact_act_orb,mo_num,n_act_orb,n_core_inact_orb, &
+  !$OMP   bielec_PQxx_no,bielec_PQxx,list_act,natorbsCI)
 
   allocate (f(n_act_orb,mo_num,n_core_inact_act_orb), &
       d(n_act_orb,mo_num,n_core_inact_act_orb))
 
+  !$OMP DO
   do l=1,n_core_inact_act_orb
+    bielec_PQxx_no(:,:,:,l) = bielec_PQxx(:,:,:,l)
 
     do k=1,n_core_inact_act_orb
       do j=1,mo_num
@@ -55,11 +61,13 @@
       end do
     end do
   end do
+  !$OMP END DO NOWAIT
 
   deallocate (f,d)
 
   allocate (f(mo_num,mo_num,n_act_orb),d(mo_num,mo_num,n_act_orb))
 
+  !$OMP DO 
   do l=1,n_core_inact_act_orb
 
     do p=1,n_act_orb
@@ -83,7 +91,11 @@
       end do
     end do
   end do
+  !$OMP END DO NOWAIT
 
+  !$OMP BARRIER 
+
+  !$OMP DO
   do l=1,n_core_inact_act_orb
     do p=1,n_act_orb
       do k=1,mo_num
@@ -106,8 +118,10 @@
       end do
     end do
   end do
+  !$OMP END DO
 
   deallocate (f,d)
+  !$OMP END PARALLEL
 
 END_PROVIDER
 
@@ -122,12 +136,18 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac
   integer                        :: i,j,k,l,t,u,p,q,pp
   double precision, allocatable  :: f(:,:,:), d(:,:,:)
 
-  bielec_PxxQ_no(:,:,:,:) = bielec_PxxQ(:,:,:,:)
+  !$OMP PARALLEL DEFAULT(NONE) &
+  !$OMP PRIVATE(j,k,l,p,pp,d,f) &
+  !$OMP SHARED(n_core_inact_act_orb,mo_num,n_act_orb,n_core_inact_orb, &
+  !$OMP   bielec_PxxQ_no,bielec_PxxQ,list_act,natorbsCI)
+
 
   allocate (f(n_act_orb,n_core_inact_act_orb,n_core_inact_act_orb), &
       d(n_act_orb,n_core_inact_act_orb,n_core_inact_act_orb))
 
+  !$OMP DO
   do j=1,mo_num
+    bielec_PxxQ_no(:,:,:,j) = bielec_PxxQ(:,:,:,j)
     do l=1,n_core_inact_act_orb
       do k=1,n_core_inact_act_orb
         do p=1,n_act_orb
@@ -149,12 +169,14 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac
       end do
     end do
   end do
+  !$OMP END DO NOWAIT
 
   deallocate (f,d)
 
   allocate (f(n_act_orb,mo_num,n_core_inact_act_orb), &
     d(n_act_orb,mo_num,n_core_inact_act_orb))
 
+  !$OMP DO
   do k=1,mo_num
     do l=1,n_core_inact_act_orb
       do j=1,mo_num
@@ -177,12 +199,14 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac
       end do
     end do
   end do
+  !$OMP END DO NOWAIT
 
   deallocate(f,d)
 
   allocate(f(mo_num,n_core_inact_act_orb,n_act_orb), &
     d(mo_num,n_core_inact_act_orb,n_act_orb) )
 
+  !$OMP DO
   do k=1,mo_num
     do p=1,n_act_orb
       do l=1,n_core_inact_act_orb
@@ -205,8 +229,11 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac
       end do
     end do
   end do
+  !$OMP END DO NOWAIT
 
+  !$OMP BARRIER
 
+  !$OMP DO 
   do l=1,n_core_inact_act_orb
     do p=1,n_act_orb
       do k=1,n_core_inact_act_orb
@@ -229,6 +256,9 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac
       end do
     end do
   end do
+  !$OMP END DO NOWAIT
+  deallocate(f,d)
+  !$OMP END PARALLEL
 
 END_PROVIDER
 
@@ -242,12 +272,17 @@ BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
   integer                        :: i,j,k,l,t,u,p,q,pp
   double precision, allocatable  :: f(:,:,:), d(:,:,:)
   
-  bielecCI_no(:,:,:,:) = bielecCI(:,:,:,:)
+  !$OMP PARALLEL DEFAULT(NONE) &
+  !$OMP PRIVATE(j,k,l,p,pp,d,f) &
+  !$OMP SHARED(n_core_inact_act_orb,mo_num,n_act_orb,n_core_inact_orb, &
+  !$OMP   bielecCI_no,bielecCI,list_act,natorbsCI)
 
   allocate (f(n_act_orb,n_act_orb,mo_num), &
       d(n_act_orb,n_act_orb,mo_num))
 
+  !$OMP DO
   do l=1,mo_num
+    bielecCI_no(:,:,:,l) = bielecCI(:,:,:,l)
     do k=1,n_act_orb
       do j=1,n_act_orb
         do p=1,n_act_orb
@@ -267,110 +302,80 @@ BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
           bielecCI_no(p,j,k,l)=d(pp,j,k)
         end do
       end do
-    end do
-  end do
 
-  do j=1,n_act_orb
-    do k=1,n_act_orb
-      do l=1,mo_num
+      do j=1,n_act_orb
         do p=1,n_act_orb
-          d(p,1,1)=0.D0
-        end do
-        do p=1,n_act_orb
-          do q=1,n_act_orb
-            d(p,1,1)+=bielecCI_no(j,q,k,l)*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          pp=n_act_orb-p+1
-          bielecCI_no(j,p,k,l)=d(pp,1,1)
+          f(p,j,k)=bielecCI_no(j,p,k,l)
         end do
       end do
     end do
-  end do
-  ! 3rd quarter
-  do j=1,n_act_orb
+    call dgemm('T','N',n_act_orb,n_act_orb*n_act_orb,n_act_orb,1.d0,  &
+          natorbsCI, n_act_orb,                                      &
+          f, n_act_orb,                                              &
+          0.d0,                                                      &
+          d, n_act_orb)
     do k=1,n_act_orb
-      do l=1,mo_num
-        do p=1,n_act_orb
-          d(p,1,1)=0.D0
+      do p=1,n_act_orb
+        pp=n_act_orb-p+1
+        do j=1,n_act_orb
+          bielecCI_no(j,p,k,l)=d(pp,j,k)
         end do
-        do p=1,n_act_orb
-          do q=1,n_act_orb
-            d(p,1,1)+=bielecCI_no(j,k,q,l)*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          pp=n_act_orb-p+1
-          bielecCI_no(j,k,p,l)=d(pp,1,1)
-        end do
-      end do 
-    end do  
-  end do  
+      end do
+    end do
 
-!  do l=1,mo_num
-!    do k=1,n_act_orb
-!      do p=1,n_act_orb
-!        do j=1,n_act_orb
-!          f(j,p,k)=bielecCI_no(j,p,k,l)
-!        end do
-!      end do
-!    end do
-!    call dgemm('T','N',n_act_orb,n_act_orb*n_act_orb,n_act_orb,1.d0,  &
-!          natorbsCI, n_act_orb,                                      &
-!          f, n_act_orb,                                              &
-!          0.d0,                                                      &
-!          d, n_act_orb)
-!    do k=1,n_act_orb
-!      do p=1,n_act_orb
-!        pp=n_act_orb-p+1
-!        do j=1,n_act_orb
-!          bielecCI_no(j,p,k,l)=d(j,pp,k)
-!        end do
-!      end do
-!    end do
-!
-!    do p=1,n_act_orb
-!      do k=1,n_act_orb
-!        do j=1,n_act_orb
-!          f(j,k,p)=bielecCI_no(j,k,p,l)
-!        end do
-!      end do
-!    end do
-!    call dgemm('N','N',n_act_orb*n_act_orb,n_act_orb,n_act_orb,1.d0,       &
-!          f, n_act_orb*n_act_orb,                                          &
-!          natorbsCI, n_act_orb,                                      &
-!          0.d0,                                                      &
-!          d, n_act_orb*n_act_orb)
-!
-!    do p=1,n_act_orb
-!      pp=n_act_orb-p+1
-!      do k=1,n_act_orb
-!        do j=1,n_act_orb
-!          bielecCI_no(j,k,p,l)=d(j,k,pp)
-!        end do
-!      end do 
-!    end do  
-!  end do  
-!
-  do l=1,n_act_orb
-    do k=1,n_act_orb
-      do j=1,n_act_orb
-        do p=1,n_act_orb
-          d(p,1,1)=0.D0
+    do p=1,n_act_orb
+      do k=1,n_act_orb
+        do j=1,n_act_orb
+          f(j,k,p)=bielecCI_no(j,k,p,l)
         end do
-        do p=1,n_act_orb
-          do q=1,n_act_orb
-            d(p,1,1)+=bielecCI_no(j,k,l,list_act(q))*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          pp=n_act_orb-p+1
-          bielecCI_no(j,k,l,list_act(p))=d(pp,1,1)
+      end do
+    end do
+    call dgemm('N','N',n_act_orb*n_act_orb,n_act_orb,n_act_orb,1.d0,       &
+          f, n_act_orb*n_act_orb,                                          &
+          natorbsCI, n_act_orb,                                      &
+          0.d0,                                                      &
+          d, n_act_orb*n_act_orb)
+
+    do p=1,n_act_orb
+      pp=n_act_orb-p+1
+      do k=1,n_act_orb
+        do j=1,n_act_orb
+          bielecCI_no(j,k,p,l)=d(j,k,pp)
         end do
       end do 
     end do  
   end do  
+  !$OMP END DO
+
+  !$OMP DO
+  do l=1,n_act_orb
+    do p=1,n_act_orb
+      do k=1,n_act_orb
+        do j=1,n_act_orb
+          f(j,k,p)=bielecCI_no(j,k,l,list_act(p))
+        end do
+      end do
+    end do
+    call dgemm('N','N',n_act_orb*n_act_orb,n_act_orb,n_act_orb,1.d0,       &
+          f, n_act_orb*n_act_orb,                                          &
+          natorbsCI, n_act_orb,                                      &
+          0.d0,                                                      &
+          d, n_act_orb*n_act_orb)
+
+    do p=1,n_act_orb
+      pp=n_act_orb-p+1
+      do k=1,n_act_orb
+        do j=1,n_act_orb
+          bielecCI_no(j,k,l,list_act(p))=d(j,k,pp)
+        end do
+      end do 
+    end do  
+  end do  
+  !$OMP END DO
+  
+  deallocate(d,f)
+  !$OMP END PARALLEL
+
 
 END_PROVIDER