From bd767bbb36991be2229f24a94a387e5e64011b0e Mon Sep 17 00:00:00 2001 From: Yann Garniron Date: Mon, 26 Feb 2018 11:33:32 +0100 Subject: [PATCH] i_h_j_s2 - optimization in mrcc_sto --- plugins/dress_zmq/dressing_vector.irp.f | 8 -- plugins/mrcc_sto/mrcc_sto.irp.f | 57 +++++++------ plugins/mrcepa0/dressing.irp.f | 19 +---- plugins/mrcepa0/dressing_vector.irp.f | 8 -- src/Determinants/slater_rules.irp.f | 101 ++++++++++++++++++++++++ 5 files changed, 136 insertions(+), 57 deletions(-) diff --git a/plugins/dress_zmq/dressing_vector.irp.f b/plugins/dress_zmq/dressing_vector.irp.f index 0d24ed94..d3e465bd 100644 --- a/plugins/dress_zmq/dressing_vector.irp.f +++ b/plugins/dress_zmq/dressing_vector.irp.f @@ -1,8 +1,3 @@ - BEGIN_PROVIDER [ integer, nalp ] -&BEGIN_PROVIDER [ integer, ninc ] - nalp = 0 - ninc = 0 - END_PROVIDER BEGIN_PROVIDER [ double precision, dressing_column_h, (N_det,N_states) ] &BEGIN_PROVIDER [ double precision, dressing_column_s, (N_det,N_states) ] @@ -30,8 +25,5 @@ tmp = u_dot_v(dressing_column_s(1,k), psi_coef(1,k), N_det) dressing_column_s(l,k) -= tmp * f enddo - print *, "NALP", nalp - print *, "NINC", ninc - print *, "DELTA_IJ", dressing_column_h(:10,1) END_PROVIDER diff --git a/plugins/mrcc_sto/mrcc_sto.irp.f b/plugins/mrcc_sto/mrcc_sto.irp.f index 09b5d626..0eb6a75e 100644 --- a/plugins/mrcc_sto/mrcc_sto.irp.f +++ b/plugins/mrcc_sto/mrcc_sto.irp.f @@ -12,7 +12,9 @@ end &BEGIN_PROVIDER [ double precision, sij_cache_, (N_det,Nproc) ] &BEGIN_PROVIDER [ double precision, dIa_hla_, (N_states,N_det,Nproc) ] &BEGIN_PROVIDER [ double precision, dIa_sla_, (N_states,N_det,Nproc) ] - BEGIN_DOC +&BEGIN_PROVIDER [ integer, excs_ , (0:2,2,2,N_det,Nproc) ] +&BEGIN_PROVIDER [ double precision, phases_, (N_det, Nproc) ] +BEGIN_DOC ! temporay arrays for dress_with_alpha_buffer. Avoids realocation. END_DOC END_PROVIDER @@ -46,8 +48,7 @@ subroutine dress_with_alpha_buffer(delta_ij_loc, minilist, det_minilist, n_minil double precision :: Delta_E_inv(N_states) double precision :: sdress, hdress logical :: ok, ok2 - integer :: old_ninc - double precision :: shdress + integer :: canbediamond PROVIDE mo_class @@ -57,28 +58,33 @@ subroutine dress_with_alpha_buffer(delta_ij_loc, minilist, det_minilist, n_minil if(idx_non_ref_rev(minilist(i)) == 0) return end do - shdress = 0d0 - old_ninc = ninc - if (perturbative_triples) then PROVIDE one_anhil fock_virt_total fock_core_inactive_total one_creat endif - + canbediamond = 0 do l_sd=1,n_minilist - !call get_excitation(det_minilist(1,1,l_sd),alpha,exc,degree1,phase,N_int) - !if(degree1 == 0 .or. degree1 > 2) stop "minilist error" - !call decode_exc(exc,degree1,h1,p1,h2,p2,s1,s2) - ! - !ok = (mo_class(h1)(1:1) == 'A' .or. mo_class(h1)(1:1) == 'I') .and. & - ! (mo_class(p1)(1:1) == 'A' .or. mo_class(p1)(1:1) == 'V') - !if(ok .and. degree1 == 2) then - ! ok = (mo_class(h2)(1:1) == 'A' .or. mo_class(h2)(1:1) == 'I') .and. & - ! (mo_class(p2)(1:1) == 'A' .or. mo_class(p2)(1:1) == 'V') - !end if - call i_h_j(alpha,det_minilist(1,1,l_sd),N_int,hij_cache_(l_sd,iproc)) - call get_s2(alpha,det_minilist(1,1,l_sd),N_int,sij_cache_(l_sd,iproc)) + call get_excitation(det_minilist(1,1,l_sd),alpha,exc,degree1,phase,N_int) + call decode_exc(exc,degree1,h1,p1,h2,p2,s1,s2) + + ok = (mo_class(h1)(1:1) == 'A' .or. mo_class(h1)(1:1) == 'I') .and. & + (mo_class(p1)(1:1) == 'A' .or. mo_class(p1)(1:1) == 'V') + if(ok .and. degree1 == 2) then + ok = (mo_class(h2)(1:1) == 'A' .or. mo_class(h2)(1:1) == 'I') .and. & + (mo_class(p2)(1:1) == 'A' .or. mo_class(p2)(1:1) == 'V') + end if + + if(ok) then + canbediamond += 1 + excs_(:,:,:,l_sd,iproc) = exc(:,:,:) + phases_(l_sd, iproc) = phase + else + phases_(l_sd, iproc) = 0d0 + end if + !call i_h_j(alpha,det_minilist(1,1,l_sd),N_int,hij_cache_(l_sd,iproc)) + !call get_s2(alpha,det_minilist(1,1,l_sd),N_int,sij_cache_(l_sd,iproc)) + call i_h_j_s2(alpha,det_minilist(1,1,l_sd),N_int,hij_cache_(l_sd,iproc), sij_cache_(l_sd,iproc)) enddo - + if(canbediamond <= 1) return do i_I=1,N_det_ref call get_excitation_degree(alpha,psi_ref(1,1,i_I),degree1,N_int) @@ -91,12 +97,16 @@ subroutine dress_with_alpha_buffer(delta_ij_loc, minilist, det_minilist, n_minil enddo do k_sd=1,n_minilist + if(phases_(k_sd,iproc) == 0d0) cycle call get_excitation_degree(psi_ref(1,1,i_I),det_minilist(1,1,k_sd),degree,N_int) if (degree > 2) then cycle endif - call get_excitation(det_minilist(1,1,k_sd),alpha,exc,degree2,phase,N_int) + !call get_excitation(det_minilist(1,1,k_sd),alpha,exc,degree2,phase,N_int) + phase = phases_(k_sd, iproc) + exc(:,:,:) = excs_(:,:,:,k_sd,iproc) + degree2 = exc(0,1,1) + exc(0,1,2) call apply_excitation(psi_ref(1,1,i_I), exc, tmp_det, ok, N_int) if((.not. ok) .and. (.not. perturbative_triples)) cycle @@ -118,6 +128,7 @@ subroutine dress_with_alpha_buffer(delta_ij_loc, minilist, det_minilist, n_minil if (ok) then phase2 = 0d0 do l_sd=k_sd+1,n_minilist + if(phases_(l_sd, iproc) == 0d0) cycle call get_excitation_degree(tmp_det,det_minilist(1,1,l_sd),degree,N_int) if (degree == 0) then do i_state=1,N_states @@ -204,11 +215,11 @@ subroutine test_minilist(minilist, n_minilist, alpha) refc = 0 testc = 0 do i=1,N_det - call get_excitation_degree(psi_det_sorted(1,1,i), alpha, deg, N_int) + call get_excitation_degree(psi_det(1,1,i), alpha, deg, N_int) if(deg <= 2) refc(i) = refc(i) + 1 end do do i=1,n_minilist - call get_excitation_degree(psi_det_sorted(1,1,minilist(i)), alpha, deg, N_int) + call get_excitation_degree(psi_det(1,1,minilist(i)), alpha, deg, N_int) if(deg <= 2) then testc(minilist(i)) += 1 else diff --git a/plugins/mrcepa0/dressing.irp.f b/plugins/mrcepa0/dressing.irp.f index 1fc024b5..a376585c 100644 --- a/plugins/mrcepa0/dressing.irp.f +++ b/plugins/mrcepa0/dressing.irp.f @@ -343,10 +343,8 @@ subroutine mrcc_part_dress(delta_ij_, delta_ij_s2_, i_generator,n_selected,det_b !double precision, external :: get_dij, get_dij_index double precision :: Delta_E_inv(N_states) double precision, intent(inout) :: contrib(N_states) - double precision :: sdress, hdress, shdress - integer :: old_ninc + double precision :: sdress, hdress - old_ninc = ninc if (perturbative_triples) then PROVIDE one_anhil fock_virt_total fock_core_inactive_total one_creat @@ -411,9 +409,6 @@ subroutine mrcc_part_dress(delta_ij_, delta_ij_s2_, i_generator,n_selected,det_b do i_alpha=1,N_tq - old_ninc = ninc - shdress = 0d0 - if(key_mask(1,1) /= 0) then call getMobiles(tq(1,1,i_alpha), key_mask, mobiles, Nint) @@ -552,8 +547,6 @@ subroutine mrcc_part_dress(delta_ij_, delta_ij_s2_, i_generator,n_selected,det_b hdress = dIa_hla(i_state,k_sd) * psi_ref_coef(i_I,i_state) sdress = dIa_sla(i_state,k_sd) * psi_ref_coef(i_I,i_state) !!$OMP ATOMIC - if(hdress /= 0d0) ninc = ninc + 1 - shdress += hdress !$OMP ATOMIC contrib(i_state) += hdress * psi_coef(dressed_column_idx(i_state), i_state) * psi_non_ref_coef(k_sd, i_state) !$OMP ATOMIC @@ -563,16 +556,6 @@ subroutine mrcc_part_dress(delta_ij_, delta_ij_s2_, i_generator,n_selected,det_b enddo enddo enddo - if(ninc /= old_ninc) then - nalp = nalp + 1 - !print "(A8,I20,I20,E15.5)", "grepme", tq(1,1,i_alpha), tq(1,2,i_alpha), shdress - !if(tq(1,1,i_alpha) == 1007 .and. tq(1,2,i_alpha) == 17301943) then - ! print *, "foinder", i_generator - ! call debug_det(psi_det_generators(1,1, i_generator), N_int) - ! call debug_det(tq(1,1,i_alpha), N_int) - ! stop - ! end if - end if enddo deallocate (dIa_hla,dIa_sla,hij_cache,sij_cache) diff --git a/plugins/mrcepa0/dressing_vector.irp.f b/plugins/mrcepa0/dressing_vector.irp.f index 22394a92..933e57b9 100644 --- a/plugins/mrcepa0/dressing_vector.irp.f +++ b/plugins/mrcepa0/dressing_vector.irp.f @@ -1,9 +1,4 @@ - BEGIN_PROVIDER [ integer, nalp ] -&BEGIN_PROVIDER [ integer, ninc ] - nalp = 0 - ninc = 0 -END_PROVIDER BEGIN_PROVIDER [ double precision, dressing_column_h, (N_det,N_states) ] &BEGIN_PROVIDER [ double precision, dressing_column_s, (N_det,N_states) ] @@ -31,9 +26,6 @@ END_PROVIDER tmp = u_dot_v(dressing_column_s(1,k), psi_coef(1,k), N_det) dressing_column_s(l,k) -= tmp * f enddo - print *, "NALP", nalp - print *, "NINC", ninc - print *, "DRESS", dressing_column_h(:10,1) ! stop END_PROVIDER diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index 75baf269..ee597720 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -474,6 +474,107 @@ subroutine bitstring_to_list_ab_old( string, list, n_elements, Nint) end +subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2) + use bitmasks + implicit none + BEGIN_DOC + ! Returns where i and j are determinants + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + double precision, intent(out) :: hij, s2 + + integer :: exc(0:2,2,2) + integer :: degree + double precision :: get_mo_bielec_integral + integer :: m,n,p,q + integer :: i,j,k + integer :: occ(Nint*bit_kind_size,2) + double precision :: diag_H_mat_elem, phase,phase_2 + integer :: n_occ_ab(2) + PROVIDE mo_bielec_integrals_in_map mo_integrals_map big_array_exchange_integrals + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num) + ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num) + + hij = 0.d0 + s2 = 0d0 + !DIR$ FORCEINLINE + call get_excitation_degree(key_i,key_j,degree,Nint) + integer :: spin + select case (degree) + case (2) + call get_double_excitation(key_i,key_j,exc,phase,Nint) + + if (exc(0,1,1) == 1) then + ! Mono alpha, mono beta + if(exc(1,1,1) == exc(1,2,2) )then + if(exc(1,1,2) == exc(1,2,1)) s2 = -phase !!!!! + hij = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1)) + else if (exc(1,2,1) ==exc(1,1,2))then + hij = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2)) + else + hij = phase*get_mo_bielec_integral( & + exc(1,1,1), & + exc(1,1,2), & + exc(1,2,1), & + exc(1,2,2) ,mo_integrals_map) + endif + else if (exc(0,1,1) == 2) then + ! Double alpha + hij = phase*(get_mo_bielec_integral( & + exc(1,1,1), & + exc(2,1,1), & + exc(1,2,1), & + exc(2,2,1) ,mo_integrals_map) - & + get_mo_bielec_integral( & + exc(1,1,1), & + exc(2,1,1), & + exc(2,2,1), & + exc(1,2,1) ,mo_integrals_map) ) + else if (exc(0,1,2) == 2) then + ! Double beta + hij = phase*(get_mo_bielec_integral( & + exc(1,1,2), & + exc(2,1,2), & + exc(1,2,2), & + exc(2,2,2) ,mo_integrals_map) - & + get_mo_bielec_integral( & + exc(1,1,2), & + exc(2,1,2), & + exc(2,2,2), & + exc(1,2,2) ,mo_integrals_map) ) + endif + case (1) + call get_mono_excitation(key_i,key_j,exc,phase,Nint) + !DIR$ FORCEINLINE + call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint) + if (exc(0,1,1) == 1) then + ! Mono alpha + m = exc(1,1,1) + p = exc(1,2,1) + spin = 1 + else + ! Mono beta + m = exc(1,1,2) + p = exc(1,2,2) + spin = 2 + endif + call get_mono_excitation_from_fock(key_i,key_j,p,m,spin,phase,hij) + + case (0) + print *," ZERO" + double precision, external :: diag_S_mat_elem + s2 = diag_S_mat_elem(key_i,Nint) + hij = diag_H_mat_elem(key_i,Nint) + end select +end + + subroutine i_H_j(key_i,key_j,Nint,hij) use bitmasks