From 6cf3dcca0b00c9ee619a267bb8bb91a84fa73db7 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 24 Nov 2015 11:40:49 +0100 Subject: [PATCH] Extended template for pt2 --- plugins/Perturbation/pert_sc2.irp.f | 214 ------------------- plugins/Perturbation/perturbation.template.f | 34 +-- plugins/Perturbation/pt2_equations.irp.f | 199 +++++++++++++++++ scripts/generate_h_apply.py | 9 +- src/Determinants/create_excitations.irp.f | 2 +- 5 files changed, 212 insertions(+), 246 deletions(-) diff --git a/plugins/Perturbation/pert_sc2.irp.f b/plugins/Perturbation/pert_sc2.irp.f index 1ae6ce73..f225e757 100644 --- a/plugins/Perturbation/pert_sc2.irp.f +++ b/plugins/Perturbation/pert_sc2.irp.f @@ -1,166 +1,3 @@ - -subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist) - use bitmasks - implicit none - integer, intent(in) :: Nint,ndet,N_st - integer(bit_kind), intent(in) :: det_pert(Nint,2) - double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st) - double precision :: i_H_psi_array(N_st) - integer :: idx_repeat(0:ndet) - - integer, intent(in) :: N_minilist - integer, intent(in) :: idx_minilist(0:N_det_selectors) - integer(bit_kind), intent(in) :: minilist(Nint,2,N_det_selectors) - - BEGIN_DOC - ! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution - ! - ! for the various N_st states, - ! - ! but with the correction in the denominator - ! - ! comming from the interaction of that determinant with all the others determinants - ! - ! that can be repeated by repeating all the double excitations - ! - ! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1) - ! - ! that could be repeated to this determinant. - ! - ! In addition, for the perturbative energetic contribution you have the standard second order - ! - ! e_2_pert = ^2/(Delta_E) - ! - ! and also the purely projected contribution - ! - ! H_pert_diag = c_pert - END_DOC - - integer :: i,j,degree,l - double precision :: diag_H_mat_elem,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) - ASSERT (Nint > 0) - - call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat) - accu_e_corr = 0.d0 - !$IVDEP - 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 - delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h) - - - c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h) - e_2_pert(1) = i_H_psi_array(1) * c_pert(1) - - do i =2,N_st - H_pert_diag(i) = h - if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then - c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h)) - e_2_pert(i) = (c_pert(i) * i_H_psi_array(i)) - else - c_pert(i) = i_H_psi_array(i) - e_2_pert(i) = -dabs(i_H_psi_array(i)) - endif - enddo - - degree = popcnt(xor( ref_bitmask(1,1), det_pert(1,1))) + & - popcnt(xor( ref_bitmask(1,2), det_pert(1,2))) - !DEC$ NOUNROLL - do l=2,Nint - degree = degree+ popcnt(xor( ref_bitmask(l,1), det_pert(l,1))) + & - popcnt(xor( ref_bitmask(l,2), det_pert(l,2))) - enddo - if(degree==4)then - ! - e_2_pert_fonda = e_2_pert(1) - H_pert_diag(1) = e_2_pert(1) * c_pert(1) * c_pert(1) - do i = 1, N_st - do j = 1, idx_repeat(0) - e_2_pert(i) += e_2_pert_fonda * psi_selectors_coef(idx_repeat(j),i) * psi_selectors_coef(idx_repeat(j),i) - enddo - enddo - endif - -end - - -subroutine pt2_epstein_nesbet_SC2_no_projected(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist) - use bitmasks - implicit none - integer, intent(in) :: Nint,ndet,N_st - integer(bit_kind), intent(in) :: det_pert(Nint,2) - double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st) - double precision :: i_H_psi_array(N_st) - integer :: idx_repeat(0:ndet) - - integer, intent(in) :: N_minilist - integer, intent(in) :: idx_minilist(0:N_det_selectors) - integer(bit_kind), intent(in) :: minilist(Nint,2,N_det_selectors) - - BEGIN_DOC - ! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution - ! - ! for the various N_st states, - ! - ! but with the correction in the denominator - ! - ! comming from the interaction of that determinant with all the others determinants - ! - ! that can be repeated by repeating all the double excitations - ! - ! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1) - ! - ! that could be repeated to this determinant. - ! - ! In addition, for the perturbative energetic contribution you have the standard second order - ! - ! e_2_pert = ^2/(Delta_E) - ! - ! and also the purely projected contribution - ! - ! H_pert_diag = c_pert - END_DOC - - integer :: i,j,degree,l - double precision :: diag_H_mat_elem,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) - ASSERT (Nint > 0) - - call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat) - accu_e_corr = 0.d0 - !$IVDEP - 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 - delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h) - - - c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h) - e_2_pert(1) = i_H_psi_array(1) * c_pert(1) - - do i =2,N_st - H_pert_diag(i) = h - if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then - c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h)) - e_2_pert(i) = (c_pert(i) * i_H_psi_array(i)) - else - c_pert(i) = i_H_psi_array(i) - e_2_pert(i) = -dabs(i_H_psi_array(i)) - endif - enddo -end - - - - - double precision function repeat_all_e_corr(key_in) implicit none integer(bit_kind), intent(in) :: key_in(N_int,2) @@ -190,54 +27,3 @@ double precision function repeat_all_e_corr(key_in) end - -subroutine pt2_epstein_nesbet_sc2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist) - use bitmasks - implicit none - integer, intent(in) :: Nint,ndet,N_st - integer(bit_kind), intent(in) :: det_pert(Nint,2) - double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st) - double precision :: i_H_psi_array(N_st) - - integer, intent(in) :: N_minilist - integer, intent(in) :: idx_minilist(0:N_det_selectors) - integer(bit_kind), intent(in) :: minilist(Nint,2,N_det_selectors) - - BEGIN_DOC - ! compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution - ! - ! for the various N_st states, but with the CISD_SC2 energies and coefficients - ! - ! c_pert(i) = /( E(i) - ) - ! - ! e_2_pert(i) = ^2/( E(i) - ) - ! - END_DOC - - integer :: i,j - double precision :: diag_H_mat_elem, h - PROVIDE selection_criterion - - ASSERT (Nint == N_int) - ASSERT (Nint > 0) - !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) - 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 - e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0 - else if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then - c_pert(i) = i_H_psi_array(i) / (CI_SC2_electronic_energy(i) - h) - H_pert_diag(i) = h*c_pert(i)*c_pert(i) - e_2_pert(i) = c_pert(i) * i_H_psi_array(i) - else - c_pert(i) = -1.d0 - e_2_pert(i) = -dabs(i_H_psi_array(i)) - H_pert_diag(i) = h - endif - enddo - -end diff --git a/plugins/Perturbation/perturbation.template.f b/plugins/Perturbation/perturbation.template.f index 05176fe6..13099bbd 100644 --- a/plugins/Perturbation/perturbation.template.f +++ b/plugins/Perturbation/perturbation.template.f @@ -43,24 +43,8 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c end if - buffer_loop : do i = 1,buffer_size + do i=1,buffer_size -! do k=1,N_minilist_gen -! ex = 0 -! do ni=1,Nint -! ex += popcnt(xor(minilist_gen(ni,1,k), buffer(ni,1,i))) + popcnt(xor(minilist_gen(ni,2,k), buffer(ni,2,i))) -! end do -! if(ex <= 4) then -! cycle buffer_loop -! end if -! end do - -! c_ref = connected_to_ref(buffer(1,1,i),miniList_gen,Nint,N_minilist_gen+1,N_minilist_gen) -! -! if (c_ref /= 0) then -! cycle -! endif - if(is_connected_to(buffer(1,1,i), miniList_gen, Nint, N_minilist_gen)) then cycle end if @@ -71,20 +55,18 @@ 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), & -! c_pert,e_2_pert,H_pert_diag,Nint,N_det_selectors,n_st,minilist,idx_minilist) call pt2_$PERT(buffer(1,1,i), & - c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist) !!!!!!!!!!!!!!!!! MAUVAISE SIGNATURE PR LES AUTRES PT2_* !!!!! + c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist) do k = 1,N_st - e_2_pert_buffer(k,i) = e_2_pert(k) - coef_pert_buffer(k,i) = c_pert(k) - sum_norm_pert(k) += c_pert(k) * c_pert(k) - sum_e_2_pert(k) += e_2_pert(k) - sum_H_pert_diag(k) += H_pert_diag(k) + e_2_pert_buffer(k,i) = e_2_pert(k) + coef_pert_buffer(k,i) = c_pert(k) + sum_norm_pert(k) = sum_norm_pert(k) + c_pert(k) * c_pert(k) + sum_e_2_pert(k) = sum_e_2_pert(k) + e_2_pert(k) + sum_H_pert_diag(k) = sum_H_pert_diag(k) + H_pert_diag(k) enddo - enddo buffer_loop + enddo end diff --git a/plugins/Perturbation/pt2_equations.irp.f b/plugins/Perturbation/pt2_equations.irp.f index 79bba2be..c2c8026c 100644 --- a/plugins/Perturbation/pt2_equations.irp.f +++ b/plugins/Perturbation/pt2_equations.irp.f @@ -143,6 +143,202 @@ subroutine pt2_moller_plesset ($arguments) enddo end + + +subroutine pt2_epstein_nesbet_SC2_projected ($arguments) + use bitmasks + implicit none + $declarations + BEGIN_DOC + ! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution + ! + ! for the various N_st states, + ! + ! but with the correction in the denominator + ! + ! comming from the interaction of that determinant with all the others determinants + ! + ! that can be repeated by repeating all the double excitations + ! + ! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1) + ! + ! that could be repeated to this determinant. + ! + ! In addition, for the perturbative energetic contribution you have the standard second order + ! + ! e_2_pert = ^2/(Delta_E) + ! + ! and also the purely projected contribution + ! + ! H_pert_diag = c_pert + END_DOC + + 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 :: repeat_all_e_corr,accu_e_corr_tmp,e_2_pert_fonda + + ASSERT (Nint == N_int) + ASSERT (Nint > 0) + + call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat) + accu_e_corr = 0.d0 + !$IVDEP + 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 + delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h) + + + c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h) + e_2_pert(1) = i_H_psi_array(1) * c_pert(1) + + do i =2,N_st + H_pert_diag(i) = h + if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then + c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h)) + e_2_pert(i) = (c_pert(i) * i_H_psi_array(i)) + else + c_pert(i) = i_H_psi_array(i) + e_2_pert(i) = -dabs(i_H_psi_array(i)) + endif + enddo + + degree = popcnt(xor( ref_bitmask(1,1), det_pert(1,1))) + & + popcnt(xor( ref_bitmask(1,2), det_pert(1,2))) + !DEC$ NOUNROLL + do l=2,Nint + degree = degree+ popcnt(xor( ref_bitmask(l,1), det_pert(l,1))) + & + popcnt(xor( ref_bitmask(l,2), det_pert(l,2))) + enddo + if(degree==4)then + ! + e_2_pert_fonda = e_2_pert(1) + H_pert_diag(1) = e_2_pert(1) * c_pert(1) * c_pert(1) + do i = 1, N_st + do j = 1, idx_repeat(0) + e_2_pert(i) += e_2_pert_fonda * psi_selectors_coef(idx_repeat(j),i) * psi_selectors_coef(idx_repeat(j),i) + enddo + enddo + endif + +end + + +subroutine pt2_epstein_nesbet_SC2_no_projected ($arguments) + use bitmasks + implicit none + $declarations + BEGIN_DOC + ! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution + ! + ! for the various N_st states, + ! + ! but with the correction in the denominator + ! + ! comming from the interaction of that determinant with all the others determinants + ! + ! that can be repeated by repeating all the double excitations + ! + ! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1) + ! + ! that could be repeated to this determinant. + ! + ! In addition, for the perturbative energetic contribution you have the standard second order + ! + ! e_2_pert = ^2/(Delta_E) + ! + ! and also the purely projected contribution + ! + ! H_pert_diag = c_pert + END_DOC + + 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 :: repeat_all_e_corr,accu_e_corr_tmp,e_2_pert_fonda + + ASSERT (Nint == N_int) + ASSERT (Nint > 0) + + call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat) + accu_e_corr = 0.d0 + !$IVDEP + 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 + delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h) + + + c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h) + e_2_pert(1) = i_H_psi_array(1) * c_pert(1) + + do i =2,N_st + H_pert_diag(i) = h + if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then + c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h)) + e_2_pert(i) = (c_pert(i) * i_H_psi_array(i)) + else + c_pert(i) = i_H_psi_array(i) + e_2_pert(i) = -dabs(i_H_psi_array(i)) + endif + enddo +end + + + + + +subroutine pt2_epstein_nesbet_sc2 ($arguments) + use bitmasks + implicit none + $declarations + BEGIN_DOC + ! compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution + ! + ! for the various N_st states, but with the CISD_SC2 energies and coefficients + ! + ! c_pert(i) = /( E(i) - ) + ! + ! e_2_pert(i) = ^2/( E(i) - ) + ! + END_DOC + + integer :: i,j + double precision :: i_H_psi_array(N_st) + double precision :: diag_H_mat_elem, h + PROVIDE selection_criterion + + ASSERT (Nint == N_int) + ASSERT (Nint > 0) + !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) + 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 + e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0 + else if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then + c_pert(i) = i_H_psi_array(i) / (CI_SC2_electronic_energy(i) - h) + H_pert_diag(i) = h*c_pert(i)*c_pert(i) + e_2_pert(i) = c_pert(i) * i_H_psi_array(i) + else + c_pert(i) = -1.d0 + e_2_pert(i) = -dabs(i_H_psi_array(i)) + H_pert_diag(i) = h + endif + enddo + +end + + + SUBST [ arguments, declarations ] det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist ; @@ -162,3 +358,6 @@ det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_mini END_TEMPLATE +! Note : If the arguments are changed here, they should also be changed accordingly in +! the perturbation.template.f file. + diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index 26f19d0d..d23b18c8 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -131,10 +131,10 @@ class H_apply(object): def filter_vvvv_excitation(self): self["filter_vvvv_excitation"] = """ key_union_hole_part = 0_bit_kind - call set_bite_to_integer(i_a,key_union_hole_part,N_int) - call set_bite_to_integer(j_a,key_union_hole_part,N_int) - call set_bite_to_integer(i_b,key_union_hole_part,N_int) - call set_bite_to_integer(j_b,key_union_hole_part,N_int) + call set_bit_to_integer(i_a,key_union_hole_part,N_int) + call set_bit_to_integer(j_a,key_union_hole_part,N_int) + call set_bit_to_integer(i_b,key_union_hole_part,N_int) + call set_bit_to_integer(j_b,key_union_hole_part,N_int) do jtest_vvvv = 1, N_int if(iand(key_union_hole_part(jtest_vvvv),virt_bitmask(jtest_vvvv,1).ne.key_union_hole_part(jtest_vvvv)))then b_cycle = .False. @@ -157,7 +157,6 @@ class H_apply(object): def set_filter_2h_2p(self): self["filter2h2p"] = """ -! ! DIR$ FORCEINLINE if (is_a_two_holes_two_particles(key)) cycle """ diff --git a/src/Determinants/create_excitations.irp.f b/src/Determinants/create_excitations.irp.f index a2acc8df..b7233beb 100644 --- a/src/Determinants/create_excitations.irp.f +++ b/src/Determinants/create_excitations.irp.f @@ -35,7 +35,7 @@ subroutine do_mono_excitation(key_in,i_hole,i_particle,ispin,i_ok) endif end -subroutine set_bite_to_integer(i_physical,key,Nint) +subroutine set_bit_to_integer(i_physical,key,Nint) use bitmasks implicit none integer, intent(in) :: i_physical,Nint