LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-11-12 17:13:54 +01:00
Code Issues Releases Wiki Activity

Added Bi-electronic integrals module

Browse Source
This commit is contained in:
Anthony Scemama 2014-04-17 23:50:51 +02:00
parent b4904d65df
commit cd418ee484
28 changed files with 4793 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

0
src/BiInts/ASSUMPTIONS.rst Normal file
View File

8
src/BiInts/Makefile Normal file
View File

@ -0,0 +1,8 @@
default: all
# Define here all new external source files and objects.Don't forget to prefix the
# object files with IRPF90_temp/
SRC=
OBJ=
include $(QPACKAGE_ROOT)/src/Makefile.common

1
src/BiInts/NEEDED_MODULES Normal file
View File

@ -0,0 +1 @@
AOs Bitmask Electrons Ezfio_files MOs Nuclei Output Utils MonoInts

27
src/BiInts/README.rst Normal file
View File

@ -0,0 +1,27 @@
=============
BiInts Module
=============
Here, all bi-electronic integrals (:math:`1/r_{12}`) are computed. As they have
4 indices and many are zero, they are stored in a map.
To fetch an AO integral, use the ``get_ao_bielec_integral(i,j,k,l,ao_integrals_map)``
function, and to fetch and MO integral, use
.
Needed Modules
==============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
* `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
* `Bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask>`_
* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
* `MOs <http://github.com/LCPQ/quantum_package/tree/master/src/MOs>`_
* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_

1016
src/BiInts/ao_bi_integrals.irp.f Normal file
View File

File diff suppressed because it is too large Load Diff

8
src/BiInts/bi_integrals.ezfio_config Normal file
View File

@ -0,0 +1,8 @@
bielec_integrals
read_ao_integrals logical
read_mo_integrals logical
write_ao_integrals logical
write_mo_integrals logical
threshold_ao double precision
threshold_mo double precision

57
src/BiInts/gauss_legendre.irp.f Normal file
View File

@ -0,0 +1,57 @@
BEGIN_PROVIDER [ double precision, gauleg_t2, (n_pt_max_integrals,n_pt_max_integrals/2) ]
&BEGIN_PROVIDER [ double precision, gauleg_w, (n_pt_max_integrals,n_pt_max_integrals/2) ]
implicit none
BEGIN_DOC
! t_w(i,1,k) = w(i)
! t_w(i,2,k) = t(i)
END_DOC
integer :: i,j,l
l=0
do i = 2,n_pt_max_integrals,2
l = l+1
call gauleg(0.d0,1.d0,gauleg_t2(1,l),gauleg_w(1,l),i)
do j=1,i
gauleg_t2(j,l) *= gauleg_t2(j,l)
enddo
enddo
END_PROVIDER
subroutine gauleg(x1,x2,x,w,n)
implicit none
BEGIN_DOC
! Gauss-Legendre
END_DOC
integer, intent(in) :: n
double precision, intent(in) :: x1, x2
double precision, intent (out) :: x(n),w(n)
double precision, parameter :: eps=3.d-14
integer :: m,i,j
double precision :: xm, xl, z, z1, p1, p2, p3, pp, dn
m=(n+1)/2
xm=0.5d0*(x2+x1)
xl=0.5d0*(x2-x1)
dn = dble(n)
do i=1,m
z=dcos(3.141592654d0*(dble(i)-.25d0)/(dble(n)+.5d0))
z1 = z+1.d0
do while (dabs(z-z1) > eps)
p1=1.d0
p2=0.d0
do j=1,n
p3=p2
p2=p1
p1=(dble(j+j-1)*z*p2-dble(j-1)*p3)/j
enddo
pp=dn*(z*p1-p2)/(z*z-1.d0)
z1=z
z=z1-p1/pp
end do
x(i)=xm-xl*z
x(n+1-i)=xm+xl*z
w(i)=(xl+xl)/((1.d0-z*z)*pp*pp)
w(n+1-i)=w(i)
enddo
end

334
src/BiInts/map_integrals.irp.f Normal file
View File

@ -0,0 +1,334 @@
use map_module
!! AO Map
!! ======
BEGIN_PROVIDER [ type(map_type), ao_integrals_map ]
implicit none
BEGIN_DOC
! AO integrals
END_DOC
integer*8 :: sze
call bielec_integrals_index(ao_num,ao_num,ao_num,ao_num,sze)
call map_init(ao_integrals_map,sze)
write(output_BiInts,*) 'AO map initialized'
END_PROVIDER
subroutine bielec_integrals_index(i,j,k,l,i1)
implicit none
integer, intent(in) :: i,j,k,l
integer*8, intent(out) :: i1
integer*8 :: p,q,r,s,i2
p = min(i,k)
r = max(i,k)
p = p+ishft(r*r-r,-1)
q = min(j,l)
s = max(j,l)
q = q+ishft(s*s-s,-1)
i1 = min(p,q)
i2 = max(p,q)
i1 = i1+ishft(i2*i2-i2,-1)
end
double precision function get_ao_bielec_integral(i,j,k,l,map)
use map_module
implicit none
BEGIN_DOC
! Gets one AO bi-electronic integral from the AO map
END_DOC
integer, intent(in) :: i,j,k,l
integer*8 :: idx
type(map_type), intent(inout) :: map
real(integral_kind) :: tmp
PROVIDE ao_bielec_integrals_in_map
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,idx)
call map_get(map,idx,tmp)
get_ao_bielec_integral = tmp
end
subroutine get_ao_bielec_integrals(j,k,l,sze,out_val)
use map_module
BEGIN_DOC
! Gets multiple AO bi-electronic integral from the AO map .
! All i are retrieved for j,k,l fixed.
END_DOC
implicit none
integer, intent(in) :: j,k,l, sze
real(integral_kind), intent(out) :: out_val(sze)
integer :: i
integer*8 :: hash
double precision :: thresh
PROVIDE ao_bielec_integrals_in_map
thresh = ao_integrals_threshold
if (ao_overlap_abs(j,l) < thresh) then
out_val = 0.d0
return
endif
do i=1,sze
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh ) then
out_val(i) = 0.d0
else
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,hash)
call map_get(ao_integrals_map, hash, out_val(i))
endif
enddo
end
subroutine get_ao_bielec_integrals_non_zero(j,k,l,sze,out_val,out_val_index,non_zero_int)
use map_module
implicit none
BEGIN_DOC
! Gets multiple AO bi-electronic integral from the AO map .
! All non-zero i are retrieved for j,k,l fixed.
END_DOC
integer, intent(in) :: j,k,l, sze
real(integral_kind), intent(out) :: out_val(sze)
integer, intent(out) :: out_val_index(sze),non_zero_int
integer :: i
integer*8 :: hash
double precision :: thresh,tmp
PROVIDE ao_bielec_integrals_in_map
thresh = ao_integrals_threshold
if (ao_overlap_abs(j,l) < thresh) then
out_val = 0.d0
return
endif
non_zero_int = 0
do i=1,sze
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,hash)
call map_get(ao_integrals_map, hash,tmp)
if (dabs(tmp) < thresh ) cycle
non_zero_int = non_zero_int+1
out_val_index(non_zero_int) = i
out_val(non_zero_int) = tmp
enddo
end
integer*8 function get_ao_map_size()
implicit none
BEGIN_DOC
! Returns the number of elements in the AO map
END_DOC
get_ao_map_size = ao_integrals_map % n_elements
end
subroutine clear_ao_map
implicit none
BEGIN_DOC
! Frees the memory of the AO map
END_DOC
call map_deinit(ao_integrals_map)
FREE ao_integrals_map
end
!! MO Map
!! ======
BEGIN_PROVIDER [ type(map_type), mo_integrals_map ]
implicit none
BEGIN_DOC
! MO integrals
END_DOC
integer*8 :: sze
call bielec_integrals_index(mo_tot_num,mo_tot_num,mo_tot_num,mo_tot_num,sze)
call map_init(mo_integrals_map,sze)
write(output_BiInts,*) 'MO map initialized'
END_PROVIDER
subroutine insert_into_ao_integrals_map(n_integrals, &
buffer_i, buffer_values)
use map_module
implicit none
BEGIN_DOC
! Create new entry into AO map
END_DOC
integer, intent(in) :: n_integrals
integer*8, intent(inout) :: buffer_i(n_integrals)
real(integral_kind), intent(inout) :: buffer_values(n_integrals)
call map_append(ao_integrals_map, buffer_i, buffer_values, n_integrals)
end
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*8, 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
double precision function get_mo_bielec_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*8 :: idx
type(map_type), intent(inout) :: map
real(integral_kind) :: tmp
PROVIDE mo_bielec_integrals_in_map
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,idx)
call map_get(map,idx,tmp)
get_mo_bielec_integral = dble(tmp)
end
double precision function mo_bielec_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_mo_bielec_integral
PROVIDE mo_bielec_integrals_in_map
mo_bielec_integral = get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
return
end
subroutine get_mo_bielec_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
real(integral_kind), intent(out) :: out_val(sze)
type(map_type), intent(inout) :: map
integer :: i
integer*8 :: hash(sze)
PROVIDE mo_bielec_integrals_in_map
do i=1,sze
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,hash(i))
enddo
call map_get_many(map, hash, out_val, sze)
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
end
subroutine clear_mo_map
implicit none
BEGIN_DOC
! Frees the memory of the MO map
END_DOC
call map_deinit(mo_integrals_map)
FREE mo_integrals_map
end
BEGIN_TEMPLATE
subroutine dump_$ao_integrals(filename)
use map_module
implicit none
BEGIN_DOC
! Save to disk the $ao integrals
END_DOC
character*(*), intent(in) :: filename
integer(cache_key_kind), pointer :: key(:)
real(integral_kind), pointer :: val(:)
integer*8 :: i,j, n
open(unit=66,file=filename,FORM='unformatted')
write(66) integral_kind, key_kind
write(66) $ao_integrals_map%sorted, $ao_integrals_map%map_size, &
$ao_integrals_map%n_elements
do i=0_8,$ao_integrals_map%map_size
write(66) $ao_integrals_map%map(i)%sorted, $ao_integrals_map%map(i)%map_size,&
$ao_integrals_map%map(i)%n_elements
enddo
do i=0_8,$ao_integrals_map%map_size
key => $ao_integrals_map%map(i)%key
val => $ao_integrals_map%map(i)%value
n = $ao_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_$ao_integrals(filename)
implicit none
BEGIN_DOC
! Read from disk the $ao integrals
END_DOC
character*(*), intent(in) :: filename
integer*8 :: i
integer(cache_key_kind), pointer :: key(:)
real(integral_kind), pointer :: val(:)
integer :: iknd, kknd, n, j
load_$ao_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) $ao_integrals_map%sorted, $ao_integrals_map%map_size,&
$ao_integrals_map%n_elements
do i=0_8, $ao_integrals_map%map_size
read(66,err=99,end=99) $ao_integrals_map%map(i)%sorted, &
$ao_integrals_map%map(i)%map_size, $ao_integrals_map%map(i)%n_elements
call cache_map_reallocate($ao_integrals_map%map(i),$ao_integrals_map%map(i)%map_size)
enddo
do i=0_8, $ao_integrals_map%map_size
key => $ao_integrals_map%map(i)%key
val => $ao_integrals_map%map(i)%value
n = $ao_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($ao_integrals_map)
load_$ao_integrals = 0
return
99 continue
call map_deinit($ao_integrals_map)
FREE $ao_integrals_map
if (.True.) then
PROVIDE $ao_integrals_map
endif
stop 'Problem reading $ao_integrals_map file in work/'
98 continue
end
SUBST [ ao_integrals_map, ao_integrals, ao_num , get_ao_bielec_integral ]
ao_integrals_map ; ao_integrals ; ao_num ; get_ao_bielec_integral ;;
mo_integrals_map ; mo_integrals ; n_act ; get_mo_bielec_integral ;;
END_TEMPLATE

356
src/BiInts/mo_bi_integrals.irp.f Normal file
View File

@ -0,0 +1,356 @@
subroutine mo_bielec_integrals_index(i,j,k,l,i1)
implicit none
BEGIN_DOC
! Computes an unique index for i,j,k,l integrals
END_DOC
integer, intent(in) :: i,j,k,l
integer*8, intent(out) :: i1
integer*8 :: p,q,r,s,i2
p = min(i,k)
r = max(i,k)
p = p+ishft(r*r-r,-1)
q = min(j,l)
s = max(j,l)
q = q+ishft(s*s-s,-1)
i1 = min(p,q)
i2 = max(p,q)
i1 = i1+ishft(i2*i2-i2,-1)
end
BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ]
implicit none
BEGIN_DOC
! If True, the map of MO bielectronic integrals is provided
END_DOC
mo_bielec_integrals_in_map = .True.
if (read_mo_integrals) then
integer :: load_mo_integrals
if (load_mo_integrals(trim(ezfio_filename)//'/work/mo_integrals.bin') == 0) then
write(output_BiInts,*) 'MO integrals provided'
return
endif
endif
call add_integrals_to_map(full_ijkl_bitmask)
END_PROVIDER
subroutine add_integrals_to_map(mask_ijkl)
use bitmasks
implicit none
BEGIN_DOC
! Adds integrals to tha MO map according to some bitmask
END_DOC
integer(bit_kind), intent(in) :: mask_ijkl(N_int,4)
integer :: i,j,k,l
integer :: i0,j0,k0,l0
double precision :: c, cpu_1, cpu_2, wall_1, wall_2
integer, allocatable :: list_ijkl(:,:)
integer :: n_i, n_j, n_k, n_l
integer, allocatable :: bielec_tmp_0_idx(:)
real(integral_kind), allocatable :: bielec_tmp_0(:,:)
double precision, allocatable :: bielec_tmp_1(:)
double precision, allocatable :: bielec_tmp_2(:,:)
double precision, allocatable :: bielec_tmp_3(:,:,:)
!DEC$ ATTRIBUTES ALIGN : 64 :: bielec_tmp_1, bielec_tmp_2, bielec_tmp_3
integer :: n_integrals
integer :: size_buffer
integer*8,allocatable :: buffer_i(:)
real(integral_kind),allocatable :: buffer_value(:)
real :: map_mb
integer :: i1,j1,k1,l1, ii1, kmax, l1_global
integer :: i2,i3,i4
double precision,parameter :: thr_coef = 0.d0
PROVIDE N_int ao_bielec_integrals_in_map ao_integrals_map mo_coef mo_coef_transp
!Get list of MOs for i,j,k and l
!-------------------------------
allocate(list_ijkl(mo_tot_num,4))
call bitstring_to_list( mask_ijkl(1,1), list_ijkl(1,1), n_i, N_int )
call bitstring_to_list( mask_ijkl(1,2), list_ijkl(1,2), n_j, N_int )
call bitstring_to_list( mask_ijkl(1,3), list_ijkl(1,3), n_k, N_int )
call bitstring_to_list( mask_ijkl(1,4), list_ijkl(1,4), n_l, N_int )
l1_global = 0
size_buffer = min(ao_num*ao_num*ao_num,16000000)
write(output_BiInts,*) 'Providing the molecular integrals '
write(output_BiInts,*) 'Buffers : ', 8.*(mo_tot_num_align*(n_j)*(n_k+1) + mo_tot_num_align +&
ao_num+ao_num*ao_num+ size_buffer*3)/(1024*1024), 'MB / core'
call wall_time(wall_1)
call cpu_time(cpu_1)
mo_integrals_threshold = 0.d0
!$OMP PARALLEL PRIVATE(l1,k1,j1,i1,i2,i3,i4,i,j,k,l,c, ii1,kmax, &
!$OMP bielec_tmp_0_idx, bielec_tmp_0, bielec_tmp_1,bielec_tmp_2,bielec_tmp_3,&
!$OMP buffer_i,buffer_value,n_integrals,wall_2,i0,j0,k0,l0) &
!$OMP DEFAULT(NONE) &
!$OMP SHARED(size_buffer,ao_num,mo_tot_num,n_i,n_j,n_k,n_l,mo_tot_num_align,&
!$OMP mo_coef_transp,output_BiInts, &
!$OMP mo_coef_transp_is_built, list_ijkl, &
!$OMP mo_coef_is_built, wall_1, &
!$OMP mo_coef,mo_integrals_threshold,l1_global,ao_integrals_map,mo_integrals_map)
n_integrals = 0
allocate(bielec_tmp_3(mo_tot_num_align, n_j, n_k), &
bielec_tmp_1(mo_tot_num_align), &
bielec_tmp_0(ao_num,ao_num), &
bielec_tmp_0_idx(ao_num), &
bielec_tmp_2(mo_tot_num_align, n_j), &
buffer_i(size_buffer), &
buffer_value(size_buffer) )
!$OMP DO SCHEDULE(guided)
do l1 = 1,ao_num
!DEC$ VECTOR ALIGNED
bielec_tmp_3 = 0.d0
do k1 = 1,ao_num
!DEC$ VECTOR ALIGNED
bielec_tmp_2 = 0.d0
do j1 = 1,ao_num
call get_ao_bielec_integrals(j1,k1,l1,ao_num,bielec_tmp_0(1,j1))
! call compute_ao_bielec_integrals(j1,k1,l1,ao_num,bielec_tmp_0(1,j1))
enddo
do j1 = 1,ao_num
kmax = 0
do i1 = 1,ao_num
c = bielec_tmp_0(i1,j1)
if (c == 0.d0) then
cycle
endif
kmax += 1
bielec_tmp_0(kmax,j1) = c
bielec_tmp_0_idx(kmax) = i1
enddo
if (kmax==0) then
cycle
endif
!DEC$ VECTOR ALIGNED
bielec_tmp_1 = 0.d0
ii1=1
do ii1 = 1,kmax-4,4
i1 = bielec_tmp_0_idx(ii1)
i2 = bielec_tmp_0_idx(ii1+1)
i3 = bielec_tmp_0_idx(ii1+2)
i4 = bielec_tmp_0_idx(ii1+3)
do i = list_ijkl(1,1), list_ijkl(n_i,1)
bielec_tmp_1(i) = bielec_tmp_1(i) + &
mo_coef_transp(i,i1) * bielec_tmp_0(ii1,j1) + &
mo_coef_transp(i,i2) * bielec_tmp_0(ii1+1,j1) + &
mo_coef_transp(i,i3) * bielec_tmp_0(ii1+2,j1) + &
mo_coef_transp(i,i4) * bielec_tmp_0(ii1+3,j1)
enddo ! i
enddo ! ii1
i2 = ii1
do ii1 = i2,kmax
i1 = bielec_tmp_0_idx(ii1)
do i = list_ijkl(1,1), list_ijkl(n_i,1)
bielec_tmp_1(i) = bielec_tmp_1(i) + mo_coef_transp(i,i1) * bielec_tmp_0(ii1,j1)
enddo ! i
enddo ! ii1
c = 0.d0
do i = list_ijkl(1,1), list_ijkl(n_i,1)
c = max(c,abs(bielec_tmp_1(i)))
if (c>mo_integrals_threshold) exit
enddo
if ( c < mo_integrals_threshold ) then
cycle
endif
do j0 = 1, n_j
j = list_ijkl(j0,2)
c = mo_coef_transp(j,j1)
if (abs(c) < thr_coef) then
cycle
endif
do i = list_ijkl(1,1), list_ijkl(n_i,1)
bielec_tmp_2(i,j0) = bielec_tmp_2(i,j0) + c * bielec_tmp_1(i)
enddo ! i
enddo ! j
enddo !j1
if ( maxval(abs(bielec_tmp_2)) < mo_integrals_threshold ) then
cycle
endif
do k0 = 1, n_k
k = list_ijkl(k0,3)
c = mo_coef_transp(k,k1)
if (abs(c) < thr_coef) then
cycle
endif
do j0 = 1, n_j
j = list_ijkl(j0,2)
do i = list_ijkl(1,1), k
bielec_tmp_3(i,j0,k0) = bielec_tmp_3(i,j0,k0) + c* bielec_tmp_2(i,j0)
enddo!i
enddo !j
enddo !k
enddo !k1
do l0 = 1,n_l
l = list_ijkl(l0,4)
c = mo_coef_transp(l,l1)
if (abs(c) < thr_coef) then
cycle
endif
j1 = ishft((l*l-l),-1)
do j0 = 1, n_j
j = list_ijkl(j0,2)
if (j > l) then
exit
endif
j1 += 1
do k0 = 1, n_k
k = list_ijkl(k0,3)
i1 = ishft((k*k-k),-1)
if (i1<=j1) then
continue
else
exit
endif
bielec_tmp_1 = 0.d0
do i0 = 1, n_i
i = list_ijkl(i0,1)
if (i>k) then
exit
endif
bielec_tmp_1(i) = c*bielec_tmp_3(i,j0,k0)
enddo
do i = 1, min(k,j1-i1)
if (abs(bielec_tmp_1(i)) < mo_integrals_threshold) then
cycle
endif
n_integrals += 1
buffer_value(n_integrals) = bielec_tmp_1(i)
!DEC$ FORCEINLINE
call mo_bielec_integrals_index(i,j,k,l,buffer_i(n_integrals))
if (n_integrals == size_buffer) then
call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,&
real(mo_integrals_threshold,integral_kind))
n_integrals = 0
endif
enddo
enddo
enddo
enddo
!$OMP ATOMIC
l1_global +=1
call wall_time(wall_2)
write(output_BiInts,*) 100.*float(l1_global)/float(ao_num), '% in ',&
wall_2-wall_1, 's', map_mb(mo_integrals_map) ,'MB'
enddo
!$OMP END DO NOWAIT
deallocate (bielec_tmp_1,bielec_tmp_2,bielec_tmp_3)
call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,&
real(mo_integrals_threshold,integral_kind))
deallocate(buffer_i, buffer_value)
!$OMP END PARALLEL
call map_unique(mo_integrals_map)
call wall_time(wall_2)
call cpu_time(cpu_2)
integer*8 :: get_mo_map_size, mo_map_size
mo_map_size = get_mo_map_size()
deallocate(list_ijkl)
write(output_BiInts,*)'Molecular integrals provided:'
write(output_BiInts,*)' Size of MO map ', map_mb(mo_integrals_map) ,'MB'
write(output_BiInts,*)' Number of MO integrals: ', mo_map_size
write(output_BiInts,*)' cpu time :',cpu_2 - cpu_1, 's'
write(output_BiInts,*)' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')'
if (write_mo_integrals) then
call dump_mo_integrals(trim(ezfio_filename)//'/work/mo_integrals.bin')
endif
end
BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj, (mo_tot_num_align,mo_tot_num)]
implicit none
BEGIN_DOC
! Bi-electronic integrals <ij|ij>
END_DOC
integer :: i,j
double precision :: get_mo_bielec_integral
PROVIDE mo_bielec_integrals_in_map
do j= 1, mo_tot_num
!DIR$ VECTOR ALIGNED
do i= 1, mo_tot_num
mo_bielec_integral_jj(i,j) = get_mo_bielec_integral(i,j,i,j,mo_integrals_map)
enddo
! Padding
do i= mo_tot_num+1, mo_tot_num_align
mo_bielec_integral_jj(i,j) = 0.d0
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_exchange, (mo_tot_num_align,mo_tot_num)]
implicit none
BEGIN_DOC
! Bi-electronic integrals <ij|ij> - <ij|ji>
END_DOC
integer :: i,j
double precision :: get_mo_bielec_integral
PROVIDE mo_bielec_integrals_in_map
do j = 1, mo_tot_num
!DIR$ VECTOR ALIGNED
do i = 1,mo_tot_num
mo_bielec_integral_jj_exchange(i,j) = get_mo_bielec_integral(i,j,j,i,mo_integrals_map)
enddo
! Padding
do i= mo_tot_num+1, mo_tot_num_align
mo_bielec_integral_jj_exchange(i,j) = 0.d0
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_anti, (mo_tot_num_align,mo_tot_num)]
implicit none
BEGIN_DOC
! Bi-electronic integrals <ij|ij> - <ij|ji>
END_DOC
integer :: i,j
PROVIDE mo_bielec_integrals_in_map
do j = 1, mo_tot_num
!DIR$ VECTOR ALIGNED
do i = 1,mo_tot_num_align
mo_bielec_integral_jj_anti(i,j) = mo_bielec_integral_jj(i,j) - mo_bielec_integral_jj_exchange(i,j)
enddo
enddo
END_PROVIDER

109
src/BiInts/options.irp.f Normal file
View File

@ -0,0 +1,109 @@
BEGIN_PROVIDER [ logical, write_mo_integrals ]
implicit none
BEGIN_DOC
! If true, write MO integrals in EZFIO
END_DOC
logical :: has
PROVIDE ezfio_filename
call ezfio_has_bielec_integrals_write_mo_integrals(has)
if (has) then
call ezfio_get_bielec_integrals_write_mo_integrals(write_mo_integrals)
else
write_mo_integrals = .False.
call ezfio_set_bielec_integrals_write_mo_integrals(write_mo_integrals)
endif
END_PROVIDER
BEGIN_PROVIDER [ logical, write_ao_integrals ]
implicit none
BEGIN_DOC
! If true, write AO integrals in EZFIO
END_DOC
logical :: has
PROVIDE ezfio_filename
call ezfio_has_bielec_integrals_write_ao_integrals(has)
if (has) then
call ezfio_get_bielec_integrals_write_ao_integrals(write_ao_integrals)
else
write_ao_integrals = .False.
call ezfio_set_bielec_integrals_write_ao_integrals(write_ao_integrals)
endif
END_PROVIDER
BEGIN_PROVIDER [ logical, read_mo_integrals ]
implicit none
BEGIN_DOC
! If true, read MO integrals in EZFIO
END_DOC
logical :: has
PROVIDE ezfio_filename
call ezfio_has_bielec_integrals_read_mo_integrals(has)
if (has) then
call ezfio_get_bielec_integrals_read_mo_integrals(read_mo_integrals)
else
read_mo_integrals = .False.
call ezfio_set_bielec_integrals_read_mo_integrals(read_mo_integrals)
endif
END_PROVIDER
BEGIN_PROVIDER [ logical, read_ao_integrals ]
implicit none
BEGIN_DOC
! If true, read AO integrals in EZFIO
END_DOC
logical :: has
PROVIDE ezfio_filename
call ezfio_has_bielec_integrals_read_ao_integrals(has)
if (has) then
call ezfio_get_bielec_integrals_read_ao_integrals(read_ao_integrals)
else
read_ao_integrals = .False.
call ezfio_set_bielec_integrals_read_ao_integrals(read_ao_integrals)
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_integrals_threshold ]
implicit none
BEGIN_DOC
! If <pq|rs> < ao_integrals_threshold, <pq|rs> = 0
END_DOC
logical :: has
PROVIDE ezfio_filename
call ezfio_has_bielec_integrals_threshold_ao(has)
if (has) then
call ezfio_get_bielec_integrals_threshold_ao(ao_integrals_threshold)
else
ao_integrals_threshold = 1.d-15
call ezfio_set_bielec_integrals_threshold_ao(ao_integrals_threshold)
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_integrals_threshold ]
implicit none
BEGIN_DOC
! If <ij|kl> < mo_integrals_threshold, <ij|kl> = 0
END_DOC
logical :: has
PROVIDE ezfio_filename
call ezfio_has_bielec_integrals_threshold_mo(has)
if (has) then
call ezfio_get_bielec_integrals_threshold_mo(mo_integrals_threshold)
else
mo_integrals_threshold = 1.d-11
call ezfio_set_bielec_integrals_threshold_mo(mo_integrals_threshold)
endif
END_PROVIDER

33
src/BiInts/tests/Makefile Normal file
View File

@ -0,0 +1,33 @@
OPENMP =1
PROFILE =0
DEBUG = 0
IRPF90+= -I tests
REF_FILES=$(subst %.irp.f, %.ref, $(wildcard *.irp.f))
.PHONY: clean executables serial_tests parallel_tests
all: clean executables serial_tests parallel_tests
parallel_tests: $(REF_FILES)
@echo ; echo " ---- Running parallel tests ----" ; echo
@OMP_NUM_THREADS=10 ${QPACKAGE_ROOT}/scripts/run_tests.py
serial_tests: $(REF_FILES)
@echo ; echo " ---- Running serial tests ----" ; echo
@OMP_NUM_THREADS=1 ${QPACKAGE_ROOT}/scripts/run_tests.py
executables: $(wildcard *.irp.f) veryclean
$(MAKE) -C ..
%.ref: $(wildcard $(QPACKAGE_ROOT)/data/inputs/*.md5) executables
$(QPACKAGE_ROOT)/scripts/create_test_ref.sh $*
clean:
$(MAKE) -C .. clean
veryclean:
$(MAKE) -C .. veryclean

424
src/BiInts/tests/test_ao.ref Normal file
View File

@ -0,0 +1,424 @@
data = {
'AlCl.ezfio' : {
},
'Al.ezfio' : {
},
'Al+.ezfio' : {
},
'AlH2.ezfio' : {
},
'AlH3.ezfio' : {
},
'AlH.ezfio' : {
},
'BCl.ezfio' : {
},
'BeCl.ezfio' : {
},
'Be.ezfio' : {
},
'Be+.ezfio' : {
},
'BeF.ezfio' : {
},
'BeH2.ezfio' : {
},
'BeH.ezfio' : {
},
'BeO.ezfio' : {
},
'BeOH.ezfio' : {
},
'BeS.ezfio' : {
},
'B.ezfio' : {
},
'B+.ezfio' : {
},
'BH2.ezfio' : {
},
'BH3.ezfio' : {
},
'BH.ezfio' : {
},
'BO.ezfio' : {
},
'BS.ezfio' : {
},
'C2.ezfio' : {
},
'C2H2.ezfio' : {
},
'C2H2+.ezfio' : {
},
'C2H3.ezfio' : {
},
'C2H3+.ezfio' : {
},
'C2H4.ezfio' : {
},
'C2H4+.ezfio' : {
},
'C2H5.ezfio' : {
},
'C2H6.ezfio' : {
},
'C2H.ezfio' : {
},
'CCl.ezfio' : {
},
'C-.ezfio' : {
},
'C.ezfio' : {
},
'C+.ezfio' : {
},
'CF.ezfio' : {
},
'CH2_1A1.ezfio' : {
},
'CH2_3B1.ezfio' : {
},
'CH2-.ezfio' : {
},
'CH3-.ezfio' : {
},
'CH3.ezfio' : {
},
'CH4.ezfio' : {
},
'CH4+.ezfio' : {
},
'CH-.ezfio' : {
},
'CH.ezfio' : {
},
'Cl2-.ezfio' : {
},
'Cl2.ezfio' : {
},
'Cl2+.ezfio' : {
},
'Cl-.ezfio' : {
},
'Cl.ezfio' : {
},
'Cl+.ezfio' : {
},
'ClH2+.ezfio' : {
},
'ClH.ezfio' : {
},
'ClH+.ezfio' : {
},
'ClS.ezfio' : {
},
'ClSiH3.ezfio' : {
},
'CN-.ezfio' : {
},
'CN.ezfio' : {
},
'CO2.ezfio' : {
},
'CO.ezfio' : {
},
'CO+.ezfio' : {
},
'COS.ezfio' : {
},
'CP.ezfio' : {
},
'CS2.ezfio' : {
},
'CS.ezfio' : {
},
'CS+.ezfio' : {
},
'CSi.ezfio' : {
},
'F2.ezfio' : {
},
'FAl.ezfio' : {
},
'FCl.ezfio' : {
},
'FCl+.ezfio' : {
},
'F-.ezfio' : {
},
'F.ezfio' : {
},
'F+.ezfio' : {
},
'FH.ezfio' : {
},
'FH+.ezfio' : {
},
'FMg.ezfio' : {
},
'FNa.ezfio' : {
},
'FP.ezfio' : {
},
'FS.ezfio' : {
},
'FSi.ezfio' : {
},
'FSiH3.ezfio' : {
},
'H2CNH.ezfio' : {
},
'H2CO.ezfio' : {
},
'H2CPH.ezfio' : {
},
'H2CS.ezfio' : {
},
'H2.ezfio' : {
},
'H2NNH2.ezfio' : {
},
'H2PPH2.ezfio' : {
},
'H3CCl.ezfio' : {
},
'H3CF.ezfio' : {
},
'H3CNH2.ezfio' : {
},
'H3COH.ezfio' : {
},
'H3CSH.ezfio' : {
},
'H3SiSiH3.ezfio' : {
},
'HBO.ezfio' : {
},
'HBS.ezfio' : {
},
'HCF.ezfio' : {
},
'HCN.ezfio' : {
},
'HCO.ezfio' : {
},
'HCP.ezfio' : {
},
'H.ezfio' : {
},
'HNO.ezfio' : {
},
'HOCl.ezfio' : {
},
'HOF.ezfio' : {
},
'HOMg.ezfio' : {
},
'HONa.ezfio' : {
},
'HOO.ezfio' : {
},
'HOOH.ezfio' : {
},
'HSSH.ezfio' : {
},
'Li2.ezfio' : {
},
'LiCl.ezfio' : {
},
'Li.ezfio' : {
},
'Li+.ezfio' : {
},
'LiF.ezfio' : {
},
'LiH.ezfio' : {
},
'LiN.ezfio' : {
},
'LiO.ezfio' : {
},
'LiOH.ezfio' : {
},
'MgCl.ezfio' : {
},
'Mg.ezfio' : {
},
'Mg+.ezfio' : {
},
'MgH.ezfio' : {
},
'MgS.ezfio' : {
},
'N2.ezfio' : {
},
'N2+.ezfio' : {
},
'Na2.ezfio' : {
},
'NaCl.ezfio' : {
},
'Na.ezfio' : {
},
'Na+.ezfio' : {
},
'NaH.ezfio' : {
},
'N.ezfio' : {
},
'N+.ezfio' : {
},
'NF.ezfio' : {
},
'NH2-.ezfio' : {
},
'NH2.ezfio' : {
},
'NH3.ezfio' : {
},
'NH3+.ezfio' : {
},
'NH4+.ezfio' : {
},
'NH-.ezfio' : {
},
'NH.ezfio' : {
},
'NO-.ezfio' : {
},
'NO.ezfio' : {
},
'NP.ezfio' : {
},
'NS.ezfio' : {
},
'NSi.ezfio' : {
},
'O2Cl.ezfio' : {
},
'O2-.ezfio' : {
},
'O2.ezfio' : {
},
'O2+.ezfio' : {
},
'O2S.ezfio' : {
},
'O2Si.ezfio' : {
},
'O3.ezfio' : {
},
'OCl.ezfio' : {
},
'O-.ezfio' : {
},
'O.ezfio' : {
},
'O+.ezfio' : {
},
'OH2.ezfio' : {
},
'OH2+.ezfio' : {
},
'OH3+.ezfio' : {
},
'OH-.ezfio' : {
},
'OH.ezfio' : {
},
'OH+.ezfio' : {
},
'OMg.ezfio' : {
},
'ONa.ezfio' : {
},
'OP-.ezfio' : {
},
'OP.ezfio' : {
},
'OPH.ezfio' : {
},
'OS.ezfio' : {
},
'OSi.ezfio' : {
},
'P2.ezfio' : {
},
'P2+.ezfio' : {
},
'PCl.ezfio' : {
},
'P-.ezfio' : {
},
'P.ezfio' : {
},
'PH2-.ezfio' : {
},
'PH2.ezfio' : {
},
'PH2+.ezfio' : {
},
'PH3.ezfio' : {
},
'PH3+.ezfio' : {
},
'PH4+.ezfio' : {
},
'PH-.ezfio' : {
},
'PH.ezfio' : {
},
'PS.ezfio' : {
},
'S2-.ezfio' : {
},
'S2.ezfio' : {
},
'S-.ezfio' : {
},
'S.ezfio' : {
},
'S+.ezfio' : {
},
'SH2.ezfio' : {
},
'SH2+.ezfio' : {
},
'SH3+.ezfio' : {
},
'SH-.ezfio' : {
},
'SH.ezfio' : {
},
'SH+.ezfio' : {
},
'Si2.ezfio' : {
},
'SiCl.ezfio' : {
},
'Si-.ezfio' : {
},
'Si.ezfio' : {
},
'SiH2_1A1.ezfio' : {
},
'SiH2_3B1.ezfio' : {
},
'SiH2-.ezfio' : {
},
'SiH3-.ezfio' : {
},
'SiH3.ezfio' : {
},
'SiH4.ezfio' : {
},
'SiH4+.ezfio' : {
},
'SiH-.ezfio' : {
},
'SiH.ezfio' : {
},
'SiS.ezfio' : {
},
}

35
src/BiInts/tests/test_ao_map.irp.f Normal file
View File

@ -0,0 +1,35 @@
program test_ao_map
implicit none
integer :: i,j,k,l
double precision :: get_ao_bielec_integral, tmp
double precision, allocatable :: buffer_value(:,:)
double precision :: ao_bielec_integral,tmp2
integer :: non_zero_int
double precision :: s
integer :: OK
allocate(buffer_value(ao_num,2))
OK = 1
do l=1,ao_num,7
s = 0.d0
do k=1, ao_num
do j=1,l
call compute_ao_bielec_integrals(j,k,l,ao_num,buffer_value(1,1))
call get_ao_bielec_integrals(j,k,l,ao_num,buffer_value(1,2))
do i=1,k
s += (buffer_value(i,1)-buffer_value(i,2))*(buffer_value(i,1)-buffer_value(i,2))
enddo
enddo
enddo
if (s > 1.d-20) then
print *, l, ' : ', s
OK=0
endif
enddo
deallocate(buffer_value)
print *, ' OK : ', OK
end

635
src/BiInts/tests/test_ao_map.ref Normal file
View File

@ -0,0 +1,635 @@
data = {
'AlCl.ezfio' : {
'OK' : 1,
},
'Al.ezfio' : {
'OK' : 1,
},
'Al+.ezfio' : {
'OK' : 1,
},
'AlH2.ezfio' : {
'OK' : 1,
},
'AlH3.ezfio' : {
'OK' : 1,
},
'AlH.ezfio' : {
'OK' : 1,
},
'BCl.ezfio' : {
'OK' : 1,
},
'BeCl.ezfio' : {
'OK' : 1,
},
'Be.ezfio' : {
'OK' : 1,
},
'Be+.ezfio' : {
'OK' : 1,
},
'BeF.ezfio' : {
'OK' : 1,
},
'BeH2.ezfio' : {
'OK' : 1,
},
'BeH.ezfio' : {
'OK' : 1,
},
'BeO.ezfio' : {
'OK' : 1,
},
'BeOH.ezfio' : {
'OK' : 1,
},
'BeS.ezfio' : {
'OK' : 1,
},
'B.ezfio' : {
'OK' : 1,
},
'B+.ezfio' : {
'OK' : 1,
},
'BH2.ezfio' : {
'OK' : 1,
},
'BH3.ezfio' : {
'OK' : 1,
},
'BH.ezfio' : {
'OK' : 1,
},
'BO.ezfio' : {
'OK' : 1,
},
'BS.ezfio' : {
'OK' : 1,
},
'C2.ezfio' : {
'OK' : 1,
},
'C2H2.ezfio' : {
'OK' : 1,
},
'C2H2+.ezfio' : {
'OK' : 1,
},
'C2H3.ezfio' : {
'OK' : 1,
},
'C2H3+.ezfio' : {
'OK' : 1,
},
'C2H4.ezfio' : {
'OK' : 1,
},
'C2H4+.ezfio' : {
'OK' : 1,
},
'C2H5.ezfio' : {
'OK' : 1,
},
'C2H6.ezfio' : {
'OK' : 1,
},
'C2H.ezfio' : {
'OK' : 1,
},
'CCl.ezfio' : {
'OK' : 1,
},
'C-.ezfio' : {
'OK' : 1,
},
'C.ezfio' : {
'OK' : 1,
},
'C+.ezfio' : {
'OK' : 1,
},
'CF.ezfio' : {
'OK' : 1,
},
'CH2_1A1.ezfio' : {
'OK' : 1,
},
'CH2_3B1.ezfio' : {
'OK' : 1,
},
'CH2-.ezfio' : {
'OK' : 1,
},
'CH3-.ezfio' : {
'OK' : 1,
},
'CH3.ezfio' : {
'OK' : 1,
},
'CH4.ezfio' : {
'OK' : 1,
},
'CH4+.ezfio' : {
'OK' : 1,
},
'CH-.ezfio' : {
'OK' : 1,
},
'CH.ezfio' : {
'OK' : 1,
},
'Cl2-.ezfio' : {
'OK' : 1,
},
'Cl2.ezfio' : {
'OK' : 1,
},
'Cl2+.ezfio' : {
'OK' : 1,
},
'Cl-.ezfio' : {
'OK' : 1,
},
'Cl.ezfio' : {
'OK' : 1,
},
'Cl+.ezfio' : {
'OK' : 1,
},
'ClH2+.ezfio' : {
'OK' : 1,
},
'ClH.ezfio' : {
'OK' : 1,
},
'ClH+.ezfio' : {
'OK' : 1,
},
'ClS.ezfio' : {
'OK' : 1,
},
'ClSiH3.ezfio' : {
'OK' : 1,
},
'CN-.ezfio' : {
'OK' : 1,
},
'CN.ezfio' : {
'OK' : 1,
},
'CO2.ezfio' : {
'OK' : 1,
},
'CO.ezfio' : {
'OK' : 1,
},
'CO+.ezfio' : {
'OK' : 1,
},
'COS.ezfio' : {
'OK' : 1,
},
'CP.ezfio' : {
'OK' : 1,
},
'CS2.ezfio' : {
'OK' : 1,
},
'CS.ezfio' : {
'OK' : 1,
},
'CS+.ezfio' : {
'OK' : 1,
},
'CSi.ezfio' : {
'OK' : 1,
},
'F2.ezfio' : {
'OK' : 1,
},
'FAl.ezfio' : {
'OK' : 1,
},
'FCl.ezfio' : {
'OK' : 1,
},
'FCl+.ezfio' : {
'OK' : 1,
},
'F-.ezfio' : {
'OK' : 1,
},
'F.ezfio' : {
'OK' : 1,
},
'F+.ezfio' : {
'OK' : 1,
},
'FH.ezfio' : {
'OK' : 1,
},
'FH+.ezfio' : {
'OK' : 1,
},
'FMg.ezfio' : {
'OK' : 1,
},
'FNa.ezfio' : {
'OK' : 1,
},
'FP.ezfio' : {
'OK' : 1,
},
'FS.ezfio' : {
'OK' : 1,
},
'FSi.ezfio' : {
'OK' : 1,
},
'FSiH3.ezfio' : {
'OK' : 1,
},
'H2CNH.ezfio' : {
'OK' : 1,
},
'H2CO.ezfio' : {
'OK' : 1,
},
'H2CPH.ezfio' : {
'OK' : 1,
},
'H2CS.ezfio' : {
'OK' : 1,
},
'H2.ezfio' : {
'OK' : 1,
},
'H2NNH2.ezfio' : {
'OK' : 1,
},
'H2PPH2.ezfio' : {
'OK' : 1,
},
'H3CCl.ezfio' : {
'OK' : 1,
},
'H3CF.ezfio' : {
'OK' : 1,
},
'H3CNH2.ezfio' : {
'OK' : 1,
},
'H3COH.ezfio' : {
'OK' : 1,
},
'H3CSH.ezfio' : {
'OK' : 1,
},
'H3SiSiH3.ezfio' : {
'OK' : 1,
},
'HBO.ezfio' : {
'OK' : 1,
},
'HBS.ezfio' : {
'OK' : 1,
},
'HCF.ezfio' : {
'OK' : 1,
},
'HCN.ezfio' : {
'OK' : 1,
},
'HCO.ezfio' : {
'OK' : 1,
},
'HCP.ezfio' : {
'OK' : 1,
},
'H.ezfio' : {
'OK' : 1,
},
'HNO.ezfio' : {
'OK' : 1,
},
'HOCl.ezfio' : {
'OK' : 1,
},
'HOF.ezfio' : {
'OK' : 1,
},
'HOMg.ezfio' : {
'OK' : 1,
},
'HONa.ezfio' : {
'OK' : 1,
},
'HOO.ezfio' : {
'OK' : 1,
},
'HOOH.ezfio' : {
'OK' : 1,
},
'HSSH.ezfio' : {
'OK' : 1,
},
'Li2.ezfio' : {
'OK' : 1,
},
'LiCl.ezfio' : {
'OK' : 1,
},
'Li.ezfio' : {
'OK' : 1,
},
'Li+.ezfio' : {
'OK' : 1,
},
'LiF.ezfio' : {
'OK' : 1,
},
'LiH.ezfio' : {
'OK' : 1,
},
'LiN.ezfio' : {
'OK' : 1,
},
'LiO.ezfio' : {
'OK' : 1,
},
'LiOH.ezfio' : {
'OK' : 1,
},
'MgCl.ezfio' : {
'OK' : 1,
},
'Mg.ezfio' : {
'OK' : 1,
},
'Mg+.ezfio' : {
'OK' : 1,
},
'MgH.ezfio' : {
'OK' : 1,
},
'MgS.ezfio' : {
'OK' : 1,
},
'N2.ezfio' : {
'OK' : 1,
},
'N2+.ezfio' : {
'OK' : 1,
},
'Na2.ezfio' : {
'OK' : 1,
},
'NaCl.ezfio' : {
'OK' : 1,
},
'Na.ezfio' : {
'OK' : 1,
},
'Na+.ezfio' : {
'OK' : 1,
},
'NaH.ezfio' : {
'OK' : 1,
},
'N.ezfio' : {
'OK' : 1,
},
'N+.ezfio' : {
'OK' : 1,
},
'NF.ezfio' : {
'OK' : 1,
},
'NH2-.ezfio' : {
'OK' : 1,
},
'NH2.ezfio' : {
'OK' : 1,
},
'NH3.ezfio' : {
'OK' : 1,
},
'NH3+.ezfio' : {
'OK' : 1,
},
'NH4+.ezfio' : {
'OK' : 1,
},
'NH-.ezfio' : {
'OK' : 1,
},
'NH.ezfio' : {
'OK' : 1,
},
'NO-.ezfio' : {
'OK' : 1,
},
'NO.ezfio' : {
'OK' : 1,
},
'NP.ezfio' : {
'OK' : 1,
},
'NS.ezfio' : {
'OK' : 1,
},
'NSi.ezfio' : {
'OK' : 1,
},
'O2Cl.ezfio' : {
'OK' : 1,
},
'O2-.ezfio' : {
'OK' : 1,
},
'O2.ezfio' : {
'OK' : 1,
},
'O2+.ezfio' : {
'OK' : 1,
},
'O2S.ezfio' : {
'OK' : 1,
},
'O2Si.ezfio' : {
'OK' : 1,
},
'O3.ezfio' : {
'OK' : 1,
},
'OCl.ezfio' : {
'OK' : 1,
},
'O-.ezfio' : {
'OK' : 1,
},
'O.ezfio' : {
'OK' : 1,
},
'O+.ezfio' : {
'OK' : 1,
},
'OH2.ezfio' : {
'OK' : 1,
},
'OH2+.ezfio' : {
'OK' : 1,
},
'OH3+.ezfio' : {
'OK' : 1,
},
'OH-.ezfio' : {
'OK' : 1,
},
'OH.ezfio' : {
'OK' : 1,
},
'OH+.ezfio' : {
'OK' : 1,
},
'OMg.ezfio' : {
'OK' : 1,
},
'ONa.ezfio' : {
'OK' : 1,
},
'OP-.ezfio' : {
'OK' : 1,
},
'OP.ezfio' : {
'OK' : 1,
},
'OPH.ezfio' : {
'OK' : 1,
},
'OS.ezfio' : {
'OK' : 1,
},
'OSi.ezfio' : {
'OK' : 1,
},
'P2.ezfio' : {
'OK' : 1,
},
'P2+.ezfio' : {
'OK' : 1,
},
'PCl.ezfio' : {
'OK' : 1,
},
'P-.ezfio' : {
'OK' : 1,
},
'P.ezfio' : {
'OK' : 1,
},
'PH2-.ezfio' : {
'OK' : 1,
},
'PH2.ezfio' : {
'OK' : 1,
},
'PH2+.ezfio' : {
'OK' : 1,
},
'PH3.ezfio' : {
'OK' : 1,
},
'PH3+.ezfio' : {
'OK' : 1,
},
'PH4+.ezfio' : {
'OK' : 1,
},
'PH-.ezfio' : {
'OK' : 1,
},
'PH.ezfio' : {
'OK' : 1,
},
'PS.ezfio' : {
'OK' : 1,
},
'S2-.ezfio' : {
'OK' : 1,
},
'S2.ezfio' : {
'OK' : 1,
},
'S-.ezfio' : {
'OK' : 1,
},
'S.ezfio' : {
'OK' : 1,
},
'S+.ezfio' : {
'OK' : 1,
},
'SH2.ezfio' : {
'OK' : 1,
},
'SH2+.ezfio' : {
'OK' : 1,
},
'SH3+.ezfio' : {
'OK' : 1,
},
'SH-.ezfio' : {
'OK' : 1,
},
'SH.ezfio' : {
'OK' : 1,
},
'SH+.ezfio' : {
'OK' : 1,
},
'Si2.ezfio' : {
'OK' : 1,
},
'SiCl.ezfio' : {
'OK' : 1,
},
'Si-.ezfio' : {
'OK' : 1,
},
'Si.ezfio' : {
'OK' : 1,
},
'SiH2_1A1.ezfio' : {
'OK' : 1,
},
'SiH2_3B1.ezfio' : {
'OK' : 1,
},
'SiH2-.ezfio' : {
'OK' : 1,
},
'SiH3-.ezfio' : {
'OK' : 1,
},
'SiH3.ezfio' : {
'OK' : 1,
},
'SiH4.ezfio' : {
'OK' : 1,
},
'SiH4+.ezfio' : {
'OK' : 1,
},
'SiH-.ezfio' : {
'OK' : 1,
},
'SiH.ezfio' : {
'OK' : 1,
},
'SiS.ezfio' : {
'OK' : 1,
},
}

15
src/BiInts/tests/test_mo.irp.f Normal file
View File

@ -0,0 +1,15 @@
program test
implicit none
integer :: i,j,k,l
double precision :: accu,accu2,get_mo_bielec_integral,accu3
accu = 0.d0
accu3 = 0.d0
do i = 1, elec_beta_num
do j = 1, elec_beta_num
accu += (2.d0 * get_mo_bielec_integral(i,j,i,j,mo_integrals_map) - get_mo_bielec_integral(i,j,j,i,mo_integrals_map) )
accu3 += (2.d0 * mo_bielec_integral_jj(i,j) - mo_bielec_integral_jj_exchange(i,j) )
enddo
enddo
print *,'closed EE 1 : ',accu
print *,'closed EE 2 : ',accu3
end

846
src/BiInts/tests/test_mo.ref Normal file
View File

@ -0,0 +1,846 @@
data = {
'AlCl.ezfio' : {
'closed EE 1' : 324.604108318567,
'closed EE 2' : 324.604108318567,
},
'Al.ezfio' : {
'closed EE 1' : 90.6007174995339,
'closed EE 2' : 90.6007174995339,
},
'Al+.ezfio' : {
'closed EE 1' : 91.1277124469289,
'closed EE 2' : 91.1277124469289,
},
'AlH2.ezfio' : {
'closed EE 1' : 98.3932867799503,
'closed EE 2' : 98.3932867799503,
},
'AlH3.ezfio' : {
'closed EE 1' : 108.322607071765,
'closed EE 2' : 108.322607071765,
},
'AlH.ezfio' : {
'closed EE 1' : 98.9361676259577,
'closed EE 2' : 98.9361676259577,
},
'BCl.ezfio' : {
'closed EE 1' : 208.669206621826,
'closed EE 2' : 208.669206621826,
},
'BeCl.ezfio' : {
'closed EE 1' : 195.478154262549,
'closed EE 2' : 195.478154262549,
},
'Be.ezfio' : {
'closed EE 1' : 4.48960803917240,
'closed EE 2' : 4.48960803917240,
},
'Be+.ezfio' : {
'closed EE 1' : 2.27296598314244,
'closed EE 2' : 2.27296598314244,
},
'BeF.ezfio' : {
'closed EE 1' : 55.2667378076567,
'closed EE 2' : 55.2667378076567,
},
'BeH2.ezfio' : {
'closed EE 1' : 7.91863331588170,
'closed EE 2' : 7.91863331588170,
},
'BeH.ezfio' : {
'closed EE 1' : 4.60494282657902,
'closed EE 2' : 4.60494282657902,
},
'BeO.ezfio' : {
'closed EE 1' : 47.0178059698967,
'closed EE 2' : 47.0178059698967,
},
'BeOH.ezfio' : {
'closed EE 1' : 47.4095700462092,
'closed EE 2' : 47.4095700462092,
},
'BeS.ezfio' : {
'closed EE 1' : 177.044582851287,
'closed EE 2' : 177.044582851287,
},
'B.ezfio' : {
'closed EE 1' : 5.87051480040891,
'closed EE 2' : 5.87051480040891,
},
'B+.ezfio' : {
'closed EE 1' : 6.08872058329241,
'closed EE 2' : 6.08872058329241,
},
'BH2.ezfio' : {
'closed EE 1' : 9.82092714954701,
'closed EE 2' : 9.82092714954701,
},
'BH3.ezfio' : {
'closed EE 1' : 15.1738039588984,
'closed EE 2' : 15.1738039588984,
},
'BH.ezfio' : {
'closed EE 1' : 9.89523045878027,
'closed EE 2' : 9.89523045878027,
},
'BO.ezfio' : {
'closed EE 1' : 49.5244841428460,
'closed EE 2' : 49.5244841428460,
},
'BS.ezfio' : {
'closed EE 1' : 179.847811135927,
'closed EE 2' : 179.847811135927,
},
'C2.ezfio' : {
'closed EE 1' : 40.1266484535661,
'closed EE 2' : 40.1266484535661,
},
'C2H2.ezfio' : {
'closed EE 1' : 50.0072453359780,
'closed EE 2' : 50.0072453359780,
},
'C2H2+.ezfio' : {
'closed EE 1' : 38.8813663363468,
'closed EE 2' : 38.8813663363468,
},
'C2H3.ezfio' : {
'closed EE 1' : 47.7960364243417,
'closed EE 2' : 47.7960364243417,
},
'C2H3+.ezfio' : {
'closed EE 1' : 49.4483058762631,
'closed EE 2' : 49.4483058762631,
},
'C2H4.ezfio' : {
'closed EE 1' : 58.5930731008834,
'closed EE 2' : 58.5930731008834,
},
'C2H4+.ezfio' : {
'closed EE 1' : 46.8050385693024,
'closed EE 2' : 46.8050385693024,
},
'C2H5.ezfio' : {
'closed EE 1' : 55.8648362479812,
'closed EE 2' : 55.8648362479812,
},
'C2H6.ezfio' : {
'closed EE 1' : 67.4231505770701,
'closed EE 2' : 67.4231505770701,
},
'C2H.ezfio' : {
'closed EE 1' : 39.7941997946625,
'closed EE 2' : 39.7941997946625,
},
'CCl.ezfio' : {
'closed EE 1' : 211.332787485635,
'closed EE 2' : 211.332787485635,
},
'C-.ezfio' : {
'closed EE 1' : 7.05668474724391,
'closed EE 2' : 7.05668474724391,
},
'C.ezfio' : {
'closed EE 1' : 7.22275234701424,
'closed EE 2' : 7.22275234701424,
},
'C+.ezfio' : {
'closed EE 1' : 7.43316493990123,
'closed EE 2' : 7.43316493990123,
},
'CF.ezfio' : {
'closed EE 1' : 67.5749607457798,
'closed EE 2' : 67.5749607457798,
},
'CH2_1A1.ezfio' : {
'closed EE 1' : 18.6032101860148,
'closed EE 2' : 18.6032101860148,
},
'CH2_3B1.ezfio' : {
'closed EE 1' : 11.7811748791077,
'closed EE 2' : 11.7811748791077,
},
'CH2-.ezfio' : {
'closed EE 1' : 17.4224231187863,
'closed EE 2' : 17.4224231187863,
},
'CH3-.ezfio' : {
'closed EE 1' : 25.0349104359763,
'closed EE 2' : 25.0349104359763,
},
'CH3.ezfio' : {
'closed EE 1' : 18.1595297280646,
'closed EE 2' : 18.1595297280646,
},
'CH4.ezfio' : {
'closed EE 1' : 26.0266830207216,
'closed EE 2' : 26.0266830207216,
},
'CH4+.ezfio' : {
'closed EE 1' : 18.4285252097398,
'closed EE 2' : 18.4285252097398,
},
'CH-.ezfio' : {
'closed EE 1' : 11.4343449775700,
'closed EE 2' : 11.4343449775700,
},
'CH.ezfio' : {
'closed EE 1' : 12.0838764958443,
'closed EE 2' : 12.0838764958443,
},
'Cl2-.ezfio' : {
'closed EE 1' : 402.671444472536,
'closed EE 2' : 402.671444472536,
},
'Cl2.ezfio' : {
'closed EE 1' : 425.050531461454,
'closed EE 2' : 425.050531461454,
},
'Cl2+.ezfio' : {
'closed EE 1' : 403.762363169362,
'closed EE 2' : 403.762363169362,
},
'Cl-.ezfio' : {
'closed EE 1' : 182.327075272228,
'closed EE 2' : 182.327075272228,
},
'Cl.ezfio' : {
'closed EE 1' : 165.731725346204,
'closed EE 2' : 165.731725346204,
},
'Cl+.ezfio' : {
'closed EE 1' : 149.193201694481,
'closed EE 2' : 149.193201694481,
},
'ClH2+.ezfio' : {
'closed EE 1' : 182.780324107490,
'closed EE 2' : 182.780324107490,
},
'ClH.ezfio' : {
'closed EE 1' : 182.691147725717,
'closed EE 2' : 182.691147725717,
},
'ClH+.ezfio' : {
'closed EE 1' : 165.325459132812,
'closed EE 2' : 165.325459132812,
},
'ClS.ezfio' : {
'closed EE 1' : 385.435861712905,
'closed EE 2' : 385.435861712905,
},
'ClSiH3.ezfio' : {
'closed EE 1' : 372.487466205289,
'closed EE 2' : 372.487466205289,
},
'CN-.ezfio' : {
'closed EE 1' : 55.3766122391237,
'closed EE 2' : 55.3766122391237,
},
'CN.ezfio' : {
'closed EE 1' : 45.1366654659749,
'closed EE 2' : 45.1366654659749,
},
'CO2.ezfio' : {
'closed EE 1' : 126.003917383945,
'closed EE 2' : 126.003917383945,
},
'CO.ezfio' : {
'closed EE 1' : 62.7345020629198,
'closed EE 2' : 62.7345020629198,
},
'CO+.ezfio' : {
'closed EE 1' : 52.0257806106742,
'closed EE 2' : 52.0257806106742,
},
'COS.ezfio' : {
'closed EE 1' : 270.518850152602,
'closed EE 2' : 270.518850152602,
},
'CP.ezfio' : {
'closed EE 1' : 165.996004819480,
'closed EE 2' : 165.996004819480,
},
'CS2.ezfio' : {
'closed EE 1' : 423.926547892897,
'closed EE 2' : 423.926547892897,
},
'CS.ezfio' : {
'closed EE 1' : 196.955019712066,
'closed EE 2' : 196.955019712066,
},
'CS+.ezfio' : {
'closed EE 1' : 183.101079356045,
'closed EE 2' : 183.101079356045,
},
'CSi.ezfio' : {
'closed EE 1' : 136.339370158432,
'closed EE 2' : 136.339370158432,
},
'F2.ezfio' : {
'closed EE 1' : 109.589586860957,
'closed EE 2' : 109.589586860957,
},
'FAl.ezfio' : {
'closed EE 1' : 172.329074560241,
'closed EE 2' : 172.329074560241,
},
'FCl.ezfio' : {
'closed EE 1' : 263.250821798940,
'closed EE 2' : 263.250821798940,
},
'FCl+.ezfio' : {
'closed EE 1' : 243.398671131879,
'closed EE 2' : 243.398671131879,
},
'F-.ezfio' : {
'closed EE 1' : 45.5127217885110,
'closed EE 2' : 45.5127217885110,
},
'F.ezfio' : {
'closed EE 1' : 32.5100994559105,
'closed EE 2' : 32.5100994559105,
},
'F+.ezfio' : {
'closed EE 1' : 21.0616333765295,
'closed EE 2' : 21.0616333765295,
},
'FH.ezfio' : {
'closed EE 1' : 45.4165968585922,
'closed EE 2' : 45.4165968585922,
},
'FH+.ezfio' : {
'closed EE 1' : 32.0966404042882,
'closed EE 2' : 32.0966404042882,
},
'FMg.ezfio' : {
'closed EE 1' : 147.723081703366,
'closed EE 2' : 147.723081703366,
},
'FNa.ezfio' : {
'closed EE 1' : 136.229809908309,
'closed EE 2' : 136.229809908309,
},
'FP.ezfio' : {
'closed EE 1' : 195.769918393858,
'closed EE 2' : 195.769918393858,
},
'FS.ezfio' : {
'closed EE 1' : 226.988766049930,
'closed EE 2' : 226.988766049930,
},
'FSi.ezfio' : {
'closed EE 1' : 184.557018579786,
'closed EE 2' : 184.557018579786,
},
'FSiH3.ezfio' : {
'closed EE 1' : 213.848410787899,
'closed EE 2' : 213.848410787899,
},
'H2CNH.ezfio' : {
'closed EE 1' : 64.8118892917432,
'closed EE 2' : 64.8118892917432,
},
'H2CO.ezfio' : {
'closed EE 1' : 72.1212876559760,
'closed EE 2' : 72.1212876559760,
},
'H2CPH.ezfio' : {
'closed EE 1' : 191.727599087985,
'closed EE 2' : 191.727599087985,
},
'H2CS.ezfio' : {
'closed EE 1' : 209.231104868214,
'closed EE 2' : 209.231104868214,
},
'H2.ezfio' : {
'closed EE 1' : 0.658003361873990,
'closed EE 2' : 0.658003361873990,
},
'H2NNH2.ezfio' : {
'closed EE 1' : 80.4177432036452,
'closed EE 2' : 80.4177432036452,
},
'H2PPH2.ezfio' : {
'closed EE 1' : 348.965444263036,
'closed EE 2' : 348.965444263036,
},
'H3CCl.ezfio' : {
'closed EE 1' : 239.272625318642,
'closed EE 2' : 239.272625318642,
},
'H3CF.ezfio' : {
'closed EE 1' : 89.9110566746969,
'closed EE 2' : 89.9110566746969,
},
'H3CNH2.ezfio' : {
'closed EE 1' : 74.0553942748957,
'closed EE 2' : 74.0553942748957,
},
'H3COH.ezfio' : {
'closed EE 1' : 81.4464317182977,
'closed EE 2' : 81.4464317182977,
},
'H3CSH.ezfio' : {
'closed EE 1' : 220.631290723648,
'closed EE 2' : 220.631290723648,
},
'H3SiSiH3.ezfio' : {
'closed EE 1' : 313.726268238800,
'closed EE 2' : 313.726268238800,
},
'HBO.ezfio' : {
'closed EE 1' : 58.0559325869888,
'closed EE 2' : 58.0559325869888,
},
'HBS.ezfio' : {
'closed EE 1' : 191.322628710909,
'closed EE 2' : 191.322628710909,
},
'HCF.ezfio' : {
'closed EE 1' : 79.0947411389797,
'closed EE 2' : 79.0947411389797,
},
'HCN.ezfio' : {
'closed EE 1' : 55.7130916664010,
'closed EE 2' : 55.7130916664010,
},
'HCO.ezfio' : {
'closed EE 1' : 60.6786870182187,
'closed EE 2' : 60.6786870182187,
},
'HCP.ezfio' : {
'closed EE 1' : 180.296919079404,
'closed EE 2' : 180.296919079404,
},
'H.ezfio' : {
'closed EE 1' : 0.000000000000000E+000,
'closed EE 2' : 0.000000000000000E+000,
},
'HNO.ezfio' : {
'closed EE 1' : 77.5573550103421,
'closed EE 2' : 77.5573550103421,
},
'HOCl.ezfio' : {
'closed EE 1' : 253.799452346657,
'closed EE 2' : 253.799452346657,
},
'HOF.ezfio' : {
'closed EE 1' : 101.663005370608,
'closed EE 2' : 101.663005370608,
},
'HOMg.ezfio' : {
'closed EE 1' : 138.443726972836,
'closed EE 2' : 138.443726972836,
},
'HONa.ezfio' : {
'closed EE 1' : 126.747456382832,
'closed EE 2' : 126.747456382832,
},
'HOO.ezfio' : {
'closed EE 1' : 79.4438064643797,
'closed EE 2' : 79.4438064643797,
},
'HOOH.ezfio' : {
'closed EE 1' : 93.6527230199755,
'closed EE 2' : 93.6527230199755,
},
'HSSH.ezfio' : {
'closed EE 1' : 387.218345293825,
'closed EE 2' : 387.218345293825,
},
'Li2.ezfio' : {
'closed EE 1' : 6.46521052077037,
'closed EE 2' : 6.46521052077037,
},
'LiCl.ezfio' : {
'closed EE 1' : 193.091759412769,
'closed EE 2' : 193.091759412769,
},
'Li.ezfio' : {
'closed EE 1' : 1.64991166351969,
'closed EE 2' : 1.64991166351969,
},
'Li+.ezfio' : {
'closed EE 1' : 1.65066035543863,
'closed EE 2' : 1.65066035543863,
},
'LiF.ezfio' : {
'closed EE 1' : 53.4448357111712,
'closed EE 2' : 53.4448357111712,
},
'LiH.ezfio' : {
'closed EE 1' : 3.45968704863106,
'closed EE 2' : 3.45968704863106,
},
'LiN.ezfio' : {
'closed EE 1' : 19.4901024446162,
'closed EE 2' : 19.4901024446162,
},
'LiO.ezfio' : {
'closed EE 1' : 33.6641656402794,
'closed EE 2' : 33.6641656402794,
},
'LiOH.ezfio' : {
'closed EE 1' : 45.3987062615096,
'closed EE 2' : 45.3987062615096,
},
'MgCl.ezfio' : {
'closed EE 1' : 297.457458920467,
'closed EE 2' : 297.457458920467,
},
'Mg.ezfio' : {
'closed EE 1' : 79.8152268844060,
'closed EE 2' : 79.8152268844060,
},
'Mg+.ezfio' : {
'closed EE 1' : 72.1282818471152,
'closed EE 2' : 72.1282818471152,
},
'MgH.ezfio' : {
'closed EE 1' : 79.3190204144503,
'closed EE 2' : 79.3190204144503,
},
'MgS.ezfio' : {
'closed EE 1' : 279.954838755823,
'closed EE 2' : 279.954838755823,
},
'N2.ezfio' : {
'closed EE 1' : 61.5022395201583,
'closed EE 2' : 61.5022395201583,
},
'N2+.ezfio' : {
'closed EE 1' : 48.5010781414573,
'closed EE 2' : 48.5010781414573,
},
'Na2.ezfio' : {
'closed EE 1' : 153.459222056420,
'closed EE 2' : 153.459222056420,
},
'NaCl.ezfio' : {
'closed EE 1' : 285.552878597068,
'closed EE 2' : 285.552878597068,
},
'Na.ezfio' : {
'closed EE 1' : 63.1438465495254,
'closed EE 2' : 63.1438465495254,
},
'Na+.ezfio' : {
'closed EE 1' : 63.1844537555369,
'closed EE 2' : 63.1844537555369,
},
'NaH.ezfio' : {
'closed EE 1' : 69.3524887715318,
'closed EE 2' : 69.3524887715318,
},
'N.ezfio' : {
'closed EE 1' : 8.56327997038596,
'closed EE 2' : 8.56327997038596,
},
'N+.ezfio' : {
'closed EE 1' : 8.76869142211657,
'closed EE 2' : 8.76869142211657,
},
'NF.ezfio' : {
'closed EE 1' : 69.0051111837754,
'closed EE 2' : 69.0051111837754,
},
'NH2-.ezfio' : {
'closed EE 1' : 30.6309164550511,
'closed EE 2' : 30.6309164550511,
},
'NH2.ezfio' : {
'closed EE 1' : 21.9260220960358,
'closed EE 2' : 21.9260220960358,
},
'NH3.ezfio' : {
'closed EE 1' : 31.4406045429538,
'closed EE 2' : 31.4406045429538,
},
'NH3+.ezfio' : {
'closed EE 1' : 22.1802907966178,
'closed EE 2' : 22.1802907966178,
},
'NH4+.ezfio' : {
'closed EE 1' : 31.8659857288766,
'closed EE 2' : 31.8659857288766,
},
'NH-.ezfio' : {
'closed EE 1' : 21.6689399304262,
'closed EE 2' : 21.6689399304262,
},
'NH.ezfio' : {
'closed EE 1' : 14.2929080048555,
'closed EE 2' : 14.2929080048555,
},
'NO-.ezfio' : {
'closed EE 1' : 61.7198631648835,
'closed EE 2' : 61.7198631648835,
},
'NO.ezfio' : {
'closed EE 1' : 64.6696341812441,
'closed EE 2' : 64.6696341812441,
},
'NP.ezfio' : {
'closed EE 1' : 187.100108926969,
'closed EE 2' : 187.100108926969,
},
'NS.ezfio' : {
'closed EE 1' : 200.128710693974,
'closed EE 2' : 200.128710693974,
},
'NSi.ezfio' : {
'closed EE 1' : 155.279906480313,
'closed EE 2' : 155.279906480313,
},
'O2Cl.ezfio' : {
'closed EE 1' : 301.747915430568,
'closed EE 2' : 301.747915430568,
},
'O2-.ezfio' : {
'closed EE 1' : 79.2328963139684,
'closed EE 2' : 79.2328963139684,
},
'O2.ezfio' : {
'closed EE 1' : 66.1914541099256,
'closed EE 2' : 66.1914541099256,
},
'O2+.ezfio' : {
'closed EE 1' : 69.4144048026520,
'closed EE 2' : 69.4144048026520,
},
'O2S.ezfio' : {
'closed EE 1' : 312.243487598150,
'closed EE 2' : 312.243487598150,
},
'O2Si.ezfio' : {
'closed EE 1' : 257.255582658461,
'closed EE 2' : 257.255582658461,
},
'O3.ezfio' : {
'closed EE 1' : 152.241108922056,
'closed EE 2' : 152.241108922056,
},
'OCl.ezfio' : {
'closed EE 1' : 236.412413075442,
'closed EE 2' : 236.412413075442,
},
'O-.ezfio' : {
'closed EE 1' : 26.9389565391700,
'closed EE 2' : 26.9389565391700,
},
'O.ezfio' : {
'closed EE 1' : 17.7360424805971,
'closed EE 2' : 17.7360424805971,
},
'O+.ezfio' : {
'closed EE 1' : 10.1017140483289,
'closed EE 2' : 10.1017140483289,
},
'OH2.ezfio' : {
'closed EE 1' : 37.8440317785205,
'closed EE 2' : 37.8440317785205,
},
'OH2+.ezfio' : {
'closed EE 1' : 26.6519594351168,
'closed EE 2' : 26.6519594351168,
},
'OH3+.ezfio' : {
'closed EE 1' : 38.1275235387040,
'closed EE 2' : 38.1275235387040,
},
'OH-.ezfio' : {
'closed EE 1' : 37.4117044437012,
'closed EE 2' : 37.4117044437012,
},
'OH.ezfio' : {
'closed EE 1' : 26.7707296794756,
'closed EE 2' : 26.7707296794756,
},
'OH+.ezfio' : {
'closed EE 1' : 17.2889916223005,
'closed EE 2' : 17.2889916223005,
},
'OMg.ezfio' : {
'closed EE 1' : 138.801601522175,
'closed EE 2' : 138.801601522175,
},
'ONa.ezfio' : {
'closed EE 1' : 112.111906727576,
'closed EE 2' : 112.111906727576,
},
'OP-.ezfio' : {
'closed EE 1' : 188.573528250670,
'closed EE 2' : 188.573528250670,
},
'OP.ezfio' : {
'closed EE 1' : 191.910766225339,
'closed EE 2' : 191.910766225339,
},
'OPH.ezfio' : {
'closed EE 1' : 208.214788904911,
'closed EE 2' : 208.214788904911,
},
'OS.ezfio' : {
'closed EE 1' : 203.101155235861,
'closed EE 2' : 203.101155235861,
},
'OSi.ezfio' : {
'closed EE 1' : 179.456084306272,
'closed EE 2' : 179.456084306272,
},
'P2.ezfio' : {
'closed EE 1' : 323.782785949801,
'closed EE 2' : 323.782785949801,
},
'P2+.ezfio' : {
'closed EE 1' : 301.882738442470,
'closed EE 2' : 301.882738442470,
},
'PCl.ezfio' : {
'closed EE 1' : 349.890715378447,
'closed EE 2' : 349.890715378447,
},
'P-.ezfio' : {
'closed EE 1' : 122.657805807297,
'closed EE 2' : 122.657805807297,
},
'P.ezfio' : {
'closed EE 1' : 111.445647659427,
'closed EE 2' : 111.445647659427,
},
'PH2-.ezfio' : {
'closed EE 1' : 147.249988541779,
'closed EE 2' : 147.249988541779,
},
'PH2.ezfio' : {
'closed EE 1' : 134.714659285095,
'closed EE 2' : 134.714659285095,
},
'PH2+.ezfio' : {
'closed EE 1' : 136.440910105336,
'closed EE 2' : 136.440910105336,
},
'PH3.ezfio' : {
'closed EE 1' : 148.202871342844,
'closed EE 2' : 148.202871342844,
},
'PH3+.ezfio' : {
'closed EE 1' : 134.975246396259,
'closed EE 2' : 134.975246396259,
},
'PH4+.ezfio' : {
'closed EE 1' : 148.989370984817,
'closed EE 2' : 148.989370984817,
},
'PH-.ezfio' : {
'closed EE 1' : 134.328964674983,
'closed EE 2' : 134.328964674983,
},
'PH.ezfio' : {
'closed EE 1' : 122.455536217569,
'closed EE 2' : 122.455536217569,
},
'PS.ezfio' : {
'closed EE 1' : 337.453574046369,
'closed EE 2' : 337.453574046369,
},
'S2-.ezfio' : {
'closed EE 1' : 366.675175225190,
'closed EE 2' : 366.675175225190,
},
'S2.ezfio' : {
'closed EE 1' : 349.448868216573,
'closed EE 2' : 349.448868216573,
},
'S-.ezfio' : {
'closed EE 1' : 149.700149385399,
'closed EE 2' : 149.700149385399,
},
'S.ezfio' : {
'closed EE 1' : 135.897852467643,
'closed EE 2' : 135.897852467643,
},
'S+.ezfio' : {
'closed EE 1' : 122.164730812357,
'closed EE 2' : 122.164730812357,
},
'SH2.ezfio' : {
'closed EE 1' : 164.906914972667,
'closed EE 2' : 164.906914972667,
},
'SH2+.ezfio' : {
'closed EE 1' : 149.596440849195,
'closed EE 2' : 149.596440849195,
},
'SH3+.ezfio' : {
'closed EE 1' : 165.246577302836,
'closed EE 2' : 165.246577302836,
},
'SH-.ezfio' : {
'closed EE 1' : 164.245517624070,
'closed EE 2' : 164.245517624070,
},
'SH.ezfio' : {
'closed EE 1' : 149.754586732578,
'closed EE 2' : 149.754586732578,
},
'SH+.ezfio' : {
'closed EE 1' : 135.225667181454,
'closed EE 2' : 135.225667181454,
},
'Si2.ezfio' : {
'closed EE 1' : 252.066473929603,
'closed EE 2' : 252.066473929603,
},
'SiCl.ezfio' : {
'closed EE 1' : 337.640041152210,
'closed EE 2' : 337.640041152210,
},
'Si-.ezfio' : {
'closed EE 1' : 100.662878826237,
'closed EE 2' : 100.662878826237,
},
'Si.ezfio' : {
'closed EE 1' : 101.098004923652,
'closed EE 2' : 101.098004923652,
},
'SiH2_1A1.ezfio' : {
'closed EE 1' : 121.587026002297,
'closed EE 2' : 121.587026002297,
},
'SiH2_3B1.ezfio' : {
'closed EE 1' : 109.761390873831,
'closed EE 2' : 109.761390873831,
},
'SiH2-.ezfio' : {
'closed EE 1' : 120.021930833955,
'closed EE 2' : 120.021930833955,
},
'SiH3-.ezfio' : {
'closed EE 1' : 131.352511652808,
'closed EE 2' : 131.352511652808,
},
'SiH3.ezfio' : {
'closed EE 1' : 120.791742854126,
'closed EE 2' : 120.791742854126,
},
'SiH4.ezfio' : {
'closed EE 1' : 132.863104887314,
'closed EE 2' : 132.863104887314,
},
'SiH4+.ezfio' : {
'closed EE 1' : 121.801965974605,
'closed EE 2' : 121.801965974605,
},
'SiH-.ezfio' : {
'closed EE 1' : 109.783147473057,
'closed EE 2' : 109.783147473057,
},
'SiH.ezfio' : {
'closed EE 1' : 110.775916738616,
'closed EE 2' : 110.775916738616,
},
'SiS.ezfio' : {
'closed EE 1' : 324.574435129198,
'closed EE 2' : 324.574435129198,
},
}

4
src/Bitmask/ASSUMPTIONS.rst
View File

@ -1,4 +1,6 @@
* ``bit_kind_shift``, ``bit_kind_size`` and ``bit_kind`` are coherent::
``bit_kind_shift``, ``bit_kind_size`` and ``bit_kind`` are coherent:
.. code_block:: fortran
2**bit_kind_shift = bit_kind_size
bit_kind = bit_kind_size / 8

4
src/Bitmask/README.rst
View File

@ -24,7 +24,9 @@ Assumptions
.. Do not edit this section. It was auto-generated from the
.. ASSUMPTIONS.rst file.
* ``bit_kind_shift``, ``bit_kind_size`` and ``bit_kind`` are coherent::
``bit_kind_shift``, ``bit_kind_size`` and ``bit_kind`` are coherent:
.. code_block:: fortran
2**bit_kind_shift = bit_kind_size
bit_kind = bit_kind_size / 8

0
src/MOGuess/ASSUMPTIONS.rst Normal file
View File

11
src/MOGuess/H_CORE_guess.irp.f Normal file
View File

@ -0,0 +1,11 @@
program H_CORE_guess
implicit none
character*(64) :: label
mo_coef = ao_ortho_lowdin_coef
TOUCH mo_coef
label = "H_CORE_GUESS"
call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),label)
call save_mos
end

8
src/MOGuess/Makefile Normal file
View File

@ -0,0 +1,8 @@
default: all
# Define here all new external source files and objects.Don't forget to prefix the
# object files with IRPF90_temp/
SRC=
OBJ=
include $(QPACKAGE_ROOT)/src/Makefile.common

4
src/MOGuess/README.rst Normal file
View File

@ -0,0 +1,4 @@
==============
MOGuess Module
==============

114
src/MOGuess/localize.irp.f Normal file
View File

@ -0,0 +1,114 @@
!TODO Ecrire un cholesky avec bitmask
subroutine localize_mos(mask, nint)
implicit none
use bitmasks
integer, intent(in) :: nint
integer(bit_kind), intent(in) :: mask(nint)
integer :: i,j,k,l
double precision, allocatable :: DM(:,:)
double precision, allocatable :: mo_coef_new(:,:), R(:,:)
integer :: n
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: DM, mo_coef_new, R
integer :: rank
integer, parameter :: n_core = 2
allocate(R(mo_tot_num,mo_tot_num))
allocate(DM(ao_num_align,ao_num))
allocate(mo_coef_new(ao_num_align,mo_tot_num))
n = ao_num
mo_coef_new = mo_coef
BEGIN_TEMPLATE
DM = 0.d0
if ($START < $END) then
do k=$START, $END
do j=1,n
!DEC$ VECTOR ALIGNED
do i=1,n
DM(i,j) += mo_coef_new(i,k)*mo_coef_new(j,k)
enddo
enddo
enddo
call cholesky_mo(n,$END-$START+1,DM,mo_coef_new(1,$START),size(mo_coef_new,1),-1.d0,rank)
endif
SUBST [ START, END ]
1 ; n_core ;;
END_TEMPLATE
deallocate(DM)
call find_rotation(mo_coef,ao_num_align,mo_coef_new,ao_num,R,mo_tot_num)
mo_coef = mo_coef_new
deallocate(mo_coef_new)
double precision,allocatable :: mo_energy_new(:)
integer, allocatable :: iorder(:)
allocate(mo_energy_new(mo_tot_num),iorder(mo_tot_num))
do i=1,mo_tot_num
iorder(i) = i
mo_energy_new(i) = 0.d0
do k=1,mo_tot_num
mo_energy_new(i) += R(k,i)*R(k,i)*mo_energy(k)
enddo
enddo
mo_energy = mo_energy_new
call dsort(mo_energy(1),iorder(1),n_core)
allocate (mo_coef_new(ao_num_align,mo_tot_num))
mo_coef_new = mo_coef
do j=1,mo_tot_num
do i=1,ao_num
mo_coef(i,j) = mo_coef_new(i,iorder(j))
enddo
enddo
deallocate (mo_coef_new,R)
deallocate(mo_energy_new,iorder)
mo_label = 'localized'
SOFT_TOUCH mo_coef mo_energy mo_label
end
subroutine cholesky_mo(n,m,P,C,LDC,tol_in,rank)
implicit none
integer, intent(in) :: n,m, LDC
double precision, intent(in) :: P(LDC,n)
double precision, intent(out) :: C(LDC,m)
double precision, intent(in) :: tol_in
integer, intent(out) :: rank
integer :: info
integer :: i,k
integer :: ipiv(n)
double precision:: tol
double precision, allocatable :: W(:,:), work(:)
!DEC$ ATTRIBUTES ALIGN: 32 :: W
!DEC$ ATTRIBUTES ALIGN: 32 :: work
!DEC$ ATTRIBUTES ALIGN: 32 :: ipiv
allocate(W(LDC,n),work(2*n))
tol=tol_in
info = 0
do i=1,n
do k=1,i
W(i,k) = P(i,k)
enddo
do k=i+1,n
W(i,k) = 0.
enddo
enddo
call DPSTRF('L', n, W, LDC, ipiv, rank, tol, work, info )
do i=1,n
do k=1,min(m,rank)
C(ipiv(i),k) = W(i,k)
enddo
enddo
deallocate(W,work)
end

46
src/MOGuess/mo_ortho_lowdin.irp.f Normal file
View File

@ -0,0 +1,46 @@
BEGIN_PROVIDER [double precision, ao_ortho_lowdin_coef, (ao_num_align,ao_num)]
BEGIN_DOC
! matrix of the coefficients of the mos generated by the
! orthonormalization by the S^{-1/2} canonical transformation of the aos
! ao_ortho_lowdin_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_lowdin orbital
END_DOC
implicit none
integer :: i,j,k,l
double precision :: tmp_matrix(ao_num_align,ao_num),accu
do i = 1, ao_num
do j = 1, ao_num
tmp_matrix(i,j) = 0.d0
enddo
enddo
do i = 1, ao_num
tmp_matrix(i,i) = 1.d0
enddo
call ortho_lowdin(ao_overlap,ao_num_align,ao_num,tmp_matrix,ao_num_align,ao_num)
do i = 1, ao_num
do j = 1, ao_num
ao_ortho_lowdin_coef(j,i) = tmp_matrix(i,j)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, ao_ortho_lowdin_overlap, (ao_num_align,ao_num)]
BEGIN_DOC
! overlap matrix of the ao_ortho_lowdin
! supposed to be the Identity
END_DOC
implicit none
integer :: i,j,k,l
double precision :: accu
do i = 1, ao_num
do j = 1, ao_num
accu = 0.d0
do k = 1, ao_num
do l = 1, ao_num
accu = accu + ao_ortho_lowdin_coef(i,k) * ao_ortho_lowdin_coef(j,l) * ao_overlap(k,l)
enddo
enddo
ao_ortho_lowdin_overlap(i,j) = accu
enddo
enddo
END_PROVIDER

25
src/MOGuess/pot_mo_ortho_lowdin_ints.irp.f Normal file
View File

@ -0,0 +1,25 @@
BEGIN_PROVIDER [double precision, ao_ortho_lowdin_nucl_elec_integral, (mo_tot_num_align,mo_tot_num)]
implicit none
integer :: i1,j1,i,j
double precision :: c_i1,c_j1
ao_ortho_lowdin_nucl_elec_integral = 0.d0
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP PRIVATE(i,j,i1,j1,c_j1,c_i1) &
!$OMP SHARED(mo_tot_num,ao_num,ao_ortho_lowdin_coef, &
!$OMP ao_ortho_lowdin_nucl_elec_integral, ao_nucl_elec_integral)
do i = 1, mo_tot_num
do j = 1, mo_tot_num
do i1 = 1,ao_num
c_i1 = ao_ortho_lowdin_coef(i1,i)
do j1 = 1,ao_num
c_j1 = c_i1*ao_ortho_lowdin_coef(j1,j)
ao_ortho_lowdin_nucl_elec_integral(j,i) = ao_ortho_lowdin_nucl_elec_integral(j,i) + &
c_j1 * ao_nucl_elec_integral(j1,i1)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER

33
src/MOGuess/tests/Makefile Normal file
View File

@ -0,0 +1,33 @@
OPENMP =1
PROFILE =0
DEBUG = 0
IRPF90+= -I tests
REF_FILES=$(subst %.irp.f, %.ref, $(wildcard *.irp.f))
.PHONY: clean executables serial_tests parallel_tests
all: clean executables serial_tests parallel_tests
parallel_tests: $(REF_FILES)
@echo ; echo " ---- Running parallel tests ----" ; echo
@OMP_NUM_THREADS=10 ${QPACKAGE_ROOT}/scripts/run_tests.py
serial_tests: $(REF_FILES)
@echo ; echo " ---- Running serial tests ----" ; echo
@OMP_NUM_THREADS=1 ${QPACKAGE_ROOT}/scripts/run_tests.py
executables: $(wildcard *.irp.f) veryclean
$(MAKE) -C ..
%.ref: $(wildcard $(QPACKAGE_ROOT)/data/inputs/*.md5) executables
$(QPACKAGE_ROOT)/scripts/create_test_ref.sh $*
clean:
$(MAKE) -C .. clean
veryclean:
$(MAKE) -C .. veryclean

4
src/Utils/Makefile
View File

@ -2,7 +2,7 @@ default: all
# Define here all new external source files and objects.Don't forget to prefix the
# object files with IRPF90_temp/
SRC=
OBJ=
SRC=map_module.f90
OBJ=IRPF90_temp/map_module.o
include $(QPACKAGE_ROOT)/src/Makefile.common

640
src/Utils/map_module.f90 Normal file
View File

@ -0,0 +1,640 @@
module map_module
use omp_lib
integer, parameter :: integral_kind = 8
integer, parameter :: cache_key_kind = 2
integer, parameter :: cache_map_size_kind = 4
integer, parameter :: key_kind = 8
integer, parameter :: map_size_kind = 8
integer, parameter :: map_shift = -15
integer*8, parameter :: map_mask = ibset(0_8,15)-1_8
type cache_map_type
integer(cache_key_kind), pointer :: key(:)
real(integral_kind), pointer :: value(:)
logical :: sorted
integer(cache_map_size_kind) :: map_size
integer(cache_map_size_kind) :: n_elements
integer(omp_lock_kind) :: lock
end type cache_map_type
type map_type
type(cache_map_type), allocatable :: map(:)
integer(map_size_kind) :: map_size
integer(map_size_kind) :: n_elements
logical :: sorted
integer(omp_lock_kind) :: lock
end type map_type
end module map_module
real function map_mb(map)
use map_module
use omp_lib
implicit none
type (map_type), intent(in) :: map
integer(map_size_kind) :: i
map_mb = 8+map_size_kind+map_size_kind+omp_lock_kind+4
do i=0,map%map_size
map_mb = map_mb + map%map(i)%map_size*(cache_key_kind+integral_kind) + &
8+8+4+cache_map_size_kind+cache_map_size_kind+omp_lock_kind
enddo
map_mb = map_mb / (1024.d0*1024.d0)
end
subroutine cache_map_init(map,sze)
use map_module
implicit none
type (cache_map_type), intent(inout) :: map
integer(cache_map_size_kind) :: sze
call omp_init_lock(map%lock)
call omp_set_lock(map%lock)
map%n_elements = 0_8
map%map_size = 0_8
map%sorted = .True.
NULLIFY(map%value, map%key)
call cache_map_reallocate(map,sze)
call omp_unset_lock(map%lock)
end
subroutine map_init(map,keymax)
use map_module
implicit none
integer*8, intent(in) :: keymax
type (map_type), intent(inout) :: map
integer(map_size_kind) :: i
integer(cache_map_size_kind) :: sze
integer :: err
call omp_init_lock(map%lock)
call omp_set_lock(map%lock)
map%n_elements = 0_8
map%map_size = ishft(keymax,map_shift)
allocate(map%map(0_8:map%map_size),stat=err)
if (err /= 0) then
print *, 'Unable to allocate map'
stop 5
endif
!sze = max(sqrt(map%map_size/16.d0),2048.d0)
sze = 2
do i=0_8,map%map_size
call cache_map_init(map%map(i),sze)
enddo
map%sorted = .True.
call omp_unset_lock(map%lock)
end
subroutine cache_map_reallocate(map,sze)
use map_module
implicit none
integer(cache_map_size_kind), intent(in) :: sze
type (cache_map_type), intent(inout) :: map
integer(cache_key_kind), pointer :: key_new(:)
real(integral_kind), pointer :: value_new(:)
integer(map_size_kind) :: i
integer :: err
!DIR$ ATTRIBUTES ALIGN : 64 :: key_new, value_new
if (sze < map%n_elements) then
print *, 'Unable to resize map : map too large'
stop 3
endif
! Resize keys
allocate( key_new(sze), stat=err )
if (err /= 0) then
print *, 'Unable to allocate map', sze
stop 1
endif
if (associated(map%key)) then
do i=1_8,min(size(map%key),map%n_elements)
key_new(i) = map%key(i)
enddo
deallocate(map%key)
endif
! Resize values
allocate( value_new(sze), stat=err )
if (err /= 0) then
print *, 'Unable to allocate map', sze
stop 2
endif
if (associated(map%value)) then
do i=1_8,min(size(map%key),map%n_elements)
value_new(i) = map%value(i)
enddo
deallocate(map%value)
endif
! Set new pointers
map%key => key_new
map%value => value_new
map%map_size = sze
end
subroutine cache_map_deinit(map)
use map_module
implicit none
type (cache_map_type), intent(inout) :: map
integer :: err
if (associated( map % value )) then
deallocate( map % value, stat=err )
if (err /= 0) then
print *, 'Unable to deallocate map'
stop 2
endif
NULLIFY(map%value)
endif
if (associated( map % key )) then
deallocate( map % key, stat=err )
if (err /= 0) then
print *, 'Unable to deallocate map'
stop 4
endif
NULLIFY(map%key)
endif
map%n_elements = 0_8
map%map_size = 0_8
call omp_destroy_lock(map%lock)
end
subroutine map_deinit(map)
use map_module
implicit none
type (map_type), intent(inout) :: map
integer :: err
integer(map_size_kind) :: i
if (allocated( map % map )) then
do i=0_8, map%map_size
call cache_map_deinit(map%map(i))
enddo
deallocate( map % map, stat=err )
if (err /= 0) then
print *, 'Unable to deallocate map'
stop 6
endif
endif
map%n_elements = 0_8
map%map_size = 0_8
call omp_destroy_lock(map%lock)
end
subroutine cache_map_sort(map)
use map_module
implicit none
type (cache_map_type), intent(inout) :: map
integer(cache_map_size_kind), allocatable :: iorder(:)
integer(cache_map_size_kind) :: i
!DIR$ ATTRIBUTES ALIGN : 64 :: iorder
if (.not.map%sorted) then
allocate(iorder(map%n_elements))
do i=1,map%n_elements
iorder(i) = i
enddo
if (cache_key_kind == 2) then
call i2radix_sort(map%key,iorder,map%n_elements,-1)
else if (cache_key_kind == 4) then
call iradix_sort(map%key,iorder,map%n_elements,-1)
else if (cache_key_kind == 8) then
call i8radix_sort(map%key,iorder,map%n_elements,-1)
endif
if (integral_kind == 4) then
call set_order(map%value,iorder,map%n_elements)
else if (integral_kind == 8) then
call dset_order(map%value,iorder,map%n_elements)
endif
deallocate(iorder)
map%sorted = .True.
endif
end
subroutine map_sort(map)
use map_module
implicit none
type (map_type), intent(inout) :: map
integer(map_size_kind) :: i
if (.not.map%sorted) then
!$OMP PARALLEL DO SCHEDULE(static,1024) DEFAULT(SHARED) PRIVATE(i)
do i=0_8,map%map_size
call omp_set_lock(map%map(i)%lock)
call cache_map_sort(map%map(i))
call omp_unset_lock(map%map(i)%lock)
enddo
!$OMP END PARALLEL DO
map%sorted = .True.
endif
end
subroutine cache_map_unique(map)
use map_module
implicit none
type (cache_map_type), intent(inout) :: map
integer(cache_key_kind) :: prev_key
integer(cache_map_size_kind) :: i, j
call cache_map_sort(map)
prev_key = -1_8
j=0
do i=1,map%n_elements
if (map%key(i) /= prev_key) then
j = j+1
map%value(j) = map%value(i)
map%key(j) = map%key(i)
prev_key = map%key(i)
else
map%value(j) = map%value(j)+map%value(i)
endif
enddo
map%n_elements = j
end
subroutine cache_map_shrink(map,thr)
use map_module
implicit none
type (cache_map_type), intent(inout) :: map
real(integral_kind) , intent(in) :: thr
integer(cache_map_size_kind) :: i,j
j=0
do i=1,map%n_elements
if (abs(map%value(i)) > thr) then
j = j+1
map%value(j) = map%value(i)
map%key(j) = map%key(i)
endif
enddo
map%n_elements = j
end
subroutine map_unique(map)
use map_module
implicit none
type (map_type), intent(inout) :: map
integer(map_size_kind) :: i
integer(map_size_kind) :: icount
icount = 0_8
!$OMP PARALLEL DO SCHEDULE(dynamic,1000) DEFAULT(SHARED) PRIVATE(i) &
!$OMP REDUCTION(+:icount)
do i=0_8,map%map_size
call omp_set_lock(map%map(i)%lock)
call cache_map_unique(map%map(i))
call omp_unset_lock(map%map(i)%lock)
icount = icount + map%map(i)%n_elements
enddo
!$OMP END PARALLEL DO
map%n_elements = icount
end
subroutine map_shrink(map,thr)
use map_module
implicit none
type (map_type), intent(inout) :: map
type (map_type), intent(in) :: thr
integer(map_size_kind) :: i
integer(map_size_kind) :: icount
icount = 0_8
!$OMP PARALLEL DO SCHEDULE(dynamic,1000) DEFAULT(SHARED) PRIVATE(i) &
!$OMP REDUCTION(+:icount)
do i=0_8,map%map_size
call omp_set_lock(map%map(i)%lock)
call cache_map_shrink(map%map(i),thr)
call omp_unset_lock(map%map(i)%lock)
icount = icount + map%map(i)%n_elements
enddo
!$OMP END PARALLEL DO
map%n_elements = icount
end
subroutine map_update(map, key, value, sze, thr)
use map_module
implicit none
type (map_type), intent(inout) :: map
integer, intent(in) :: sze
integer(key_kind), intent(inout) :: key(sze)
real(integral_kind), intent(inout):: value(sze)
real(integral_kind), intent(in) :: thr
integer :: i
integer(map_size_kind) :: idx_cache, idx_cache_new
integer(cache_map_size_kind) :: idx
integer :: sze2
integer(cache_key_kind) :: cache_key
integer(map_size_kind) :: n_elements_temp
type (cache_map_type) :: local_map
logical :: map_sorted
sze2 = sze
map_sorted = .True.
n_elements_temp = 0_8
n_elements_temp = n_elements_temp + 1_8
do while (sze2>0)
i=1
do while (i<=sze)
if (key(i) /= 0_8) then
idx_cache = ishft(key(i),map_shift)
if (omp_test_lock(map%map(idx_cache)%lock)) then
local_map%key => map%map(idx_cache)%key
local_map%value => map%map(idx_cache)%value
local_map%sorted = map%map(idx_cache)%sorted
local_map%map_size = map%map(idx_cache)%map_size
local_map%n_elements = map%map(idx_cache)%n_elements
do
!DIR$ FORCEINLINE
call search_key_big_interval(key(i),local_map%key, local_map%n_elements, idx, 1_8, local_map%n_elements)
if (idx > 0_8) then
local_map%value(idx) = local_map%value(idx) + value(i)
else
! Assert that the map has a proper size
if (local_map%n_elements == local_map%map_size) then
call cache_map_unique(local_map)
call cache_map_reallocate(local_map, local_map%n_elements + local_map%n_elements)
call cache_map_shrink(local_map,thr)
endif
cache_key = iand(key(i),map_mask)
local_map%n_elements = local_map%n_elements + 1_8
local_map%value(local_map%n_elements) = value(i)
local_map%key(local_map%n_elements) = cache_key
local_map%sorted = .False.
n_elements_temp = n_elements_temp + 1_8
endif ! idx > 0
key(i) = 0_8
i = i+1
sze2 = sze2-1
if (i>sze) then
i=1
endif
if ( (ishft(key(i),map_shift) /= idx_cache).or.(key(i)==0_8)) then
exit
endif
enddo
map%map(idx_cache)%key => local_map%key
map%map(idx_cache)%value => local_map%value
map%map(idx_cache)%sorted = local_map%sorted
map%map(idx_cache)%n_elements = local_map%n_elements
map%map(idx_cache)%map_size = local_map%map_size
map_sorted = map_sorted .and. local_map%sorted
call omp_unset_lock(map%map(idx_cache)%lock)
endif ! omp_test_lock
else
i=i+1
endif ! key = 0
enddo ! i
enddo ! sze2 > 0
call omp_set_lock(map%lock)
map%n_elements = map%n_elements + n_elements_temp
map%sorted = map%sorted .and. map_sorted
call omp_unset_lock(map%lock)
end
subroutine map_append(map, key, value, sze)
use map_module
implicit none
type (map_type), intent(inout) :: map
integer, intent(in) :: sze
integer(key_kind), intent(inout) :: key(sze)
real(integral_kind), intent(inout):: value(sze)
integer :: i
integer(cache_map_size_kind) :: n_elements
integer(map_size_kind) :: idx_cache
integer(cache_key_kind) :: cache_key
do i=1,sze
idx_cache = ishft(key(i),map_shift)
call omp_set_lock(map%map(idx_cache)%lock)
n_elements = map%map(idx_cache)%n_elements + 1
! Assert that the map has a proper size
if (n_elements == map%map(idx_cache)%map_size) then
call cache_map_reallocate(map%map(idx_cache), n_elements+ ishft(n_elements,-1))
endif
cache_key = iand(key(i),map_mask)
map%map(idx_cache)%value(n_elements) = value(i)
map%map(idx_cache)%key(n_elements) = cache_key
map%map(idx_cache)%n_elements = n_elements
if (map%map(idx_cache)%sorted.and.n_elements > 1) then
map%map(idx_cache)%sorted = (map%map(idx_cache)%key(n_elements-1) <= cache_key)
map%sorted = map%sorted .and. map%map(idx_cache)%sorted
endif
call omp_unset_lock(map%map(idx_cache)%lock)
enddo
call omp_set_lock(map%lock)
map%n_elements = map%n_elements + sze
call omp_unset_lock(map%lock)
end
subroutine map_get(map, key, value)
use map_module
implicit none
type (map_type), intent(in) :: map
integer(key_kind), intent(in) :: key
real(integral_kind), intent(out) :: value
integer(map_size_kind) :: idx_cache
integer(cache_map_size_kind) :: idx
idx_cache = ishft(key,map_shift)
!DIR$ FORCEINLINE
call cache_map_get_interval(map%map(idx_cache), key, value, 1, map%map(idx_cache)%n_elements,idx)
end
subroutine cache_map_get_interval(map, key, value, ibegin, iend, idx)
use map_module
implicit none
type (cache_map_type), intent(inout) :: map
integer(key_kind), intent(in) :: key
integer(cache_map_size_kind), intent(in) :: ibegin, iend
real(integral_kind), intent(out) :: value
integer(cache_map_size_kind), intent(in) :: idx
call search_key_big_interval(key,map%key, map%n_elements, idx, ibegin, iend)
if (idx > 0) then
value = map%value(idx)
else
value = 0.
endif
end
subroutine map_get_many(map, key, value, sze)
use map_module
implicit none
type (map_type), intent(inout) :: map
integer, intent(in) :: sze
integer(key_kind), intent(in) :: key(sze)
real(integral_kind), intent(out) :: value(sze)
integer :: i
integer(map_size_kind) :: idx_cache, idx_cache_prev
integer(cache_map_size_kind) :: ibegin, iend
integer(cache_map_size_kind), allocatable :: idx(:)
!DIR$ ATTRIBUTES ALIGN : 64 :: idx
allocate(idx(sze))
do i=1,sze
idx_cache = ishft(key(i),map_shift)
iend = map%map(idx_cache)%n_elements
!DIR$ FORCEINLINE
call search_key_big_interval(key(i),map%map(idx_cache)%key, iend, idx(i), 1, iend)
enddo
do i=1,sze
idx_cache = ishft(key(i),map_shift)
if (idx(i) > 0) then
value(i) = map%map(idx_cache)%value(idx(i))
else
value(i) = 0.
endif
enddo
deallocate(idx)
end
subroutine search_key_big(key,X,sze,idx)
use map_module
implicit none
integer(cache_map_size_kind), intent(in) :: sze
integer(key_kind) , intent(in) :: key
integer(cache_key_kind) , intent(in) :: X(sze)
integer(cache_map_size_kind), intent(out) :: idx
call search_key_big_interval(key,X,sze,idx,1,sze)
end
subroutine search_key_big_interval(key,X,sze,idx,ibegin_in,iend_in)
use map_module
implicit none
integer(cache_map_size_kind), intent(in) :: sze
integer(key_kind) , intent(in) :: key
integer(cache_key_kind) , intent(in) :: X(sze)
integer(cache_map_size_kind), intent(in) :: ibegin_in, iend_in
integer(cache_map_size_kind), intent(out) :: idx
integer(cache_map_size_kind) :: istep, ibegin, iend, i
integer(cache_key_kind) :: cache_key
if (sze /= 0) then
continue
else
idx = -1
return
endif
cache_key = iand(key,map_mask)
ibegin = min(ibegin_in,sze)
iend = min(iend_in,sze)
if ((cache_key > X(ibegin)) .and. (cache_key < X(iend))) then
istep = ishft(iend-ibegin,-1)
idx = ibegin + istep
do while (istep > 32)
idx = ibegin + istep
if (cache_key < X(idx)) then
iend = idx
istep = ishft(idx-ibegin,-1)
idx = ibegin + istep
if (cache_key < X(idx)) then
iend = idx
istep = ishft(idx-ibegin,-1)
cycle
else if (cache_key > X(idx)) then
ibegin = idx
istep = ishft(iend-idx,-1)
cycle
else
return
endif
else if (cache_key > X(idx)) then
ibegin = idx
istep = ishft(iend-idx,-1)
idx = idx + istep
if (cache_key < X(idx)) then
iend = idx
istep = ishft(idx-ibegin,-1)
cycle
else if (cache_key > X(idx)) then
ibegin = idx
istep = ishft(iend-idx,-1)
cycle
else
return
endif
else
return
endif
enddo
idx = ibegin
if (min(iend_in,sze) > ibegin+64) then
iend = ibegin+64
!DIR$ VECTOR ALIGNED
do while (cache_key > X(idx))
idx = idx+1
end do
else
!DIR$ VECTOR ALIGNED
do while (cache_key > X(idx))
idx = idx+1
if (idx /= iend) then
cycle
else
exit
endif
end do
endif
if (cache_key /= X(idx)) then
idx = 1-idx
endif
return
else
if (cache_key < X(ibegin)) then
idx = -ibegin
return
endif
if (cache_key > X(iend)) then
idx = -iend
return
endif
if (cache_key == X(ibegin)) then
idx = ibegin
return
endif
if (cache_key == X(iend)) then
idx = iend
return
endif
endif
end