diff --git a/src/det.irp.f b/src/det.irp.f
index 2b5e54c..6bd995b 100644
--- a/src/det.irp.f
+++ b/src/det.irp.f
@@ -266,17 +266,22 @@ END_PROVIDER
  integer :: ik,il,jk,jl, idx(4)
  real :: phase
  integer :: exc(4), nact, nact2
- real :: det_kl
+ real :: det_kl, det_k
  integer :: det_exc
 
  two_e_density_num = 0
 
  PROVIDE det
 
+ print *,  'Computing two-electron DM'
  do k=1,det_num
+  det_k = det_coef(k)
+  if ( abs(det_k) < 1.e-5) then
+    cycle
+  endif
   do l=k,det_num
-   det_kl = det_coef(k)*det_coef(l)
-   if ( (k /= l).and.(abs(det_kl) < 1.e-5) ) then
+   det_kl = det_k*det_coef(l)
+   if ( (k /= l).or.(abs(det_kl) < 1.e-5) ) then
      cycle
    endif
 
@@ -298,9 +303,15 @@ END_PROVIDER
 BEGIN_SHELL [ /usr/bin/python ]
 code = """
         notfound = .True.
-        idx = (/ min(%(I)s,%(J)s), max(%(I)s,%(J)s), min(%(K)s,%(L)s), max(%(K)s,%(L)s) /)
+        idx(1) = min(%(I)s,%(J)s)
+        idx(2) = max(%(I)s,%(J)s)
+        idx(3) = min(%(K)s,%(L)s)
+        idx(4) = max(%(K)s,%(L)s)
         do q=1,two_e_density_num
-         if (sum(abs(two_e_density_indice(:,q)-idx))) then
+         if ( (two_e_density_indice(1,q)==idx(1)) .and. &
+              (two_e_density_indice(2,q)==idx(2)) .and. &
+              (two_e_density_indice(3,q)==idx(3)) .and. &
+              (two_e_density_indice(4,q)==idx(4)) ) then
           two_e_density_value(1,q) += det_kl
           two_e_density_value(2,q) += det_kl
           notfound = .False.
@@ -309,15 +320,24 @@ code = """
         enddo
         if (notfound) then
          two_e_density_num += 1
-         two_e_density_indice(:,two_e_density_num)=idx
+         two_e_density_indice(1,two_e_density_num)=idx(1)
+         two_e_density_indice(2,two_e_density_num)=idx(2)
+         two_e_density_indice(3,two_e_density_num)=idx(3)
+         two_e_density_indice(4,two_e_density_num)=idx(4)
          two_e_density_value(1,two_e_density_num) = det_kl
          two_e_density_value(2,two_e_density_num) = det_kl
         endif
 
         notfound = .True.
-        idx = (/ min(%(I)s,%(K)s), max(%(I)s,%(K)s), min(%(J)s,%(L)s), max(%(J)s,%(L)s) /)
+        idx(1) = min(%(I)s,%(K)s)
+        idx(2) = max(%(I)s,%(K)s)
+        idx(3) = min(%(J)s,%(L)s)
+        idx(4) = max(%(J)s,%(L)s)
         do q=1,two_e_density_num
-         if (sum(abs(two_e_density_indice(:,q)-idx))) then
+         if ( (two_e_density_indice(1,q)==idx(1)) .and. &
+              (two_e_density_indice(2,q)==idx(2)) .and. &
+              (two_e_density_indice(3,q)==idx(3)) .and. &
+              (two_e_density_indice(4,q)==idx(4)) ) then
           two_e_density_value(1,q) -= det_kl
           notfound = .False.
           exit
@@ -325,17 +345,26 @@ code = """
         enddo
         if (notfound) then
          two_e_density_num += 1
-         two_e_density_indice(:,two_e_density_num)=idx
+         two_e_density_indice(1,two_e_density_num)=idx(1)
+         two_e_density_indice(2,two_e_density_num)=idx(2)
+         two_e_density_indice(3,two_e_density_num)=idx(3)
+         two_e_density_indice(4,two_e_density_num)=idx(4)
          two_e_density_value(1,two_e_density_num) = -det_kl
          two_e_density_value(2,two_e_density_num) = 0.
         endif
 """
 
 code1 = """
-        idx = (/ min(%(I)s,%(J)s), max(%(I)s,%(J)s), min(%(K)s,%(L)s), max(%(K)s,%(L)s) /)
+        idx(1) = min(%(I)s,%(J)s)
+        idx(2) = max(%(I)s,%(J)s)
+        idx(3) = min(%(K)s,%(L)s)
+        idx(4) = max(%(K)s,%(L)s)
         notfound = .True.
         do q=1,two_e_density_num
-         if (sum(abs(two_e_density_indice(:,q)-idx))) then
+         if ( (two_e_density_indice(1,q)==idx(1)) .and. &
+              (two_e_density_indice(2,q)==idx(2)) .and. &
+              (two_e_density_indice(3,q)==idx(3)) .and. &
+              (two_e_density_indice(4,q)==idx(4)) ) then
           two_e_density_value(1,q) += det_kl
           notfound = .False.
           exit
@@ -343,15 +372,24 @@ code1 = """
         enddo
         if (notfound) then
          two_e_density_num += 1
-         two_e_density_indice(:,two_e_density_num)=idx
+         two_e_density_indice(1,two_e_density_num)=idx(1)
+         two_e_density_indice(2,two_e_density_num)=idx(2)
+         two_e_density_indice(3,two_e_density_num)=idx(3)
+         two_e_density_indice(4,two_e_density_num)=idx(4)
          two_e_density_value(1,two_e_density_num) = det_kl
          two_e_density_value(2,two_e_density_num) = 0.
         endif
 
         notfound = .True.
-        idx = (/ min(%(I)s,%(K)s), max(%(I)s,%(K)s), min(%(J)s,%(L)s), max(%(J)s,%(L)s) /)
+        idx(1) = min(%(I)s,%(K)s)
+        idx(2) = max(%(I)s,%(K)s)
+        idx(3) = min(%(J)s,%(L)s)
+        idx(4) = max(%(J)s,%(L)s)
         do q=1,two_e_density_num
-         if (sum(abs(two_e_density_indice(:,q)-idx))) then
+         if ( (two_e_density_indice(1,q)==idx(1)) .and. &
+              (two_e_density_indice(2,q)==idx(2)) .and. &
+              (two_e_density_indice(3,q)==idx(3)) .and. &
+              (two_e_density_indice(4,q)==idx(4)) ) then
           two_e_density_value(1,q) -= det_kl
           notfound = .False.
           exit
@@ -359,7 +397,10 @@ code1 = """
         enddo
         if (notfound) then
          two_e_density_num += 1
-         two_e_density_indice(:,two_e_density_num)=idx
+         two_e_density_indice(1,two_e_density_num)=idx(1)
+         two_e_density_indice(2,two_e_density_num)=idx(2)
+         two_e_density_indice(3,two_e_density_num)=idx(3)
+         two_e_density_indice(4,two_e_density_num)=idx(4)
          two_e_density_value(1,two_e_density_num) = -det_kl
          two_e_density_value(2,two_e_density_num) = 0.
         endif
@@ -367,9 +408,15 @@ code1 = """
 
 code2 = """
         notfound = .True.
-        idx = (/ min(%(I)s,%(J)s), max(%(I)s,%(J)s), min(%(K)s,%(L)s), max(%(K)s,%(L)s) /)
+        idx(1) = min(%(I)s,%(J)s)
+        idx(2) = max(%(I)s,%(J)s)
+        idx(3) = min(%(K)s,%(L)s)
+        idx(4) = max(%(K)s,%(L)s)
         do q=1,two_e_density_num
-         if (sum(abs(two_e_density_indice(:,q)-idx))) then
+         if ( (two_e_density_indice(1,q)==idx(1)) .and. &
+              (two_e_density_indice(2,q)==idx(2)) .and. &
+              (two_e_density_indice(3,q)==idx(3)) .and. &
+              (two_e_density_indice(4,q)==idx(4)) ) then
           two_e_density_value(2,q) += det_kl
           notfound = .False.
           exit
@@ -377,7 +424,10 @@ code2 = """
         enddo
         if (notfound) then
          two_e_density_num += 1
-         two_e_density_indice(:,two_e_density_num)=idx
+         two_e_density_indice(1,two_e_density_num)=idx(1)
+         two_e_density_indice(2,two_e_density_num)=idx(2)
+         two_e_density_indice(3,two_e_density_num)=idx(3)
+         two_e_density_indice(4,two_e_density_num)=idx(4)
          two_e_density_value(1,two_e_density_num) = 0.
          two_e_density_value(2,two_e_density_num) = det_kl
         endif
@@ -471,6 +521,7 @@ END_SHELL
  call set_density_matrix_two_num(two_e_density_num)
  call set_density_matrix_two_indice(two_e_density_indice)
  call set_density_matrix_two_value(two_e_density_value)
+ print *,  'Done'
 
 END_PROVIDER
 
@@ -490,6 +541,7 @@ BEGIN_PROVIDER [ real, one_e_density_mo, (mo_active_num,mo_active_num,2) ]
    return
  endif
 
+ print *,  'Computing one-electron DM'
  do p=1,2
   do i=1,mo_active_num
    do j=1,mo_active_num
@@ -532,5 +584,6 @@ BEGIN_PROVIDER [ real, one_e_density_mo, (mo_active_num,mo_active_num,2) ]
  enddo
 
  call set_density_matrix_one(one_e_density_mo)
+ print *,  'Done'
 
 END_PROVIDER