diff --git a/src/ao_two_e_ints/map_integrals.irp.f b/src/ao_two_e_ints/map_integrals.irp.f index bf284841..57abe5e2 100644 --- a/src/ao_two_e_ints/map_integrals.irp.f +++ b/src/ao_two_e_ints/map_integrals.irp.f @@ -518,8 +518,12 @@ complex*16 function get_ao_two_e_integral_periodic_simple(i,j,k,l,map,map2) resu endif else call map_get(map2,idx,tmp_re) - call map_get(map2,idx+1,tmp_im) - tmp_im *= sign + if (sign/=0.d0) then + call map_get(map2,idx+1,tmp_im) + tmp_im *= sign + else + tmp_im=0.d0 + endif endif tmp = dcmplx(tmp_re,tmp_im) result = tmp diff --git a/src/mo_two_e_ints/map_integrals.irp.f b/src/mo_two_e_ints/map_integrals.irp.f index 41534372..9374ea80 100644 --- a/src/mo_two_e_ints/map_integrals.irp.f +++ b/src/mo_two_e_ints/map_integrals.irp.f @@ -41,22 +41,6 @@ subroutine insert_into_mo_integrals_map(n_integrals, & call map_update(mo_integrals_map, buffer_i, buffer_values, n_integrals, thr) end -subroutine insert_into_mo_integrals_map_2(n_integrals, & - buffer_i, buffer_values, thr) - use map_module - implicit none - - BEGIN_DOC - ! Create new entry into MO map, or accumulate in an existing entry - END_DOC - - integer, intent(in) :: n_integrals - integer(key_kind), intent(inout) :: buffer_i(n_integrals) - real(integral_kind), intent(inout) :: buffer_values(n_integrals) - real(integral_kind), intent(in) :: thr - call map_update(mo_integrals_map_2, buffer_i, buffer_values, n_integrals, thr) -end - BEGIN_PROVIDER [ integer*4, mo_integrals_cache_min ] &BEGIN_PROVIDER [ integer*4, mo_integrals_cache_max ] &BEGIN_PROVIDER [ integer*8, mo_integrals_cache_min_8 ] @@ -110,45 +94,6 @@ BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0_8:128_8*128_8*128_8*12 END_PROVIDER -BEGIN_PROVIDER [ complex*16, mo_integrals_cache_periodic, (0_8:128_8*128_8*128_8*128_8) ] - implicit none - BEGIN_DOC - ! Cache of MO integrals for fast access - END_DOC - PROVIDE mo_two_e_integrals_in_map - integer*8 :: i,j,k,l - integer*4 :: i4,j4,k4,l4 - integer*8 :: ii - integer(key_kind) :: idx - complex(integral_kind) :: integral - complex*16 :: get_mo_two_e_integrals_periodic_simple - FREE ao_integrals_cache - !$OMP PARALLEL DO PRIVATE (i,j,k,l,i4,j4,k4,l4,idx,ii,integral) - do l=mo_integrals_cache_min_8,mo_integrals_cache_max_8 - l4 = int(l,4) - do k=mo_integrals_cache_min_8,mo_integrals_cache_max_8 - k4 = int(k,4) - do j=mo_integrals_cache_min_8,mo_integrals_cache_max_8 - j4 = int(j,4) - do i=mo_integrals_cache_min_8,mo_integrals_cache_max_8 - i4 = int(i,4) - !DIR$ FORCEINLINE - integral = get_mo_two_e_integrals_periodic_simple(i,j,k,l,& - mo_integrals_map,mo_integrals_map_2) - ii = l-mo_integrals_cache_min_8 - ii = ior( shiftl(ii,7), k-mo_integrals_cache_min_8) - ii = ior( shiftl(ii,7), j-mo_integrals_cache_min_8) - ii = ior( shiftl(ii,7), i-mo_integrals_cache_min_8) - mo_integrals_cache_periodic(ii) = integral - enddo - enddo - enddo - enddo - !$OMP END PARALLEL DO - -END_PROVIDER - - double precision function get_two_e_integral(i,j,k,l,map) use map_module implicit none @@ -181,40 +126,6 @@ double precision function get_two_e_integral(i,j,k,l,map) endif end -complex*16 function get_two_e_integral_periodic(i,j,k,l,map,map2) - use map_module - implicit none -! BEGIN_DOC -! ! Returns one integral in the MO basis -! ! TODO: finish this -! END_DOC -! integer, intent(in) :: i,j,k,l -! integer(key_kind) :: idx -! integer :: ii -! integer*8 :: ii_8 -! type(map_type), intent(inout) :: map -! real(integral_kind) :: tmp -! PROVIDE mo_two_e_integrals_in_map mo_integrals_cache_periodic -! ii = l-mo_integrals_cache_min -! ii = ior(ii, k-mo_integrals_cache_min) -! ii = ior(ii, j-mo_integrals_cache_min) -! ii = ior(ii, i-mo_integrals_cache_min) -! if (iand(ii, -128) /= 0) then -! !DIR$ FORCEINLINE -! call two_e_integrals_index(i,j,k,l,idx) -! !DIR$ FORCEINLINE -! call map_get(map,idx,tmp) -! get_two_e_integral_periodic = dble(tmp) -! else -! ii_8 = int(l,8)-mo_integrals_cache_min_8 -! ii_8 = ior( shiftl(ii_8,7), int(k,8)-mo_integrals_cache_min_8) -! ii_8 = ior( shiftl(ii_8,7), int(j,8)-mo_integrals_cache_min_8) -! ii_8 = ior( shiftl(ii_8,7), int(i,8)-mo_integrals_cache_min_8) -! get_two_e_integral = mo_integrals_cache(ii_8) -! endif -end - - double precision function mo_two_e_integral(i,j,k,l) implicit none BEGIN_DOC @@ -462,13 +373,15 @@ subroutine get_mo_two_e_integrals_exch_ii(k,l,sze,out_val,map) endif end - integer*8 function get_mo_map_size() implicit none BEGIN_DOC ! Return the number of elements in the MO map END_DOC get_mo_map_size = mo_integrals_map % n_elements + if (is_periodic) then + get_mo_map_size += mo_integrals_map_2 % n_elements + endif end diff --git a/src/mo_two_e_ints/map_integrals_complex.irp.f b/src/mo_two_e_ints/map_integrals_complex.irp.f new file mode 100644 index 00000000..7cf4ba0d --- /dev/null +++ b/src/mo_two_e_ints/map_integrals_complex.irp.f @@ -0,0 +1,498 @@ +use map_module + +subroutine insert_into_mo_integrals_map_2(n_integrals, & + buffer_i, buffer_values, thr) + use map_module + implicit none + + BEGIN_DOC + ! Create new entry into MO map, or accumulate in an existing entry + END_DOC + + integer, intent(in) :: n_integrals + integer(key_kind), intent(inout) :: buffer_i(n_integrals) + real(integral_kind), intent(inout) :: buffer_values(n_integrals) + real(integral_kind), intent(in) :: thr + call map_update(mo_integrals_map_2, buffer_i, buffer_values, n_integrals, thr) +end + +BEGIN_PROVIDER [ complex*16, mo_integrals_cache_periodic, (0_8:128_8*128_8*128_8*128_8) ] + implicit none + BEGIN_DOC + ! Cache of MO integrals for fast access + END_DOC + PROVIDE mo_two_e_integrals_in_map + integer*8 :: i,j,k,l + integer*4 :: i4,j4,k4,l4 + integer*8 :: ii + integer(key_kind) :: idx + complex(integral_kind) :: integral + complex*16 :: get_two_e_integral_periodic_simple + FREE ao_integrals_cache + !$OMP PARALLEL DO PRIVATE (i,j,k,l,i4,j4,k4,l4,idx,ii,integral) + do l=mo_integrals_cache_min_8,mo_integrals_cache_max_8 + l4 = int(l,4) + do k=mo_integrals_cache_min_8,mo_integrals_cache_max_8 + k4 = int(k,4) + do j=mo_integrals_cache_min_8,mo_integrals_cache_max_8 + j4 = int(j,4) + do i=mo_integrals_cache_min_8,mo_integrals_cache_max_8 + i4 = int(i,4) + !DIR$ FORCEINLINE + integral = get_two_e_integral_periodic_simple(i,j,k,l,& + mo_integrals_map,mo_integrals_map_2) + ii = l-mo_integrals_cache_min_8 + ii = ior( shiftl(ii,7), k-mo_integrals_cache_min_8) + ii = ior( shiftl(ii,7), j-mo_integrals_cache_min_8) + ii = ior( shiftl(ii,7), i-mo_integrals_cache_min_8) + mo_integrals_cache_periodic(ii) = integral + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO + +END_PROVIDER + + +complex*16 function get_two_e_integral_periodic_simple(i,j,k,l,map,map2) result(result) + use map_module + implicit none + BEGIN_DOC + ! Gets one MO bi-electronic integral from the MO map + ! reuse ao map/idx/sign function + END_DOC + integer, intent(in) :: i,j,k,l + integer(key_kind) :: idx + real(integral_kind) :: tmp_re, tmp_im + type(map_type), intent(inout) :: map,map2 + complex(integral_kind) :: tmp + logical :: use_map1 + double precision :: sign + PROVIDE mo_two_e_integrals_in_map + call ao_two_e_integral_periodic_map_idx_sign(i,j,k,l,use_map1,idx,sign) + if (use_map1) then + call map_get(map,idx,tmp_re) + if (sign/=0.d0) then + call map_get(map,idx+1,tmp_im) + tmp_im *= sign + else + tmp_im=0.d0 + endif + else + call map_get(map2,idx,tmp_re) + if (sign/=0.d0) then + call map_get(map2,idx+1,tmp_im) + tmp_im *= sign + else + tmp_im=0.d0 + endif + endif + tmp = dcmplx(tmp_re,tmp_im) + result = tmp +end + +complex*16 function get_two_e_integral_periodic(i,j,k,l,map,map2) + use map_module + implicit none + BEGIN_DOC + ! Returns one integral in the MO basis + ! TODO: finish this + END_DOC + integer, intent(in) :: i,j,k,l + integer(key_kind) :: idx + integer :: ii + integer*8 :: ii_8 + type(map_type), intent(inout) :: map,map2 + complex(integral_kind) :: tmp + complex(integral_kind) :: get_two_e_integral_periodic_simple + PROVIDE mo_two_e_integrals_in_map mo_integrals_cache_periodic + ii = l-mo_integrals_cache_min + ii = ior(ii, k-mo_integrals_cache_min) + ii = ior(ii, j-mo_integrals_cache_min) + ii = ior(ii, i-mo_integrals_cache_min) + if (iand(ii, -128) /= 0) then + tmp = get_two_e_integral_periodic_simple(i,j,k,l,map,map2) + else + ii_8 = int(l,8)-mo_integrals_cache_min_8 + ii_8 = ior( shiftl(ii_8,7), int(k,8)-mo_integrals_cache_min_8) + ii_8 = ior( shiftl(ii_8,7), int(j,8)-mo_integrals_cache_min_8) + ii_8 = ior( shiftl(ii_8,7), int(i,8)-mo_integrals_cache_min_8) + tmp = mo_integrals_cache_periodic(ii_8) + endif + get_two_e_integral_periodic = tmp +end + +complex*16 function mo_two_e_integral_periodic(i,j,k,l) + implicit none + BEGIN_DOC + ! Returns one integral in the MO basis + END_DOC + integer, intent(in) :: i,j,k,l + complex*16 :: get_two_e_integral_periodic + PROVIDE mo_two_e_integrals_in_map mo_integrals_cache_periodic + PROVIDE mo_two_e_integrals_in_map + !DIR$ FORCEINLINE + mo_two_e_integral_periodic = get_two_e_integral_periodic(i,j,k,l,mo_integrals_map,mo_integrals_map_2) + return +end + +subroutine get_mo_two_e_integrals_periodic(j,k,l,sze,out_val,map,map2) + use map_module + implicit none + BEGIN_DOC + ! Returns multiple integrals in the MO basis, all + ! i for j,k,l fixed. + END_DOC + integer, intent(in) :: j,k,l, sze + complex*16, intent(out) :: out_val(sze) + type(map_type), intent(inout) :: map,map2 + integer :: i + complex*16, external :: get_two_e_integral_periodic_simple + + integer :: ii, ii0 + integer*8 :: ii_8, ii0_8 + complex(integral_kind) :: tmp + integer(key_kind) :: i1, idx + integer(key_kind) :: p,q,r,s,i2 + PROVIDE mo_two_e_integrals_in_map mo_integrals_cache_periodic + +!DEBUG +! do i=1,sze +! out_val(i) = get_two_e_integral_periodic(i,j,k,l,map,map2) +! enddo +! return +!DEBUG + + ii0 = l-mo_integrals_cache_min + ii0 = ior(ii0, k-mo_integrals_cache_min) + ii0 = ior(ii0, j-mo_integrals_cache_min) + + ii0_8 = int(l,8)-mo_integrals_cache_min_8 + ii0_8 = ior( shiftl(ii0_8,7), int(k,8)-mo_integrals_cache_min_8) + ii0_8 = ior( shiftl(ii0_8,7), int(j,8)-mo_integrals_cache_min_8) + + do i=1,sze + ii = ior(ii0, i-mo_integrals_cache_min) + if (iand(ii, -128) == 0) then + ii_8 = ior( shiftl(ii0_8,7), int(i,8)-mo_integrals_cache_min_8) + out_val(i) = mo_integrals_cache_periodic(ii_8) + else + out_val(i) = get_two_e_integral_periodic_simple(i,j,k,l,map,map2) + endif + enddo +end + +!subroutine get_mo_two_e_integrals_ij_periodic(k,l,sze,out_array,map) +! use map_module +! implicit none +! BEGIN_DOC +! ! Returns multiple integrals in the MO basis, all +! ! i(1)j(2) 1/r12 k(1)l(2) +! ! i, j for k,l fixed. +! END_DOC +! integer, intent(in) :: k,l, sze +! double precision, intent(out) :: out_array(sze,sze) +! type(map_type), intent(inout) :: map +! integer :: i,j,kk,ll,m +! integer(key_kind),allocatable :: hash(:) +! integer ,allocatable :: pairs(:,:), iorder(:) +! real(integral_kind), allocatable :: tmp_val(:) +! +! PROVIDE mo_two_e_integrals_in_map +! allocate (hash(sze*sze), pairs(2,sze*sze),iorder(sze*sze), & +! tmp_val(sze*sze)) +! +! kk=0 +! out_array = 0.d0 +! do j=1,sze +! do i=1,sze +! kk += 1 +! !DIR$ FORCEINLINE +! call two_e_integrals_index(i,j,k,l,hash(kk)) +! pairs(1,kk) = i +! pairs(2,kk) = j +! iorder(kk) = kk +! enddo +! enddo +! +! logical :: integral_is_in_map +! if (key_kind == 8) then +! call i8radix_sort(hash,iorder,kk,-1) +! else if (key_kind == 4) then +! call iradix_sort(hash,iorder,kk,-1) +! else if (key_kind == 2) then +! call i2radix_sort(hash,iorder,kk,-1) +! endif +! +! call map_get_many(mo_integrals_map, hash, tmp_val, kk) +! +! do ll=1,kk +! m = iorder(ll) +! i=pairs(1,m) +! j=pairs(2,m) +! out_array(i,j) = tmp_val(ll) +! enddo +! +! deallocate(pairs,hash,iorder,tmp_val) +!end + +!subroutine get_mo_two_e_integrals_i1j1_periodic(k,l,sze,out_array,map) +! use map_module +! implicit none +! BEGIN_DOC +! ! Returns multiple integrals in the MO basis, all +! ! i(1)j(1) 1/r12 k(2)l(2) +! ! i, j for k,l fixed. +! END_DOC +! integer, intent(in) :: k,l, sze +! double precision, intent(out) :: out_array(sze,sze) +! type(map_type), intent(inout) :: map +! integer :: i,j,kk,ll,m +! integer(key_kind),allocatable :: hash(:) +! integer ,allocatable :: pairs(:,:), iorder(:) +! real(integral_kind), allocatable :: tmp_val(:) +! +! PROVIDE mo_two_e_integrals_in_map +! allocate (hash(sze*sze), pairs(2,sze*sze),iorder(sze*sze), & +! tmp_val(sze*sze)) +! +! kk=0 +! out_array = 0.d0 +! do j=1,sze +! do i=1,sze +! kk += 1 +! !DIR$ FORCEINLINE +! call two_e_integrals_index(i,k,j,l,hash(kk)) +! pairs(1,kk) = i +! pairs(2,kk) = j +! iorder(kk) = kk +! enddo +! enddo +! +! logical :: integral_is_in_map +! if (key_kind == 8) then +! call i8radix_sort(hash,iorder,kk,-1) +! else if (key_kind == 4) then +! call iradix_sort(hash,iorder,kk,-1) +! else if (key_kind == 2) then +! call i2radix_sort(hash,iorder,kk,-1) +! endif +! +! call map_get_many(mo_integrals_map, hash, tmp_val, kk) +! +! do ll=1,kk +! m = iorder(ll) +! i=pairs(1,m) +! j=pairs(2,m) +! out_array(i,j) = tmp_val(ll) +! enddo +! +! deallocate(pairs,hash,iorder,tmp_val) +!end + +subroutine get_mo_two_e_integrals_coulomb_ii_periodic(k,l,sze,out_val,map) + use map_module + implicit none + BEGIN_DOC + ! Returns multiple integrals + ! k(1)i(2) 1/r12 l(1)i(2) :: out_val(i1) + ! for k,l fixed. + ! always in map1, take conjugate if k>l, real if k==l + ! TODO: determine best way to structure code + ! to account for single/double integral_kind, real/complex, and +/- imag part + END_DOC + integer, intent(in) :: k,l, sze + complex*16, intent(out) :: out_val(sze) + type(map_type), intent(inout) :: map + integer :: i + integer(key_kind) :: hash(sze),hash_re(sze),hash_im(sze) + real(integral_kind) :: tmp_re(sze),tmp_im(sze) + complex*16 :: out_re(sze),out_im(sze) + double precision :: sign + PROVIDE mo_two_e_integrals_in_map + + if (k.eq.l) then ! real, call other function + call get_mo_two_e_integrals_coulomb_ijij_periodic(k,sze,out_re,map) + do i=1,sze + out_val(i) = dcmplx(out_re(i),0.d0) + enddo + else ! complex + if (k.gt.l) then + sign = -1.d0 + else + sign = 1.d0 + endif + + do i=1,sze + !DIR$ FORCEINLINE + call two_e_integrals_index(k,i,l,i,hash(i)) + !hash_im(i) = hash(i)*2 + hash_im(i) = shiftl(hash(i),1) + hash_re(i) = hash_im(i)-1 + enddo + + if (integral_kind == 8) then + call map_get_many(map, hash_re, out_re, sze) + call map_get_many(map, hash_im, out_im, sze) + do i=1,sze + out_val(i) = dcmplx(out_re(i),sign*out_im(i)) + enddo + else + call map_get_many(map, hash_re, tmp_re, sze) + call map_get_many(map, hash_im, tmp_im, sze) + ! Conversion to double complex + do i=1,sze + out_val(i) = dcmplx(tmp_re(i),sign*tmp_im(i)) + enddo + endif + endif +end + +subroutine get_mo_two_e_integrals_coulomb_ijij_periodic(j,sze,out_val,map) + use map_module + implicit none + BEGIN_DOC + ! Returns multiple integrals + ! i*(1)j*(2) 1/r12 i(1)j(2) :: out_val(i) + ! for j fixed. + ! always in map1, always real + END_DOC + integer, intent(in) :: j, sze + double precision, intent(out) :: out_val(sze) + type(map_type), intent(inout) :: map + integer :: i + integer(key_kind) :: hash(sze),hash_re(sze) + real(integral_kind) :: tmp_re(sze) + PROVIDE mo_two_e_integrals_in_map + + do i=1,sze + !DIR$ FORCEINLINE + call two_e_integrals_index(i,j,i,j,hash(i)) + !hash_re(i) = hash(i)*2 - 1 + hash_re(i) = shiftl(hash(i),1) - 1 + enddo + + if (integral_kind == 8) then + call map_get_many(map, hash_re, out_val, sze) + else + call map_get_many(map, hash_re, tmp_re, sze) + ! Conversion to double complex + do i=1,sze + out_val(i) = dble(tmp_re(i)) + enddo + endif +end + +subroutine get_mo_two_e_integrals_exch_ii_periodic(k,l,sze,out_val,map,map2) + use map_module + implicit none + BEGIN_DOC + ! Returns multiple integrals + ! k*(1)i*(2) 1/r12 i(1)l(2) :: out_val(i1) + ! for k,l fixed. + ! + ! if k + ! i*(1)j*(2) 1/r12 j(1)i(2) :: out_val(i) + ! for j fixed. + ! always real, always in map2 (except when j==i, then in map1) + END_DOC + integer, intent(in) :: j, sze + double precision, intent(out) :: out_val(sze) + type(map_type), intent(inout) :: map,map2 + integer :: i + integer(key_kind) :: hash(sze),hash_re(sze) + real(integral_kind) :: tmp_val(sze) + PROVIDE mo_two_e_integrals_in_map + + do i=1,sze + !DIR$ FORCEINLINE + call two_e_integrals_index(i,j,j,i,hash(i)) + !hash_re(i) = 2*hash(i) - 1 + hash_re(i) = shiftl(hash(i),1) - 1 + enddo + + if (integral_kind == 8) then + call map_get_many(map2, hash_re, out_val, sze) + call map_get(map,hash_re(j), out_val(j)) + else + call map_get_many(map2, hash_re, tmp_val, sze) + call map_get(map, hash_re(j), tmp_val(j)) + ! Conversion to double precision + do i=1,sze + out_val(i) = dble(tmp_val(i)) + enddo + endif +end +