Fixed bugs of DDCI

scripts/generate_h_apply.py

@@ -24,6 +24,10 @@ skip
 init_main
 filter_integrals
 filter2h2p
+filterhole
+filterparticle
+do_double_excitations
+check_double_excitation
 """.split()
 
 class H_apply(object):
@@ -116,6 +120,21 @@ class H_apply(object):
       buffer = buffer.replace('$'+key, value)
     return buffer
 
+  def unset_double_excitations(self):
+    self["do_double_excitations"] = ".False."
+    self["check_double_excitation"] = """
+     check_double_excitation = .False.
+    """
+  def set_filter_holes(self):
+    self["filterhole"] = """
+     if(iand(ibset(0_bit_kind,j),hole(k,other_spin)).eq.0_bit_kind )cycle
+    """
+  def set_filter_particl(self):
+    self["filterparticle"] = """
+     if(iand(ibset(0_bit_kind,j_a),hole(k_a,other_spin)).eq.0_bit_kind )cycle
+    """
+
+
   def set_filter_2h_2p(self):
     self["filter2h2p"] = """
 !    ! DIR$ FORCEINLINE

src/Dets/H_apply_template.f

@@ -36,6 +36,11 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene
     ifirst=1
   endif
   
+  logical :: check_double_excitation 
+  check_double_excitation = .True.
+
+
+
   $initialization
   
   $omp_parallel
@@ -276,6 +281,12 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc $param
   logical, allocatable           :: array_pairs(:,:)
   double precision               :: diag_H_mat_elem
   integer(omp_lock_kind), save   :: lck, ifirst=0
+
+  logical :: check_double_excitation 
+  check_double_excitation = .True.
+  $check_double_excitation
+
+
   if (ifirst == 0) then
     ifirst=1
 !$    call omp_init_lock(lck)
@@ -333,9 +344,11 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc $param
       hole = key_in
       k = ishft(i_a-1,-bit_kind_shift)+1
       j = i_a-ishft(k-1,bit_kind_shift)-1
+  $filterhole
       hole(k,ispin) = ibclr(hole(k,ispin),j)
       k_a = ishft(j_a-1,-bit_kind_shift)+1
       l_a = j_a-ishft(k_a-1,bit_kind_shift)-1
+  $filterparticle
       hole(k_a,ispin) = ibset(hole(k_a,ispin),l_a)
       $filter2h2p
       key_idx += 1

src/Full_CI/H_apply.irp.f

@@ -11,22 +11,21 @@ s.set_perturbation("epstein_nesbet_2x2")
 print s
 
 
-if False:
-  s = H_apply("FCI_mono")
-  s.set_selection_pt2("epstein_nesbet_2x2")
-  s.unset_double_excitations()
-  print s
+s = H_apply("FCI_mono")
+s.set_selection_pt2("epstein_nesbet_2x2")
+s.unset_double_excitations()
+print s
 
 
-  s = H_apply("select_mono_delta_rho")
-  s.unset_double_excitations()
-  s.set_selection_pt2("delta_rho_one_point")
-  print s
+s = H_apply("select_mono_delta_rho")
+s.unset_double_excitations()
+s.set_selection_pt2("delta_rho_one_point")
+print s
 
-  s = H_apply("pt2_mono_delta_rho")
-  s.unset_double_excitations()
-  s.set_perturbation("delta_rho_one_point")
-  print s
+s = H_apply("pt2_mono_delta_rho")
+s.unset_double_excitations()
+s.set_perturbation("delta_rho_one_point")
+print s
 
 s = H_apply("select_mono_di_delta_rho")
 s.set_selection_pt2("delta_rho_one_point")

src/Utils/map_module.f90

@@ -11,6 +11,10 @@ module map_module
 ! as integer*2 and is found by applying the map_mask
 ! to the initial key. The element are found in the
 ! cache_map using a binary search
+!
+! When using the map_update subroutine to build the map,
+! the map_unique subroutine
+! should be called before getting data from the map.
 
  use omp_lib
  
@@ -433,6 +437,97 @@ call omp_unset_lock(map%lock)
 
 end
 
+subroutine map_update_verbose(map, key, value, sze, thr)
+  use map_module
+  implicit none
+  type (map_type), intent(inout) :: map
+  integer, intent(in)            :: sze
+  integer(key_kind), intent(inout) :: key(sze)
+  real(integral_kind), intent(inout) :: value(sze)
+  real(integral_kind), intent(in) :: thr
+  
+  integer                        :: i
+  integer(map_size_kind)         :: idx_cache, idx_cache_new
+  integer(cache_map_size_kind)   :: idx
+  integer                        :: sze2
+  integer(cache_key_kind)        :: cache_key
+  integer(map_size_kind)         :: n_elements_temp
+  type (cache_map_type)          :: local_map
+  logical                        :: map_sorted
+! do i = 1, sze
+!  print*,'value in map = ',value(i)
+! enddo
+  
+  sze2 = sze
+  map_sorted = .True.
+  
+  n_elements_temp = 0_8
+  n_elements_temp = n_elements_temp + 1_8
+  do while (sze2>0)
+    i=1
+    do while (i<=sze)
+      if (key(i) /= 0_8) then
+        idx_cache = ishft(key(i),map_shift)
+        if (omp_test_lock(map%map(idx_cache)%lock)) then
+          local_map%key => map%map(idx_cache)%key
+          local_map%value => map%map(idx_cache)%value
+          local_map%sorted = map%map(idx_cache)%sorted
+          local_map%map_size = map%map(idx_cache)%map_size
+          local_map%n_elements = map%map(idx_cache)%n_elements
+          do
+          !DIR$ FORCEINLINE
+          call search_key_big_interval(key(i),local_map%key, local_map%n_elements, idx, 1, local_map%n_elements)
+          if (idx > 0_8) then
+!           print*,'AHAAH'
+!           print*,'local_map%value(idx) = ',local_map%value(idx)
+            local_map%value(idx) = local_map%value(idx) + value(i)
+!           print*,'not a new value !'
+!           print*,'local_map%value(idx) = ',local_map%value(idx)
+          else
+            ! Assert that the map has a proper size
+            if (local_map%n_elements == local_map%map_size) then
+              call cache_map_unique(local_map)
+              call cache_map_reallocate(local_map, local_map%n_elements + local_map%n_elements)
+              call cache_map_shrink(local_map,thr)
+            endif
+            cache_key = iand(key(i),map_mask)
+            local_map%n_elements = local_map%n_elements + 1_8
+            local_map%value(local_map%n_elements) = value(i)
+!           print*,'new value !'
+            local_map%key(local_map%n_elements) = cache_key
+            local_map%sorted = .False.
+            n_elements_temp = n_elements_temp + 1_8
+          endif  ! idx > 0
+          key(i) = 0_8
+          i = i+1
+          sze2 = sze2-1
+          if (i>sze) then
+            i=1
+          endif
+          if ( (ishft(key(i),map_shift) /= idx_cache).or.(key(i)==0_8)) then
+            exit
+          endif
+        enddo
+        map%map(idx_cache)%key => local_map%key
+        map%map(idx_cache)%value => local_map%value
+        map%map(idx_cache)%sorted = local_map%sorted
+        map%map(idx_cache)%n_elements = local_map%n_elements
+        map%map(idx_cache)%map_size = local_map%map_size
+        map_sorted = map_sorted .and. local_map%sorted
+        call omp_unset_lock(map%map(idx_cache)%lock)
+      endif  ! omp_test_lock
+    else
+      i=i+1
+    endif  ! key = 0
+  enddo  ! i
+enddo  ! sze2 > 0
+call omp_set_lock(map%lock)
+map%n_elements = map%n_elements + n_elements_temp
+map%sorted = map%sorted .and. map_sorted
+call omp_unset_lock(map%lock)
+
+end
+
 subroutine map_append(map, key, value, sze)
   use map_module
   implicit none