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

scripts/generate_h_apply.py

@@ -20,6 +20,8 @@ init_thread
 printout_now
 printout_always
 deinit_thread
+skip
+init_main
 """.split()
 
 class H_apply(object):
@@ -201,6 +203,11 @@ class H_apply(object):
     self.set_perturbation(pert)
     self.selection_pt2 = pert
     if pert is not None:
+      self.data["parameters"] += ",select_max_out"
+      self.data["declarations"] += """
+      double precision, intent(inout) :: select_max_out"""
+
+      self.data["params_post"] += ", select_max(i_generator)"
       self.data["size_max"] = str(1024*128) 
       self.data["copy_buffer"] = """
       call copy_h_apply_buffer_to_wf
@@ -212,7 +219,24 @@ class H_apply(object):
       coef_pert_buffer = 0.d0
       """ + self.data["keys_work"]
       self.data["keys_work"] += """
-      call fill_H_apply_buffer_selection(key_idx,keys_out,e_2_pert_buffer,coef_pert_buffer,N_st,N_int,iproc)
+      call fill_H_apply_buffer_selection(key_idx,keys_out,e_2_pert_buffer, &
+        coef_pert_buffer,N_st,N_int,iproc,select_max_out) 
+      """
+      self.data["omp_parallel"]    += """&
+ !$OMP REDUCTION (max:select_max_out)"""
+      self.data["skip"] = """
+      if ((i_generator < size(select_max)).and. &
+        (select_max(i_generator) < selection_criterion_min*selection_criterion_factor)) then
+        !$ call omp_set_lock(lck)
+        do k=1,N_st
+          norm_psi(k) = norm_psi(k) + psi_coef(i_generator,k)*psi_coef(i_generator,k)
+          delta_pt2(k) = 0.d0
+          pt2_old(k) = 0.d0
+        enddo
+        !$ call omp_unset_lock(lck)
+        cycle
+      endif
+      select_max(i_generator) = 0.d0
       """
 
 

src/BiInts/ao_bi_integrals.irp.f

@@ -26,7 +26,6 @@ double precision function ao_bielec_integral(i,j,k,l)
   ao_bielec_integral = 0.d0
   double precision               :: thresh
   thresh = ao_integrals_threshold
-  thresh = 0.d0
   
   if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k)then
     do p = 1, 3

src/BiInts/mo_bi_integrals.irp.f

@@ -88,7 +88,6 @@ subroutine add_integrals_to_map(mask_ijkl)
   
   call wall_time(wall_1)
   call cpu_time(cpu_1)
-  mo_integrals_threshold = 0.d0
   
   !$OMP PARALLEL PRIVATE(l1,k1,j1,i1,i2,i3,i4,i,j,k,l,c, ii1,kmax,   &
       !$OMP  bielec_tmp_0_idx, bielec_tmp_0, bielec_tmp_1,bielec_tmp_2,bielec_tmp_3,&

src/Bitmask/README.rst

@@ -57,15 +57,40 @@ Documentation
 `full_ijkl_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L12>`_
   Bitmask to include all possible MOs
 
+`generators_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L91>`_
+  Bitmasks for generator determinants. (N_int, alpha/beta, hole/particle, generator).
+  3rd index is :
+  * 1 : hole     for single exc
+  * 1 : particle for single exc
+  * 3 : hole     for 1st exc of double
+  * 4 : particle for 1st exc of double
+  * 5 : hole     for 2dn exc of double
+  * 6 : particle for 2dn exc of double
+
 `hf_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L32>`_
   Hartree Fock bit mask
 
+`i_bitmask_gen <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L182>`_
+  Current bitmask for the generators
+
+`i_bitmask_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L190>`_
+  Current bitmask for the reference
+
+`n_generators_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L58>`_
+  Number of bitmasks for generators
+
 `n_int <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L3>`_
   Number of 64-bit integers needed to represent determinants as binary strings
 
+`n_reference_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L126>`_
+  Number of bitmasks for reference
+
 `ref_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L50>`_
   Reference bit mask, used in Slater rules, chosen as Hartree-Fock bitmask
 
+`reference_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L159>`_
+  Bitmasks for reference determinants. (N_int, alpha/beta, hole/particle, reference)
+
 `bitstring_to_hexa <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks_routines.irp.f#L95>`_
   Transform a bit string to a string in hexadecimal format for printing
 

src/Bitmask/bitmasks.irp.f

@@ -107,12 +107,17 @@ BEGIN_PROVIDER [ integer(bit_kind), generators_bitmask, (N_int,2,6,N_generators_
  if (exists) then
    call ezfio_get_bitmasks_generators(generators_bitmask)
  else
-   generators_bitmask(:,:,s_hole ,1) = HF_bitmask
+   integer :: k, ispin
-   generators_bitmask(:,:,s_part ,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:))
+   do k=1,N_generators_bitmask
-   generators_bitmask(:,:,d_hole1,1) = HF_bitmask
+     do ispin=1,2
-   generators_bitmask(:,:,d_part1,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:))
+       generators_bitmask(:,ispin,s_hole ,k) = full_ijkl_bitmask(:,d_hole1)
-   generators_bitmask(:,:,d_hole2,1) = HF_bitmask
+       generators_bitmask(:,ispin,s_part ,k) = full_ijkl_bitmask(:,d_part1)
-   generators_bitmask(:,:,d_part2,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:))
+       generators_bitmask(:,ispin,d_hole1,k) = full_ijkl_bitmask(:,d_hole1)
+       generators_bitmask(:,ispin,d_part1,k) = full_ijkl_bitmask(:,d_part1)
+       generators_bitmask(:,ispin,d_hole2,k) = full_ijkl_bitmask(:,d_hole2)
+       generators_bitmask(:,ispin,d_part2,k) = full_ijkl_bitmask(:,d_part2)
+     enddo
+   enddo
    call ezfio_set_bitmasks_generators(generators_bitmask)
  endif
 

src/Bitmask/bitmasks_module.f90

@@ -2,10 +2,10 @@ module bitmasks
   integer, parameter             :: bit_kind_shift = 6   ! 5: 32 bits, 6: 64 bits
   integer, parameter             :: bit_kind_size = 64
   integer, parameter             :: bit_kind = 64/8
-  integer, parameter             :: s_hole  = 1
+  integer, parameter             :: d_hole1 = 1
-  integer, parameter             :: s_part  = 2
+  integer, parameter             :: d_part1 = 2
-  integer, parameter             :: d_hole1 = 3
+  integer, parameter             :: d_hole2 = 3
-  integer, parameter             :: d_part1 = 4
+  integer, parameter             :: d_part2 = 4
-  integer, parameter             :: d_hole2 = 5
+  integer, parameter             :: s_hole  = 5
-  integer, parameter             :: d_part2 = 6
+  integer, parameter             :: s_part  = 6
 end module

src/Dets/H_apply.irp.f

@@ -19,7 +19,7 @@ BEGIN_PROVIDER [ logical, H_apply_buffer_allocated ]
   ! Uninitialized. Filled by H_apply subroutines.
   END_DOC
   integer                        :: iproc, sze
-  sze = 100
+  sze = 10000
   if (.not.associated(H_apply_buffer)) then
     allocate(H_apply_buffer(0:nproc-1))
     iproc = 0
@@ -237,4 +237,3 @@ subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc)
 end
 
 
-

src/Dets/H_apply_template.f

@@ -47,14 +47,6 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene
       occ_hole_tmp(N_int*bit_kind_size,2))
   $init_thread
   
-  !print*,'key_in !!'
-  !call print_key(key_in)
-  !print*,'hole_1, particl_1'
-  !call print_key(hole_1)
-  !call print_key(particl_1)
-  !print*,'hole_2, particl_2'
-  !call print_key(hole_2)
-  !call print_key(particl_2)
   
   
   !!!! First couple hole particle
@@ -352,9 +344,11 @@ subroutine $subroutine($params_main)
   
   $decls_main
   
-  integer                        :: i_generator, k, nmax
+  integer                        :: i_generator, nmax
-  double precision               :: wall_0, wall_1, wall_2, d
+  double precision               :: wall_0, wall_1, wall_2
   integer(omp_lock_kind)         :: lck
+  integer(bit_kind), allocatable :: mask(:,:,:)
+  integer                        :: ispin, k
 
   PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map N_det_selectors psi_generators
   PROVIDE psi_det_sorted_bit coef_hf_selector
@@ -363,21 +357,45 @@ subroutine $subroutine($params_main)
   call wall_time(wall_1)
   !$ call omp_init_lock(lck)
   !$OMP PARALLEL DEFAULT(SHARED) &
-  !$OMP PRIVATE(i_generator,wall_2) 
+  !$OMP PRIVATE(i_generator,wall_2,ispin,k,mask) 
+  allocate( mask(N_int,2,6) )
   !$OMP DO SCHEDULE(guided)
   do i_generator=1,nmax
     if (abort_here) then
       cycle
     endif
+    $skip
+
+    ! Create bit masks for holes and particles
+    do ispin=1,2
+      do k=1,N_int
+        mask(k,ispin,s_hole) =                                      &
+            iand(generators_bitmask(k,ispin,s_hole,i_bitmask_gen),  &
+            psi_generators(k,ispin,i_generator) )
+        mask(k,ispin,s_part) =                                      &
+            iand(generators_bitmask(k,ispin,s_part,i_bitmask_gen),  &
+            not(psi_generators(k,ispin,i_generator)) )
+        mask(k,ispin,d_hole1) =                                      &
+            iand(generators_bitmask(k,ispin,d_hole1,i_bitmask_gen),  &
+            psi_generators(k,ispin,i_generator) )
+        mask(k,ispin,d_part1) =                                      &
+            iand(generators_bitmask(k,ispin,d_part1,i_bitmask_gen),  &
+            not(psi_generators(k,ispin,i_generator)) )
+        mask(k,ispin,d_hole2) =                                      &
+            iand(generators_bitmask(k,ispin,d_hole2,i_bitmask_gen),  &
+            psi_generators(k,ispin,i_generator) )
+        mask(k,ispin,d_part2) =                                      &
+            iand(generators_bitmask(k,ispin,d_part2,i_bitmask_gen),  &
+            not (psi_generators(k,ispin,i_generator)) )
+      enddo
+    enddo
+
     call $subroutine_diexc(psi_generators(1,1,i_generator),          &
-        generators_bitmask(1,1,d_hole1,i_bitmask_gen),               &
+        mask(1,1,d_hole1), mask(1,1,d_part1),                        &
-        generators_bitmask(1,1,d_part1,i_bitmask_gen),               &
+        mask(1,1,d_hole2), mask(1,1,d_part2),                        &
-        generators_bitmask(1,1,d_hole2,i_bitmask_gen),               &
-        generators_bitmask(1,1,d_part2,i_bitmask_gen),               &
         i_generator $params_post)
     call $subroutine_monoexc(psi_generators(1,1,i_generator),        &
-        generators_bitmask(1,1,s_hole ,i_bitmask_gen),               &
+        mask(1,1,s_hole ), mask(1,1,s_part ),                        &
-        generators_bitmask(1,1,s_part ,i_bitmask_gen),               &
         i_generator $params_post)
     !$ call omp_set_lock(lck)
     call wall_time(wall_2)
@@ -388,23 +406,47 @@ subroutine $subroutine($params_main)
     endif
     !$ call omp_unset_lock(lck)
   enddo
-  !$OMP END DO NOWAIT
+  !$OMP END DO 
+  deallocate( mask )
   !$OMP END PARALLEL
   !$ call omp_destroy_lock(lck)
 
+  allocate( mask(N_int,2,6) )
   do i_generator=nmax+1,N_det_generators
     if (abort_here) then
       exit
     endif
+    $skip
+
+    ! Create bit masks for holes and particles
+    do ispin=1,2
+      do k=1,N_int
+        mask(k,ispin,s_hole) =                                      &
+            iand(generators_bitmask(k,ispin,s_hole,i_bitmask_gen),  &
+            psi_generators(k,ispin,i_generator) )
+        mask(k,ispin,s_part) =                                      &
+            iand(generators_bitmask(k,ispin,s_part,i_bitmask_gen),  &
+            not(psi_generators(k,ispin,i_generator)) )
+        mask(k,ispin,d_hole1) =                                      &
+            iand(generators_bitmask(k,ispin,d_hole1,i_bitmask_gen),  &
+            psi_generators(k,ispin,i_generator) )
+        mask(k,ispin,d_part1) =                                      &
+            iand(generators_bitmask(k,ispin,d_part1,i_bitmask_gen),  &
+            not(psi_generators(k,ispin,i_generator)) )
+        mask(k,ispin,d_hole2) =                                      &
+            iand(generators_bitmask(k,ispin,d_hole2,i_bitmask_gen),  &
+            psi_generators(k,ispin,i_generator) )
+        mask(k,ispin,d_part2) =                                      &
+            iand(generators_bitmask(k,ispin,d_part2,i_bitmask_gen),  &
+            not(psi_generators(k,ispin,i_generator)) )
+      enddo
+    enddo
     call $subroutine_diexc(psi_generators(1,1,i_generator),          &
-        generators_bitmask(1,1,d_hole1,i_bitmask_gen),               &
+        mask(1,1,d_hole1), mask(1,1,d_part1),                        &
-        generators_bitmask(1,1,d_part1,i_bitmask_gen),               &
+        mask(1,1,d_hole2), mask(1,1,d_part2),                        &
-        generators_bitmask(1,1,d_hole2,i_bitmask_gen),               &
-        generators_bitmask(1,1,d_part2,i_bitmask_gen),               &
         i_generator $params_post)
     call $subroutine_monoexc(psi_generators(1,1,i_generator),        &
-        generators_bitmask(1,1,s_hole ,i_bitmask_gen),               &
+        mask(1,1,s_hole ), mask(1,1,s_part ),                        &
-        generators_bitmask(1,1,s_part ,i_bitmask_gen),               &
         i_generator $params_post)
     call wall_time(wall_2)
     $printout_always
@@ -417,6 +459,7 @@ subroutine $subroutine($params_main)
   $copy_buffer
   $generate_psi_guess
   abort_here = abort_all
+  deallocate( mask )
   
 end
 

src/Dets/connected_to_ref.irp.f

@@ -11,8 +11,6 @@ logical function is_in_wavefunction(key,Nint,Ndet)
   is_in_wavefunction = .False.
   ibegin = 1
   iend   = N_det+1
-  ASSERT (N_past > 0)
-  ASSERT (N_det >= N_past)
   
   det_ref = det_search_key(key,Nint)
   det_search = det_search_key(psi_det_sorted_bit(1,1,1),Nint)
@@ -32,10 +30,6 @@ logical function is_in_wavefunction(key,Nint,Ndet)
     i = ibegin + istep 
   end do
 
-! if (det_search /= det_ref) then
-!   return
-! endif
-
   do while (det_search_key(psi_det_sorted_bit(1,1,i),Nint) == det_ref)
     i = i-1
     if (i == 0) then
@@ -43,6 +37,10 @@ logical function is_in_wavefunction(key,Nint,Ndet)
     endif
   enddo
   i += 1
+  if (i > N_det) then
+    return
+  endif
+
   do while (det_search_key(psi_det_sorted_bit(1,1,i),Nint) == det_ref)
     if ( (key(1,1) /= psi_det_sorted_bit(1,1,i)).or.                               &
           (key(1,2) /= psi_det_sorted_bit(1,2,i)) ) then
@@ -58,16 +56,40 @@ logical function is_in_wavefunction(key,Nint,Ndet)
         endif
       enddo
       if (is_in_wavefunction) then
-        return
+        exit
       endif
     endif
     i += 1
     if (i > N_det) then
-      exit
+      return
+ !    exit
     endif
     
   enddo
   
+
+! DEBUG is_in_wf
+! if (is_in_wavefunction) then
+!   degree = 1
+!   do i=1,N_det
+!     integer                        :: degree
+!     call get_excitation_degree(key,psi_det(1,1,i),degree,N_int)
+!     if (degree == 0) then
+!       exit
+!     endif
+!   enddo
+!   if (degree /=0) then
+!     stop 'pouet 1'
+!   endif
+! else
+!   do i=1,N_det
+!     call get_excitation_degree(key,psi_det(1,1,i),degree,N_int)
+!     if (degree == 0) then
+!       stop 'pouet 2'
+!     endif
+!   enddo
+! endif
+! END DEBUG is_in_wf
 end
 
 integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet)
@@ -94,14 +116,10 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet)
   N_past = max(1,N_past_in)
   if (Nint == 1) then
     
-    do i=N_past-1,1,-1
+    do i=1,N_past-1
       degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) +              &
           popcnt(xor( key(1,2), keys(1,2,i)))
-      if(degree_x2 == 0)then
+      if (degree_x2 > 4) then
-        connected_to_ref = -i
-        return
-      endif
-      if (degree_x2 > 5) then
         cycle
       else
         connected_to_ref = i
@@ -114,16 +132,12 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet)
     
   else if (Nint==2) then
     
-    do i=N_past-1,1,-1
+    do i=1,N_past-1
       degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) +              &
           popcnt(xor( key(1,2), keys(1,2,i))) +                      &
           popcnt(xor( key(2,1), keys(2,1,i))) +                      &
           popcnt(xor( key(2,2), keys(2,2,i)))
-      if(degree_x2 == 0)then
+      if (degree_x2 > 4) then
-        connected_to_ref = -i
-        return
-      endif
-      if (degree_x2 > 5) then
         cycle
       else
         connected_to_ref = i
@@ -135,18 +149,14 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet)
     
   else if (Nint==3) then
     
-    do i=N_past-1,1,-1
+    do i=1,N_past-1
       degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) +              &
           popcnt(xor( key(1,2), keys(1,2,i))) +                      &
           popcnt(xor( key(2,1), keys(2,1,i))) +                      &
           popcnt(xor( key(2,2), keys(2,2,i))) +                      &
           popcnt(xor( key(3,1), keys(3,1,i))) +                      &
           popcnt(xor( key(3,2), keys(3,2,i)))
-      if(degree_x2 == 0)then
+      if (degree_x2 > 4) then
-        connected_to_ref = -i
-        return
-      endif
-      if (degree_x2 > 5) then
         cycle
       else
         connected_to_ref = i
@@ -158,7 +168,7 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet)
     
   else
     
-    do i=N_past-1,1,-1
+    do i=1,N_past-1
       degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) +              &
           popcnt(xor( key(1,2), keys(1,2,i)))
       !DEC$ LOOP COUNT MIN(3)
@@ -166,11 +176,7 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet)
         degree_x2 = degree_x2 + popcnt(xor( key(l,1), keys(l,1,i))) +&
             popcnt(xor( key(l,2), keys(l,2,i)))
       enddo
-      if(degree_x2 == 0)then
+      if (degree_x2 > 4) then
-        connected_to_ref = -i
-        return
-      endif
-      if (degree_x2 > 5) then
         cycle
       else
         connected_to_ref = i

src/Full_CI/H_apply.irp.f

@@ -6,5 +6,9 @@ s = H_apply("FCI",openmp=True)
 s.set_selection_pt2("epstein_nesbet_2x2")
 print s
 
+s = H_apply("FCI_PT2",openmp=True)
+s.set_perturbation("epstein_nesbet_2x2")
+print s
+
 END_SHELL
 

src/Full_CI/full_ci.irp.f

@@ -5,13 +5,13 @@ program cisd
   
   double precision, allocatable  :: pt2(:), norm_pert(:), H_pert_diag(:)
   integer                        :: N_st, degree
-  character*(64)                 :: perturbation
   N_st = N_states
   allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st))
+  character*(64)                 :: perturbation
   
   pt2 = 1.d0
   diag_algorithm = "Lapack"
-  do while (maxval(abs(pt2(1:N_st))) > 1.d-2)
+  do while (maxval(abs(pt2(1:N_st))) > 1.d-3)
     call H_apply_FCI(pt2, norm_pert, H_pert_diag,  N_st)
     call diagonalize_CI
     print *,  'N_det    = ', N_det

src/Generators_full/generators.irp.f

@@ -33,8 +33,8 @@ BEGIN_PROVIDER [ integer, N_det_generators ]
  N_det_generators = N_det
  do i=1,N_det
    norm = norm + psi_average_norm_contrib_sorted(i)
-   if (norm > threshold_generators) then
+   if (norm >= threshold_generators) then
-     N_det_generators = i-1
+     N_det_generators = i
      exit
    endif
  enddo
@@ -52,4 +52,11 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ]
 
 END_PROVIDER
 
+ BEGIN_PROVIDER [ double precision, select_max, (3000) ]
+ implicit none
+ BEGIN_DOC
+ ! Memo to skip useless selectors
+ END_DOC
+ select_max = huge(1.d0)
+END_PROVIDER
 

src/Perturbation/epstein_nesbet.irp.f

@@ -66,14 +66,24 @@ subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet
   call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
   h = diag_H_mat_elem(det_pert,Nint)
   do i =1,N_st
-    delta_e = h - CI_electronic_energy(i)
+    if (i_H_psi_array(i) /= 0.d0) then
-    e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
+      delta_e = h - CI_electronic_energy(i)
-    if (dabs(i_H_psi_array(i)) > 1.d-6) then
+      if (delta_e > 0.d0) then
-      c_pert(i) = e_2_pert(i)/i_H_psi_array(i)
+        e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
+      else
+        e_2_pert(i) = 0.5d0 * (delta_e + dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
+      endif
+      if (dabs(i_H_psi_array(i)) > 1.d-6) then
+        c_pert(i) = e_2_pert(i)/i_H_psi_array(i)
+      else
+        c_pert(i) = 0.d0
+      endif
+      H_pert_diag(i) = h*c_pert(i)*c_pert(i)
     else
-      c_pert(i) = -1.d0
+      e_2_pert(i) = 0.d0
+      c_pert(i) = 0.d0
+      H_pert_diag(i) = 0.d0
     endif
-    H_pert_diag(i) = h*c_pert(i)*c_pert(i)
   enddo
 
 end

src/Perturbation/selection.irp.f

@@ -1,4 +1,5 @@
-subroutine fill_H_apply_buffer_selection(n_selected,det_buffer,e_2_pert_buffer,coef_pert_buffer,N_st,Nint,iproc)
+subroutine fill_H_apply_buffer_selection(n_selected,det_buffer,e_2_pert_buffer,coef_pert_buffer, &
+   N_st,Nint,iproc,select_max_out)
   use bitmasks
   implicit none
   BEGIN_DOC
@@ -9,6 +10,7 @@ subroutine fill_H_apply_buffer_selection(n_selected,det_buffer,e_2_pert_buffer,c
   integer(bit_kind), intent(in)  :: det_buffer(Nint,2,n_selected)
   double precision, intent(in)   :: e_2_pert_buffer(N_st,n_selected)
   double precision, intent(in)   :: coef_pert_buffer(N_st,n_selected)
+  double precision, intent(inout):: select_max_out
   integer                        :: i,j,k,l
   integer                        :: new_size
   double precision               :: s, smin, smax
@@ -35,6 +37,7 @@ subroutine fill_H_apply_buffer_selection(n_selected,det_buffer,e_2_pert_buffer,c
     do j=1,N_st
       s = -e_2_pert_buffer(j,i)
       is_selected = s > selection_criterion*selection_criterion_factor .or. is_selected
+      select_max_out = max(select_max_out,s)
     enddo
     
      
@@ -72,7 +75,7 @@ end
  BEGIN_DOC
  ! Threshold to select determinants. Set by selection routines.
  END_DOC
- selection_criterion =  10.d0 
+ selection_criterion =  .1d0 
  selection_criterion_factor = 0.01d0
  selection_criterion_min = selection_criterion
 
@@ -132,5 +135,3 @@ subroutine remove_small_contributions
 end
 
 
-
-

src/Selectors_full/e_corr_selectors.irp.f

@@ -51,7 +51,11 @@ END_PROVIDER
    E_corr_per_selectors(i) = -1000.d0
   endif
  enddo
- inv_selectors_coef_hf = 1.d0/coef_hf_selector
+ if (dabs(coef_hf_selector) > 1.d-8) then
+   inv_selectors_coef_hf = 1.d0/coef_hf_selector
+ else
+   inv_selectors_coef_hf = 0.d0
+ endif
  do i = 1,n_double_selectors
   E_corr_per_selectors(double_index_selectors(i)) *=inv_selectors_coef_hf
  enddo

src/SingleRefMethod/generators.irp.f

@@ -40,4 +40,10 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ]
 
 END_PROVIDER
 
-
+BEGIN_PROVIDER [ double precision, select_max, (1) ]
+ implicit none
+ BEGIN_DOC
+ ! Memo to skip useless selectors
+ END_DOC
+ select_max(1) = huge(1.d0)
+END_PROVIDER