Merge branch 'master' of github.com:LCPQ/quantum_package

scripts/generate_h_apply.py

@@ -22,17 +22,17 @@ printout_always
 deinit_thread
 skip
 init_main
+filter_integrals
 """.split()
 
 class H_apply(object):
 
-  def __init__(self,sub,openmp=True):
+  def __init__(self,sub,SingleRef=False):
     s = {}
     for k in keywords:
       s[k] = ""
     s["subroutine"] = "H_apply_%s"%(sub)
     s["params_post"] = ""
-    self.openmp = openmp
 
     self.selection_pt2 = None
     self.perturbation = None
@@ -41,7 +41,7 @@ class H_apply(object):
     s["omp_parallel"]     = """!$OMP PARALLEL DEFAULT(SHARED)        &
         !$OMP PRIVATE(i,j,k,l,keys_out,hole,particle,                &
         !$OMP  occ_particle,occ_hole,j_a,k_a,other_spin,             &
-        !$OMP  hole_save,ispin,jj,l_a,                               &
+        !$OMP  hole_save,ispin,jj,l_a,ib_jb_pairs,array_pairs,       &
         !$OMP  accu,i_a,hole_tmp,particle_tmp,occ_particle_tmp,      &
         !$OMP  occ_hole_tmp,key_idx,i_b,j_b,key,N_elec_in_key_part_1,&
         !$OMP  N_elec_in_key_hole_1,N_elec_in_key_part_2,            &
@@ -58,6 +58,12 @@ class H_apply(object):
 
     s["keys_work"]  += "call fill_H_apply_buffer_no_selection(key_idx,keys_out,N_int,iproc)"
 
+    s["filter_integrals"] = "array_pairs = .True."
+    if SingleRef:
+      s["filter_integrals"] = """
+      call get_mo_bielec_integrals_existing_ik(i_a,j_a,mo_tot_num,array_pairs,mo_integrals_map)
+      """
+
     s["generate_psi_guess"]  = """
   ! Sort H_jj to find the N_states lowest states
   integer                        :: i
@@ -83,9 +89,6 @@ class H_apply(object):
   deallocate(H_jj,iorder)
     """
 
-    if not openmp:
-      for k in s:
-        s[k] = ""
     s["size_max"] = str(1024*128) 
     s["copy_buffer"] = "call copy_h_apply_buffer_to_wf"
     s["printout_now"]   = """write(output_Dets,*)  &
@@ -192,8 +195,7 @@ class H_apply(object):
          pt2_old(k) = pt2(k)
       enddo
       """
-      if self.openmp:
+      self.data["omp_parallel"]    += """&
-        self.data["omp_parallel"]    += """&
  !$OMP SHARED(N_st) PRIVATE(e_2_pert_buffer,coef_pert_buffer) &
  !$OMP PRIVATE(sum_e_2_pert, sum_norm_pert, sum_H_pert_diag)"""
 

src/BiInts/map_integrals.irp.f

@@ -232,6 +232,48 @@ subroutine get_mo_bielec_integrals(j,k,l,sze,out_val,map)
   call map_get_many(map, hash, out_val, sze)
 end
 
+subroutine get_mo_bielec_integrals_existing_ik(j,l,sze,out_array,map)
+  use map_module
+  implicit none
+  BEGIN_DOC
+  ! Returns multiple integrals <ij|kl> in the MO basis, all
+  ! i for j,k,l fixed.
+  END_DOC
+  integer, intent(in)            :: j,l, sze
+  logical, intent(out)           :: out_array(sze,sze)
+  type(map_type), intent(inout)  :: map
+  integer                        :: i,k,kk,ll,m
+  integer*8,allocatable          :: hash(:)
+  integer  ,allocatable          :: pairs(:,:), iorder(:)
+  PROVIDE mo_bielec_integrals_in_map
+  allocate (hash(sze*sze), pairs(2,sze*sze),iorder(sze*sze))
+  
+  kk=0
+  do k=1,sze
+   do i=1,sze
+    kk += 1
+    !DIR$ FORCEINLINE
+    call bielec_integrals_index(i,j,k,l,hash(kk))
+    pairs(1,kk) = i
+    pairs(2,kk) = k
+    iorder(kk) = kk
+   enddo
+  enddo
+
+  logical :: integral_is_in_map
+  call i8radix_sort(hash,iorder,kk,-1)
+  call map_exists_many(mo_integrals_map, hash, kk)
+
+  do ll=1,kk
+    m = iorder(ll)
+    i=pairs(1,m)
+    k=pairs(2,m)
+    out_array(i,k) = (hash(ll) /= 0_8)
+  enddo  
+
+  deallocate(pairs,hash,iorder)
+end
+
 integer*8 function get_mo_map_size()
   implicit none
   BEGIN_DOC

src/CISD/H_apply.irp.f

@@ -3,7 +3,7 @@
 
 BEGIN_SHELL [ /usr/bin/env python ]
 from generate_h_apply import H_apply
-H = H_apply("cisd",openmp=True)
+H = H_apply("cisd")
 print H
 END_SHELL
 

src/CISD_SC2_selected/H_apply.irp.f

@@ -3,7 +3,7 @@ BEGIN_SHELL [ /usr/bin/env python ]
 from generate_h_apply import *
 from perturbation import perturbations
 
-s = H_apply("PT2",openmp=True)
+s = H_apply("PT2",SingleRef=True)
 s.set_perturbation("epstein_nesbet_sc2_projected")
 print s
 END_SHELL

src/CISD_selected/H_apply.irp.f

@@ -4,7 +4,7 @@ from generate_h_apply import *
 from perturbation import perturbations
 
 for perturbation in perturbations:
-  s = H_apply("cisd_selection_"+perturbation,openmp=True)
+  s = H_apply("cisd_selection_"+perturbation)
   s.set_selection_pt2(perturbation)
   print s
 END_SHELL

src/Dets/H_apply_template.f

@@ -26,6 +26,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene
   
   double precision               :: mo_bielec_integral
   integer, allocatable           :: ia_ja_pairs(:,:,:)
+  integer, allocatable           :: ib_jb_pairs(:,:)
   double precision               :: diag_H_mat_elem
   integer                        :: iproc
   integer(omp_lock_kind), save   :: lck, ifirst=0
@@ -60,7 +61,8 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene
   call bitstring_to_list(particle(1,2),occ_particle(1,2),N_elec_in_key_part_1(2),N_int)
   call bitstring_to_list(hole(1,1),occ_hole(1,1),N_elec_in_key_hole_1(1),N_int)
   call bitstring_to_list(hole(1,2),occ_hole(1,2),N_elec_in_key_hole_1(2),N_int)
-  allocate (ia_ja_pairs(2,0:(elec_alpha_num)*mo_tot_num,2))
+  allocate (ia_ja_pairs(2,0:(elec_alpha_num)*mo_tot_num,2),          &
+            ib_jb_pairs(2,0:(elec_alpha_num)*mo_tot_num))
   
   do ispin=1,2
     i=0
@@ -86,6 +88,8 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene
   integer                        :: i_a,j_a,i_b,j_b,k_a,l_a,k_b,l_b
   integer(bit_kind)              :: test(N_int,2)
   double precision               :: accu
+  logical, allocatable           :: array_pairs(:,:)
+  allocate(array_pairs(mo_tot_num,mo_tot_num))
   accu = 0.d0
   do ispin=1,2
     other_spin = iand(ispin,1)+1
@@ -126,63 +130,41 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene
       
       !   hole = a^(+)_j_a(ispin) a_i_a(ispin)|key_in> : mono exc :: orb(i_a,ispin) --> orb(j_a,ispin)
       hole_save = hole
+
+      ! Build array of the non-zero integrals of second excitation
+      $filter_integrals
       if (ispin == 1) then
         integer                        :: jjj
         
+        i=0
         do kk = 1,N_elec_in_key_hole_2(other_spin)
-          hole = hole_save
           i_b = occ_hole_tmp(kk,other_spin)
           ASSERT (i_b > 0)
           ASSERT (i_b <= mo_tot_num)
-          k = ishft(i_b-1,-bit_kind_shift)+1
-          j = i_b-ishft(k-1,bit_kind_shift)-1
-          hole(k,other_spin) = ibclr(hole(k,other_spin),j)
           do jjj=1,N_elec_in_key_part_2(other_spin)     ! particule
             j_b = occ_particle_tmp(jjj,other_spin)
             ASSERT (j_b > 0)
             ASSERT (j_b <= mo_tot_num)
-            key = hole
+            if (array_pairs(i_b,j_b)) then
-            k = ishft(j_b-1,-bit_kind_shift)+1
+              i+= 1
-            l = j_b-ishft(k-1,bit_kind_shift)-1
+              ib_jb_pairs(1,i) = i_b
-            key(k,other_spin) = ibset(key(k,other_spin),l)
+              ib_jb_pairs(2,i) = j_b
-            key_idx += 1
-            do k=1,N_int
-              keys_out(k,1,key_idx) = key(k,1)
-              keys_out(k,2,key_idx) = key(k,2)
-            enddo
-            ASSERT (key_idx <= size_max)
-            if (key_idx == size_max) then
-              $keys_work
-              key_idx = 0
             endif
-            !           endif
           enddo
-          if (abort_here) then
-            exit
-          endif
         enddo
-      endif
+        ib_jb_pairs(1,0) = i
-      !   does all the mono excitations of the same spin
+
-      do kk = 1,N_elec_in_key_hole_2(ispin)
+        do kk = 1,ib_jb_pairs(1,0)
-        i_b = occ_hole_tmp(kk,ispin)
+          hole = hole_save
-        ASSERT (i_b > 0)
+          i_b = ib_jb_pairs(1,kk)
-        ASSERT (i_b <= mo_tot_num)
+          j_b = ib_jb_pairs(2,kk)
-        if (i_b <= i_a.or.i_b == j_a) cycle
+          k = ishft(i_b-1,-bit_kind_shift)+1
-        hole = hole_save
+          j = i_b-ishft(k-1,bit_kind_shift)-1
-        k = ishft(i_b-1,-bit_kind_shift)+1
+          hole(k,other_spin) = ibclr(hole(k,other_spin),j)
-        j = i_b-ishft(k-1,bit_kind_shift)-1
-        hole(k,ispin) = ibclr(hole(k,ispin),j)
-        do jjj=1,N_elec_in_key_part_2(ispin)
-          j_b = occ_particle_tmp(jjj,ispin)
-          ASSERT (j_b > 0)
-          ASSERT (j_b <= mo_tot_num)
-          if (j_b <= j_a) cycle
           key = hole
           k = ishft(j_b-1,-bit_kind_shift)+1
           l = j_b-ishft(k-1,bit_kind_shift)-1
-          key(k,ispin) = ibset(key(k,ispin),l)
+          key(k,other_spin) = ibset(key(k,other_spin),l)
-          !!   a^((+)_j_b(ispin) a_i_b(ispin) : mono exc :: orb(i_b,ispin) --> orb(j_b,ispin)
-          
           key_idx += 1
           do k=1,N_int
             keys_out(k,1,key_idx) = key(k,1)
@@ -197,18 +179,66 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene
             exit
           endif
         enddo
-      enddo! kk
+      endif
+
+      !   does all the mono excitations of the same spin
+      i=0
+      do kk = 1,N_elec_in_key_hole_2(ispin)
+        i_b = occ_hole_tmp(kk,ispin)
+        if (i_b <= i_a.or.i_b == j_a) cycle
+        ASSERT (i_b > 0)
+        ASSERT (i_b <= mo_tot_num)
+        do jjj=1,N_elec_in_key_part_2(ispin)     ! particule
+          j_b = occ_particle_tmp(jjj,ispin)
+          ASSERT (j_b > 0)
+          ASSERT (j_b <= mo_tot_num)
+          if (j_b <= j_a) cycle
+          if (array_pairs(i_b,j_b)) then
+            i+= 1
+            ib_jb_pairs(1,i) = i_b
+            ib_jb_pairs(2,i) = j_b
+          endif
+        enddo
+      enddo
+      ib_jb_pairs(1,0) = i
+
+      do kk = 1,ib_jb_pairs(1,0)
+        hole = hole_save
+        i_b = ib_jb_pairs(1,kk)
+        j_b = ib_jb_pairs(2,kk)
+        k = ishft(i_b-1,-bit_kind_shift)+1
+        j = i_b-ishft(k-1,bit_kind_shift)-1
+        hole(k,ispin) = ibclr(hole(k,ispin),j)
+        key = hole
+        k = ishft(j_b-1,-bit_kind_shift)+1
+        l = j_b-ishft(k-1,bit_kind_shift)-1
+        key(k,ispin) = ibset(key(k,ispin),l)
+        key_idx += 1
+        do k=1,N_int
+          keys_out(k,1,key_idx) = key(k,1)
+          keys_out(k,2,key_idx) = key(k,2)
+        enddo
+        ASSERT (key_idx <= size_max)
+        if (key_idx == size_max) then
+          $keys_work
+          key_idx = 0
+        endif
+        if (abort_here) then
+          exit
+        endif
+      enddo ! kk
+
     enddo  ! ii
     $omp_enddo
   enddo   ! ispin
   $keys_work
   $deinit_thread
-  deallocate (ia_ja_pairs, &
+  deallocate (ia_ja_pairs, ib_jb_pairs, &
       keys_out, hole_save,          &
       key,hole, particle, hole_tmp,&
       particle_tmp, occ_particle,    &
       occ_hole, occ_particle_tmp,&
-      occ_hole_tmp)
+      occ_hole_tmp,array_pairs)
   $omp_end_parallel
   $finalization
 end
@@ -235,11 +265,13 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator $parameters
   integer                        :: ii,i,jj,j,k,ispin,l
   integer,allocatable            :: occ_particle(:,:), occ_hole(:,:)
   integer,allocatable            :: occ_particle_tmp(:,:), occ_hole_tmp(:,:)
+  integer,allocatable            :: ib_jb_pairs(:,:)
   integer                        :: kk,pp,other_spin,key_idx
   integer                        :: N_elec_in_key_hole_1(2),N_elec_in_key_part_1(2)
   integer                        :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2)
   
   integer, allocatable           :: ia_ja_pairs(:,:,:)
+  logical, allocatable           :: array_pairs(:,:)
   double precision               :: diag_H_mat_elem
   integer                        :: iproc
   integer(omp_lock_kind), save   :: lck, ifirst=0

src/Dets/README.rst

@@ -83,7 +83,7 @@ Documentation
   .br
   Initial guess vectors are not necessarily orthonormal
 
-`repeat_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/SC2.irp.f#L215>`_
+`repeat_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/SC2.irp.f#L220>`_
   Undocumented
 
 `connected_to_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/connected_to_ref.irp.f#L95>`_
@@ -152,7 +152,7 @@ Documentation
 `davidson_threshold <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L374>`_
   Can be : [  energy  | residual | both | wall_time | cpu_time | iterations ]
 
-`det_search_key <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L180>`_
+`det_search_key <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L195>`_
   Return an integer*8 corresponding to a determinant index for searching
 
 `n_det <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L20>`_
@@ -164,35 +164,35 @@ Documentation
 `n_states <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L3>`_
   Number of states to consider
 
-`psi_average_norm_contrib <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L93>`_
+`psi_average_norm_contrib <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L108>`_
   Contribution of determinants to the state-averaged density
 
-`psi_average_norm_contrib_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L114>`_
+`psi_average_norm_contrib_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L129>`_
   Wave function sorted by determinants contribution to the norm (state-averaged)
 
-`psi_coef <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L65>`_
+`psi_coef <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L86>`_
-  The wave function. Initialized with Hartree-Fock if the EZFIO file
+  The wave function coefficients. Initialized with Hartree-Fock if the EZFIO file
   is empty
 
-`psi_coef_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L113>`_
+`psi_coef_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L128>`_
   Wave function sorted by determinants contribution to the norm (state-averaged)
 
-`psi_coef_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L144>`_
+`psi_coef_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L159>`_
   Determinants on which we apply <i|H|psi> for perturbation.
   o They are sorted by determinants interpreted as integers. Useful
   to accelerate the search of a determinant
 
 `psi_det <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L64>`_
-  The wave function. Initialized with Hartree-Fock if the EZFIO file
+  The wave function determinants. Initialized with Hartree-Fock if the EZFIO file
   is empty
 
 `psi_det_size <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L56>`_
   Size of the psi_det/psi_coef arrays
 
-`psi_det_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L112>`_
+`psi_det_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L127>`_
   Wave function sorted by determinants contribution to the norm (state-averaged)
 
-`psi_det_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L143>`_
+`psi_det_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L158>`_
   Determinants on which we apply <i|H|psi> for perturbation.
   o They are sorted by determinants interpreted as integers. Useful
   to accelerate the search of a determinant

src/Full_CI/H_apply.irp.f

@@ -2,11 +2,11 @@ use bitmasks
 BEGIN_SHELL [ /usr/bin/env python ]
 from generate_h_apply import *
 
-s = H_apply("FCI",openmp=True)
+s = H_apply("FCI")
 s.set_selection_pt2("epstein_nesbet_2x2")
 print s
 
-s = H_apply("FCI_PT2",openmp=True)
+s = H_apply("FCI_PT2")
 s.set_perturbation("epstein_nesbet_2x2")
 print s
 

src/MP2/H_apply.irp.f

@@ -3,7 +3,7 @@ BEGIN_SHELL [ /usr/bin/env python ]
 from generate_h_apply import *
 from perturbation import perturbations
 
-s = H_apply("mp2",openmp=True)
+s = H_apply("mp2")
 s.set_perturbation("Moller_plesset")
 print s
 END_SHELL

src/Utils/map_module.f90

@@ -495,7 +495,7 @@ subroutine map_get_many(map, key, value, sze)
  integer(key_kind), intent(in)     :: key(sze)
  real(integral_kind), intent(out)  :: value(sze)
  integer                           :: i
- integer(map_size_kind)            :: idx_cache, idx_cache_prev
+ integer(map_size_kind)            :: idx_cache
  integer(cache_map_size_kind)      :: ibegin, iend
  integer(cache_map_size_kind), allocatable :: idx(:)
  !DIR$ ATTRIBUTES ALIGN : 64 :: idx
@@ -518,6 +518,43 @@ subroutine map_get_many(map, key, value, sze)
  deallocate(idx)
 end
 
+subroutine map_exists_many(map, key, sze)
+ use map_module
+ implicit none
+ type (map_type), intent(inout) :: map
+ integer, intent(in)            :: sze
+ integer(key_kind), intent(inout):: key(sze)
+ integer                        :: i
+ integer(map_size_kind)         :: idx_cache, idx_cache_prev
+ integer(cache_map_size_kind)   :: ibegin, iend
+ integer(cache_map_size_kind), allocatable :: idx(:)
+ !DIR$ ATTRIBUTES ALIGN : 64 :: idx
+
+ idx_cache_prev = -1_map_size_kind
+ allocate(idx(sze))
+ do i=1,sze
+   idx_cache = ishft(key(i),map_shift)
+   iend = map%map(idx_cache)%n_elements
+   if (idx_cache == idx_cache_prev) then
+      if ((idx(i-1) > 0_cache_map_size_kind).and.(idx(i-1) < iend)) then
+        if ((key(i) == key(i-1)+1).and.(map%map(idx_cache)%key(idx(i-1))+1) == key(i)) then
+          idx(i) = idx(i-1)+1
+          cycle
+        endif
+      endif
+   endif
+   !DIR$ FORCEINLINE
+   call search_key_big_interval(key(i),map%map(idx_cache)%key, iend, idx(i), 1, iend)
+   idx_cache_prev = idx_cache
+ enddo
+ do i=1,sze
+   idx_cache = ishft(key(i),map_shift)
+   if (idx(i) <= 0) then
+     key(i) = 0_key_kind
+   endif
+ enddo
+ deallocate(idx)
+end
 
 subroutine search_key_big(key,X,sze,idx)
  use map_module