Accelerated diagonal element calculations in PT2

plugins/Perturbation/perturbation.template.f

@@ -2,7 +2,7 @@ BEGIN_SHELL [ /usr/bin/env python ]
 import perturbation
 END_SHELL
 
-subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint,key_mask)
+subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint,key_mask,fock_diag_tmp)
   implicit none
   BEGIN_DOC
   !  Applly pertubration ``$PERT`` to the buffer of determinants generated in the H_apply
@@ -12,6 +12,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
   integer, intent(in)            :: Nint, N_st, buffer_size, i_generator
   integer(bit_kind), intent(in)  :: buffer(Nint,2,buffer_size)
   integer(bit_kind),intent(in)    :: key_mask(Nint,2)
+  double precision, intent(in)    :: fock_diag_tmp(2,0:mo_tot_num)
   double precision, intent(inout) :: sum_norm_pert(N_st),sum_e_2_pert(N_st)
   double precision, intent(inout) :: coef_pert_buffer(N_st,buffer_size),e_2_pert_buffer(N_st,buffer_size),sum_H_pert_diag(N_st)
   double precision               :: c_pert(N_st), e_2_pert(N_st),  H_pert_diag(N_st)
@@ -55,7 +56,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
     
     integer :: degree
     call get_excitation_degree(HF_bitmask,buffer(1,1,i),degree,N_int)
-    call pt2_$PERT(buffer(1,1,i),         &
+    call pt2_$PERT(psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp,        &
          c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist) 
 
     do k = 1,N_st

plugins/Perturbation/pt2_equations.irp.f

@@ -17,7 +17,7 @@ subroutine pt2_epstein_nesbet ($arguments)
   END_DOC
   
   integer                        :: i,j
-  double precision               :: diag_H_mat_elem, h
+  double precision               :: diag_H_mat_elem_fock, h
   double precision               :: i_H_psi_array(N_st)
   PROVIDE  selection_criterion
 
@@ -27,7 +27,7 @@ subroutine pt2_epstein_nesbet ($arguments)
   call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array)
   
   
-  h = diag_H_mat_elem(det_pert,Nint)
+  h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
   do i =1,N_st
     if(CI_electronic_energy(i)>h.and.CI_electronic_energy(i).ne.0.d0)then
       c_pert(i) = -1.d0
@@ -62,7 +62,7 @@ subroutine pt2_epstein_nesbet_2x2 ($arguments)
   END_DOC
   
   integer                        :: i,j
-  double precision               :: diag_H_mat_elem,delta_e, h
+  double precision               :: diag_H_mat_elem_fock,delta_e, h
   double precision               :: i_H_psi_array(N_st)
   ASSERT (Nint == N_int)
   ASSERT (Nint > 0)
@@ -71,7 +71,7 @@ subroutine pt2_epstein_nesbet_2x2 ($arguments)
   !call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
   call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array)
   
-  h = diag_H_mat_elem(det_pert,Nint)
+  h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
   do i =1,N_st
     if (i_H_psi_array(i) /= 0.d0) then
       delta_e = h - CI_electronic_energy(i)
@@ -112,7 +112,7 @@ subroutine pt2_moller_plesset ($arguments)
   END_DOC
   
   integer                        :: i,j
-  double precision               :: diag_H_mat_elem
+  double precision               :: diag_H_mat_elem_fock
   integer                        :: exc(0:2,2,2)
   integer                        :: degree
   double precision               :: phase,delta_e,h
@@ -135,7 +135,7 @@ subroutine pt2_moller_plesset ($arguments)
   endif
 
   call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det,psi_selectors_size,n_st,i_H_psi_array)
-  h = diag_H_mat_elem(det_pert,Nint)
+  h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
   do i =1,n_st
     H_pert_diag(i) = h
     c_pert(i) = i_H_psi_array(i) *delta_e
@@ -176,7 +176,7 @@ subroutine pt2_epstein_nesbet_SC2_projected ($arguments)
   double precision               :: i_H_psi_array(N_st)
   integer                        :: idx_repeat(0:ndet)
   integer                        :: i,j,degree,l
-  double precision               :: diag_H_mat_elem,accu_e_corr,hij,h0j,h,delta_E
+  double precision               :: diag_H_mat_elem_fock,accu_e_corr,hij,h0j,h,delta_E
   double precision               :: repeat_all_e_corr,accu_e_corr_tmp,e_2_pert_fonda
 
   ASSERT (Nint == N_int)
@@ -188,7 +188,8 @@ subroutine pt2_epstein_nesbet_SC2_projected ($arguments)
   do i = 1, idx_repeat(0)
    accu_e_corr = accu_e_corr + E_corr_per_selectors(idx_repeat(i))
   enddo
-  h =  diag_H_mat_elem(det_pert,Nint) + accu_e_corr
+  h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
+  h = h + accu_e_corr
   delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h)
 
 
@@ -258,7 +259,7 @@ subroutine pt2_epstein_nesbet_SC2_no_projected ($arguments)
   double precision               :: i_H_psi_array(N_st)
   integer                        :: idx_repeat(0:ndet)
   integer                        :: i,j,degree,l
-  double precision               :: diag_H_mat_elem,accu_e_corr,hij,h0j,h,delta_E
+  double precision               :: diag_H_mat_elem_fock,accu_e_corr,hij,h0j,h,delta_E
   double precision               :: repeat_all_e_corr,accu_e_corr_tmp,e_2_pert_fonda
 
   ASSERT (Nint == N_int)
@@ -270,7 +271,8 @@ subroutine pt2_epstein_nesbet_SC2_no_projected ($arguments)
   do i = 1, idx_repeat(0)
    accu_e_corr = accu_e_corr + E_corr_per_selectors(idx_repeat(i))
   enddo
-  h =  diag_H_mat_elem(det_pert,Nint) + accu_e_corr
+  h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
+  h = h + accu_e_corr
   delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h)
 
 
@@ -310,7 +312,7 @@ subroutine pt2_epstein_nesbet_sc2 ($arguments)
   
   integer                        :: i,j
   double precision               :: i_H_psi_array(N_st)
-  double precision               :: diag_H_mat_elem, h
+  double precision               :: diag_H_mat_elem_fock, h
   PROVIDE  selection_criterion
 
   ASSERT (Nint == N_int)
@@ -319,7 +321,7 @@ subroutine pt2_epstein_nesbet_sc2 ($arguments)
   call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array)
 
   
-  h = diag_H_mat_elem(det_pert,Nint)
+  h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
   do i =1,N_st
     if(CI_SC2_electronic_energy(i)>h.and.CI_SC2_electronic_energy(i).ne.0.d0)then
       c_pert(i) = -1.d0
@@ -341,13 +343,15 @@ end
 
 SUBST [ arguments, declarations ]
 
-det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist ;
+det_ref,det_pert,fock_diag_tmp,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist ;
 
     integer, intent(in)             :: Nint
     integer, intent(in)             :: ndet
     integer, intent(in)             :: N_st
     integer, intent(in)             :: N_minilist
+    integer(bit_kind), intent(in)   :: det_ref (Nint,2)
     integer(bit_kind), intent(in)   :: det_pert(Nint,2)
+    double precision , intent(in)   :: fock_diag_tmp(2,mo_tot_num+1)
     double precision , intent(out)  :: c_pert(N_st)
     double precision , intent(out)  :: e_2_pert(N_st)
     double precision, intent(out)   :: H_pert_diag(N_st)

scripts/generate_h_apply.py

@@ -204,7 +204,7 @@ class H_apply(object):
       """
       self.data["keys_work"] = """
       call perturb_buffer_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, &
-       sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask)
+       sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask,fock_diag_tmp)
       """%(pert,)
       self.data["finalization"] = """
       """

src/Determinants/Fock_diag.irp.f

@@ -0,0 +1,84 @@
+subroutine build_fock_tmp(fock_diag_tmp,det_ref,Nint)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+! Build the diagonal of the Fock matrix corresponding to a generator
+! determinant. F_00 is <i|H|i> = E0.
+  END_DOC
+  integer, intent(in)            :: Nint
+  integer(bit_kind), intent(in)  :: det_ref(Nint,2)
+  double precision, intent(out)  :: fock_diag_tmp(2,mo_tot_num+1)
+  
+  integer                        :: occ(Nint*bit_kind_size,2)
+  integer                        :: ne(2), i, j, ii, jj
+  double precision               :: E0
+  
+  ! Compute Fock matrix diagonal elements
+  call bitstring_to_list_ab(det_ref,occ,Ne,Nint)
+
+  fock_diag_tmp = 0.d0
+  E0 = 0.d0
+
+  ! Occupied MOs
+  do ii=1,elec_alpha_num
+    i = occ(ii,1)
+    fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_mono_elec_integral(i,i)
+    E0 = E0 + mo_mono_elec_integral(i,i)
+    do jj=1,elec_alpha_num
+      j = occ(jj,1)
+      if (i==j) cycle
+      fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bielec_integral_jj_anti(i,j)
+      E0 = E0 + 0.5d0*mo_bielec_integral_jj_anti(i,j)
+    enddo
+    do jj=1,elec_beta_num
+      j = occ(jj,2)
+      fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bielec_integral_jj(i,j)
+      E0 = E0 + mo_bielec_integral_jj(i,j)
+    enddo
+  enddo
+  do ii=1,elec_beta_num
+    i = occ(ii,2)
+    fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_mono_elec_integral(i,i)
+    E0 = E0 + mo_mono_elec_integral(i,i)
+    do jj=1,elec_beta_num
+      j = occ(jj,2)
+      if (i==j) cycle
+      fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bielec_integral_jj_anti(i,j)
+      E0 = E0 + 0.5d0*mo_bielec_integral_jj_anti(i,j)
+    enddo
+    do jj=1,elec_alpha_num
+      j = occ(jj,1)
+      fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bielec_integral_jj(i,j)
+    enddo
+  enddo
+
+  ! Virtual MOs
+  do i=1,mo_tot_num
+    if (fock_diag_tmp(1,i) /= 0.d0) cycle
+    fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_mono_elec_integral(i,i)
+    do jj=1,elec_alpha_num
+      j = occ(jj,1)
+      fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bielec_integral_jj_anti(i,j)
+    enddo
+    do jj=1,elec_beta_num
+      j = occ(jj,2)
+      fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bielec_integral_jj(i,j)
+    enddo
+  enddo
+  do i=1,mo_tot_num
+    if (fock_diag_tmp(2,i) /= 0.d0) cycle
+    fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_mono_elec_integral(i,i)
+    do jj=1,elec_beta_num
+      j = occ(jj,2)
+      fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bielec_integral_jj_anti(i,j)
+    enddo
+    do jj=1,elec_alpha_num
+      j = occ(jj,1)
+      fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bielec_integral_jj(i,j)
+    enddo
+  enddo
+
+  fock_diag_tmp(1,mo_tot_num+1) = E0
+  fock_diag_tmp(2,mo_tot_num+1) = E0
+
+end

src/Determinants/H_apply.template.f

@@ -1,6 +1,6 @@
 
 
-subroutine $subroutine_diexc(key_in, key_prev, hole_1,particl_1, hole_2, particl_2, i_generator, iproc_in $parameters )
+subroutine $subroutine_diexc(key_in, key_prev, hole_1,particl_1, hole_2, particl_2, fock_diag_tmp, i_generator, iproc_in $parameters )
   
   integer(bit_kind), intent(in)         :: key_in(N_int, 2), hole_1(N_int, 2), hole_2(N_int, 2)
   integer(bit_kind), intent(in)         :: particl_1(N_int, 2), particl_2(N_int, 2)
@@ -8,8 +8,8 @@ subroutine $subroutine_diexc(key_in, key_prev, hole_1,particl_1, hole_2, particl
   integer,intent(in)                    :: i_generator,iproc_in
   integer(bit_kind)                     :: status(N_int*bit_kind_size, 2)
   integer                               :: highest, p1,p2,sp,ni,i,mi,nt,ns
-  
+  double precision, intent(in)          :: fock_diag_tmp(2,mo_tot_num+1)
-  integer(bit_kind), intent(in)          :: key_prev(N_int, 2, *)
+  integer(bit_kind), intent(in)         :: key_prev(N_int, 2, *)
   PROVIDE N_int
   PROVIDE N_det
  
@@ -72,7 +72,7 @@ subroutine $subroutine_diexc(key_in, key_prev, hole_1,particl_1, hole_2, particl
         if((status(p1, sp) == 1 .and. status(p2, sp) > 1) .or. &
             (status(p1, sp) == 2 .and. status(p2, sp) == 3) .or. &
             (status(p1, sp) == 3 .and. status(p2, sp) == 3 .and. p2 > p1)) then
-          call $subroutine_diexcP(key_in, sp, p1, particl_1, sp, p2, particl_2, i_generator, iproc_in $parameters )
+          call $subroutine_diexcP(key_in, sp, p1, particl_1, sp, p2, particl_2, fock_diag_tmp, i_generator, iproc_in $parameters )
         end if
       end do
     end do
@@ -89,16 +89,17 @@ subroutine $subroutine_diexc(key_in, key_prev, hole_1,particl_1, hole_2, particl
           (status(p1, 1) == 1 .and. status(p2, 2) >= 2) .or. &
           (status(p1, 1) == 2 .and. status(p2, 2) /= 2)) then
           
-          call $subroutine_diexcP(key_in, 1, p1, particl_1, 2, p2, particl_2, i_generator, iproc_in $parameters )
+          call $subroutine_diexcP(key_in, 1, p1, particl_1, 2, p2, particl_2, fock_diag_tmp, i_generator, iproc_in $parameters )
       end if
     end do
   end do
 end subroutine
 
 
-subroutine $subroutine_diexcP(key_in, fs1, fh1, particl_1, fs2, fh2, particl_2, i_generator, iproc_in $parameters )
+subroutine $subroutine_diexcP(key_in, fs1, fh1, particl_1, fs2, fh2, particl_2, fock_diag_tmp, i_generator, iproc_in $parameters )
   
   integer(bit_kind), intent(in)         :: key_in(N_int, 2), particl_1(N_int, 2), particl_2(N_int, 2)
+  double precision, intent(in)          :: fock_diag_tmp(2,mo_tot_num+1)
   integer(bit_kind)                     :: p1_mask(N_int, 2), p2_mask(N_int, 2), key_mask(N_int, 2)
   integer,intent(in)                    :: fh1,fh2,fs1,fs2,i_generator,iproc_in
   integer(bit_kind)                     :: miniList(N_int, 2, N_det)
@@ -115,11 +116,11 @@ subroutine $subroutine_diexcP(key_in, fs1, fh1, particl_1, fs2, fh2, particl_2,
   key_mask(ishft(fh1,-bit_kind_shift) + 1, fs1) -= ishft(1,iand(fh1-1,bit_kind_size-1))
   key_mask(ishft(fh2,-bit_kind_shift) + 1, fs2) -= ishft(1,iand(fh2-1,bit_kind_size-1))
 
-  call $subroutine_diexcOrg(key_in, key_mask, p1_mask, particl_1, p2_mask, particl_2, i_generator, iproc_in $parameters )
+  call $subroutine_diexcOrg(key_in, key_mask, p1_mask, particl_1, p2_mask, particl_2, fock_diag_tmp, i_generator, iproc_in $parameters )
 end subroutine
 
 
-subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl_2, i_generator, iproc_in $parameters )
+subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl_2, fock_diag_tmp, i_generator, iproc_in $parameters )
   use omp_lib
   use bitmasks
   implicit none
@@ -136,6 +137,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
   integer(bit_kind), intent(in)  :: hole_1(N_int,2), particl_1(N_int,2)
   integer(bit_kind), intent(in)  :: hole_2(N_int,2), particl_2(N_int,2)
   integer, intent(in)            :: iproc_in
+  double precision, intent(in)   :: fock_diag_tmp(2,mo_tot_num+1)
   integer(bit_kind), allocatable :: hole_save(:,:)
   integer(bit_kind), allocatable :: key(:,:),hole(:,:), particle(:,:)
   integer(bit_kind), allocatable :: hole_tmp(:,:), particle_tmp(:,:)
@@ -175,6 +177,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
       particle_tmp(N_int,2), occ_particle(N_int*bit_kind_size,2),    &
       occ_hole(N_int*bit_kind_size,2), occ_particle_tmp(N_int*bit_kind_size,2),&
       occ_hole_tmp(N_int*bit_kind_size,2),key_union_hole_part(N_int))
+
   $init_thread
   
   
@@ -362,17 +365,17 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
   enddo   ! ispin
   $keys_work
   $deinit_thread
-  deallocate (ia_ja_pairs, ib_jb_pairs, &
+  deallocate (ia_ja_pairs, ib_jb_pairs,                              &
-      keys_out, hole_save,          &
+      keys_out, hole_save,                                           &
-      key,hole, particle, hole_tmp,&
+      key,hole, particle, hole_tmp,                                  &
-      particle_tmp, occ_particle,    &
+      particle_tmp, occ_particle,                                    &
-      occ_hole, occ_particle_tmp,&
+      occ_hole, occ_particle_tmp,                                    &
       occ_hole_tmp,array_pairs,key_union_hole_part)
   $omp_end_parallel
   $finalization
 end
 
-subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc_in $parameters )
+subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generator,iproc_in $parameters )
   use omp_lib
   use bitmasks
   implicit none
@@ -387,6 +390,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc_in $pa
   integer(bit_kind),intent(in)   :: key_in(N_int,2)
   integer(bit_kind),intent(in)   :: hole_1(N_int,2), particl_1(N_int,2)
   integer, intent(in)            :: iproc_in
+  double precision, intent(in)   :: fock_diag_tmp(2,mo_tot_num+1)
   integer(bit_kind),allocatable  :: keys_out(:,:,:)
   integer(bit_kind),allocatable  :: hole_save(:,:)
   integer(bit_kind),allocatable  :: key(:,:),hole(:,:), particle(:,:)
@@ -526,6 +530,7 @@ subroutine $subroutine($params_main)
   integer(bit_kind), allocatable :: mask(:,:,:)
   integer                        :: ispin, k
   integer                        :: iproc
+  double precision, allocatable  :: fock_diag_tmp(:,:)
 
   $initialization
   PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map psi_det_generators psi_coef_generators
@@ -539,7 +544,7 @@ subroutine $subroutine($params_main)
   call wall_time(wall_0)
 
   iproc = 0
-  allocate( mask(N_int,2,6) )
+  allocate( mask(N_int,2,6), fock_diag_tmp(2,mo_tot_num+1) )
   do i_generator=1,nmax
 
     progress_bar(1) = i_generator
@@ -549,6 +554,9 @@ subroutine $subroutine($params_main)
     endif
     $skip
 
+    ! Compute diagonal of the Fock matrix
+    call build_fock_tmp(fock_diag_tmp,psi_det_generators(1,1,i_generator),N_int)
+
     ! Create bit masks for holes and particles
     do ispin=1,2
       do k=1,N_int
@@ -577,12 +585,12 @@ subroutine $subroutine($params_main)
          psi_det_generators(1,1,1),                                   &
          mask(1,1,d_hole1), mask(1,1,d_part1),                        &
          mask(1,1,d_hole2), mask(1,1,d_part2),                        &
-         i_generator, iproc $params_post)
+         fock_diag_tmp, i_generator, iproc $params_post)
     endif
     if($do_mono_excitations)then
      call $subroutine_monoexc(psi_det_generators(1,1,i_generator),    &
          mask(1,1,s_hole ), mask(1,1,s_part ),                        &
-         i_generator, iproc $params_post)
+         fock_diag_tmp, i_generator, iproc $params_post)
     endif
     call wall_time(wall_1)
     $printout_always
@@ -592,13 +600,13 @@ subroutine $subroutine($params_main)
     endif
   enddo
 
-  deallocate( mask )
+  deallocate( mask, fock_diag_tmp )
 
   !$OMP PARALLEL DEFAULT(SHARED) &
-  !$OMP PRIVATE(i_generator,wall_1,wall_0,ispin,k,mask,iproc) 
+  !$OMP PRIVATE(i_generator,wall_1,wall_0,ispin,k,mask,iproc,fock_diag_tmp) 
   call wall_time(wall_0)
   !$ iproc = omp_get_thread_num()
-  allocate( mask(N_int,2,6) )
+  allocate( mask(N_int,2,6), fock_diag_tmp(2,mo_tot_num+1) )
   !$OMP DO SCHEDULE(dynamic,1)
   do i_generator=nmax+1,N_det_generators
     if (iproc == 0) then
@@ -609,6 +617,9 @@ subroutine $subroutine($params_main)
     endif
     $skip
 
+    ! Compute diagonal of the Fock matrix
+    call build_fock_tmp(fock_diag_tmp,psi_det_generators(1,1,i_generator),N_int)
+
     ! Create bit masks for holes and particles
     do ispin=1,2
       do k=1,N_int
@@ -638,12 +649,12 @@ subroutine $subroutine($params_main)
         psi_det_generators(1,1,1),                                   &
         mask(1,1,d_hole1), mask(1,1,d_part1),                        &
         mask(1,1,d_hole2), mask(1,1,d_part2),                        &
-        i_generator, iproc $params_post)
+        fock_diag_tmp, i_generator, iproc $params_post)
     endif
     if($do_mono_excitations)then
       call $subroutine_monoexc(psi_det_generators(1,1,i_generator),  &
         mask(1,1,s_hole ), mask(1,1,s_part ),                        &
-        i_generator, iproc $params_post)
+        fock_diag_tmp, i_generator, iproc $params_post)
     endif
     !$ call omp_set_lock(lck)
     call wall_time(wall_1)
@@ -655,7 +666,7 @@ subroutine $subroutine($params_main)
     !$ call omp_unset_lock(lck)
   enddo
   !$OMP END DO 
-  deallocate( mask )
+  deallocate( mask, fock_diag_tmp )
   !$OMP END PARALLEL
   !$ call omp_destroy_lock(lck)
 

src/Determinants/slater_rules.irp.f

@@ -1264,6 +1264,75 @@ end
 
 
 
+double precision function diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Computes <i|H|i> when i is at most a double excitation from
+  ! a reference.
+  END_DOC
+  integer,intent(in)             :: Nint
+  integer(bit_kind),intent(in)   :: det_ref(Nint,2), det_pert(Nint,2)
+  double precision, intent(in)   :: fock_diag_tmp(2,mo_tot_num+1)
+
+  integer                        :: degree
+  double precision               :: phase, E0
+  integer                        :: exc(0:2,2,2)
+  integer                        :: h1, p1, h2, p2, s1, s2
+
+  call get_excitation_degree(det_ref,det_pert,degree,Nint)
+  E0 = fock_diag_tmp(1,mo_tot_num+1)
+  if (degree == 2) then
+    call get_double_excitation(det_ref,det_pert,exc,phase,Nint)
+    call decode_exc(exc,2,h1,p1,h2,p2,s1,s2)
+
+    if ( (s1 == 1).and.(s2 == 1) ) then      ! alpha/alpha
+      diag_H_mat_elem_fock = E0 &
+        - fock_diag_tmp(1,h1) &
+        + ( fock_diag_tmp(1,p1) - mo_bielec_integral_jj_anti(h1,p1) ) &
+        - ( fock_diag_tmp(1,h2) - mo_bielec_integral_jj_anti(h1,h2)   &
+            + mo_bielec_integral_jj_anti(p1,h2) )                     &
+        + ( fock_diag_tmp(1,p2) - mo_bielec_integral_jj_anti(h1,p2)   &
+            + mo_bielec_integral_jj_anti(p1,p2) - mo_bielec_integral_jj_anti(h2,p2) )
+
+    else if ( (s1 == 2).and.(s2 == 2) ) then ! beta/beta
+      diag_H_mat_elem_fock = E0 &
+        - fock_diag_tmp(2,h1) &
+        + ( fock_diag_tmp(2,p1) - mo_bielec_integral_jj_anti(h1,p1) ) &
+        - ( fock_diag_tmp(2,h2) - mo_bielec_integral_jj_anti(h1,h2)   &
+            + mo_bielec_integral_jj_anti(p1,h2) )                     &
+        + ( fock_diag_tmp(2,p2) - mo_bielec_integral_jj_anti(h1,p2)   &
+            + mo_bielec_integral_jj_anti(p1,p2) - mo_bielec_integral_jj_anti(h2,p2) )
+
+    else                                    ! alpha/beta
+      diag_H_mat_elem_fock = E0 &
+        - fock_diag_tmp(1,h1) &
+        + ( fock_diag_tmp(1,p1) - mo_bielec_integral_jj_anti(h1,p1) ) &
+        - ( fock_diag_tmp(2,h2) - mo_bielec_integral_jj(h1,h2)        &
+            + mo_bielec_integral_jj(p1,h2) )                          &
+        + ( fock_diag_tmp(2,p2) - mo_bielec_integral_jj(h1,p2)        &
+            + mo_bielec_integral_jj(p1,p2) - mo_bielec_integral_jj_anti(h2,p2) )
+
+    endif
+
+  else if (degree == 1) then
+    call get_mono_excitation(det_ref,det_pert,exc,phase,Nint)
+    call decode_exc(exc,1,h1,p1,h2,p2,s1,s2)
+    if (s1 == 1) then
+      diag_H_mat_elem_fock = E0 - fock_diag_tmp(1,h1) &
+        + ( fock_diag_tmp(1,p1) - mo_bielec_integral_jj_anti(h1,p1) ) 
+    else 
+      diag_H_mat_elem_fock = E0 - fock_diag_tmp(2,h1) &
+        + ( fock_diag_tmp(2,p1) - mo_bielec_integral_jj_anti(h1,p1) ) 
+    endif
+
+  else if (degree == 0) then
+    diag_H_mat_elem_fock = E0
+  else
+    STOP 'Bug in diag_H_mat_elem_fock'
+  endif
+end
+
 double precision function diag_H_mat_elem(det_in,Nint)
   implicit none
   BEGIN_DOC
@@ -1541,8 +1610,8 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
   !$OMP DO SCHEDULE(dynamic)
   do sh=1,shortcut(0)
     do i=shortcut(sh),shortcut(sh+1)-1
-      local_threshold = threshold_davidson - dabs(u_0(org_i))
       org_i = sort_idx(i)
+      local_threshold = threshold_davidson - dabs(u_0(org_i))
       do j=shortcut(sh),i-1
         org_j = sort_idx(j)
         if ( dabs(u_0(org_j)) > local_threshold ) then