From eb25d52e22391ad7a3a40cc9faad72e6ad1ec2de Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sat, 20 Sep 2014 23:45:18 +0200 Subject: [PATCH] Not working --- src/BiInts/map_integrals.irp.f | 63 ++++++++++++++++++++--- src/Hartree_Fock/Fock_matrix.irp.f | 80 +++++++++++++----------------- 2 files changed, 91 insertions(+), 52 deletions(-) diff --git a/src/BiInts/map_integrals.irp.f b/src/BiInts/map_integrals.irp.f index a638bc39..d15bb8e3 100644 --- a/src/BiInts/map_integrals.irp.f +++ b/src/BiInts/map_integrals.irp.f @@ -32,21 +32,72 @@ end subroutine bielec_integrals_index_reverse(i,j,k,l,i1) implicit none - integer, intent(out) :: i,j,k,l + integer, intent(out) :: i(8),j(8),k(8),l(8) integer*8, intent(in) :: i1 integer*8 :: i2,i3 real :: x + i = 0 x = 0.5*(sqrt(8.*real(i1)+1.)-1.) i2 = ceiling(x) i3 = i1 - ishft(i2*i2-i2,-1) - l = 0.5*(sqrt(8.*real(i2)+1.)-1.) - l = ceiling(x) - j = i2 - ishft(l*l-l,-1) + x = 0.5*(sqrt(8.*real(i2)+1.)-1.) + l(1) = ceiling(x) + j(1) = i2 - ishft(l(1)*l(1)-l(1),-1) x = 0.5*(sqrt(8.*real(i3)+1.)-1.) - k = ceiling(x) - i = i3 - ishft(k*k-k,-1) + k(1) = ceiling(x) + i(1) = i3 - ishft(k(1)*k(1)-k(1),-1) + + !ijkl + i(2) = k(1) !kjil + j(2) = j(1) + k(2) = i(1) + l(2) = l(1) + + i(3) = i(1) !ilkj + j(3) = l(1) + k(3) = k(1) + l(3) = j(1) + + i(4) = k(1) !klij + j(4) = l(1) + k(4) = i(1) + l(4) = j(1) + + i(5) = j(1) !jilk + j(5) = i(1) + k(5) = l(1) + l(5) = k(1) + + i(6) = l(1) !lijk + j(6) = i(1) + k(6) = j(1) + l(6) = k(1) + + i(7) = j(1) !jkli + j(7) = k(1) + k(7) = l(1) + l(7) = i(1) + + i(8) = l(1) !lkji + j(8) = k(1) + k(8) = j(1) + l(8) = i(1) + + integer :: ii, jj + do ii=2,8 + do jj=1,ii-1 + if ( (i(ii) == i(jj)).and. & + (j(ii) == j(jj)).and. & + (k(ii) == k(jj)).and. & + (l(ii) == l(jj)) ) then + i(jj) = 0 + exit + endif + enddo + enddo + end diff --git a/src/Hartree_Fock/Fock_matrix.irp.f b/src/Hartree_Fock/Fock_matrix.irp.f index 05148937..2fd9e841 100644 --- a/src/Hartree_Fock/Fock_matrix.irp.f +++ b/src/Hartree_Fock/Fock_matrix.irp.f @@ -101,6 +101,7 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, ao_bi_elec_integral_alpha, (ao_num_align, ao_num) ] &BEGIN_PROVIDER [ double precision, ao_bi_elec_integral_beta , (ao_num_align, ao_num) ] + use map_module implicit none BEGIN_DOC ! Alpha Fock matrix in AO basis set @@ -159,7 +160,6 @@ END_PROVIDER !$OMP SHARED(ao_num,ao_bi_elec_integral_alpha,ao_mono_elec_integral,& !$OMP ao_num_align,ao_bi_elec_integral_beta,HF_density_matrix_ao_alpha, & !$OMP HF_density_matrix_ao_beta,lck) - allocate(ao_ints_idx(ao_num_align),ao_ints_val(ao_num_align)) !$OMP DO do j=1,ao_num !DIR$ VECTOR ALIGNED @@ -169,54 +169,42 @@ END_PROVIDER enddo enddo !$OMP END DO - !$OMP DO SCHEDULE(GUIDED) - do j=1,ao_num - do l=1,ao_num - do i=1,j-1 - call get_ao_bielec_integrals_non_zero(i,l,j,ao_num,ao_ints_val,ao_ints_idx,kmax) - call OMP_SET_LOCK(lck(i)) - call OMP_SET_LOCK(lck(j)) - do k1=1,kmax - k = ao_ints_idx(k1) - integral = (HF_density_matrix_ao_alpha(k,l)+HF_density_matrix_ao_beta(k,l)) * ao_ints_val(k1) - ao_bi_elec_integral_alpha(i,j) += integral - ao_bi_elec_integral_beta (i,j) += integral - ao_bi_elec_integral_alpha(j,i) += integral - ao_bi_elec_integral_beta (j,i) += integral - enddo - do k1=1,kmax - k = ao_ints_idx(k1) - integral = ao_ints_val(k1) - ao_bi_elec_integral_alpha(l,j) -= HF_density_matrix_ao_alpha(k,i) * integral - ao_bi_elec_integral_beta (l,j) -= HF_density_matrix_ao_beta (k,i) * integral - ao_bi_elec_integral_alpha(l,i) -= HF_density_matrix_ao_alpha(k,j) * integral - ao_bi_elec_integral_beta (l,i) -= HF_density_matrix_ao_beta (k,j) * integral - enddo - call OMP_UNSET_LOCK(lck(i)) - call OMP_UNSET_LOCK(lck(j)) - enddo - do i=j,j - call get_ao_bielec_integrals_non_zero(i,l,j,ao_num,ao_ints_val,ao_ints_idx,kmax) - call OMP_SET_LOCK(lck(i)) - do k1=1,kmax - k = ao_ints_idx(k1) - integral = (HF_density_matrix_ao_alpha(k,l)+HF_density_matrix_ao_beta(k,l)) * ao_ints_val(k1) - ao_bi_elec_integral_alpha(i,j) += integral - ao_bi_elec_integral_beta (i,j) += integral - enddo - do k1=1,kmax - k = ao_ints_idx(k1) - integral = ao_ints_val(k1) - ao_bi_elec_integral_alpha(l,j) -= HF_density_matrix_ao_alpha(k,i) * integral - ao_bi_elec_integral_beta (l,j) -= HF_density_matrix_ao_beta (k,i) * integral - enddo - call OMP_UNSET_LOCK(lck(i)) + !$OMP END PARALLEL + + integer(cache_map_size_kind) :: n_elements_max, n_elements + integer(key_kind), allocatable :: keys(:) + double precision, allocatable :: values(:) + call get_cache_map_n_elements_max(ao_integrals_map,n_elements_max) + allocate(keys(n_elements_max), values(n_elements_max)) + + integer*8 :: i8 + integer :: ii(8), jj(8), kk(8), ll(8), k2 + + do i8=1,ao_integrals_map%map_size + n_elements = n_elements_max + call get_cache_map(ao_integrals_map,i8,keys,values,n_elements) + if (n_elements == 0) then + cycle + endif + do k1=1,n_elements + call bielec_integrals_index_reverse(kk,ii,ll,jj,keys(k1)) + do k2=1,8 + if (kk(k2) == 0) cycle + i = ii(k2) + j = jj(k2) + k = kk(k2) + l = ll(k2) + integral = (HF_density_matrix_ao_alpha(k,l)+HF_density_matrix_ao_beta(k,l)) * values(k1) + ao_bi_elec_integral_alpha(i,j) += integral + ao_bi_elec_integral_beta (i,j) += integral +! integral = values(k1) +! ao_bi_elec_integral_alpha(i,l) -= HF_density_matrix_ao_alpha(k,l) * integral +! ao_bi_elec_integral_beta (i,l) -= HF_density_matrix_ao_beta (k,l) * integral enddo enddo enddo - !$OMP END DO - deallocate(ao_ints_val,ao_ints_idx) - !$OMP END PARALLEL + deallocate(keys,values) + do i=1,ao_num call omp_destroy_lock(lck(i)) enddo