From ab799d5acb878c77edeb9b644d6b48549a7c8e22 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 29 May 2014 03:40:06 +0200
Subject: [PATCH] Multistate corrected

---
 scripts/generate_h_apply.py      |  2 +-
 src/Full_CI/full_ci.irp.f        |  2 +-
 src/Perturbation/selection.irp.f | 44 ++++++++++++++++++++++----------
 3 files changed, 33 insertions(+), 15 deletions(-)

diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py
index 73848488..eeac4669 100755
--- a/scripts/generate_h_apply.py
+++ b/scripts/generate_h_apply.py
@@ -177,7 +177,7 @@ class H_apply(object):
       call copy_h_apply_buffer_to_wf
       selection_criterion_min = selection_criterion_min*0.1d0
       selection_criterion = selection_criterion_min
-      call remove_small_contributions
+      !call remove_small_contributions
       """
       self.data["keys_work"] = """
       e_2_pert_buffer = 0.d0
diff --git a/src/Full_CI/full_ci.irp.f b/src/Full_CI/full_ci.irp.f
index 7d563420..72eec40e 100644
--- a/src/Full_CI/full_ci.irp.f
+++ b/src/Full_CI/full_ci.irp.f
@@ -11,7 +11,7 @@ program cisd
   
   pt2 = 1.d0
   diag_algorithm = "Lapack"
-  do while (maxval(abs(pt2(1:N_st))) > 1.d-6)
+  do while (maxval(abs(pt2(1:N_st))) > 1.d-4)
     call H_apply_FCI(pt2, norm_pert, H_pert_diag,  N_st)
     call diagonalize_CI
     print *,  'N_det    = ', N_det
diff --git a/src/Perturbation/selection.irp.f b/src/Perturbation/selection.irp.f
index 400d6976..203dccec 100644
--- a/src/Perturbation/selection.irp.f
+++ b/src/Perturbation/selection.irp.f
@@ -33,10 +33,9 @@ subroutine fill_H_apply_buffer_selection(n_selected,det_buffer,e_2_pert_buffer,c
 
     is_selected = .False.
     do j=1,N_st
-      s = dabs(e_2_pert_buffer(j,i))
+      s = -e_2_pert_buffer(j,i)
       is_selected = s > selection_criterion*selection_criterion_factor .or. is_selected
     enddo
-    ASSERT (s>=-1.d-8)
     
      
     if (is_selected) then
@@ -60,8 +59,10 @@ subroutine fill_H_apply_buffer_selection(n_selected,det_buffer,e_2_pert_buffer,c
     ASSERT (sum(popcnt(h_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num)
     ASSERT (sum(popcnt(h_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num)
   enddo
-  selection_criterion = smax
-  selection_criterion_min = smin
+  !$OMP CRITICAL
+  selection_criterion = max(selection_criterion,smax)
+  selection_criterion_min = min(selection_criterion_min,smin)
+  !$OMP END CRITICAL
 end
 
  BEGIN_PROVIDER [ double precision, selection_criterion ]
@@ -84,17 +85,32 @@ subroutine remove_small_contributions
 !  provided.
   END_DOC
   integer :: i,j,k, N_removed
-  logical keep
+  logical, allocatable :: keep(:)
   double precision :: i_H_psi_array(N_states)
-  k = 0
-  N_removed = 0
-  do i=N_det, 50
-    call i_H_psi(psi_det_sorted(1,1,i),psi_det_sorted,psi_coef_sorted,N_int,N_det,psi_det_size,N_states,i_H_psi_array)
-    keep = .False.
+  allocate (keep(N_det))
+  call diagonalize_CI
+  do i=1,N_det
+    keep(i) = .True.
+  enddo
+  !$OMP PARALLEL DEFAULT(NONE) &
+  !$OMP PRIVATE(i,j,i_H_psi_array) &
+  !$OMP SHARED(k,psi_det_sorted,psi_coef_sorted,N_int,N_det,psi_det_size,N_states, &
+  !$OMP   selection_criterion_min,keep,N_det_generators) &
+  !$OMP REDUCTION(+:N_removed)
+  !$OMP DO
+  do i=2*N_det_generators+1, N_det
+    call i_H_psi(psi_det_sorted(1,1,i),psi_det_sorted,psi_coef_sorted,N_int,min(N_det,2*N_det_generators),psi_det_size,N_states,i_H_psi_array)
+    keep(i) = .False.
     do j=1,N_states
-      keep = keep .or. (dabs(psi_coef_sorted(i,j)*i_H_psi_array(j)) > selection_criterion_min)
+      keep(i) = keep(i) .or. (-(psi_coef_sorted(i,j)*i_H_psi_array(j)) > selection_criterion_min)
     enddo
-    if (keep) then
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+  N_removed = 0
+  k = 0
+  do i=1, N_det
+    if (keep(i)) then
       k += 1
       do j=1,N_int
          psi_det(j,1,k) = psi_det_sorted(j,1,i)
@@ -107,10 +123,12 @@ subroutine remove_small_contributions
       N_removed += 1
     endif
   enddo
+  deallocate(keep)
   if (N_removed > 0) then
+    N_det = N_det - N_removed
+    SOFT_TOUCH N_det psi_det psi_coef 
     call write_int(output_dets,N_removed, 'Removed determinants')
   endif
-  SOFT_TOUCH N_det psi_det psi_coef
 end