From 3d7687e3c3c0946f110e9346973f937cc7b774db Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 24 Nov 2015 17:01:09 +0100 Subject: [PATCH] Accelerated diagonal element calculations in PT2 --- plugins/Perturbation/perturbation.template.f | 5 +- plugins/Perturbation/pt2_equations.irp.f | 30 ++++--- scripts/generate_h_apply.py | 2 +- src/Determinants/Fock_diag.irp.f | 84 ++++++++++++++++++++ src/Determinants/H_apply.template.f | 57 +++++++------ src/Determinants/slater_rules.irp.f | 71 ++++++++++++++++- 6 files changed, 209 insertions(+), 40 deletions(-) create mode 100644 src/Determinants/Fock_diag.irp.f diff --git a/plugins/Perturbation/perturbation.template.f b/plugins/Perturbation/perturbation.template.f index 13099bbd..fa01cc99 100644 --- a/plugins/Perturbation/perturbation.template.f +++ b/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 diff --git a/plugins/Perturbation/pt2_equations.irp.f b/plugins/Perturbation/pt2_equations.irp.f index c2c8026c..f49ee2ff 100644 --- a/plugins/Perturbation/pt2_equations.irp.f +++ b/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) diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index d23b18c8..e1c915bc 100755 --- a/scripts/generate_h_apply.py +++ b/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"] = """ """ diff --git a/src/Determinants/Fock_diag.irp.f b/src/Determinants/Fock_diag.irp.f new file mode 100644 index 00000000..a99bbcad --- /dev/null +++ b/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 = 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 diff --git a/src/Determinants/H_apply.template.f b/src/Determinants/H_apply.template.f index 7ee88e28..58ae8b08 100644 --- a/src/Determinants/H_apply.template.f +++ b/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 - - integer(bit_kind), intent(in) :: key_prev(N_int, 2, *) + double precision, intent(in) :: fock_diag_tmp(2,mo_tot_num+1) + 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, & - keys_out, hole_save, & - key,hole, particle, hole_tmp,& - particle_tmp, occ_particle, & - occ_hole, occ_particle_tmp,& + 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,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) diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index 46be7c18..32e84532 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/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 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