From 3946c710feb8454d53f74ba8dc16390c2de01429 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 28 Oct 2016 22:26:20 +0200 Subject: [PATCH] Accelerated mono-excitations (mipi miip) --- src/Determinants/slater_rules.irp.f | 131 +++------------------ src/Integrals_Bielec/map_integrals.irp.f | 1 + src/Integrals_Bielec/mo_bi_integrals.irp.f | 25 ++++ 3 files changed, 40 insertions(+), 117 deletions(-) diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index 67463088..7df6e79e 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -515,8 +515,6 @@ subroutine i_H_j(key_i,key_j,Nint,hij) integer :: occ(Nint*bit_kind_size,2) double precision :: diag_H_mat_elem, phase,phase_2 integer :: n_occ_ab(2) - logical :: has_mipi(Nint*bit_kind_size) - double precision :: mipi(Nint*bit_kind_size), miip(Nint*bit_kind_size) PROVIDE mo_bielec_integrals_in_map mo_integrals_map ASSERT (Nint > 0) @@ -568,59 +566,27 @@ subroutine i_H_j(key_i,key_j,Nint,hij) 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) - has_mipi = .False. if (exc(0,1,1) == 1) then ! Mono alpha m = exc(1,1,1) p = exc(1,2,1) do k = 1, elec_alpha_num - i = occ(k,1) - if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) - has_mipi(i) = .True. - endif + hij = hij + mo_bielec_integral_mipi_anti(occ(k,1),m,p) enddo do k = 1, elec_beta_num - i = occ(k,2) - if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - has_mipi(i) = .True. - endif - enddo - - do k = 1, elec_alpha_num - hij = hij + mipi(occ(k,1)) - miip(occ(k,1)) - enddo - do k = 1, elec_beta_num - hij = hij + mipi(occ(k,2)) + hij = hij + mo_bielec_integral_mipi(occ(k,2),m,p) enddo else ! Mono beta m = exc(1,1,2) p = exc(1,2,2) - do k = 1, elec_beta_num - i = occ(k,2) - if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) - has_mipi(i) = .True. - endif - enddo + do k = 1, elec_alpha_num - i = occ(k,1) - if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - has_mipi(i) = .True. - endif - enddo - - do k = 1, elec_alpha_num - hij = hij + mipi(occ(k,1)) + hij = hij + mo_bielec_integral_mipi(occ(k,1),m,p) enddo do k = 1, elec_beta_num - hij = hij + mipi(occ(k,2)) - miip(occ(k,2)) + hij = hij + mo_bielec_integral_mipi_anti(occ(k,2),m,p) enddo endif @@ -651,8 +617,6 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree) integer :: occ(Nint*bit_kind_size,2) double precision :: diag_H_mat_elem integer :: n_occ_ab(2) - logical :: has_mipi(Nint*bit_kind_size) - double precision :: mipi(Nint*bit_kind_size), miip(Nint*bit_kind_size) PROVIDE mo_bielec_integrals_in_map mo_integrals_map ASSERT (Nint > 0) @@ -704,59 +668,27 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree) 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) - has_mipi = .False. if (exc(0,1,1) == 1) then ! Mono alpha m = exc(1,1,1) p = exc(1,2,1) do k = 1, elec_alpha_num - i = occ(k,1) - if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) - has_mipi(i) = .True. - endif + hij = hij + mo_bielec_integral_mipi_anti(occ(k,1),m,p) enddo do k = 1, elec_beta_num - i = occ(k,2) - if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - has_mipi(i) = .True. - endif + hij = hij + mo_bielec_integral_mipi(occ(k,2),m,p) enddo - - do k = 1, elec_alpha_num - hij = hij + mipi(occ(k,1)) - miip(occ(k,1)) - enddo - do k = 1, elec_beta_num - hij = hij + mipi(occ(k,2)) - enddo - + else ! Mono beta m = exc(1,1,2) p = exc(1,2,2) - do k = 1, elec_beta_num - i = occ(k,2) - if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) - has_mipi(i) = .True. - endif - enddo - do k = 1, elec_alpha_num - i = occ(k,1) - if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - has_mipi(i) = .True. - endif - enddo do k = 1, elec_alpha_num - hij = hij + mipi(occ(k,1)) + hij = hij + mo_bielec_integral_mipi(occ(k,1),m,p) enddo do k = 1, elec_beta_num - hij = hij + mipi(occ(k,2)) - miip(occ(k,2)) + hij = hij + mo_bielec_integral_mipi_anti(occ(k,2),m,p) enddo endif @@ -787,8 +719,6 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble) integer :: occ(Nint*bit_kind_size,2) double precision :: diag_H_mat_elem, phase,phase_2 integer :: n_occ_ab(2) - logical :: has_mipi(Nint*bit_kind_size) - double precision :: mipi(Nint*bit_kind_size), miip(Nint*bit_kind_size) PROVIDE mo_bielec_integrals_in_map mo_integrals_map ASSERT (Nint > 0) @@ -842,59 +772,26 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble) 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) - has_mipi = .False. if (exc(0,1,1) == 1) then ! Mono alpha m = exc(1,1,1) p = exc(1,2,1) do k = 1, elec_alpha_num - i = occ(k,1) - if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) - has_mipi(i) = .True. - endif + hdouble = hdouble + mo_bielec_integral_mipi_anti(occ(k,1),m,p) enddo do k = 1, elec_beta_num - i = occ(k,2) - if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - has_mipi(i) = .True. - endif - enddo - - do k = 1, elec_alpha_num - hdouble = hdouble + mipi(occ(k,1)) - miip(occ(k,1)) - enddo - do k = 1, elec_beta_num - hdouble = hdouble + mipi(occ(k,2)) + hdouble = hdouble + mo_bielec_integral_mipi(occ(k,2),m,p) enddo else ! Mono beta m = exc(1,1,2) p = exc(1,2,2) - do k = 1, elec_beta_num - i = occ(k,2) - if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) - has_mipi(i) = .True. - endif - enddo do k = 1, elec_alpha_num - i = occ(k,1) - if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - has_mipi(i) = .True. - endif - enddo - - do k = 1, elec_alpha_num - hdouble = hdouble + mipi(occ(k,1)) + hdouble = hdouble + mo_bielec_integral_mipi(occ(k,1),m,p) enddo do k = 1, elec_beta_num - hdouble = hdouble + mipi(occ(k,2)) - miip(occ(k,2)) + hdouble = hdouble + mo_bielec_integral_mipi_anti(occ(k,2),m,p) enddo endif diff --git a/src/Integrals_Bielec/map_integrals.irp.f b/src/Integrals_Bielec/map_integrals.irp.f index 22fd48a6..b41a3177 100644 --- a/src/Integrals_Bielec/map_integrals.irp.f +++ b/src/Integrals_Bielec/map_integrals.irp.f @@ -370,6 +370,7 @@ BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0:64*64*64*64) ] END_PROVIDER + double precision function get_mo_bielec_integral(i,j,k,l,map) use map_module implicit none diff --git a/src/Integrals_Bielec/mo_bi_integrals.irp.f b/src/Integrals_Bielec/mo_bi_integrals.irp.f index 0a468c24..e581b536 100644 --- a/src/Integrals_Bielec/mo_bi_integrals.irp.f +++ b/src/Integrals_Bielec/mo_bi_integrals.irp.f @@ -467,6 +467,31 @@ END_PROVIDER enddo enddo +END_PROVIDER + + BEGIN_PROVIDER [ double precision, mo_bielec_integral_mipi, (mo_tot_num_align,mo_tot_num,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, mo_bielec_integral_mipi_anti, (mo_tot_num_align,mo_tot_num,mo_tot_num) ] + implicit none + BEGIN_DOC + ! and - . Indices are (i,m,p) + END_DOC + + integer :: m,i,p + double precision :: get_mo_bielec_integral + + PROVIDE mo_bielec_integrals_in_map + + mo_bielec_integral_mipi = 0.d0 + mo_bielec_integral_mipi_anti = 0.d0 + do p=1,mo_tot_num + do m=1,mo_tot_num + do i=1,mo_tot_num + mo_bielec_integral_mipi(i,m,p) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + mo_bielec_integral_mipi_anti(i,m,p) = mo_bielec_integral_mipi(i,m,p) - get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + enddo + enddo + enddo + END_PROVIDER BEGIN_PROVIDER [ double precision, mo_bielec_integral_schwartz,(mo_tot_num,mo_tot_num) ]