From 29143cbe0a1c7091ae05646762b097a6bafc3dba Mon Sep 17 00:00:00 2001 From: Manu Date: Wed, 25 Mar 2015 12:06:50 +0100 Subject: [PATCH] Fixed bugs of DDCI --- scripts/generate_h_apply.py | 19 ++++++++ src/Dets/H_apply_template.f | 13 +++++ src/Full_CI/H_apply.irp.f | 25 +++++----- src/Utils/map_module.f90 | 95 +++++++++++++++++++++++++++++++++++++ 4 files changed, 139 insertions(+), 13 deletions(-) diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index 8bde1d11..0ac11353 100755 --- a/scripts/generate_h_apply.py +++ b/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 diff --git a/src/Dets/H_apply_template.f b/src/Dets/H_apply_template.f index 557a6325..f4486c36 100644 --- a/src/Dets/H_apply_template.f +++ b/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 diff --git a/src/Full_CI/H_apply.irp.f b/src/Full_CI/H_apply.irp.f index 213cc9fd..5c6fcdc7 100644 --- a/src/Full_CI/H_apply.irp.f +++ b/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") diff --git a/src/Utils/map_module.f90 b/src/Utils/map_module.f90 index 820d3aaf..ecff478f 100644 --- a/src/Utils/map_module.f90 +++ b/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