diff --git a/src/det.irp.f b/src/det.irp.f
index 85494f2..17410f0 100644
--- a/src/det.irp.f
+++ b/src/det.irp.f
@@ -1543,53 +1543,64 @@ END_PROVIDER
   DaC = 0.d0
   CDb = 0.d0
 
-  det_num4 = iand(det_num,not(3))
-  !DIR$ VECTOR ALIGNED
-  do k=1,det_num4,4
-    i1 = det_coef_matrix_rows(k  )
-    i2 = det_coef_matrix_rows(k+1)
-    i3 = det_coef_matrix_rows(k+2)
-    i4 = det_coef_matrix_rows(k+3)
-    j1 = det_coef_matrix_columns(k  )
-    j2 = det_coef_matrix_columns(k+1)
-    j3 = det_coef_matrix_columns(k+2)
-    j4 = det_coef_matrix_columns(k+3)
-    if ( (j1 == j2).and.(j1 == j3).and.(j1 == j4) ) then
-      f = det_beta_value (j1)
-      CDb(i1) = CDb(i1) + det_coef_matrix_values(k  )*f
-      CDb(i2) = CDb(i2) + det_coef_matrix_values(k+1)*f
-      CDb(i3) = CDb(i3) + det_coef_matrix_values(k+2)*f
-      CDb(i4) = CDb(i4) + det_coef_matrix_values(k+3)*f
+  if (det_num < ishft(det_alpha_num*det_beta_num,2)) then
 
-      if ( ((i2-i1) == 1).and.((i3-i1) == 2).and.((i4-i1) == 3) ) then
-        DaC(j1) = DaC(j1) + det_coef_matrix_values(k)*det_alpha_value(i1) &
-         + det_coef_matrix_values(k+1)*det_alpha_value(i1+1) &
-         + det_coef_matrix_values(k+2)*det_alpha_value(i1+2) &
-         + det_coef_matrix_values(k+3)*det_alpha_value(i1+3)
+    det_num4 = iand(det_num,not(3))
+    !DIR$ VECTOR ALIGNED
+    do k=1,det_num4,4
+      i1 = det_coef_matrix_rows(k  )
+      i2 = det_coef_matrix_rows(k+1)
+      i3 = det_coef_matrix_rows(k+2)
+      i4 = det_coef_matrix_rows(k+3)
+      j1 = det_coef_matrix_columns(k  )
+      j2 = det_coef_matrix_columns(k+1)
+      j3 = det_coef_matrix_columns(k+2)
+      j4 = det_coef_matrix_columns(k+3)
+      if ( (j1 == j2).and.(j1 == j3).and.(j1 == j4) ) then
+        f = det_beta_value (j1)
+        CDb(i1) = CDb(i1) + det_coef_matrix_values(k  )*f
+        CDb(i2) = CDb(i2) + det_coef_matrix_values(k+1)*f
+        CDb(i3) = CDb(i3) + det_coef_matrix_values(k+2)*f
+        CDb(i4) = CDb(i4) + det_coef_matrix_values(k+3)*f
+
+        if ( ((i2-i1) == 1).and.((i3-i1) == 2).and.((i4-i1) == 3) ) then
+          DaC(j1) = DaC(j1) + det_coef_matrix_values(k)*det_alpha_value(i1) &
+          + det_coef_matrix_values(k+1)*det_alpha_value(i1+1) &
+          + det_coef_matrix_values(k+2)*det_alpha_value(i1+2) &
+          + det_coef_matrix_values(k+3)*det_alpha_value(i1+3)
+        else
+          DaC(j1) = DaC(j1) + det_coef_matrix_values(k)*det_alpha_value(i1) &
+          + det_coef_matrix_values(k+1)*det_alpha_value(i2) &
+          + det_coef_matrix_values(k+2)*det_alpha_value(i3) &
+          + det_coef_matrix_values(k+3)*det_alpha_value(i4)
+        endif
       else
-        DaC(j1) = DaC(j1) + det_coef_matrix_values(k)*det_alpha_value(i1) &
-         + det_coef_matrix_values(k+1)*det_alpha_value(i2) &
-         + det_coef_matrix_values(k+2)*det_alpha_value(i3) &
-         + det_coef_matrix_values(k+3)*det_alpha_value(i4)
+        DaC(j1) = DaC(j1) + det_coef_matrix_values(k  )*det_alpha_value(i1)
+        DaC(j2) = DaC(j2) + det_coef_matrix_values(k+1)*det_alpha_value(i2)
+        DaC(j3) = DaC(j3) + det_coef_matrix_values(k+2)*det_alpha_value(i3)
+        DaC(j4) = DaC(j4) + det_coef_matrix_values(k+3)*det_alpha_value(i4)
+        CDb(i1) = CDb(i1) + det_coef_matrix_values(k  )*det_beta_value (j1)
+        CDb(i2) = CDb(i2) + det_coef_matrix_values(k+1)*det_beta_value (j2)
+        CDb(i3) = CDb(i3) + det_coef_matrix_values(k+2)*det_beta_value (j3)
+        CDb(i4) = CDb(i4) + det_coef_matrix_values(k+3)*det_beta_value (j4)
       endif
-    else
-      DaC(j1) = DaC(j1) + det_coef_matrix_values(k  )*det_alpha_value(i1)
-      DaC(j2) = DaC(j2) + det_coef_matrix_values(k+1)*det_alpha_value(i2)
-      DaC(j3) = DaC(j3) + det_coef_matrix_values(k+2)*det_alpha_value(i3)
-      DaC(j4) = DaC(j4) + det_coef_matrix_values(k+3)*det_alpha_value(i4)
-      CDb(i1) = CDb(i1) + det_coef_matrix_values(k  )*det_beta_value (j1)
-      CDb(i2) = CDb(i2) + det_coef_matrix_values(k+1)*det_beta_value (j2)
-      CDb(i3) = CDb(i3) + det_coef_matrix_values(k+2)*det_beta_value (j3)
-      CDb(i4) = CDb(i4) + det_coef_matrix_values(k+3)*det_beta_value (j4)
-    endif
-  enddo
+    enddo
 
-  do k=det_num4+1,det_num
-    i = det_coef_matrix_rows(k)
-    j = det_coef_matrix_columns(k)
-    DaC(j) = DaC(j) + det_coef_matrix_values(k)*det_alpha_value(i)
-    CDb(i) = CDb(i) + det_coef_matrix_values(k)*det_beta_value (j)
-  enddo
+    do k=det_num4+1,det_num
+      i = det_coef_matrix_rows(k)
+      j = det_coef_matrix_columns(k)
+      DaC(j) = DaC(j) + det_coef_matrix_values(k)*det_alpha_value(i)
+      CDb(i) = CDb(i) + det_coef_matrix_values(k)*det_beta_value (j)
+    enddo
+
+  else
+
+    call dgemv('T',det_alpha_num,det_beta_num,1.d0,det_coef_matrix_dense, &
+      size(det_coef_matrix_dense,1), det_alpha_value, 1, 0.d0, DaC, 1)
+    call dgemv('N',det_alpha_num,det_beta_num,1.d0,det_coef_matrix_dense, &
+      size(det_coef_matrix_dense,1), det_beta_value, 1, 0.d0, CDb, 1)
+
+  endif
 
   ! Value
   ! -----
diff --git a/src/wf.irp.f b/src/wf.irp.f
index cee557b..5e117f5 100644
--- a/src/wf.irp.f
+++ b/src/wf.irp.f
@@ -80,6 +80,21 @@ END_PROVIDER
   deallocate(buffer)
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, det_coef_matrix_dense, (det_alpha_num, det_beta_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Dense version of det_coef_matrix
+ END_DOC
+ integer :: i,j,k
+ det_coef_matrix_dense = 0.d0
+ do k=1,det_num
+   i = det_coef_matrix_rows(k)
+   j = det_coef_matrix_columns(k)
+   det_coef_matrix_dense(i,j) = det_coef_matrix_values(k)
+ enddo
+END_PROVIDER
+
+
 BEGIN_PROVIDER [ integer, det_num ]
   implicit none
   BEGIN_DOC