From 90c3db31036b3aeffa3b94445960f4d092c6f929 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 11 Jun 2024 14:38:50 +0200 Subject: [PATCH] Accelerated cache --- src/mo_two_e_ints/cholesky.irp.f | 2 +- src/mo_two_e_ints/four_idx_novvvv.irp.f | 180 ------------------------ src/mo_two_e_ints/map_integrals.irp.f | 163 ++++++++++++++++----- src/mo_two_e_ints/mo_bi_integrals.irp.f | 4 - 4 files changed, 128 insertions(+), 221 deletions(-) delete mode 100644 src/mo_two_e_ints/four_idx_novvvv.irp.f diff --git a/src/mo_two_e_ints/cholesky.irp.f b/src/mo_two_e_ints/cholesky.irp.f index 773d240a..7e2c8b37 100644 --- a/src/mo_two_e_ints/cholesky.irp.f +++ b/src/mo_two_e_ints/cholesky.irp.f @@ -4,7 +4,7 @@ BEGIN_PROVIDER [ logical, do_mo_cholesky ] ! If True, use Cholesky vectors for MO integrals END_DOC do_mo_cholesky = do_ao_cholesky - do_mo_cholesky = .False. +! do_mo_cholesky = .False. END_PROVIDER BEGIN_PROVIDER [ integer, cholesky_mo_num ] diff --git a/src/mo_two_e_ints/four_idx_novvvv.irp.f b/src/mo_two_e_ints/four_idx_novvvv.irp.f deleted file mode 100644 index 80af35dc..00000000 --- a/src/mo_two_e_ints/four_idx_novvvv.irp.f +++ /dev/null @@ -1,180 +0,0 @@ -BEGIN_PROVIDER [ double precision, mo_coef_novirt, (ao_num,n_core_inact_act_orb) ] - implicit none - BEGIN_DOC - ! MO coefficients without virtual MOs - END_DOC - integer :: j,jj - - do j=1,n_core_inact_act_orb - jj = list_core_inact_act(j) - mo_coef_novirt(:,j) = mo_coef(:,jj) - enddo - -END_PROVIDER - -subroutine ao_to_mo_novirt(A_ao,LDA_ao,A_mo,LDA_mo) - implicit none - BEGIN_DOC - ! Transform A from the |AO| basis to the |MO| basis excluding virtuals - ! - ! $C^\dagger.A_{ao}.C$ - END_DOC - integer, intent(in) :: LDA_ao,LDA_mo - double precision, intent(in) :: A_ao(LDA_ao,ao_num) - double precision, intent(out) :: A_mo(LDA_mo,n_core_inact_act_orb) - double precision, allocatable :: T(:,:) - - allocate ( T(ao_num,n_core_inact_act_orb) ) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T - - call dgemm('N','N', ao_num, n_core_inact_act_orb, ao_num, & - 1.d0, A_ao,LDA_ao, & - mo_coef_novirt, size(mo_coef_novirt,1), & - 0.d0, T, size(T,1)) - - call dgemm('T','N', n_core_inact_act_orb, n_core_inact_act_orb, ao_num,& - 1.d0, mo_coef_novirt,size(mo_coef_novirt,1), & - T, ao_num, & - 0.d0, A_mo, size(A_mo,1)) - - deallocate(T) -end - - -subroutine four_idx_novvvv - print*,'********' - print*,'********' - print*,'********' - print*,'WARNING :: Using four_idx_novvvv, and we are not sure that this routine is not bugged ...' - print*,'********' - print*,'********' - print*,'********' - use map_module - implicit none - BEGIN_DOC - ! Retransform MO integrals for next CAS-SCF step - END_DOC - print*,'Using partial transformation' - print*,'It will not transform all integrals with at least 3 indices within the virtuals' - integer :: i,j,k,l,n_integrals - double precision, allocatable :: f(:,:,:), f2(:,:,:), d(:,:), T(:,:,:,:), T2(:,:,:,:) - double precision, external :: get_ao_two_e_integral - integer(key_kind), allocatable :: idx(:) - real(integral_kind), allocatable :: values(:) - - integer :: p,q,r,s - double precision :: c - allocate( T(n_core_inact_act_orb,n_core_inact_act_orb,ao_num,ao_num) , & - T2(n_core_inact_act_orb,n_core_inact_act_orb,ao_num,ao_num) ) - - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP SHARED(mo_num,ao_num,T,n_core_inact_act_orb, mo_coef_transp, & - !$OMP mo_integrals_threshold,mo_coef,mo_integrals_map, & - !$OMP list_core_inact_act,T2,ao_integrals_map) & - !$OMP PRIVATE(i,j,k,l,p,q,r,s,idx,values,n_integrals, & - !$OMP f,f2,d,c) - allocate(f(ao_num,ao_num,ao_num), f2(ao_num,ao_num,ao_num), d(mo_num,mo_num), & - idx(mo_num*mo_num), values(mo_num*mo_num) ) - - ! - !$OMP DO - do s=1,ao_num - do r=1,ao_num - do q=1,ao_num - do p=1,r - f (p,q,r) = get_ao_two_e_integral(p,q,r,s,ao_integrals_map) - f (r,q,p) = f(p,q,r) - enddo - enddo - enddo - do r=1,ao_num - do q=1,ao_num - do p=1,ao_num - f2(p,q,r) = f(p,r,q) - enddo - enddo - enddo - ! f (p,q,r) = - ! f2(p,q,r) = - - do r=1,ao_num - call ao_to_mo_novirt(f (1,1,r),size(f ,1),T (1,1,r,s),size(T,1)) - call ao_to_mo_novirt(f2(1,1,r),size(f2,1),T2(1,1,r,s),size(T,1)) - enddo - ! T (i,j,p,q) = - ! T2(i,j,p,q) = - - enddo - !$OMP END DO - - !$OMP DO - do j=1,n_core_inact_act_orb - do i=1,n_core_inact_act_orb - do s=1,ao_num - do r=1,ao_num - f (r,s,1) = T (i,j,r,s) - f2(r,s,1) = T2(i,j,r,s) - enddo - enddo - call ao_to_mo(f ,size(f ,1),d,size(d,1)) - n_integrals = 0 - do l=1,mo_num - do k=1,mo_num - n_integrals+=1 - call two_e_integrals_index(list_core_inact_act(i),list_core_inact_act(j),k,l,idx(n_integrals)) - values(n_integrals) = d(k,l) - enddo - enddo - call map_append(mo_integrals_map, idx, values, n_integrals) - - call ao_to_mo(f2,size(f2,1),d,size(d,1)) - n_integrals = 0 - do l=1,mo_num - do k=1,mo_num - n_integrals+=1 - call two_e_integrals_index(list_core_inact_act(i),k,list_core_inact_act(j),l,idx(n_integrals)) - values(n_integrals) = d(k,l) - enddo - enddo - call map_append(mo_integrals_map, idx, values, n_integrals) - enddo - enddo - !$OMP END DO - deallocate(f,f2,d,idx,values) - - !$OMP END PARALLEL - - deallocate(T,T2) - - - call map_sort(mo_integrals_map) - call map_unique(mo_integrals_map) - call map_shrink(mo_integrals_map,real(mo_integrals_threshold,integral_kind)) - -end - -subroutine four_idx_novvvv2 - use bitmasks - implicit none - integer :: i - integer(bit_kind) :: mask_ijkl(N_int,4) - - print*, '' - do i = 1,N_int - mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) - mask_ijkl(i,2) = full_ijkl_bitmask_4(i,1) - mask_ijkl(i,3) = core_inact_act_bitmask_4(i,1) - mask_ijkl(i,4) = full_ijkl_bitmask_4(i,1) - enddo - call add_integrals_to_map(mask_ijkl) - - print*, '' - do i = 1,N_int - mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) - mask_ijkl(i,2) = core_inact_act_bitmask_4(i,1) - mask_ijkl(i,3) = virt_bitmask(i,1) - mask_ijkl(i,4) = virt_bitmask(i,1) - enddo - call add_integrals_to_map(mask_ijkl) - -end diff --git a/src/mo_two_e_ints/map_integrals.irp.f b/src/mo_two_e_ints/map_integrals.irp.f index 571de573..90257bbd 100644 --- a/src/mo_two_e_ints/map_integrals.irp.f +++ b/src/mo_two_e_ints/map_integrals.irp.f @@ -48,6 +48,8 @@ end mo_integrals_cache_shift = 7 ! 7 = log(128). Max 7 endif +!mo_integrals_cache_shift = 2 ! 5 = log(32). + mo_integrals_cache_size = 2**mo_integrals_cache_shift mo_integrals_cache_min = max(1,elec_alpha_num - (mo_integrals_cache_size/2 - 1) ) @@ -112,6 +114,24 @@ BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0:mo_integrals_cache_siz END_PROVIDER +double precision function get_two_e_integral_cache(i,j,k,l) + use map_module + implicit none + BEGIN_DOC + ! Returns one integral in the MO basis taken from the cache + END_DOC + integer, intent(in) :: i,j,k,l + integer :: ii + + ii = l-mo_integrals_cache_min + ii = ior( shiftl(ii,mo_integrals_cache_shift), k-mo_integrals_cache_min) + ii = ior( shiftl(ii,mo_integrals_cache_shift), j-mo_integrals_cache_min) + ii = ior( shiftl(ii,mo_integrals_cache_shift), i-mo_integrals_cache_min) + get_two_e_integral_cache = mo_integrals_cache(ii) + +end + + double precision function get_two_e_integral(i,j,k,l,map) use map_module implicit none @@ -123,7 +143,6 @@ double precision function get_two_e_integral(i,j,k,l,map) integer :: ii type(map_type), intent(inout) :: map real(integral_kind) :: tmp - integer :: kk PROVIDE mo_two_e_integrals_in_map mo_integrals_cache do_mo_cholesky @@ -145,13 +164,9 @@ double precision function get_two_e_integral(i,j,k,l,map) ii = ior(ii, i-mo_integrals_cache_min) if (iand(ii, -mo_integrals_cache_size) == 0) then - ! Integrals is in the cache - ii = l-mo_integrals_cache_min - ii = ior( shiftl(ii,mo_integrals_cache_shift), k-mo_integrals_cache_min) - ii = ior( shiftl(ii,mo_integrals_cache_shift), j-mo_integrals_cache_min) - ii = ior( shiftl(ii,mo_integrals_cache_shift), i-mo_integrals_cache_min) - get_two_e_integral = mo_integrals_cache(ii) + double precision, external :: get_two_e_integral_cache + get_two_e_integral = get_two_e_integral_cache(i,j,k,l) else @@ -199,7 +214,6 @@ subroutine get_mo_two_e_integrals(j,k,l,sze,out_val,map) double precision, intent(out) :: out_val(sze) type(map_type), intent(inout) :: map integer :: i - double precision, external :: get_two_e_integral integer :: ii real(integral_kind) :: tmp @@ -256,13 +270,7 @@ subroutine get_mo_two_e_integrals(j,k,l,sze,out_val,map) endif - - ii = l-mo_integrals_cache_min - ii = ior( shiftl(ii, mo_integrals_cache_shift), k-mo_integrals_cache_min) - ii = ior( shiftl(ii, mo_integrals_cache_shift), j-mo_integrals_cache_min) - ii = shiftl(ii, mo_integrals_cache_shift) - out_val(mo_integrals_cache_min:mo_integrals_cache_max) = & - mo_integrals_cache(ii:ii+mo_integrals_cache_max-mo_integrals_cache_min) + call get_mo_two_e_integrals_cache(j,k,l,sze,out_val) if (mo_integrals_cache_max < mo_num) then @@ -333,6 +341,26 @@ subroutine get_mo_two_e_integrals(j,k,l,sze,out_val,map) end +subroutine get_mo_two_e_integrals_cache(j,k,l,sze,out_val) + use map_module + implicit none + BEGIN_DOC + ! Returns multiple integrals in the MO basis, all + ! i for j,k,l fixed, all integrals from the cache + END_DOC + integer, intent(in) :: j,k,l, sze + double precision, intent(out) :: out_val(sze) + integer :: ii + + ii = l-mo_integrals_cache_min + ii = ior( shiftl(ii, mo_integrals_cache_shift), k-mo_integrals_cache_min) + ii = ior( shiftl(ii, mo_integrals_cache_shift), j-mo_integrals_cache_min) + ii = shiftl(ii, mo_integrals_cache_shift) + out_val(mo_integrals_cache_min:mo_integrals_cache_max) = & + mo_integrals_cache(ii:ii+mo_integrals_cache_max-mo_integrals_cache_min) + +end + subroutine get_mo_two_e_integrals_ij(k,l,sze,out_array,map) use map_module implicit none @@ -347,16 +375,32 @@ subroutine get_mo_two_e_integrals_ij(k,l,sze,out_array,map) integer :: j real(integral_kind), allocatable :: tmp_val(:) - if (do_mo_cholesky) then - call dgemm('T', 'N', mo_num, mo_num, cholesky_mo_num, 1.d0, & - cholesky_mo_transp(1,1,k), cholesky_mo_num, & - cholesky_mo_transp(1,1,l), cholesky_mo_num, 0.d0, & - out_array, sze) + if ( (mo_integrals_cache_min>1).or.(mo_integrals_cache_max1).or.(mo_integrals_cache_max1).or.(mo_integrals_cache_max1).or.(mo_integrals_cache_max