diff --git a/src/casscf/bielec.irp.f b/src/casscf/bielec.irp.f
index 73c4cea7..2fca1a8c 100644
--- a/src/casscf/bielec.irp.f
+++ b/src/casscf/bielec.irp.f
@@ -59,9 +59,9 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_inact_orb+n_act_orb,n_core_i
     ii=list_core_inact(i)
     do j=i,n_core_inact_orb
       jj=list_core_inact(j)
-      call get_mo_two_e_integrals_ij  (ii,jj,mo_num,integrals_array,mo_integrals_map)
-      do p=1,mo_num
-        do q=1,mo_num
+      call get_mo_two_e_integrals_ij(ii,jj,mo_num,integrals_array,mo_integrals_map)
+      do q=1,mo_num
+        do p=1,mo_num
           bielec_PxxQ(p,i,j,q)=integrals_array(p,q)
           bielec_PxxQ(p,j,i,q)=integrals_array(q,p)
         end do
@@ -70,9 +70,9 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_inact_orb+n_act_orb,n_core_i
     do j=1,n_act_orb
       jj=list_act(j)
       j3=j+n_core_inact_orb
-      call get_mo_two_e_integrals_ij  (ii,jj,mo_num,integrals_array,mo_integrals_map)
-      do p=1,mo_num
-        do q=1,mo_num
+      call get_mo_two_e_integrals_ij(ii,jj,mo_num,integrals_array,mo_integrals_map)
+      do q=1,mo_num
+        do p=1,mo_num
           bielec_PxxQ(p,i,j3,q)=integrals_array(p,q)
           bielec_PxxQ(p,j3,i,q)=integrals_array(q,p)
         end do
@@ -88,9 +88,9 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_inact_orb+n_act_orb,n_core_i
     do j=i,n_act_orb
       jj=list_act(j)
       j3=j+n_core_inact_orb
-      call get_mo_two_e_integrals_ij  (ii,jj,mo_num,integrals_array,mo_integrals_map)
-      do p=1,mo_num
-        do q=1,mo_num
+      call get_mo_two_e_integrals_ij(ii,jj,mo_num,integrals_array,mo_integrals_map)
+      do q=1,mo_num
+        do p=1,mo_num
           bielec_PxxQ(p,i3,j3,q)=integrals_array(p,q)
           bielec_PxxQ(p,j3,i3,q)=integrals_array(q,p)
         end do
@@ -107,24 +107,19 @@ BEGIN_PROVIDER [real*8, bielecCI, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
   END_DOC
   implicit none
   integer                        :: i,j,k,p,t,u,v
-  double precision, allocatable  :: integrals_array(:)
-  real*8                         :: mo_two_e_integral
+  double precision, external     :: mo_two_e_integral
   
-  allocate(integrals_array(mo_num))
-  
-  do i=1,n_act_orb
-    t=list_act(i)
+  do p=1,mo_num
     do j=1,n_act_orb
       u=list_act(j)
       do k=1,n_act_orb
         v=list_act(k)
-        ! (tu|vp)
-        call get_mo_two_e_integrals(t,u,v,mo_num,integrals_array,mo_integrals_map)
-        do p=1,mo_num
-          bielecCI(i,k,j,p)=integrals_array(p)
+        do i=1,n_act_orb
+          t=list_act(i)
+          bielecCI(i,k,j,p) = mo_two_e_integral(t,u,v,p)
         end do
       end do
     end do
   end do
 END_PROVIDER
-
+  
diff --git a/src/casscf/bielec_natorb.irp.f b/src/casscf/bielec_natorb.irp.f
index 53d74e14..9826d80c 100644
--- a/src/casscf/bielec_natorb.irp.f
+++ b/src/casscf/bielec_natorb.irp.f
@@ -1,90 +1,114 @@
- BEGIN_PROVIDER [real*8, bielec_PQxx_no, (mo_num, mo_num,n_core_inact_orb+n_act_orb,n_core_inact_orb+n_act_orb)]
+ BEGIN_PROVIDER [real*8, bielec_PQxx_no, (mo_num, mo_num,n_core_inact_act_orb,n_core_inact_act_orb)]
   BEGIN_DOC
   ! integral (pq|xx) in the basis of natural MOs
   ! indices are unshifted orbital numbers
   END_DOC
   implicit none
   integer                        :: i,j,k,l,t,u,p,q,pp
-  real*8                         :: d(n_act_orb)
+  double precision, allocatable  :: f(:,:,:), d(:,:,:)
+
 
   bielec_PQxx_no(:,:,:,:) = bielec_PQxx(:,:,:,:)
 
-  do j=1,mo_num
-    do k=1,n_core_inact_orb+n_act_orb
-      do l=1,n_core_inact_orb+n_act_orb
+  allocate (f(n_act_orb,mo_num,n_core_inact_act_orb), &
+      d(n_act_orb,mo_num,n_core_inact_act_orb))
+
+  do l=1,n_core_inact_act_orb
+
+    do k=1,n_core_inact_act_orb
+      do j=1,mo_num
         do p=1,n_act_orb
-          d(p)=0.D0
+          f(p,j,k)=bielec_PQxx_no(list_act(p),j,k,l)
         end do
+      end do
+    end do
+    call dgemm('T','N',n_act_orb,mo_num*n_core_inact_act_orb,n_act_orb,1.d0,  &
+          natorbsCI, size(natorbsCI,1),                              &
+          f, n_act_orb,                                              &
+          0.d0,                                                      &
+          d, n_act_orb)
+    do k=1,n_core_inact_act_orb
+      do j=1,mo_num
         do p=1,n_act_orb
           pp=n_act_orb-p+1
-          do q=1,n_act_orb
-            d(pp)+=bielec_PQxx_no(list_act(q),j,k,l)*natorbsCI(q,p)
-          end do
+          bielec_PQxx_no(list_act(p),j,k,l)=d(pp,j,k)
         end do
+      end do
+
+      do j=1,mo_num
         do p=1,n_act_orb
-          bielec_PQxx_no(list_act(p),j,k,l)=d(p)
+          f(p,j,k)=bielec_PQxx_no(j,list_act(p),k,l)
+        end do
+      end do
+    end do
+    call dgemm('T','N',n_act_orb,mo_num*n_core_inact_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_core_inact_act_orb
+      do p=1,n_act_orb
+        pp=n_act_orb-p+1
+        do j=1,mo_num
+          bielec_PQxx_no(j,list_act(p),k,l)=d(pp,j,k)
         end do
       end do
     end do
   end do
-  ! 2nd quarter
-  do j=1,mo_num
-    do k=1,n_core_inact_orb+n_act_orb
-      do l=1,n_core_inact_orb+n_act_orb
-        do p=1,n_act_orb
-          d(p)=0.D0
+
+  deallocate (f,d)
+
+  allocate (f(mo_num,mo_num,n_act_orb),d(mo_num,mo_num,n_act_orb))
+
+  do l=1,n_core_inact_act_orb
+
+    do p=1,n_act_orb
+      do k=1,mo_num
+        do j=1,mo_num
+          f(j,k,p) = bielec_PQxx_no(j,k,n_core_inact_orb+p,l)
         end do
-        do p=1,n_act_orb
-          pp=n_act_orb-p+1
-          do q=1,n_act_orb
-            d(pp)+=bielec_PQxx_no(j,list_act(q),k,l)*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          bielec_PQxx_no(j,list_act(p),k,l)=d(p)
+      end do
+    end do
+    call dgemm('N','N',mo_num*mo_num,n_act_orb,n_act_orb,1.d0,       &
+          f, mo_num*mo_num,                                          &
+          natorbsCI, n_act_orb,                                      &
+          0.d0,                                                      &
+          d, mo_num*mo_num)
+    do p=1,n_act_orb
+      pp=n_act_orb-p+1
+      do k=1,mo_num
+        do j=1,mo_num
+          bielec_PQxx_no(j,k,n_core_inact_orb+p,l)=d(j,k,pp)
         end do
       end do
     end do
   end do
-  ! 3rd quarter
-  do j=1,mo_num
-    do k=1,mo_num
-      do l=1,n_core_inact_orb+n_act_orb
-        do p=1,n_act_orb
-          d(p)=0.D0
+
+  do l=1,n_core_inact_act_orb
+    do p=1,n_act_orb
+      do k=1,mo_num
+        do j=1,mo_num
+          f(j,k,p) = bielec_PQxx_no(j,k,l,n_core_inact_orb+p)
         end do
-        do p=1,n_act_orb
-          pp=n_act_orb-p+1
-          do q=1,n_act_orb
-            d(pp)+=bielec_PQxx_no(j,k,n_core_inact_orb+q,l)*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          bielec_PQxx_no(j,k,n_core_inact_orb+p,l)=d(p)
-        end do
-      end do
-    end do
-  end do
-  ! 4th quarter
-  do j=1,mo_num
-    do k=1,mo_num
-      do l=1,n_core_inact_orb+n_act_orb
-        do p=1,n_act_orb
-          d(p)=0.D0
-        end do
-        do p=1,n_act_orb
-          pp=n_act_orb-p+1
-          do q=1,n_act_orb
-            d(pp)+=bielec_PQxx_no(j,k,l,n_core_inact_orb+q)*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          bielec_PQxx_no(j,k,l,n_core_inact_orb+p)=d(p)
+      end do
+    end do
+    call dgemm('N','N',mo_num*mo_num,n_act_orb,n_act_orb,1.d0,       &
+          f, mo_num*mo_num,                                          &
+          natorbsCI, n_act_orb,                                      &
+          0.d0,                                                      &
+          d, mo_num*mo_num)
+    do p=1,n_act_orb
+      pp=n_act_orb-p+1
+      do k=1,mo_num
+        do j=1,mo_num
+          bielec_PQxx_no(j,k,l,n_core_inact_orb+p)=d(j,k,pp)
         end do
       end do
     end do
   end do
 
+  deallocate (f,d)
+
 END_PROVIDER