10
0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 20:34:58 +01:00

Cleaned mo_get for multiple integrals

This commit is contained in:
Anthony Scemama 2020-05-28 11:57:12 +02:00
parent b6ebd8fd6d
commit 4ccb17a5dd
2 changed files with 9 additions and 103 deletions

View File

@ -127,7 +127,6 @@ double precision function mo_two_e_integral(i,j,k,l)
integer, intent(in) :: i,j,k,l integer, intent(in) :: i,j,k,l
double precision :: get_two_e_integral double precision :: get_two_e_integral
PROVIDE mo_two_e_integrals_in_map mo_integrals_cache PROVIDE mo_two_e_integrals_in_map mo_integrals_cache
PROVIDE mo_two_e_integrals_in_map
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
mo_two_e_integral = get_two_e_integral(i,j,k,l,mo_integrals_map) mo_two_e_integral = get_two_e_integral(i,j,k,l,mo_integrals_map)
return return
@ -202,47 +201,12 @@ subroutine get_mo_two_e_integrals_ij(k,l,sze,out_array,map)
integer, intent(in) :: k,l, sze integer, intent(in) :: k,l, sze
double precision, intent(out) :: out_array(sze,sze) double precision, intent(out) :: out_array(sze,sze)
type(map_type), intent(inout) :: map type(map_type), intent(inout) :: map
integer :: i,j,kk,ll,m integer :: j
integer(key_kind),allocatable :: hash(:)
integer ,allocatable :: pairs(:,:), iorder(:)
real(integral_kind), allocatable :: tmp_val(:) 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 j=1,sze
do i=1,sze call get_mo_two_e_integrals(j,k,l,sze,out_array(1,j),map)
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
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 end
subroutine get_mo_two_e_integrals_i1j1(k,l,sze,out_array,map) subroutine get_mo_two_e_integrals_i1j1(k,l,sze,out_array,map)
@ -256,47 +220,13 @@ subroutine get_mo_two_e_integrals_i1j1(k,l,sze,out_array,map)
integer, intent(in) :: k,l, sze integer, intent(in) :: k,l, sze
double precision, intent(out) :: out_array(sze,sze) double precision, intent(out) :: out_array(sze,sze)
type(map_type), intent(inout) :: map type(map_type), intent(inout) :: map
integer :: i,j,kk,ll,m integer :: j
integer(key_kind),allocatable :: hash(:)
integer ,allocatable :: pairs(:,:), iorder(:)
real(integral_kind), allocatable :: tmp_val(:)
PROVIDE mo_two_e_integrals_in_map 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 j=1,sze
do i=1,sze call get_mo_two_e_integrals(k,j,l,sze,out_array(1,j),map)
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 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 end
@ -312,25 +242,13 @@ subroutine get_mo_two_e_integrals_coulomb_ii(k,l,sze,out_val,map)
double precision, intent(out) :: out_val(sze) double precision, intent(out) :: out_val(sze)
type(map_type), intent(inout) :: map type(map_type), intent(inout) :: map
integer :: i integer :: i
integer(key_kind) :: hash(sze) double precision, external :: get_two_e_integral
real(integral_kind) :: tmp_val(sze)
PROVIDE mo_two_e_integrals_in_map PROVIDE mo_two_e_integrals_in_map
integer :: kk
do i=1,sze do i=1,sze
!DIR$ FORCEINLINE out_val(i) = get_two_e_integral(k,i,l,i,map)
call two_e_integrals_index(k,i,l,i,hash(i))
enddo 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 end
subroutine get_mo_two_e_integrals_exch_ii(k,l,sze,out_val,map) subroutine get_mo_two_e_integrals_exch_ii(k,l,sze,out_val,map)
@ -345,25 +263,13 @@ subroutine get_mo_two_e_integrals_exch_ii(k,l,sze,out_val,map)
double precision, intent(out) :: out_val(sze) double precision, intent(out) :: out_val(sze)
type(map_type), intent(inout) :: map type(map_type), intent(inout) :: map
integer :: i integer :: i
integer(key_kind) :: hash(sze) double precision, external :: get_two_e_integral
real(integral_kind) :: tmp_val(sze)
PROVIDE mo_two_e_integrals_in_map PROVIDE mo_two_e_integrals_in_map
integer :: kk
do i=1,sze do i=1,sze
!DIR$ FORCEINLINE out_val(i) = get_two_e_integral(k,i,i,l,map)
call two_e_integrals_index(k,i,i,l,hash(i))
enddo 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 end

View File

@ -17,7 +17,7 @@ program molden
write(i_unit_output,'(A)') '[Molden Format]' write(i_unit_output,'(A)') '[Molden Format]'
write(i_unit_output,'(A)') '[Atoms] AU' write(i_unit_output,'(A)') '[Atoms] ANGSTROM'
do i = 1, nucl_num do i = 1, nucl_num
write(i_unit_output,'(A2,2X,I4,2X,I4,3(2X,F15.10))') & write(i_unit_output,'(A2,2X,I4,2X,I4,3(2X,F15.10))') &
trim(element_name(int(nucl_charge(i)))), & trim(element_name(int(nucl_charge(i)))), &