LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-02 11:25:26 +02:00
Code Issues Releases Wiki Activity

Renamed integrals

Browse Source
This commit is contained in:
Anthony Scemama 2018-12-20 16:45:21 +01:00
parent a4344fc41f
commit 36f69cdb43
40 changed files with 910 additions and 927 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
data/module_gitignore
View File

@ -6,6 +6,7 @@ tags
Makefile
ao_basis
ao_one_e_integrals
ao_two_e_integrals
becke_numerical_grid
bitmask
cis
@ -19,7 +20,6 @@ dummy
electrons
ezfio_files
fci
four_idx
generators_cas
generators_full
hartree_fock
@ -28,6 +28,7 @@ iterations
mo_basis
mo_guess
mo_one_e_integrals
mo_two_e_integrals
mpi
mrpt_utils
nuclei

33
src/ao_two_e_integrals/EZFIO.cfg Normal file
View File

@ -0,0 +1,33 @@
[disk_access_ao_integrals]
type: Disk_access
doc: Read/Write |AO| integrals from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml
default: None
[disk_access_ao_integrals_erf]
type: Disk_access
doc: Read/Write |AO| integrals with the long range interaction from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml
default: None
[ao_integrals_threshold]
type: Threshold
doc: If | (pq|rs) | < `ao_integrals_threshold` then (pq|rs) is zero
interface: ezfio,provider,ocaml
default: 1.e-15
ezfio_name: threshold_ao
[do_direct_integrals]
type: logical
doc: Compute integrals on the fly (very slow, only for debugging)
interface: ezfio,provider,ocaml
default: False
ezfio_name: direct
[mu_erf]
type: double precision
doc: cutting of the interaction in the range separated model
interface: ezfio,provider,ocaml
default: 0.5
ezfio_name: mu_erf

2
src/integrals_bielec/NEED → src/ao_two_e_integrals/NEED
View File

@ -1,7 +1,5 @@
ao_one_e_integrals
mo_one_e_integrals
pseudo
bitmask
zmq
ao_basis
mo_basis

0
src/integrals_bielec/gauss_legendre.irp.f → src/ao_two_e_integrals/gauss_legendre.irp.f
View File

0
src/integrals_bielec_erf/ao_bielec_integrals_erf_in_map_slave.irp.f → src/ao_two_e_integrals/integrals_erf_in_map_slave.irp.f
View File

0
src/integrals_bielec/ao_bielec_integrals_in_map_slave.irp.f → src/ao_two_e_integrals/integrals_in_map_slave.irp.f
View File

375
src/ao_two_e_integrals/map_integrals.irp.f Normal file
View File

@ -0,0 +1,375 @@
use map_module
!! AO Map
!! ======
BEGIN_PROVIDER [ type(map_type), ao_integrals_map ]
implicit none
BEGIN_DOC
! AO integrals
END_DOC
integer(key_kind) :: key_max
integer(map_size_kind) :: sze
call bielec_integrals_index(ao_num,ao_num,ao_num,ao_num,key_max)
sze = key_max
call map_init(ao_integrals_map,sze)
print*, 'AO map initialized : ', sze
END_PROVIDER
subroutine bielec_integrals_index(i,j,k,l,i1)
use map_module
implicit none
integer, intent(in) :: i,j,k,l
integer(key_kind), intent(out) :: i1
integer(key_kind) :: p,q,r,s,i2
p = min(i,k)
r = max(i,k)
p = p+shiftr(r*r-r,1)
q = min(j,l)
s = max(j,l)
q = q+shiftr(s*s-s,1)
i1 = min(p,q)
i2 = max(p,q)
i1 = i1+shiftr(i2*i2-i2,1)
end
subroutine bielec_integrals_index_reverse(i,j,k,l,i1)
use map_module
implicit none
integer, intent(out) :: i(8),j(8),k(8),l(8)
integer(key_kind), intent(in) :: i1
integer(key_kind) :: i2,i3
i = 0
i2 = ceiling(0.5d0*(dsqrt(8.d0*dble(i1)+1.d0)-1.d0))
l(1) = ceiling(0.5d0*(dsqrt(8.d0*dble(i2)+1.d0)-1.d0))
i3 = i1 - shiftr(i2*i2-i2,1)
k(1) = ceiling(0.5d0*(dsqrt(8.d0*dble(i3)+1.d0)-1.d0))
j(1) = int(i2 - shiftr(l(1)*l(1)-l(1),1),4)
i(1) = int(i3 - shiftr(k(1)*k(1)-k(1),1),4)
!ijkl
i(2) = i(1) !ilkj
j(2) = l(1)
k(2) = k(1)
l(2) = j(1)
i(3) = k(1) !kjil
j(3) = j(1)
k(3) = i(1)
l(3) = l(1)
i(4) = k(1) !klij
j(4) = l(1)
k(4) = i(1)
l(4) = j(1)
i(5) = j(1) !jilk
j(5) = i(1)
k(5) = l(1)
l(5) = k(1)
i(6) = j(1) !jkli
j(6) = k(1)
k(6) = l(1)
l(6) = i(1)
i(7) = l(1) !lijk
j(7) = i(1)
k(7) = j(1)
l(7) = k(1)
i(8) = l(1) !lkji
j(8) = k(1)
k(8) = j(1)
l(8) = i(1)
integer :: ii, jj
do ii=2,8
do jj=1,ii-1
if ( (i(ii) == i(jj)).and. &
(j(ii) == j(jj)).and. &
(k(ii) == k(jj)).and. &
(l(ii) == l(jj)) ) then
i(ii) = 0
exit
endif
enddo
enddo
do ii=1,8
if (i(ii) /= 0) then
call bielec_integrals_index(i(ii),j(ii),k(ii),l(ii),i2)
if (i1 /= i2) then
print *, i1, i2
print *, i(ii), j(ii), k(ii), l(ii)
stop 'bielec_integrals_index_reverse failed'
endif
endif
enddo
end
BEGIN_PROVIDER [ integer, ao_integrals_cache_min ]
&BEGIN_PROVIDER [ integer, ao_integrals_cache_max ]
implicit none
BEGIN_DOC
! Min and max values of the AOs for which the integrals are in the cache
END_DOC
ao_integrals_cache_min = max(1,ao_num - 63)
ao_integrals_cache_max = ao_num
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_integrals_cache, (0:64*64*64*64) ]
implicit none
BEGIN_DOC
! Cache of AO integrals for fast access
END_DOC
PROVIDE ao_bielec_integrals_in_map
integer :: i,j,k,l,ii
integer(key_kind) :: idx
real(integral_kind) :: integral
!$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral)
do l=ao_integrals_cache_min,ao_integrals_cache_max
do k=ao_integrals_cache_min,ao_integrals_cache_max
do j=ao_integrals_cache_min,ao_integrals_cache_max
do i=ao_integrals_cache_min,ao_integrals_cache_max
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,idx)
!DIR$ FORCEINLINE
call map_get(ao_integrals_map,idx,integral)
ii = l-ao_integrals_cache_min
ii = ior( shiftl(ii,6), k-ao_integrals_cache_min)
ii = ior( shiftl(ii,6), j-ao_integrals_cache_min)
ii = ior( shiftl(ii,6), i-ao_integrals_cache_min)
ao_integrals_cache(ii) = integral
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
double precision function get_ao_bielec_integral(i,j,k,l,map) result(result)
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(key_kind) :: idx
type(map_type), intent(inout) :: map
integer :: ii
real(integral_kind) :: tmp
PROVIDE ao_bielec_integrals_in_map ao_integrals_cache ao_integrals_cache_min
!DIR$ FORCEINLINE
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < ao_integrals_threshold ) then
tmp = 0.d0
else if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < ao_integrals_threshold) then
tmp = 0.d0
else
ii = l-ao_integrals_cache_min
ii = ior(ii, k-ao_integrals_cache_min)
ii = ior(ii, j-ao_integrals_cache_min)
ii = ior(ii, i-ao_integrals_cache_min)
if (iand(ii, -64) /= 0) then
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,idx)
!DIR$ FORCEINLINE
call map_get(map,idx,tmp)
else
ii = l-ao_integrals_cache_min
ii = ior( shiftl(ii,6), k-ao_integrals_cache_min)
ii = ior( shiftl(ii,6), j-ao_integrals_cache_min)
ii = ior( shiftl(ii,6), i-ao_integrals_cache_min)
tmp = ao_integrals_cache(ii)
endif
endif
result = 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(key_kind) :: hash
double precision :: thresh
PROVIDE ao_bielec_integrals_in_map ao_integrals_map
thresh = ao_integrals_threshold
if (ao_overlap_abs(j,l) < thresh) then
out_val = 0.d0
return
endif
double precision :: get_ao_bielec_integral
do i=1,sze
out_val(i) = get_ao_bielec_integral(i,j,k,l,ao_integrals_map)
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(key_kind) :: hash
double precision :: thresh,tmp
PROVIDE ao_bielec_integrals_in_map
thresh = ao_integrals_threshold
non_zero_int = 0
if (ao_overlap_abs(j,l) < thresh) then
out_val = 0.d0
return
endif
non_zero_int = 0
do i=1,sze
integer, external :: ao_l4
double precision, external :: ao_bielec_integral
!DIR$ FORCEINLINE
if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < thresh) then
cycle
endif
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
function get_ao_map_size()
implicit none
integer (map_size_kind) :: get_ao_map_size
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
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(key_kind), 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 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
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) 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
integer*8 :: 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)
98 continue
stop 'Problem reading ao_integrals_map file in work/'
end

285
src/ao_two_e_integrals/map_integrals_erf.irp.f Normal file
View File

@ -0,0 +1,285 @@
use map_module
!! AO Map
!! ======
BEGIN_PROVIDER [ type(map_type), ao_integrals_erf_map ]
implicit none
BEGIN_DOC
! AO integrals
END_DOC
integer(key_kind) :: key_max
integer(map_size_kind) :: sze
call bielec_integrals_index(ao_num,ao_num,ao_num,ao_num,key_max)
sze = key_max
call map_init(ao_integrals_erf_map,sze)
print*, 'AO map initialized : ', sze
END_PROVIDER
BEGIN_PROVIDER [ integer, ao_integrals_erf_cache_min ]
&BEGIN_PROVIDER [ integer, ao_integrals_erf_cache_max ]
implicit none
BEGIN_DOC
! Min and max values of the AOs for which the integrals are in the cache
END_DOC
ao_integrals_erf_cache_min = max(1,ao_num - 63)
ao_integrals_erf_cache_max = ao_num
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_integrals_erf_cache, (0:64*64*64*64) ]
use map_module
implicit none
BEGIN_DOC
! Cache of AO integrals for fast access
END_DOC
PROVIDE ao_bielec_integrals_erf_in_map
integer :: i,j,k,l,ii
integer(key_kind) :: idx
real(integral_kind) :: integral
!$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral)
do l=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max
do k=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max
do j=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max
do i=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,idx)
!DIR$ FORCEINLINE
call map_get(ao_integrals_erf_map,idx,integral)
ii = l-ao_integrals_erf_cache_min
ii = ior( ishft(ii,6), k-ao_integrals_erf_cache_min)
ii = ior( ishft(ii,6), j-ao_integrals_erf_cache_min)
ii = ior( ishft(ii,6), i-ao_integrals_erf_cache_min)
ao_integrals_erf_cache(ii) = integral
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
subroutine insert_into_ao_integrals_erf_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(key_kind), intent(inout) :: buffer_i(n_integrals)
real(integral_kind), intent(inout) :: buffer_values(n_integrals)
call map_append(ao_integrals_erf_map, buffer_i, buffer_values, n_integrals)
end
double precision function get_ao_bielec_integral_erf(i,j,k,l,map) result(result)
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(key_kind) :: idx
type(map_type), intent(inout) :: map
integer :: ii
real(integral_kind) :: tmp
PROVIDE ao_bielec_integrals_erf_in_map ao_integrals_erf_cache ao_integrals_erf_cache_min
!DIR$ FORCEINLINE
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < ao_integrals_threshold ) then
tmp = 0.d0
else if (ao_bielec_integral_erf_schwartz(i,k)*ao_bielec_integral_erf_schwartz(j,l) < ao_integrals_threshold) then
tmp = 0.d0
else
ii = l-ao_integrals_erf_cache_min
ii = ior(ii, k-ao_integrals_erf_cache_min)
ii = ior(ii, j-ao_integrals_erf_cache_min)
ii = ior(ii, i-ao_integrals_erf_cache_min)
if (iand(ii, -64) /= 0) then
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,idx)
!DIR$ FORCEINLINE
call map_get(map,idx,tmp)
tmp = tmp
else
ii = l-ao_integrals_erf_cache_min
ii = ior( ishft(ii,6), k-ao_integrals_erf_cache_min)
ii = ior( ishft(ii,6), j-ao_integrals_erf_cache_min)
ii = ior( ishft(ii,6), i-ao_integrals_erf_cache_min)
tmp = ao_integrals_erf_cache(ii)
endif
endif
result = tmp
end
subroutine get_ao_bielec_integrals_erf(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(key_kind) :: hash
double precision :: thresh
PROVIDE ao_bielec_integrals_erf_in_map ao_integrals_erf_map
thresh = ao_integrals_threshold
if (ao_overlap_abs(j,l) < thresh) then
out_val = 0.d0
return
endif
double precision :: get_ao_bielec_integral_erf
do i=1,sze
out_val(i) = get_ao_bielec_integral_erf(i,j,k,l,ao_integrals_erf_map)
enddo
end
subroutine get_ao_bielec_integrals_erf_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(key_kind) :: hash
double precision :: thresh,tmp
PROVIDE ao_bielec_integrals_erf_in_map
thresh = ao_integrals_threshold
non_zero_int = 0
if (ao_overlap_abs(j,l) < thresh) then
out_val = 0.d0
return
endif
non_zero_int = 0
do i=1,sze
integer, external :: ao_l4
double precision, external :: ao_bielec_integral_erf
!DIR$ FORCEINLINE
if (ao_bielec_integral_erf_schwartz(i,k)*ao_bielec_integral_erf_schwartz(j,l) < thresh) then
cycle
endif
call bielec_integrals_index(i,j,k,l,hash)
call map_get(ao_integrals_erf_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
function get_ao_erf_map_size()
implicit none
integer (map_size_kind) :: get_ao_erf_map_size
BEGIN_DOC
! Returns the number of elements in the AO map
END_DOC
get_ao_erf_map_size = ao_integrals_erf_map % n_elements
end
subroutine clear_ao_erf_map
implicit none
BEGIN_DOC
! Frees the memory of the AO map
END_DOC
call map_deinit(ao_integrals_erf_map)
FREE ao_integrals_erf_map
end
subroutine dump_ao_integrals_erf(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
call ezfio_set_work_empty(.False.)
open(unit=66,file=filename,FORM='unformatted')
write(66) integral_kind, key_kind
write(66) ao_integrals_erf_map%sorted, ao_integrals_erf_map%map_size, &
ao_integrals_erf_map%n_elements
do i=0_8,ao_integrals_erf_map%map_size
write(66) ao_integrals_erf_map%map(i)%sorted, ao_integrals_erf_map%map(i)%map_size,&
ao_integrals_erf_map%map(i)%n_elements
enddo
do i=0_8,ao_integrals_erf_map%map_size
key => ao_integrals_erf_map%map(i)%key
val => ao_integrals_erf_map%map(i)%value
n = ao_integrals_erf_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_erf(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
integer*8 :: n, j
load_ao_integrals_erf = 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_erf_map%sorted, ao_integrals_erf_map%map_size,&
ao_integrals_erf_map%n_elements
do i=0_8, ao_integrals_erf_map%map_size
read(66,err=99,end=99) ao_integrals_erf_map%map(i)%sorted, &
ao_integrals_erf_map%map(i)%map_size, ao_integrals_erf_map%map(i)%n_elements
call cache_map_reallocate(ao_integrals_erf_map%map(i),ao_integrals_erf_map%map(i)%map_size)
enddo
do i=0_8, ao_integrals_erf_map%map_size
key => ao_integrals_erf_map%map(i)%key
val => ao_integrals_erf_map%map(i)%value
n = ao_integrals_erf_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_erf_map)
load_ao_integrals_erf = 0
return
99 continue
call map_deinit(ao_integrals_erf_map)
98 continue
stop 'Problem reading ao_integrals_erf_map file in work/'
end

2
src/integrals_bielec_erf/providers_ao_erf.irp.f → src/ao_two_e_integrals/providers_ao_erf.irp.f
View File

@ -90,7 +90,7 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_erf_in_map ]
if (write_ao_integrals_erf) then
call ezfio_set_work_empty(.False.)
call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_erf',ao_integrals_erf_map)
call ezfio_set_integrals_bielec_erf_disk_access_ao_integrals_erf("Read")
call ezfio_set_ao_two_e_integrals_disk_access_ao_integrals_erf("Read")
endif
END_PROVIDER

55
src/ao_two_e_integrals/read_write.irp.f Normal file
View File

@ -0,0 +1,55 @@
BEGIN_PROVIDER [ logical, read_ao_integrals ]
&BEGIN_PROVIDER [ logical, write_ao_integrals ]
BEGIN_DOC
! Flag to read or write the |AO| integrals
END_DOC
implicit none
if (disk_access_ao_integrals.EQ.'Read') then
read_ao_integrals = .True.
write_ao_integrals = .False.
else if (disk_access_ao_integrals.EQ.'Write') then
read_ao_integrals = .False.
write_ao_integrals = .True.
else if (disk_access_ao_integrals.EQ.'None') then
read_ao_integrals = .False.
write_ao_integrals = .False.
else
print *, 'disk_access_ao_integrals has a wrong type'
stop 1
endif
END_PROVIDER
BEGIN_PROVIDER [ logical, read_ao_integrals_erf ]
&BEGIN_PROVIDER [ logical, write_ao_integrals_erf ]
implicit none
BEGIN_DOC
! Flag to read or write the |AO| erf integrals
END_DOC
if (disk_access_ao_integrals_erf.EQ.'Read') then
read_ao_integrals_erf = .True.
write_ao_integrals_erf = .False.
else if (disk_access_ao_integrals_erf.EQ.'Write') then
read_ao_integrals_erf = .False.
write_ao_integrals_erf = .True.
else if (disk_access_ao_integrals_erf.EQ.'None') then
read_ao_integrals_erf = .False.
write_ao_integrals_erf = .False.
else
print *, 'disk_access_ao_integrals_erf has a wrong type'
stop 1
endif
END_PROVIDER

18
src/ao_two_e_integrals/routines_save_integrals_erf.irp.f Normal file
View File

@ -0,0 +1,18 @@
subroutine save_erf_bi_elec_integrals_ao
implicit none
integer :: i,j,k,l
PROVIDE ao_bielec_integrals_erf_in_map
call ezfio_set_work_empty(.False.)
call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_erf',ao_integrals_erf_map)
call ezfio_set_integrals_bielec_disk_access_ao_integrals("Read")
end
subroutine save_erf_bielec_ints_ao_into_ints_ao
implicit none
integer :: i,j,k,l
PROVIDE ao_bielec_integrals_erf_in_map
call ezfio_set_work_empty(.False.)
call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_erf_map)
call ezfio_set_integrals_bielec_disk_access_ao_integrals("Read")
end

0
src/integrals_bielec/ao_bi_integrals.irp.f → src/ao_two_e_integrals/two_e_integrals.irp.f
View File

2
src/integrals_bielec_erf/ao_bi_integrals_erf.irp.f → src/ao_two_e_integrals/two_e_integrals_erf.irp.f
View File

@ -527,7 +527,7 @@ subroutine integrale_new_erf(I_f,a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z
include 'utils/constants.include.F'
double precision :: p,q
integer :: a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z
integer :: i, n_iter, n_pt, j
integer :: i, n_pt, j
double precision :: I_f, pq_inv, p10_1, p10_2, p01_1, p01_2,rho,pq_inv_2
integer :: ix,iy,iz, jx,jy,jz, sx,sy,sz

2
src/determinants/NEED
View File

@ -1,3 +1,3 @@
mo_basis
mo_one_e_integrals
integrals_bielec
mo_two_e_integrals

8
src/dft_keywords/EZFIO.cfg
View File

@ -16,14 +16,6 @@ doc: Percentage of HF exchange in the DFT model
interface: ezfio,provider,ocaml
default: 0.
[mu_erf]
type: double precision
doc: cutting of the interaction in the range separated model
interface: ezfio,provider,ocaml
default: 0.5
ezfio_name: mu_erf
[density_for_dft]
type: character*(32)
doc: Type of density used for DFT calculation. If WFT it uses the density of the WFT stored in terms of determinants. If input_density it uses the one-body dm stored in data_.../ . If damping_rs_dft it uses the damping density between WFT and input_density

4
src/dft_utils_one_body/NEED
View File

@ -1,4 +1,4 @@
density_for_dft
dft_utils_on_grid
integrals_bielec
integrals_bielec_erf
mo_one_e_integrals
mo_two_e_integrals

4
src/dummy/NEED
View File

@ -1,5 +1,6 @@
ao_basis
ao_one_e_integrals
ao_two_e_integrals
becke_numerical_grid
bitmask
cis
@ -13,14 +14,13 @@ dressing
electrons
ezfio_files
fci
four_idx
generators_cas
generators_full
hartree_fock
integrals_bielec
mo_basis
mo_guess
mo_one_e_integrals
mo_two_e_integrals
mpi
mrpt_utils
nuclei

1
src/four_idx/NEED
View File

@ -1 +0,0 @@
integrals_bielec

6
src/four_idx/README.rst
View File

@ -1,6 +0,0 @@
=======
FourIdx
=======
Four-index transformation.

2
src/hartree_fock/NEED
View File

@ -1,4 +1,4 @@
integrals_bielec
ao_one_e_integrals
ao_two_e_integrals
mo_guess
bitmask

22
src/integrals_bielec/integrals_3_index.irp.f
View File

@ -1,22 +0,0 @@
BEGIN_PROVIDER [double precision, big_array_coulomb_integrals, (mo_tot_num,mo_tot_num, mo_tot_num)]
&BEGIN_PROVIDER [double precision, big_array_exchange_integrals,(mo_tot_num,mo_tot_num, mo_tot_num)]
implicit none
integer :: i,j,k,l
double precision :: get_mo_bielec_integral
double precision :: integral
do k = 1, mo_tot_num
do i = 1, mo_tot_num
do j = 1, mo_tot_num
l = j
integral = get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
big_array_coulomb_integrals(j,i,k) = integral
l = j
integral = get_mo_bielec_integral(i,j,l,k,mo_integrals_map)
big_array_exchange_integrals(j,i,k) = integral
enddo
enddo
enddo
END_PROVIDER

47
src/integrals_bielec/read_write.irp.f
View File

@ -1,47 +0,0 @@
BEGIN_PROVIDER [ logical, read_ao_integrals ]
&BEGIN_PROVIDER [ logical, read_mo_integrals ]
&BEGIN_PROVIDER [ logical, write_ao_integrals ]
&BEGIN_PROVIDER [ logical, write_mo_integrals ]
BEGIN_DOC
! One level of abstraction for disk_access_ao_integrals and disk_access_mo_integrals
END_DOC
implicit none
if (disk_access_ao_integrals.EQ.'Read') then
read_ao_integrals = .True.
write_ao_integrals = .False.
else if (disk_access_ao_integrals.EQ.'Write') then
read_ao_integrals = .False.
write_ao_integrals = .True.
else if (disk_access_ao_integrals.EQ.'None') then
read_ao_integrals = .False.
write_ao_integrals = .False.
else
print *, 'bielec_integrals/disk_access_ao_integrals has a wrong type'
stop 1
endif
if (disk_access_mo_integrals.EQ.'Read') then
read_mo_integrals = .True.
write_mo_integrals = .False.
else if (disk_access_mo_integrals.EQ.'Write') then
read_mo_integrals = .False.
write_mo_integrals = .True.
else if (disk_access_mo_integrals.EQ.'None') then
read_mo_integrals = .False.
write_mo_integrals = .False.
else
print *, 'bielec_integrals/disk_access_mo_integrals has a wrong type'
stop 1
endif
END_PROVIDER

25
src/integrals_bielec_erf/EZFIO.cfg
View File

@ -1,25 +0,0 @@
[disk_access_ao_integrals_erf]
type: Disk_access
doc: Read/Write AO integrals with the long range interaction from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml
default: None
[disk_access_mo_integrals_erf]
type: Disk_access
doc: Read/Write MO integrals with the long range interaction from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml
default: None
[disk_access_mo_integrals_sr]
type: Disk_access
doc: Read/Write MO integrals with the short range interaction from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml
default: None
[disk_access_ao_integrals_sr]
type: Disk_access
doc: Read/Write AO integrals with the short range interaction from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml
default: None

5
src/integrals_bielec_erf/NEED
View File

@ -1,5 +0,0 @@
pseudo
bitmask
zmq
integrals_bielec
dft_keywords

32
src/integrals_bielec_erf/qp_ao_erf_ints.irp.f
View File

@ -1,32 +0,0 @@
program qp_ao_ints
use omp_lib
implicit none
BEGIN_DOC
! Increments a running calculation to compute AO integral_erfs
END_DOC
integer :: i
call switch_qp_run_to_master
zmq_context = f77_zmq_ctx_new ()
! Set the state of the ZMQ
zmq_state = 'ao_integral_erfs'
! Provide everything needed
double precision :: integral_erf, ao_bielec_integral_erf
integral_erf = ao_bielec_integral_erf(1,1,1,1)
character*(64) :: state
call wait_for_state(zmq_state,state)
do while (state /= 'Stopped')
!$OMP PARALLEL DEFAULT(PRIVATE) PRIVATE(i)
i = omp_get_thread_num()
call ao_bielec_integrals_erf_in_map_slave_tcp(i)
!$OMP END PARALLEL
call wait_for_state(zmq_state,state)
enddo
print *, 'Done'
end

47
src/integrals_bielec_erf/read_write.irp.f
View File

@ -1,47 +0,0 @@
BEGIN_PROVIDER [ logical, read_ao_integrals_erf ]
&BEGIN_PROVIDER [ logical, read_mo_integrals_erf ]
&BEGIN_PROVIDER [ logical, write_ao_integrals_erf ]
&BEGIN_PROVIDER [ logical, write_mo_integrals_erf ]
BEGIN_DOC
! One level of abstraction for disk_access_ao_integrals_erf and disk_access_mo_integrals_erf
END_DOC
implicit none
if (disk_access_ao_integrals_erf.EQ.'Read') then
read_ao_integrals_erf = .True.
write_ao_integrals_erf = .False.
else if (disk_access_ao_integrals_erf.EQ.'Write') then
read_ao_integrals_erf = .False.
write_ao_integrals_erf = .True.
else if (disk_access_ao_integrals_erf.EQ.'None') then
read_ao_integrals_erf = .False.
write_ao_integrals_erf = .False.
else
print *, 'bielec_integrals_erf/disk_access_ao_integrals_erf has a wrong type'
stop 1
endif
if (disk_access_mo_integrals_erf.EQ.'Read') then
read_mo_integrals_erf = .True.
write_mo_integrals_erf = .False.
else if (disk_access_mo_integrals_erf.EQ.'Write') then
read_mo_integrals_erf = .False.
write_mo_integrals_erf = .True.
else if (disk_access_mo_integrals_erf.EQ.'None') then
read_mo_integrals_erf = .False.
write_mo_integrals_erf = .False.
else
print *, 'bielec_integrals_erf/disk_access_mo_integrals_erf has a wrong type'
stop 1
endif
END_PROVIDER

21
src/integrals_bielec_erf/write_ints_erf_to_regular_ints.irp.f
View File

@ -1,21 +0,0 @@
program write_integrals_erf_to_regular_integrals
implicit none
BEGIN_DOC
! This program saves the bielec erf integrals into the EZFIO folder but at the regular bielec integrals place.
! Threfore, if the user runs a future calculation with a reading of the integrals, the calculation will be performed with the erf bielec integrals instead of the regular bielec integrals
END_DOC
disk_access_mo_integrals = "None"
touch disk_access_mo_integrals
disk_access_ao_integrals = "None"
touch disk_access_ao_integrals
call routine
end
subroutine routine
implicit none
call save_erf_bielec_ints_mo_into_ints_ao
call save_erf_bielec_ints_mo_into_ints_mo
end

19
src/integrals_bielec/EZFIO.cfg → src/mo_two_e_integrals/EZFIO.cfg
View File

@ -4,19 +4,12 @@ doc: Read/Write |MO| integrals from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml
default: None
[disk_access_ao_integrals]
[disk_access_mo_integrals_erf]
type: Disk_access
doc: Read/Write |AO| integrals from/to disk [ Write | Read | None ]
doc: Read/Write MO integrals with the long range interaction from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml
default: None
[ao_integrals_threshold]
type: Threshold
doc: If | (pq|rs) | < `ao_integrals_threshold` then (pq|rs) is zero
interface: ezfio,provider,ocaml
default: 1.e-15
ezfio_name: threshold_ao
[mo_integrals_threshold]
type: Threshold
doc: If | <ij|kl> | < `mo_integrals_threshold` then <ij|kl> is zero
@ -24,7 +17,6 @@ interface: ezfio,provider,ocaml
default: 1.e-15
ezfio_name: threshold_mo
[no_vvvv_integrals]
type: logical
doc: If `True`, computes all integrals except for the integrals having 4 virtual indices
@ -46,10 +38,3 @@ interface: ezfio,provider,ocaml
default: False
ezfio_name: no_vvv_integrals
[do_direct_integrals]
type: logical
doc: Compute integrals on the fly (very slow, only for debugging)
interface: ezfio,provider,ocaml
default: False
ezfio_name: direct

6
src/mo_two_e_integrals/NEED Normal file
View File

@ -0,0 +1,6 @@
ao_one_e_integrals
ao_two_e_integrals
pseudo
bitmask
zmq
mo_basis

0
src/integrals_bielec/README.rst → src/mo_two_e_integrals/README.rst
View File

0
src/integrals_bielec/core_quantities.irp.f → src/mo_two_e_integrals/core_quantities.irp.f
View File

47
src/integrals_bielec_erf/integrals_3_index_erf.irp.f → src/mo_two_e_integrals/integrals_3_index.irp.f
View File

@ -1,3 +1,26 @@
BEGIN_PROVIDER [double precision, big_array_coulomb_integrals, (mo_tot_num,mo_tot_num, mo_tot_num)]
&BEGIN_PROVIDER [double precision, big_array_exchange_integrals,(mo_tot_num,mo_tot_num, mo_tot_num)]
implicit none
integer :: i,j,k,l
double precision :: get_mo_bielec_integral
double precision :: integral
do k = 1, mo_tot_num
do i = 1, mo_tot_num
do j = 1, mo_tot_num
l = j
integral = get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
big_array_coulomb_integrals(j,i,k) = integral
l = j
integral = get_mo_bielec_integral(i,j,l,k,mo_integrals_map)
big_array_exchange_integrals(j,i,k) = integral
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, big_array_coulomb_integrals_erf, (mo_tot_num,mo_tot_num, mo_tot_num)]
&BEGIN_PROVIDER [double precision, big_array_exchange_integrals_erf,(mo_tot_num,mo_tot_num, mo_tot_num)]
implicit none
@ -21,27 +44,3 @@
END_PROVIDER
!BEGIN_PROVIDER [double precision, big_array_coulomb_integrals_sr, (mo_tot_num,mo_tot_num, mo_tot_num)]
!&BEGIN_PROVIDER [double precision, big_array_exchange_integrals_sr,(mo_tot_num,mo_tot_num, mo_tot_num)]
!implicit none
!integer :: i,j,k,l
!double precision :: get_mo_bielec_integral_sr
!double precision :: integral
!do k = 1, mo_tot_num
! do i = 1, mo_tot_num
! do j = 1, mo_tot_num
! l = j
! integral = get_mo_bielec_integral_sr(i,j,k,l,mo_integrals_sr_map)
! big_array_coulomb_integrals_sr(j,i,k) = integral
! l = j
! integral = get_mo_bielec_integral_sr(i,j,l,k,mo_integrals_sr_map)
! big_array_exchange_integrals_sr(j,i,k) = integral
! enddo
! enddo
!enddo
!ND_PROVIDER

355
src/integrals_bielec/map_integrals.irp.f → src/mo_two_e_integrals/map_integrals.irp.f
View File

@ -1,285 +1,5 @@
use map_module
!! AO Map
!! ======
BEGIN_PROVIDER [ type(map_type), ao_integrals_map ]
implicit none
BEGIN_DOC
! AO integrals
END_DOC
integer(key_kind) :: key_max
integer(map_size_kind) :: sze
call bielec_integrals_index(ao_num,ao_num,ao_num,ao_num,key_max)
sze = key_max
call map_init(ao_integrals_map,sze)
print*, 'AO map initialized : ', sze
END_PROVIDER
subroutine bielec_integrals_index(i,j,k,l,i1)
use map_module
implicit none
integer, intent(in) :: i,j,k,l
integer(key_kind), intent(out) :: i1
integer(key_kind) :: p,q,r,s,i2
p = min(i,k)
r = max(i,k)
p = p+shiftr(r*r-r,1)
q = min(j,l)
s = max(j,l)
q = q+shiftr(s*s-s,1)
i1 = min(p,q)
i2 = max(p,q)
i1 = i1+shiftr(i2*i2-i2,1)
end
subroutine bielec_integrals_index_reverse(i,j,k,l,i1)
use map_module
implicit none
integer, intent(out) :: i(8),j(8),k(8),l(8)
integer(key_kind), intent(in) :: i1
integer(key_kind) :: i2,i3
i = 0
i2 = ceiling(0.5d0*(dsqrt(8.d0*dble(i1)+1.d0)-1.d0))
l(1) = ceiling(0.5d0*(dsqrt(8.d0*dble(i2)+1.d0)-1.d0))
i3 = i1 - shiftr(i2*i2-i2,1)
k(1) = ceiling(0.5d0*(dsqrt(8.d0*dble(i3)+1.d0)-1.d0))
j(1) = int(i2 - shiftr(l(1)*l(1)-l(1),1),4)
i(1) = int(i3 - shiftr(k(1)*k(1)-k(1),1),4)
!ijkl
i(2) = i(1) !ilkj
j(2) = l(1)
k(2) = k(1)
l(2) = j(1)
i(3) = k(1) !kjil
j(3) = j(1)
k(3) = i(1)
l(3) = l(1)
i(4) = k(1) !klij
j(4) = l(1)
k(4) = i(1)
l(4) = j(1)
i(5) = j(1) !jilk
j(5) = i(1)
k(5) = l(1)
l(5) = k(1)
i(6) = j(1) !jkli
j(6) = k(1)
k(6) = l(1)
l(6) = i(1)
i(7) = l(1) !lijk
j(7) = i(1)
k(7) = j(1)
l(7) = k(1)
i(8) = l(1) !lkji
j(8) = k(1)
k(8) = j(1)
l(8) = i(1)
integer :: ii, jj
do ii=2,8
do jj=1,ii-1
if ( (i(ii) == i(jj)).and. &
(j(ii) == j(jj)).and. &
(k(ii) == k(jj)).and. &
(l(ii) == l(jj)) ) then
i(ii) = 0
exit
endif
enddo
enddo
do ii=1,8
if (i(ii) /= 0) then
call bielec_integrals_index(i(ii),j(ii),k(ii),l(ii),i2)
if (i1 /= i2) then
print *, i1, i2
print *, i(ii), j(ii), k(ii), l(ii)
stop 'bielec_integrals_index_reverse failed'
endif
endif
enddo
end
BEGIN_PROVIDER [ integer, ao_integrals_cache_min ]
&BEGIN_PROVIDER [ integer, ao_integrals_cache_max ]
implicit none
BEGIN_DOC
! Min and max values of the AOs for which the integrals are in the cache
END_DOC
ao_integrals_cache_min = max(1,ao_num - 63)
ao_integrals_cache_max = ao_num
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_integrals_cache, (0:64*64*64*64) ]
implicit none
BEGIN_DOC
! Cache of AO integrals for fast access
END_DOC
PROVIDE ao_bielec_integrals_in_map
integer :: i,j,k,l,ii
integer(key_kind) :: idx
real(integral_kind) :: integral
!$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral)
do l=ao_integrals_cache_min,ao_integrals_cache_max
do k=ao_integrals_cache_min,ao_integrals_cache_max
do j=ao_integrals_cache_min,ao_integrals_cache_max
do i=ao_integrals_cache_min,ao_integrals_cache_max
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,idx)
!DIR$ FORCEINLINE
call map_get(ao_integrals_map,idx,integral)
ii = l-ao_integrals_cache_min
ii = ior( shiftl(ii,6), k-ao_integrals_cache_min)
ii = ior( shiftl(ii,6), j-ao_integrals_cache_min)
ii = ior( shiftl(ii,6), i-ao_integrals_cache_min)
ao_integrals_cache(ii) = integral
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
double precision function get_ao_bielec_integral(i,j,k,l,map) result(result)
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(key_kind) :: idx
type(map_type), intent(inout) :: map
integer :: ii
real(integral_kind) :: tmp
PROVIDE ao_bielec_integrals_in_map ao_integrals_cache ao_integrals_cache_min
!DIR$ FORCEINLINE
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < ao_integrals_threshold ) then
tmp = 0.d0
else if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < ao_integrals_threshold) then
tmp = 0.d0
else
ii = l-ao_integrals_cache_min
ii = ior(ii, k-ao_integrals_cache_min)
ii = ior(ii, j-ao_integrals_cache_min)
ii = ior(ii, i-ao_integrals_cache_min)
if (iand(ii, -64) /= 0) then
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,idx)
!DIR$ FORCEINLINE
call map_get(map,idx,tmp)
else
ii = l-ao_integrals_cache_min
ii = ior( shiftl(ii,6), k-ao_integrals_cache_min)
ii = ior( shiftl(ii,6), j-ao_integrals_cache_min)
ii = ior( shiftl(ii,6), i-ao_integrals_cache_min)
tmp = ao_integrals_cache(ii)
endif
endif
result = 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(key_kind) :: hash
double precision :: thresh
PROVIDE ao_bielec_integrals_in_map ao_integrals_map
thresh = ao_integrals_threshold
if (ao_overlap_abs(j,l) < thresh) then
out_val = 0.d0
return
endif
double precision :: get_ao_bielec_integral
do i=1,sze
out_val(i) = get_ao_bielec_integral(i,j,k,l,ao_integrals_map)
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(key_kind) :: hash
double precision :: thresh,tmp
PROVIDE ao_bielec_integrals_in_map
thresh = ao_integrals_threshold
non_zero_int = 0
if (ao_overlap_abs(j,l) < thresh) then
out_val = 0.d0
return
endif
non_zero_int = 0
do i=1,sze
integer, external :: ao_l4
double precision, external :: ao_bielec_integral
!DIR$ FORCEINLINE
if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < thresh) then
cycle
endif
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
function get_ao_map_size()
implicit none
integer (map_size_kind) :: get_ao_map_size
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
!! ======
@ -296,20 +16,6 @@ BEGIN_PROVIDER [ type(map_type), mo_integrals_map ]
print*, 'MO map initialized: ', sze
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(key_kind), 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
@ -608,13 +314,12 @@ integer*8 function get_mo_map_size()
get_mo_map_size = mo_integrals_map % n_elements
end
BEGIN_TEMPLATE
subroutine dump_$ao_integrals(filename)
subroutine dump_mo_integrals(filename)
use map_module
implicit none
BEGIN_DOC
! Save to disk the $ao integrals
! Save to disk the |MO| integrals
END_DOC
character*(*), intent(in) :: filename
integer(cache_key_kind), pointer :: key(:)
@ -626,16 +331,16 @@ subroutine dump_$ao_integrals(filename)
call ezfio_set_work_empty(.False.)
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
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,$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
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)
@ -643,10 +348,10 @@ subroutine dump_$ao_integrals(filename)
end
integer function load_$ao_integrals(filename)
integer function load_mo_integrals(filename)
implicit none
BEGIN_DOC
! Read from disk the $ao integrals
! Read from disk the |MO| integrals
END_DOC
character*(*), intent(in) :: filename
integer*8 :: i
@ -654,7 +359,7 @@ integer function load_$ao_integrals(filename)
real(integral_kind), pointer :: val(:)
integer :: iknd, kknd
integer*8 :: n, j
load_$ao_integrals = 1
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
@ -665,30 +370,26 @@ integer function load_$ao_integrals(filename)
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)
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, $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
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($ao_integrals_map)
load_$ao_integrals = 0
call map_sort(mo_integrals_map)
load_mo_integrals = 0
return
99 continue
call map_deinit($ao_integrals_map)
call map_deinit(mo_integrals_map)
98 continue
stop 'Problem reading $ao_integrals_map file in work/'
stop 'Problem reading mo_integrals_map file in work/'
end
SUBST [ ao_integrals_map, ao_integrals, ao_num ]
ao_integrals_map ; ao_integrals ; ao_num ;;
mo_integrals_map ; mo_integrals ; mo_tot_num ;;
END_TEMPLATE

328
src/integrals_bielec_erf/map_integrals_erf.irp.f → src/mo_two_e_integrals/map_integrals_erf.irp.f
View File

@ -1,284 +1,7 @@
use map_module
!! AO Map
!! ======
BEGIN_PROVIDER [ type(map_type), ao_integrals_erf_map ]
implicit none
BEGIN_DOC
! AO integrals
END_DOC
integer(key_kind) :: key_max
integer(map_size_kind) :: sze
call bielec_integrals_index(ao_num,ao_num,ao_num,ao_num,key_max)
sze = key_max
call map_init(ao_integrals_erf_map,sze)
print*, 'AO map initialized : ', sze
END_PROVIDER
BEGIN_PROVIDER [ integer, ao_integrals_erf_cache_min ]
&BEGIN_PROVIDER [ integer, ao_integrals_erf_cache_max ]
implicit none
BEGIN_DOC
! Min and max values of the AOs for which the integrals are in the cache
END_DOC
ao_integrals_erf_cache_min = max(1,ao_num - 63)
ao_integrals_erf_cache_max = ao_num
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_integrals_erf_cache, (0:64*64*64*64) ]
use map_module
implicit none
BEGIN_DOC
! Cache of AO integrals for fast access
END_DOC
PROVIDE ao_bielec_integrals_erf_in_map
integer :: i,j,k,l,ii
integer(key_kind) :: idx
real(integral_kind) :: integral
!$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral)
do l=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max
do k=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max
do j=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max
do i=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,idx)
!DIR$ FORCEINLINE
call map_get(ao_integrals_erf_map,idx,integral)
ii = l-ao_integrals_erf_cache_min
ii = ior( ishft(ii,6), k-ao_integrals_erf_cache_min)
ii = ior( ishft(ii,6), j-ao_integrals_erf_cache_min)
ii = ior( ishft(ii,6), i-ao_integrals_erf_cache_min)
ao_integrals_erf_cache(ii) = integral
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
double precision function get_ao_bielec_integral_erf(i,j,k,l,map) result(result)
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(key_kind) :: idx
type(map_type), intent(inout) :: map
integer :: ii
real(integral_kind) :: tmp
PROVIDE ao_bielec_integrals_erf_in_map ao_integrals_erf_cache ao_integrals_erf_cache_min
!DIR$ FORCEINLINE
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < ao_integrals_threshold ) then
tmp = 0.d0
else if (ao_bielec_integral_erf_schwartz(i,k)*ao_bielec_integral_erf_schwartz(j,l) < ao_integrals_threshold) then
tmp = 0.d0
else
ii = l-ao_integrals_erf_cache_min
ii = ior(ii, k-ao_integrals_erf_cache_min)
ii = ior(ii, j-ao_integrals_erf_cache_min)
ii = ior(ii, i-ao_integrals_erf_cache_min)
if (iand(ii, -64) /= 0) then
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,idx)
!DIR$ FORCEINLINE
call map_get(map,idx,tmp)
tmp = tmp
else
ii = l-ao_integrals_erf_cache_min
ii = ior( ishft(ii,6), k-ao_integrals_erf_cache_min)
ii = ior( ishft(ii,6), j-ao_integrals_erf_cache_min)
ii = ior( ishft(ii,6), i-ao_integrals_erf_cache_min)
tmp = ao_integrals_erf_cache(ii)
endif
endif
result = tmp
end
subroutine get_ao_bielec_integrals_erf(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(key_kind) :: hash
double precision :: thresh
PROVIDE ao_bielec_integrals_erf_in_map ao_integrals_erf_map
thresh = ao_integrals_threshold
if (ao_overlap_abs(j,l) < thresh) then
out_val = 0.d0
return
endif
double precision :: get_ao_bielec_integral_erf
do i=1,sze
out_val(i) = get_ao_bielec_integral_erf(i,j,k,l,ao_integrals_erf_map)
enddo
end
subroutine get_ao_bielec_integrals_erf_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(key_kind) :: hash
double precision :: thresh,tmp
PROVIDE ao_bielec_integrals_erf_in_map
thresh = ao_integrals_threshold
non_zero_int = 0
if (ao_overlap_abs(j,l) < thresh) then
out_val = 0.d0
return
endif
non_zero_int = 0
do i=1,sze
integer, external :: ao_l4
double precision, external :: ao_bielec_integral_erf
!DIR$ FORCEINLINE
if (ao_bielec_integral_erf_schwartz(i,k)*ao_bielec_integral_erf_schwartz(j,l) < thresh) then
cycle
endif
call bielec_integrals_index(i,j,k,l,hash)
call map_get(ao_integrals_erf_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
function get_ao_erf_map_size()
implicit none
integer (map_size_kind) :: get_ao_erf_map_size
BEGIN_DOC
! Returns the number of elements in the AO map
END_DOC
get_ao_erf_map_size = ao_integrals_erf_map % n_elements
end
subroutine clear_ao_erf_map
implicit none
BEGIN_DOC
! Frees the memory of the AO map
END_DOC
call map_deinit(ao_integrals_erf_map)
FREE ao_integrals_erf_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
call ezfio_set_work_empty(.False.)
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
IRP_IF COARRAY
subroutine communicate_$ao_integrals()
use map_module
implicit none
BEGIN_DOC
! Communicate the $ao integrals with co-array
END_DOC
integer(cache_key_kind), pointer :: key(:)
real(integral_kind), pointer :: val(:)
integer*8 :: i,j, k, nmax
integer*8, save :: n[*]
integer :: copy_n
real(integral_kind), allocatable :: buffer_val(:)[:]
integer(cache_key_kind), allocatable :: buffer_key(:)[:]
real(integral_kind), allocatable :: copy_val(:)
integer(key_kind), allocatable :: copy_key(:)
n = 0_8
do i=0_8,$ao_integrals_map%map_size
n = max(n,$ao_integrals_map%map(i)%n_elements)
enddo
sync all
nmax = 0_8
do j=1,num_images()
nmax = max(nmax,n[j])
enddo
allocate( buffer_key(nmax)[*], buffer_val(nmax)[*])
allocate( copy_key(nmax), copy_val(nmax))
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
do j=1,n
buffer_key(j) = key(j)
buffer_val(j) = val(j)
enddo
sync all
do j=1,num_images()
if (j /= this_image()) then
copy_n = n[j]
do k=1,copy_n
copy_val(k) = buffer_val(k)[j]
copy_key(k) = buffer_key(k)[j]
copy_key(k) = copy_key(k)+ishft(i,-map_shift)
enddo
call map_append($ao_integrals_map, copy_key, copy_val, copy_n )
endif
enddo
sync all
enddo
deallocate( buffer_key, buffer_val, copy_val, copy_key)
end
IRP_ENDIF
integer function load_$ao_integrals(filename)
integer function load_mo_integrals_erf(filename)
implicit none
BEGIN_DOC
! Read from disk the $ao integrals
@ -289,7 +12,7 @@ integer function load_$ao_integrals(filename)
real(integral_kind), pointer :: val(:)
integer :: iknd, kknd
integer*8 :: n, j
load_$ao_integrals = 1
load_mo_integrals_erf = 1
open(unit=66,file=filename,FORM='unformatted',STATUS='UNKNOWN')
read(66,err=98,end=98) iknd, kknd
if (iknd /= integral_kind) then
@ -300,34 +23,29 @@ integer function load_$ao_integrals(filename)
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)
read(66,err=98,end=98) mo_integrals_erf_map%sorted, mo_integrals_erf_map%map_size,&
mo_integrals_erf_map%n_elements
do i=0_8, mo_integrals_erf_map%map_size
read(66,err=99,end=99) mo_integrals_erf_map%map(i)%sorted, &
mo_integrals_erf_map%map(i)%map_size, mo_integrals_erf_map%map(i)%n_elements
call cache_map_reallocate(mo_integrals_erf_map%map(i),mo_integrals_erf_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
do i=0_8, mo_integrals_erf_map%map_size
key => mo_integrals_erf_map%map(i)%key
val => mo_integrals_erf_map%map(i)%value
n = mo_integrals_erf_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
call map_sort(mo_integrals_erf_map)
load_mo_integrals_erf = 0
return
99 continue
call map_deinit($ao_integrals_map)
call map_deinit(mo_integrals_erf_map)
98 continue
stop 'Problem reading $ao_integrals_map file in work/'
stop 'Problem reading mo_integrals_erf_map file in work/'
end
SUBST [ ao_integrals_map, ao_integrals, ao_num ]
ao_integrals_erf_map ; ao_integrals_erf ; ao_num ;;
mo_integrals_erf_map ; mo_integrals_erf ; mo_tot_num;;
END_TEMPLATE
@ -344,20 +62,6 @@ BEGIN_PROVIDER [ type(map_type), mo_integrals_erf_map ]
print*, 'MO ERF map initialized'
END_PROVIDER
subroutine insert_into_ao_integrals_erf_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(key_kind), intent(inout) :: buffer_i(n_integrals)
real(integral_kind), intent(inout) :: buffer_values(n_integrals)
call map_append(ao_integrals_erf_map, buffer_i, buffer_values, n_integrals)
end
subroutine insert_into_mo_integrals_erf_map(n_integrals, &
buffer_i, buffer_values, thr)
use map_module

2
src/integrals_bielec/mo_bi_integrals.irp.f → src/mo_two_e_integrals/mo_bi_integrals.irp.f
View File

@ -152,7 +152,7 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ]
if (write_mo_integrals.and.mpi_master) then
call ezfio_set_work_empty(.False.)
call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_map)
call ezfio_set_integrals_bielec_disk_access_mo_integrals("Read")
call ezfio_set_mo_two_e_integrals_disk_access_mo_integrals_erf('Read')
endif
END_PROVIDER

0
src/integrals_bielec_erf/mo_bi_integrals_erf.irp.f → src/mo_two_e_integrals/mo_bi_integrals_erf.irp.f
View File

55
src/mo_two_e_integrals/read_write.irp.f Normal file
View File

@ -0,0 +1,55 @@
BEGIN_PROVIDER [ logical, read_mo_integrals ]
&BEGIN_PROVIDER [ logical, write_mo_integrals ]
BEGIN_DOC
! Flag to read or write the |MO| integrals
END_DOC
implicit none
if (disk_access_mo_integrals.EQ.'Read') then
read_mo_integrals = .True.
write_mo_integrals = .False.
else if (disk_access_mo_integrals.EQ.'Write') then
read_mo_integrals = .False.
write_mo_integrals = .True.
else if (disk_access_mo_integrals.EQ.'None') then
read_mo_integrals = .False.
write_mo_integrals = .False.
else
print *, 'disk_access_mo_integrals has a wrong type'
stop 1
endif
END_PROVIDER
BEGIN_PROVIDER [ logical, read_mo_integrals_erf ]
&BEGIN_PROVIDER [ logical, write_mo_integrals_erf ]
implicit none
BEGIN_DOC
! Flag to read or write the |MO| erf integrals
END_DOC
if (disk_access_mo_integrals_erf.EQ.'Read') then
read_mo_integrals_erf = .True.
write_mo_integrals_erf = .False.
else if (disk_access_mo_integrals_erf.EQ.'Write') then
read_mo_integrals_erf = .False.
write_mo_integrals_erf = .True.
else if (disk_access_mo_integrals_erf.EQ.'None') then
read_mo_integrals_erf = .False.
write_mo_integrals_erf = .False.
else
print *, 'disk_access_mo_integrals_erf has a wrong type'
stop 1
endif
END_PROVIDER

18
src/integrals_bielec_erf/routines_save_integrals_erf.irp.f → src/mo_two_e_integrals/routines_save_integrals_erf.irp.f
View File

@ -7,15 +7,6 @@ subroutine save_erf_bi_elec_integrals_mo
call ezfio_set_integrals_bielec_disk_access_mo_integrals("Read")
end
subroutine save_erf_bi_elec_integrals_ao
implicit none
integer :: i,j,k,l
PROVIDE ao_bielec_integrals_erf_in_map
call ezfio_set_work_empty(.False.)
call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_erf',ao_integrals_erf_map)
call ezfio_set_integrals_bielec_disk_access_ao_integrals("Read")
end
subroutine save_erf_bielec_ints_mo_into_ints_mo
implicit none
integer :: i,j,k,l
@ -25,12 +16,3 @@ subroutine save_erf_bielec_ints_mo_into_ints_mo
call ezfio_set_integrals_bielec_disk_access_mo_integrals("Read")
end
subroutine save_erf_bielec_ints_mo_into_ints_ao
implicit none
integer :: i,j,k,l
PROVIDE ao_bielec_integrals_erf_in_map
call ezfio_set_work_empty(.False.)
call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_erf_map)
call ezfio_set_integrals_bielec_disk_access_ao_integrals("Read")
end

0
src/four_idx/four_idx_transform.irp.f → src/tools/four_idx_transform.irp.f
View File

6
src/integrals_bielec_erf/write_integrals_erf.irp.f → src/tools/write_integrals_erf.irp.f
View File

@ -1,11 +1,11 @@
program write_integrals
implicit none
BEGIN_DOC
! This program saves the bielec erf integrals into the EZFIO folder
! Saves the bielec erf integrals into the EZFIO
END_DOC
disk_access_mo_integrals = "None"
disk_access_mo_integrals = 'None'
touch disk_access_mo_integrals
disk_access_ao_integrals = "None"
disk_access_ao_integrals = 'None'
touch disk_access_ao_integrals
call routine