From f6ad2c1a29b11ecd4fe9d0915643d19f363c9317 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@users.sourceforge.net>
Date: Wed, 23 Mar 2011 16:25:48 +0100
Subject: [PATCH] Bug

---
 eplf.config               |  2 +-
 src/det.irp.f             | 20 +++++++++++++++++---
 src/ezfio_interface.irp.f |  2 +-
 3 files changed, 19 insertions(+), 5 deletions(-)

diff --git a/eplf.config b/eplf.config
index 3ee4734..5ff470c 100644
--- a/eplf.config
+++ b/eplf.config
@@ -44,7 +44,7 @@ grid
 density_matrix
  one                 real (mo_basis_mo_active_num,mo_basis_mo_active_num,2)
  two_num             integer
- two_indice          real (4,density_matrix_two_num)
+ two_indice          integer (4,density_matrix_two_num)
  two_value           real (2,density_matrix_two_num)
 
 grid_data
diff --git a/src/det.irp.f b/src/det.irp.f
index 83b2662..a2bcd2c 100644
--- a/src/det.irp.f
+++ b/src/det.irp.f
@@ -217,6 +217,9 @@ BEGIN_PROVIDER [ integer, two_e_density_num_max ]
  integer :: det_exc
  do k=1,det_num
   do l=k,det_num
+   if ( (k /= l).and.(abs(det_coef(k)*det_coef(l)) < 1.e-5) ) then
+     cycle
+   endif
    exc(1) = abs(det_exc(k,l,1))-1
    exc(2) = abs(det_exc(k,l,2))-1
    exc(3) = exc(1)+exc(2)
@@ -272,6 +275,10 @@ END_PROVIDER
 
  do k=1,det_num
   do l=k,det_num
+   det_kl = det_coef(k)*det_coef(l)
+   if ( (k /= l).and.(abs(det_kl) < 1.e-5) ) then
+     cycle
+   endif
 
    exc(1) = det_exc(k,l,1)
    exc(2) = det_exc(k,l,2)
@@ -282,7 +289,7 @@ END_PROVIDER
    exc(4) = exc(4)/abs(exc(4))
    phase = dble(exc(4))
 
-   det_kl = phase*det_coef(k)*det_coef(l)
+   det_kl *= phase
    if (k /= l) then
     det_kl += det_kl
    endif
@@ -461,9 +468,9 @@ END_SHELL
   enddo
  enddo
  
+ 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)
- call set_density_matrix_two_num(two_e_density_num)
 
 END_PROVIDER
 
@@ -474,6 +481,9 @@ BEGIN_PROVIDER [ real, one_e_density_mo, (mo_active_num,mo_active_num,2) ]
  END_DOC
  integer :: i,j,k,l,p, il, jl
 
+ if (mo_active_num == 0) then
+   return
+ endif
  one_e_density_mo(1,1,1) = -1.
  call get_density_matrix_one(one_e_density_mo)
  if (one_e_density_mo(1,1,1) /= -1.) then
@@ -491,6 +501,10 @@ BEGIN_PROVIDER [ real, one_e_density_mo, (mo_active_num,mo_active_num,2) ]
  integer :: exc(4), det_exc
  do k=1,det_num
   do l=k,det_num
+   ckl = det_coef(k)*det_coef(l)
+   if ( (k /= l).and.(abs(ckl) < 1.e-5) ) then
+     cycle
+   endif
    exc(1) = det_exc(k,l,1)
    exc(2) = det_exc(k,l,2)
    exc(4) = exc(1)*exc(2)
@@ -499,7 +513,7 @@ BEGIN_PROVIDER [ real, one_e_density_mo, (mo_active_num,mo_active_num,2) ]
    exc(3) = exc(1)+exc(2)
    exc(4) = exc(4)/abs(exc(4))
    phase = dble(exc(4))
-   ckl = det_coef(k)*det_coef(l)*phase
+   ckl *= phase
    do p=1,2
     if (exc(3) == 0) then
       do i=1,elec_num_2(p)-mo_closed_num
diff --git a/src/ezfio_interface.irp.f b/src/ezfio_interface.irp.f
index 05905d6..8f46c67 100644
--- a/src/ezfio_interface.irp.f
+++ b/src/ezfio_interface.irp.f
@@ -47,7 +47,7 @@ data = [ \
 ("grid_data_density_partition"     , "real"          , "(grid_x_num,grid_y_num,grid_z_num)" ),
 ("density_matrix_one"              , "real"          , "(mo_active_num,mo_active_num,2)"    ),
 ("density_matrix_two_num"          , "integer"       , ""                                   ),
-("density_matrix_two_indice"       , "real"          , "(4,two_e_density_num_max)"         ),
+("density_matrix_two_indice"       , "integer"       , "(4,two_e_density_num_max)"         ),
 ("density_matrix_two_value"        , "real"          , "(2,two_e_density_num_max)"         ),
 ]