mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-08-15 08:48:30 +02:00
![Kevin Gasperich](/assets/img/avatar_default.png)
added functions to get MO 2e ints still need routines to get multiple ints reused some functions from AO 2e ints
466 lines
13 KiB
Fortran
466 lines
13 KiB
Fortran
use map_module
|
|
|
|
!! MO Map
|
|
!! ======
|
|
|
|
BEGIN_PROVIDER [ type(map_type), mo_integrals_map ]
|
|
&BEGIN_PROVIDER [ type(map_type), mo_integrals_map_2 ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! MO integrals
|
|
END_DOC
|
|
integer(key_kind) :: key_max
|
|
integer(map_size_kind) :: sze
|
|
call two_e_integrals_index(mo_num,mo_num,mo_num,mo_num,key_max)
|
|
if (is_periodic) then
|
|
sze = key_max*2
|
|
call map_init(mo_integrals_map,sze)
|
|
call map_init(mo_integrals_map_2,sze)
|
|
print*, 'MO maps initialized (complex): ', 2*sze
|
|
else
|
|
sze = key_max
|
|
call map_init(mo_integrals_map,sze)
|
|
call map_init(mo_integrals_map_2,1_map_size_kind)
|
|
print*, 'MO map initialized: ', sze
|
|
endif
|
|
END_PROVIDER
|
|
|
|
subroutine insert_into_mo_integrals_map(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, 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 ]
|
|
&BEGIN_PROVIDER [ integer*8, mo_integrals_cache_max_8 ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! Min and max values of the MOs for which the integrals are in the cache
|
|
END_DOC
|
|
mo_integrals_cache_min_8 = max(1_8,elec_alpha_num - 63_8)
|
|
mo_integrals_cache_max_8 = min(int(mo_num,8),mo_integrals_cache_min_8+127_8)
|
|
mo_integrals_cache_min = max(1,elec_alpha_num - 63)
|
|
mo_integrals_cache_max = min(mo_num,mo_integrals_cache_min+127)
|
|
|
|
END_PROVIDER
|
|
|
|
BEGIN_PROVIDER [ double precision, mo_integrals_cache, (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
|
|
real(integral_kind) :: integral
|
|
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
|
|
call two_e_integrals_index(i4,j4,k4,l4,idx)
|
|
!DIR$ FORCEINLINE
|
|
call map_get(mo_integrals_map,idx,integral)
|
|
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(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
|
|
BEGIN_DOC
|
|
! Returns one integral <ij|kl> in the MO basis
|
|
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
|
|
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 = 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
|
|
! Returns one integral <ij|kl> in the MO basis
|
|
END_DOC
|
|
integer, intent(in) :: i,j,k,l
|
|
double precision :: get_two_e_integral
|
|
PROVIDE mo_two_e_integrals_in_map mo_integrals_cache
|
|
PROVIDE mo_two_e_integrals_in_map
|
|
!DIR$ FORCEINLINE
|
|
mo_two_e_integral = get_two_e_integral(i,j,k,l,mo_integrals_map)
|
|
return
|
|
end
|
|
|
|
subroutine get_mo_two_e_integrals(j,k,l,sze,out_val,map)
|
|
use map_module
|
|
implicit none
|
|
BEGIN_DOC
|
|
! Returns multiple integrals <ij|kl> in the MO basis, all
|
|
! i for j,k,l fixed.
|
|
END_DOC
|
|
integer, intent(in) :: j,k,l, sze
|
|
double precision, intent(out) :: out_val(sze)
|
|
type(map_type), intent(inout) :: map
|
|
integer :: i
|
|
double precision, external :: get_two_e_integral
|
|
|
|
integer :: ii, ii0
|
|
integer*8 :: ii_8, ii0_8
|
|
real(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
|
|
|
|
!DEBUG
|
|
! do i=1,sze
|
|
! out_val(i) = get_two_e_integral(i,j,k,l,map)
|
|
! 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)
|
|
|
|
q = min(j,l)
|
|
s = max(j,l)
|
|
q = q+shiftr(s*s-s,1)
|
|
|
|
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(ii_8)
|
|
else
|
|
p = min(i,k)
|
|
r = max(i,k)
|
|
p = p+shiftr(r*r-r,1)
|
|
i1 = min(p,q)
|
|
i2 = max(p,q)
|
|
idx = i1+shiftr(i2*i2-i2,1)
|
|
!DIR$ FORCEINLINE
|
|
call map_get(map,idx,tmp)
|
|
out_val(i) = dble(tmp)
|
|
endif
|
|
enddo
|
|
end
|
|
|
|
subroutine get_mo_two_e_integrals_ij(k,l,sze,out_array,map)
|
|
use map_module
|
|
implicit none
|
|
BEGIN_DOC
|
|
! Returns multiple integrals <ij|kl> 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(k,l,sze,out_array,map)
|
|
use map_module
|
|
implicit none
|
|
BEGIN_DOC
|
|
! Returns multiple integrals <ik|jl> 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(k,l,sze,out_val,map)
|
|
use map_module
|
|
implicit none
|
|
BEGIN_DOC
|
|
! Returns multiple integrals <ki|li>
|
|
! k(1)i(2) 1/r12 l(1)i(2) :: out_val(i1)
|
|
! for k,l fixed.
|
|
END_DOC
|
|
integer, intent(in) :: k,l, sze
|
|
double precision, intent(out) :: out_val(sze)
|
|
type(map_type), intent(inout) :: map
|
|
integer :: i
|
|
integer(key_kind) :: hash(sze)
|
|
real(integral_kind) :: tmp_val(sze)
|
|
PROVIDE mo_two_e_integrals_in_map
|
|
|
|
integer :: kk
|
|
do i=1,sze
|
|
!DIR$ FORCEINLINE
|
|
call two_e_integrals_index(k,i,l,i,hash(i))
|
|
enddo
|
|
|
|
if (integral_kind == 8) then
|
|
call map_get_many(map, hash, out_val, sze)
|
|
else
|
|
call map_get_many(map, hash, tmp_val, sze)
|
|
! Conversion to double precision
|
|
do i=1,sze
|
|
out_val(i) = dble(tmp_val(i))
|
|
enddo
|
|
endif
|
|
end
|
|
|
|
subroutine get_mo_two_e_integrals_exch_ii(k,l,sze,out_val,map)
|
|
use map_module
|
|
implicit none
|
|
BEGIN_DOC
|
|
! Returns multiple integrals <ki|il>
|
|
! k(1)i(2) 1/r12 i(1)l(2) :: out_val(i1)
|
|
! for k,l fixed.
|
|
END_DOC
|
|
integer, intent(in) :: k,l, sze
|
|
double precision, intent(out) :: out_val(sze)
|
|
type(map_type), intent(inout) :: map
|
|
integer :: i
|
|
integer(key_kind) :: hash(sze)
|
|
real(integral_kind) :: tmp_val(sze)
|
|
PROVIDE mo_two_e_integrals_in_map
|
|
|
|
integer :: kk
|
|
do i=1,sze
|
|
!DIR$ FORCEINLINE
|
|
call two_e_integrals_index(k,i,i,l,hash(i))
|
|
enddo
|
|
|
|
if (integral_kind == 8) then
|
|
call map_get_many(map, hash, out_val, sze)
|
|
else
|
|
call map_get_many(map, hash, tmp_val, sze)
|
|
! Conversion to double precision
|
|
do i=1,sze
|
|
out_val(i) = dble(tmp_val(i))
|
|
enddo
|
|
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
|
|
|
|
|
|
subroutine dump_mo_integrals(filename)
|
|
use map_module
|
|
implicit none
|
|
BEGIN_DOC
|
|
! Save to disk the |MO| integrals
|
|
END_DOC
|
|
character*(*), intent(in) :: filename
|
|
integer(cache_key_kind), pointer :: key(:)
|
|
real(integral_kind), pointer :: val(:)
|
|
integer*8 :: i,j, n
|
|
if (.not.mpi_master) then
|
|
return
|
|
endif
|
|
call ezfio_set_work_empty(.False.)
|
|
open(unit=66,file=filename,FORM='unformatted')
|
|
write(66) integral_kind, key_kind
|
|
write(66) mo_integrals_map%sorted, mo_integrals_map%map_size, &
|
|
mo_integrals_map%n_elements
|
|
do i=0_8,mo_integrals_map%map_size
|
|
write(66) mo_integrals_map%map(i)%sorted, mo_integrals_map%map(i)%map_size,&
|
|
mo_integrals_map%map(i)%n_elements
|
|
enddo
|
|
do i=0_8,mo_integrals_map%map_size
|
|
key => mo_integrals_map%map(i)%key
|
|
val => mo_integrals_map%map(i)%value
|
|
n = mo_integrals_map%map(i)%n_elements
|
|
write(66) (key(j), j=1,n), (val(j), j=1,n)
|
|
enddo
|
|
close(66)
|
|
|
|
end
|
|
|
|
|
|
integer function load_mo_integrals(filename)
|
|
implicit none
|
|
BEGIN_DOC
|
|
! Read from disk the |MO| integrals
|
|
END_DOC
|
|
character*(*), intent(in) :: filename
|
|
integer*8 :: i
|
|
integer(cache_key_kind), pointer :: key(:)
|
|
real(integral_kind), pointer :: val(:)
|
|
integer :: iknd, kknd
|
|
integer*8 :: n, j
|
|
load_mo_integrals = 1
|
|
open(unit=66,file=filename,FORM='unformatted',STATUS='UNKNOWN')
|
|
read(66,err=98,end=98) iknd, kknd
|
|
if (iknd /= integral_kind) then
|
|
print *, 'Wrong integrals kind in file :', iknd
|
|
stop 1
|
|
endif
|
|
if (kknd /= key_kind) then
|
|
print *, 'Wrong key kind in file :', kknd
|
|
stop 1
|
|
endif
|
|
read(66,err=98,end=98) mo_integrals_map%sorted, mo_integrals_map%map_size,&
|
|
mo_integrals_map%n_elements
|
|
do i=0_8, mo_integrals_map%map_size
|
|
read(66,err=99,end=99) mo_integrals_map%map(i)%sorted, &
|
|
mo_integrals_map%map(i)%map_size, mo_integrals_map%map(i)%n_elements
|
|
call cache_map_reallocate(mo_integrals_map%map(i),mo_integrals_map%map(i)%map_size)
|
|
enddo
|
|
do i=0_8, mo_integrals_map%map_size
|
|
key => mo_integrals_map%map(i)%key
|
|
val => mo_integrals_map%map(i)%value
|
|
n = mo_integrals_map%map(i)%n_elements
|
|
read(66,err=99,end=99) (key(j), j=1,n), (val(j), j=1,n)
|
|
enddo
|
|
call map_sort(mo_integrals_map)
|
|
load_mo_integrals = 0
|
|
return
|
|
99 continue
|
|
call map_deinit(mo_integrals_map)
|
|
98 continue
|
|
stop 'Problem reading mo_integrals_map file in work/'
|
|
|
|
end
|
|
|