mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-12-23 21:03:49 +01:00
working on two e ints
This commit is contained in:
Kevin Gasperich
parent
949ff3ce3a
commit
4e93390632
@ -75,7 +75,16 @@ BEGIN_PROVIDER [ double precision, ao_overlap_imag, (ao_num, ao_num) ]
|
||||
BEGIN_DOC
|
||||
! Imaginary part of the overlap
|
||||
END_DOC
|
||||
ao_overlap_imag = 0.d0
|
||||
if (read_ao_integrals_overlap) then
|
||||
call ezfio_get_ao_one_e_ints_ao_integrals_overlap_imag(ao_overlap_imag(1:ao_num, 1:ao_num))
|
||||
print *, 'AO overlap integrals read from disk'
|
||||
else
|
||||
ao_overlap_imag = 0.d0
|
||||
endif
|
||||
if (write_ao_integrals_overlap) then
|
||||
call ezfio_set_ao_one_e_ints_ao_integrals_overlap_imag(ao_overlap_imag(1:ao_num, 1:ao_num))
|
||||
print *, 'AO overlap integrals written to disk'
|
||||
endif
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ complex*16, ao_overlap_complex, (ao_num, ao_num) ]
|
||||
|
||||
@ -160,13 +160,13 @@ BEGIN_PROVIDER [double precision, ao_kinetic_integrals_imag, (ao_num,ao_num)]
|
||||
integer :: i,j,k,l
|
||||
|
||||
if (read_ao_integrals_kinetic) then
|
||||
call ezfio_get_ao_one_e_ints_ao_integrals_kinetic(ao_kinetic_integrals_imag)
|
||||
call ezfio_get_ao_one_e_ints_ao_integrals_kinetic_imag(ao_kinetic_integrals_imag)
|
||||
print *, 'AO kinetic integrals read from disk'
|
||||
else
|
||||
print *, irp_here, ': Not yet implemented'
|
||||
endif
|
||||
if (write_ao_integrals_kinetic) then
|
||||
call ezfio_set_ao_one_e_ints_ao_integrals_kinetic(ao_kinetic_integrals_imag)
|
||||
call ezfio_set_ao_one_e_ints_ao_integrals_kinetic_imag(ao_kinetic_integrals_imag)
|
||||
print *, 'AO kinetic integrals written to disk'
|
||||
endif
|
||||
END_PROVIDER
|
||||
|
||||
@ -27,6 +27,36 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals, (ao_num,ao_num)]
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_imag, (ao_num, ao_num) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Imaginary part of the pseudo_integrals
|
||||
END_DOC
|
||||
if (read_ao_integrals_pseudo) then
|
||||
call ezfio_get_ao_one_e_ints_ao_integrals_pseudo_imag(ao_pseudo_integrals_imag(1:ao_num, 1:ao_num))
|
||||
print *, 'AO pseudo_integrals integrals read from disk'
|
||||
else
|
||||
ao_pseudo_integrals_imag = 0.d0
|
||||
endif
|
||||
if (write_ao_integrals_pseudo) then
|
||||
call ezfio_set_ao_one_e_ints_ao_integrals_pseudo_imag(ao_pseudo_integrals_imag(1:ao_num, 1:ao_num))
|
||||
print *, 'AO pseudo_integrals integrals written to disk'
|
||||
endif
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ complex*16, ao_pseudo_integrals_complex, (ao_num, ao_num) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Overlap for complex AOs
|
||||
END_DOC
|
||||
integer :: i,j
|
||||
do j=1,ao_num
|
||||
do i=1,ao_num
|
||||
ao_pseudo_integrals_complex(i,j) = dcmplx( ao_pseudo_integrals(i,j), ao_pseudo_integrals_imag(i,j) )
|
||||
enddo
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
|
||||
@ -4,6 +4,7 @@ use map_module
|
||||
!! ======
|
||||
|
||||
BEGIN_PROVIDER [ type(map_type), ao_integrals_map ]
|
||||
&BEGIN_PROVIDER [ type(map_type), ao_integrals_map_2 ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! AO integrals
|
||||
@ -11,22 +12,17 @@ BEGIN_PROVIDER [ type(map_type), ao_integrals_map ]
|
||||
integer(key_kind) :: key_max
|
||||
integer(map_size_kind) :: sze
|
||||
call two_e_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
|
||||
|
||||
BEGIN_PROVIDER [ type(map_type), ao_integrals_map_periodic ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! AO integrals
|
||||
END_DOC
|
||||
integer(key_kind) :: key_max
|
||||
integer(map_size_kind) :: sze
|
||||
call two_e_integrals_index_2fold(ao_num,ao_num,ao_num,ao_num,key_max)
|
||||
sze = key_max
|
||||
call map_init(ao_integrals_map_periodic,sze)
|
||||
print*, 'complex AO map initialized : ', sze
|
||||
if (is_periodic) then
|
||||
sze = key_max*2
|
||||
call map_init(ao_integrals_map,sze)
|
||||
call map_init(ao_integrals_map_2,sze)
|
||||
print*, 'AO maps initialized (complex): ', 2*sze
|
||||
else
|
||||
sze = key_max
|
||||
call map_init(ao_integrals_map,sze)
|
||||
call map_init(ao_integrals_map_2,1_map_size_kind)
|
||||
print*, 'AO map initialized : ', sze
|
||||
endif
|
||||
END_PROVIDER
|
||||
|
||||
subroutine two_e_integrals_index(i,j,k,l,i1)
|
||||
@ -34,7 +30,7 @@ subroutine two_e_integrals_index(i,j,k,l,i1)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Gives a unique index for i,j,k,l using permtuation symmetry.
|
||||
! i <-> k, j <-> l, and (i,k) <-> (j,l) for non-periodic systems
|
||||
! i <-> k, j <-> l, and (i,k) <-> (j,l)
|
||||
END_DOC
|
||||
integer, intent(in) :: i,j,k,l
|
||||
integer(key_kind), intent(out) :: i1
|
||||
@ -50,6 +46,28 @@ subroutine two_e_integrals_index(i,j,k,l,i1)
|
||||
i1 = i1+shiftr(i2*i2-i2,1)
|
||||
end
|
||||
|
||||
subroutine two_e_integrals_index_periodic(i,j,k,l,i1,p,q)
|
||||
use map_module
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Gives a unique index for i,j,k,l using permtuation symmetry.
|
||||
! i <-> k, j <-> l, and (i,k) <-> (j,l)
|
||||
END_DOC
|
||||
integer, intent(in) :: i,j,k,l
|
||||
integer(key_kind), intent(out) :: i1
|
||||
integer(key_kind) :: r,s,i2
|
||||
integer(key_kind),intent(out) :: p,q
|
||||
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 two_e_integrals_index_reverse(i,j,k,l,i1)
|
||||
@ -143,152 +161,6 @@ end
|
||||
|
||||
|
||||
|
||||
subroutine ao_idx2_sq(i,j,ij)
|
||||
implicit none
|
||||
integer, intent(in) :: i,j
|
||||
integer, intent(out) :: ij
|
||||
if (i<j) then
|
||||
ij=(j-1)*(j-1)+2*i-mod(j+1,2)
|
||||
else if (i>j) then
|
||||
ij=(i-1)*(i-1)+2*j-mod(i,2)
|
||||
else
|
||||
ij=i*i
|
||||
endif
|
||||
end
|
||||
|
||||
subroutine idx2_tri_int(i,j,ij)
|
||||
implicit none
|
||||
integer, intent(in) :: i,j
|
||||
integer, intent(out) :: ij
|
||||
integer :: p,q
|
||||
p = max(i,j)
|
||||
q = min(i,j)
|
||||
ij = q+ishft(p*p-p,-1)
|
||||
end
|
||||
|
||||
subroutine ao_idx2_tri_key(i,j,ij)
|
||||
use map_module
|
||||
implicit none
|
||||
integer, intent(in) :: i,j
|
||||
integer(key_kind), intent(out) :: ij
|
||||
integer(key_kind) :: p,q
|
||||
p = max(i,j)
|
||||
q = min(i,j)
|
||||
ij = q+ishft(p*p-p,-1)
|
||||
end
|
||||
|
||||
subroutine two_e_integrals_index_2fold(i,j,k,l,i1)
|
||||
use map_module
|
||||
implicit none
|
||||
integer, intent(in) :: i,j,k,l
|
||||
integer(key_kind), intent(out) :: i1
|
||||
integer :: ik,jl
|
||||
|
||||
call ao_idx2_sq(i,k,ik)
|
||||
call ao_idx2_sq(j,l,jl)
|
||||
call ao_idx2_tri_key(ik,jl,i1)
|
||||
end
|
||||
|
||||
subroutine ao_idx2_sq_rev(i,k,ik)
|
||||
BEGIN_DOC
|
||||
! reverse square compound index
|
||||
END_DOC
|
||||
! p = ceiling(dsqrt(dble(ik)))
|
||||
! q = ceiling(0.5d0*(dble(ik)-dble((p-1)*(p-1))))
|
||||
! if (mod(ik,2)==0) then
|
||||
! k=p
|
||||
! i=q
|
||||
! else
|
||||
! i=p
|
||||
! k=q
|
||||
! endif
|
||||
integer, intent(in) :: ik
|
||||
integer, intent(out) :: i,k
|
||||
integer :: pq(0:1),i1,i2
|
||||
pq(0) = ceiling(dsqrt(dble(ik)))
|
||||
pq(1) = ceiling(0.5d0*(dble(ik)-dble((pq(0)-1)*(pq(0)-1))))
|
||||
i1=mod(ik,2)
|
||||
i2=mod(ik+1,2)
|
||||
|
||||
k=pq(i1)
|
||||
i=pq(i2)
|
||||
end
|
||||
|
||||
subroutine ao_idx2_tri_rev_key(i,k,ik)
|
||||
use map_module
|
||||
BEGIN_DOC
|
||||
!return i<=k
|
||||
END_DOC
|
||||
integer(key_kind), intent(in) :: ik
|
||||
integer, intent(out) :: i,k
|
||||
integer(key_kind) :: tmp_k
|
||||
k = ceiling(0.5d0*(dsqrt(8.d0*dble(ik)+1.d0)-1.d0))
|
||||
tmp_k = k
|
||||
i = int(ik - ishft(tmp_k*tmp_k-tmp_k,-1))
|
||||
end
|
||||
|
||||
subroutine idx2_tri_rev_int(i,k,ik)
|
||||
BEGIN_DOC
|
||||
!return i<=k
|
||||
END_DOC
|
||||
integer, intent(in) :: ik
|
||||
integer, intent(out) :: i,k
|
||||
k = ceiling(0.5d0*(dsqrt(8.d0*dble(ik)+1.d0)-1.d0))
|
||||
i = int(ik - ishft(k*k-k,-1))
|
||||
end
|
||||
|
||||
subroutine two_e_integrals_index_reverse_2fold(i,j,k,l,i1)
|
||||
use map_module
|
||||
implicit none
|
||||
integer, intent(out) :: i(2),j(2),k(2),l(2)
|
||||
integer(key_kind), intent(in) :: i1
|
||||
integer(key_kind) :: i0
|
||||
integer :: i2,i3
|
||||
i = 0
|
||||
call ao_idx2_tri_rev_key(i3,i2,i1)
|
||||
|
||||
call ao_idx2_sq_rev(j(1),l(1),i2)
|
||||
call ao_idx2_sq_rev(i(1),k(1),i3)
|
||||
|
||||
!ijkl
|
||||
i(2) = j(1) !jilk
|
||||
j(2) = i(1)
|
||||
k(2) = l(1)
|
||||
l(2) = k(1)
|
||||
|
||||
! i(3) = k(1) !klij complex conjugate
|
||||
! j(3) = l(1)
|
||||
! k(3) = i(1)
|
||||
! l(3) = j(1)
|
||||
!
|
||||
! i(4) = l(1) !lkji complex conjugate
|
||||
! j(4) = k(1)
|
||||
! k(4) = j(1)
|
||||
! l(4) = i(1)
|
||||
|
||||
integer :: ii
|
||||
if ( (i(1)==i(2)).and. &
|
||||
(j(1)==j(2)).and. &
|
||||
(k(1)==k(2)).and. &
|
||||
(l(1)==l(2)) ) then
|
||||
i(2) = 0
|
||||
endif
|
||||
! This has been tested with up to 1000 AOs, and all the reverse indices are
|
||||
! correct ! We can remove the test
|
||||
! do ii=1,2
|
||||
! if (i(ii) /= 0) then
|
||||
! call two_e_integrals_index_2fold(i(ii),j(ii),k(ii),l(ii),i0)
|
||||
! if (i1 /= i0) then
|
||||
! print *, i1, i0
|
||||
! print *, i(ii), j(ii), k(ii), l(ii)
|
||||
! stop 'two_e_integrals_index_reverse_2fold failed'
|
||||
! endif
|
||||
! endif
|
||||
! enddo
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ integer, ao_integrals_cache_min ]
|
||||
&BEGIN_PROVIDER [ integer, ao_integrals_cache_max ]
|
||||
@ -384,32 +256,47 @@ BEGIN_PROVIDER [ complex*16, ao_integrals_cache_periodic, (0:64*64*64*64) ]
|
||||
real(integral_kind) :: tmp_re, tmp_im
|
||||
integer(key_kind) :: idx_re,idx_im
|
||||
complex(integral_kind) :: integral
|
||||
integer(key_kind) :: p,q,r,s,ik,jl
|
||||
logical :: ilek, jlel, iklejl
|
||||
|
||||
|
||||
!$OMP PARALLEL DO PRIVATE (i,j,k,l,idx1,idx2,tmp_re,tmp_im,idx_re,idx_im,ii,integral)
|
||||
!$OMP PARALLEL DO PRIVATE (ilek,jlel,p,q,r,s, ik,jl,iklejl, &
|
||||
!$OMP i,j,k,l,idx1,idx2,tmp_re,tmp_im,idx_re,idx_im,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 two_e_integrals_index_2fold(i,j,k,l,idx1)
|
||||
!DIR$ FORCEINLINE
|
||||
call two_e_integrals_index_2fold(k,l,i,j,idx2)
|
||||
idx_re = min(idx1,idx2)
|
||||
idx_im = max(idx1,idx2)
|
||||
!DIR$ FORCEINLINE
|
||||
call map_get(ao_integrals_map,idx_re,tmp_re)
|
||||
if (idx_re /= idx_im) then
|
||||
call map_get(ao_integrals_map,idx_im,tmp_im)
|
||||
if (idx1 < idx2) then
|
||||
call two_e_integrals_index(i,j,k,l,idx1)
|
||||
ilek = (i.le.k)
|
||||
jlel = (j.le.l)
|
||||
idx1 = 2*idx1 - 1
|
||||
if (ilek.eqv.jlel) then !map1
|
||||
!TODO: merge these calls using map_get_2
|
||||
call map_get(ao_integrals_map,idx1,tmp_re)
|
||||
call map_get(ao_integrals_map,idx1+1,tmp_im)
|
||||
if (ilek) then
|
||||
integral = dcmplx(tmp_re,tmp_im)
|
||||
else
|
||||
integral = dcmplx(tmp_re,-tmp_im)
|
||||
endif
|
||||
else
|
||||
tmp_im = 0.d0
|
||||
integral = dcmplx(tmp_re,tmp_im)
|
||||
endif
|
||||
else !map2
|
||||
!TODO: merge these calls using map_get_2
|
||||
call map_get(ao_integrals_map_2,idx1,tmp_re)
|
||||
call map_get(ao_integrals_map_2,idx1+1,tmp_im)
|
||||
p = min(i,k)
|
||||
r = max(i,k)
|
||||
ik = p+shiftr(r*r-r,1)
|
||||
q = min(j,l)
|
||||
s = max(j,l)
|
||||
jl = q+shiftr(s*s-s,1)
|
||||
iklejl = (ik.le.jl)
|
||||
if (ilek.eqv.iklejl) then
|
||||
integral = dcmplx(tmp_re,tmp_im)
|
||||
else
|
||||
integral = dcmplx(tmp_re,-tmp_im)
|
||||
endif
|
||||
endif
|
||||
|
||||
ii = l-ao_integrals_cache_min
|
||||
ii = ior( shiftl(ii,6), k-ao_integrals_cache_min)
|
||||
@ -424,8 +311,53 @@ BEGIN_PROVIDER [ complex*16, ao_integrals_cache_periodic, (0:64*64*64*64) ]
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
subroutine ao_two_e_integral_periodic_map_idx_sign(i,j,k,l,use_map1,idx,sign)
|
||||
use map_module
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! get position of periodic AO integral <ij|kl>
|
||||
! use_map1: true if integral is in first ao map, false if integral is in second ao map
|
||||
! idx: position of real part of integral in map (imag part is at idx+1)
|
||||
! sign: sign of imaginary part
|
||||
END_DOC
|
||||
integer, intent(in) :: i,j,k,l
|
||||
integer(key_kind), intent(out) :: idx
|
||||
logical, intent(out) :: use_map1
|
||||
double precision, intent(out) :: sign
|
||||
integer(key_kind) :: p,q,r,s,ik,jl
|
||||
logical :: ilek, jlel, iklejl, ikeqjl
|
||||
! i.le.k, j.le.l, tri(i,k).le.tri(j,l)
|
||||
!DIR$ FORCEINLINE
|
||||
call two_e_integrals_index_periodic(i,j,k,l,idx,ik,jl)
|
||||
ilek = (i.le.k)
|
||||
jlel = (j.le.l)
|
||||
idx = 2*idx-1
|
||||
ikeqjl = (ik.eq.jl)
|
||||
if (ilek.eqv.jlel) then !map1
|
||||
use_map1=.True.
|
||||
if (ikeqjl) then
|
||||
sign=0.d0
|
||||
else if (ilek) then
|
||||
sign=1.d0
|
||||
else
|
||||
sign=-1.d0
|
||||
endif
|
||||
else !map2
|
||||
use_map1=.False.
|
||||
if (ikeqjl) then
|
||||
sign=0.d0
|
||||
else
|
||||
iklejl = (ik.le.jl)
|
||||
if (ilek.eqv.iklejl) then
|
||||
sign=1.d0
|
||||
else
|
||||
sign=-1.d0
|
||||
endif
|
||||
endif
|
||||
endif
|
||||
end
|
||||
|
||||
complex*16 function get_ao_two_e_integral_periodic(i,j,k,l,map) result(result)
|
||||
complex*16 function get_ao_two_e_integral_periodic_simple(i,j,k,l,map,map2) result(result)
|
||||
use map_module
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
@ -435,10 +367,66 @@ complex*16 function get_ao_two_e_integral_periodic(i,j,k,l,map) result(result)
|
||||
integer(key_kind) :: idx1,idx2
|
||||
real(integral_kind) :: tmp_re, tmp_im
|
||||
integer(key_kind) :: idx_re,idx_im
|
||||
type(map_type), intent(inout) :: map
|
||||
type(map_type), intent(inout) :: map,map2
|
||||
integer :: ii
|
||||
complex(integral_kind) :: tmp
|
||||
PROVIDE ao_two_e_integrals_in_map_periodic ao_integrals_cache_periodic ao_integrals_cache_min
|
||||
integer(key_kind) :: p,q,r,s,ik,jl
|
||||
logical :: ilek, jlel, iklejl
|
||||
! a.le.c, b.le.d, tri(a,c).le.tri(b,d)
|
||||
PROVIDE ao_two_e_integrals_in_map
|
||||
!DIR$ FORCEINLINE
|
||||
call two_e_integrals_index(i,j,k,l,idx1)
|
||||
ilek = (i.le.k)
|
||||
jlel = (j.le.l)
|
||||
idx1 = idx1*2-1
|
||||
if (ilek.eqv.jlel) then !map1
|
||||
!TODO: merge these calls using map_get_2
|
||||
call map_get(map,idx1,tmp_re)
|
||||
call map_get(map,idx1+1,tmp_im)
|
||||
if (ilek) then
|
||||
tmp = dcmplx(tmp_re,tmp_im)
|
||||
else
|
||||
tmp = dcmplx(tmp_re,-tmp_im)
|
||||
endif
|
||||
else !map2
|
||||
!TODO: merge these calls using map_get_2
|
||||
call map_get(map2,idx1,tmp_re)
|
||||
call map_get(map2,idx1+1,tmp_im)
|
||||
p = min(i,k)
|
||||
r = max(i,k)
|
||||
ik = p+shiftr(r*r-r,1)
|
||||
q = min(j,l)
|
||||
s = max(j,l)
|
||||
jl = q+shiftr(s*s-s,1)
|
||||
iklejl = (ik.le.jl)
|
||||
if (ilek.eqv.iklejl) then
|
||||
tmp = dcmplx(tmp_re,tmp_im)
|
||||
else
|
||||
tmp = dcmplx(tmp_re,-tmp_im)
|
||||
endif
|
||||
endif
|
||||
result = tmp
|
||||
end
|
||||
|
||||
|
||||
complex*16 function get_ao_two_e_integral_periodic(i,j,k,l,map,map2) 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) :: idx1,idx2
|
||||
real(integral_kind) :: tmp_re, tmp_im
|
||||
integer(key_kind) :: idx_re,idx_im
|
||||
type(map_type), intent(inout) :: map,map2
|
||||
integer :: ii
|
||||
complex(integral_kind) :: tmp
|
||||
complex(integral_kind) :: get_ao_two_e_integral_periodic_simple
|
||||
integer(key_kind) :: p,q,r,s,ik,jl
|
||||
logical :: ilek, jlel, iklejl
|
||||
! a.le.c, b.le.d, tri(a,c).le.tri(b,d)
|
||||
PROVIDE ao_two_e_integrals_in_map ao_integrals_cache_periodic ao_integrals_cache_min
|
||||
!DIR$ FORCEINLINE
|
||||
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < ao_integrals_threshold ) then
|
||||
tmp = (0.d0,0.d0)
|
||||
@ -450,25 +438,7 @@ complex*16 function get_ao_two_e_integral_periodic(i,j,k,l,map) result(result)
|
||||
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 two_e_integrals_index_2fold(i,j,k,l,idx1)
|
||||
!DIR$ FORCEINLINE
|
||||
call two_e_integrals_index_2fold(k,l,i,j,idx2)
|
||||
idx_re = min(idx1,idx2)
|
||||
idx_im = max(idx1,idx2)
|
||||
!DIR$ FORCEINLINE
|
||||
call map_get(ao_integrals_map,idx_re,tmp_re)
|
||||
if (idx_re /= idx_im) then
|
||||
call map_get(ao_integrals_map,idx_im,tmp_im)
|
||||
if (idx1 < idx2) then
|
||||
tmp = dcmplx(tmp_re,tmp_im)
|
||||
else
|
||||
tmp = dcmplx(tmp_re,-tmp_im)
|
||||
endif
|
||||
else
|
||||
tmp_im = 0.d0
|
||||
tmp = dcmplx(tmp_re,tmp_im)
|
||||
endif
|
||||
tmp = get_ao_two_e_integral_periodic_simple(i,j,k,l,map,map2)
|
||||
else
|
||||
ii = l-ao_integrals_cache_min
|
||||
ii = ior( shiftl(ii,6), k-ao_integrals_cache_min)
|
||||
@ -528,14 +498,14 @@ subroutine get_ao_two_e_integrals_periodic(j,k,l,sze,out_val)
|
||||
PROVIDE ao_two_e_integrals_in_map ao_integrals_map
|
||||
thresh = ao_integrals_threshold
|
||||
|
||||
if (ao_overlap_abs_periodic(j,l) < thresh) then
|
||||
if (ao_overlap_abs(j,l) < thresh) then
|
||||
out_val = (0.d0,0.d0)
|
||||
return
|
||||
endif
|
||||
|
||||
complex*16 :: get_ao_two_e_integral_periodic
|
||||
do i=1,sze
|
||||
out_val(i) = get_ao_two_e_integral_periodic(i,j,k,l,ao_integrals_map)
|
||||
out_val(i) = get_ao_two_e_integral_periodic(i,j,k,l,ao_integrals_map,ao_integrals_map_2)
|
||||
enddo
|
||||
|
||||
end
|
||||
@ -554,6 +524,10 @@ subroutine get_ao_two_e_integrals_non_zero(j,k,l,sze,out_val,out_val_index,non_z
|
||||
integer :: i
|
||||
integer(key_kind) :: hash
|
||||
double precision :: thresh,tmp
|
||||
if(is_periodic) then
|
||||
print*,'not implemented for periodic:',irp_here
|
||||
stop -1
|
||||
endif
|
||||
PROVIDE ao_two_e_integrals_in_map
|
||||
thresh = ao_integrals_threshold
|
||||
|
||||
@ -598,6 +572,10 @@ subroutine get_ao_two_e_integrals_non_zero_jl(j,l,thresh,sze_max,sze,out_val,out
|
||||
integer(key_kind) :: hash
|
||||
double precision :: tmp
|
||||
|
||||
if(is_periodic) then
|
||||
print*,'not implemented for periodic:',irp_here
|
||||
stop -1
|
||||
endif
|
||||
PROVIDE ao_two_e_integrals_in_map
|
||||
non_zero_int = 0
|
||||
if (ao_overlap_abs(j,l) < thresh) then
|
||||
@ -644,6 +622,10 @@ subroutine get_ao_two_e_integrals_non_zero_jl_from_list(j,l,thresh,list,n_list,s
|
||||
integer(key_kind) :: hash
|
||||
double precision :: tmp
|
||||
|
||||
if(is_periodic) then
|
||||
print*,'not implemented for periodic:',irp_here
|
||||
stop -1
|
||||
endif
|
||||
PROVIDE ao_two_e_integrals_in_map
|
||||
non_zero_int = 0
|
||||
if (ao_overlap_abs(j,l) < thresh) then
|
||||
@ -682,7 +664,7 @@ function 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
|
||||
get_ao_map_size = ao_integrals_map % n_elements + ao_integrals_map_2 % n_elements
|
||||
end
|
||||
|
||||
subroutine clear_ao_map
|
||||
@ -692,6 +674,9 @@ subroutine clear_ao_map
|
||||
END_DOC
|
||||
call map_deinit(ao_integrals_map)
|
||||
FREE ao_integrals_map
|
||||
call map_deinit(ao_integrals_map_2)
|
||||
FREE ao_integrals_map_2
|
||||
|
||||
end
|
||||
|
||||
|
||||
@ -709,82 +694,18 @@ subroutine insert_into_ao_integrals_map(n_integrals,buffer_i, buffer_values)
|
||||
call map_append(ao_integrals_map, buffer_i, buffer_values, n_integrals)
|
||||
end
|
||||
|
||||
subroutine insert_into_ao_integrals_map_2(n_integrals,buffer_i, buffer_values)
|
||||
use map_module
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Create new entry into AO map
|
||||
END_DOC
|
||||
|
||||
!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, 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_2, buffer_i, buffer_values, n_integrals)
|
||||
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
|
||||
!
|
||||
|
||||
@ -348,77 +348,91 @@ BEGIN_PROVIDER [ logical, ao_two_e_integrals_in_map ]
|
||||
integer :: kk, m, j1, i1, lmax
|
||||
character*(64) :: fmt
|
||||
|
||||
integral = ao_two_e_integral(1,1,1,1)
|
||||
|
||||
double precision :: map_mb
|
||||
PROVIDE read_ao_two_e_integrals io_ao_two_e_integrals
|
||||
if (read_ao_two_e_integrals) then
|
||||
print*,'Reading the AO integrals'
|
||||
call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map)
|
||||
print*, 'AO integrals provided'
|
||||
ao_two_e_integrals_in_map = .True.
|
||||
return
|
||||
endif
|
||||
|
||||
print*, 'Providing the AO integrals'
|
||||
call wall_time(wall_0)
|
||||
call wall_time(wall_1)
|
||||
call cpu_time(cpu_1)
|
||||
|
||||
integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_socket_pull
|
||||
call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'ao_integrals')
|
||||
|
||||
character(len=:), allocatable :: task
|
||||
allocate(character(len=ao_num*12) :: task)
|
||||
write(fmt,*) '(', ao_num, '(I5,X,I5,''|''))'
|
||||
do l=1,ao_num
|
||||
write(task,fmt) (i,l, i=1,l)
|
||||
integer, external :: add_task_to_taskserver
|
||||
if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task)) == -1) then
|
||||
stop 'Unable to add task to server'
|
||||
if (is_periodic) then
|
||||
if (read_ao_two_e_integrals) then
|
||||
print*,'Reading the AO integrals (periodic)'
|
||||
call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints_periodic_1',ao_integrals_map)
|
||||
call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints_periodic_2',ao_integrals_map_2)
|
||||
print*, 'AO integrals provided (periodic)'
|
||||
ao_two_e_integrals_in_map = .True.
|
||||
return
|
||||
else
|
||||
print*,'calculation of periodic AOs not implemented'
|
||||
stop -1
|
||||
endif
|
||||
enddo
|
||||
deallocate(task)
|
||||
|
||||
integer, external :: zmq_set_running
|
||||
if (zmq_set_running(zmq_to_qp_run_socket) == -1) then
|
||||
print *, irp_here, ': Failed in zmq_set_running'
|
||||
endif
|
||||
else
|
||||
if (read_ao_two_e_integrals) then
|
||||
print*,'Reading the AO integrals'
|
||||
call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map)
|
||||
print*, 'AO integrals provided'
|
||||
ao_two_e_integrals_in_map = .True.
|
||||
return
|
||||
endif
|
||||
|
||||
PROVIDE nproc
|
||||
!$OMP PARALLEL DEFAULT(shared) private(i) num_threads(nproc+1)
|
||||
i = omp_get_thread_num()
|
||||
if (i==0) then
|
||||
call ao_two_e_integrals_in_map_collector(zmq_socket_pull)
|
||||
else
|
||||
call ao_two_e_integrals_in_map_slave_inproc(i)
|
||||
integral = ao_two_e_integral(1,1,1,1)
|
||||
print*, 'Providing the AO integrals'
|
||||
call wall_time(wall_0)
|
||||
call wall_time(wall_1)
|
||||
call cpu_time(cpu_1)
|
||||
|
||||
integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_socket_pull
|
||||
call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'ao_integrals')
|
||||
|
||||
character(len=:), allocatable :: task
|
||||
allocate(character(len=ao_num*12) :: task)
|
||||
write(fmt,*) '(', ao_num, '(I5,X,I5,''|''))'
|
||||
do l=1,ao_num
|
||||
write(task,fmt) (i,l, i=1,l)
|
||||
integer, external :: add_task_to_taskserver
|
||||
if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task)) == -1) then
|
||||
stop 'Unable to add task to server'
|
||||
endif
|
||||
!$OMP END PARALLEL
|
||||
enddo
|
||||
deallocate(task)
|
||||
|
||||
call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'ao_integrals')
|
||||
integer, external :: zmq_set_running
|
||||
if (zmq_set_running(zmq_to_qp_run_socket) == -1) then
|
||||
print *, irp_here, ': Failed in zmq_set_running'
|
||||
endif
|
||||
|
||||
PROVIDE nproc
|
||||
!$OMP PARALLEL DEFAULT(shared) private(i) num_threads(nproc+1)
|
||||
i = omp_get_thread_num()
|
||||
if (i==0) then
|
||||
call ao_two_e_integrals_in_map_collector(zmq_socket_pull)
|
||||
else
|
||||
call ao_two_e_integrals_in_map_slave_inproc(i)
|
||||
endif
|
||||
!$OMP END PARALLEL
|
||||
|
||||
call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'ao_integrals')
|
||||
|
||||
|
||||
print*, 'Sorting the map'
|
||||
call map_sort(ao_integrals_map)
|
||||
call cpu_time(cpu_2)
|
||||
call wall_time(wall_2)
|
||||
integer(map_size_kind) :: get_ao_map_size, ao_map_size
|
||||
ao_map_size = get_ao_map_size()
|
||||
print*, 'Sorting the map'
|
||||
call map_sort(ao_integrals_map)
|
||||
call cpu_time(cpu_2)
|
||||
call wall_time(wall_2)
|
||||
integer(map_size_kind) :: get_ao_map_size, ao_map_size
|
||||
ao_map_size = get_ao_map_size()
|
||||
|
||||
print*, 'AO integrals provided:'
|
||||
print*, ' Size of AO map : ', map_mb(ao_integrals_map) ,'MB'
|
||||
print*, ' Number of AO integrals :', ao_map_size
|
||||
print*, ' cpu time :',cpu_2 - cpu_1, 's'
|
||||
print*, ' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1+tiny(1.d0)), ' )'
|
||||
print*, 'AO integrals provided:'
|
||||
print*, ' Size of AO map : ', map_mb(ao_integrals_map) ,'MB'
|
||||
print*, ' Number of AO integrals :', ao_map_size
|
||||
print*, ' cpu time :',cpu_2 - cpu_1, 's'
|
||||
print*, ' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1+tiny(1.d0)), ' )'
|
||||
|
||||
ao_two_e_integrals_in_map = .True.
|
||||
ao_two_e_integrals_in_map = .True.
|
||||
|
||||
if (write_ao_two_e_integrals.and.mpi_master) then
|
||||
call ezfio_set_work_empty(.False.)
|
||||
call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map)
|
||||
call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read')
|
||||
if (write_ao_two_e_integrals.and.mpi_master) then
|
||||
call ezfio_set_work_empty(.False.)
|
||||
call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map)
|
||||
call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read')
|
||||
endif
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, ao_two_e_integral_schwartz,(ao_num,ao_num) ]
|
||||
|
||||
@ -531,13 +531,30 @@ subroutine map_get(map, key, value)
|
||||
real(integral_kind), intent(out) :: value
|
||||
integer(map_size_kind) :: idx_cache
|
||||
integer(cache_map_size_kind) :: idx
|
||||
|
||||
idx=1
|
||||
! index in tha pointers array
|
||||
idx_cache = shiftr(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 map_get_2(map, key, value1, value2)
|
||||
use map_module
|
||||
implicit none
|
||||
type (map_type), intent(inout) :: map
|
||||
integer(key_kind), intent(in) :: key
|
||||
real(integral_kind), intent(out) :: value1, value2
|
||||
integer(map_size_kind) :: idx_cache
|
||||
integer(cache_map_size_kind) :: idx
|
||||
|
||||
idx=1
|
||||
! index in tha pointers array
|
||||
idx_cache = shiftr(key,map_shift)
|
||||
!DIR$ FORCEINLINE
|
||||
call cache_map_get_interval(map%map(idx_cache), key, value1, 1, map%map(idx_cache)%n_elements,idx)
|
||||
call cache_map_get_interval(map%map(idx_cache), key+1, value2, idx+1, idx+2, idx)
|
||||
end
|
||||
|
||||
subroutine cache_map_get_interval(map, key, value, ibegin, iend, idx)
|
||||
use map_module
|
||||
implicit none
|
||||
|
||||
212
src/utils_periodic/export_integrals_ao_periodic.irp.f
Normal file
212
src/utils_periodic/export_integrals_ao_periodic.irp.f
Normal file
@ -0,0 +1,212 @@
|
||||
program print_integrals
|
||||
call run
|
||||
end
|
||||
|
||||
subroutine run
|
||||
use map_module
|
||||
implicit none
|
||||
|
||||
integer :: iunit
|
||||
integer :: getunitandopen
|
||||
|
||||
integer ::i,j,k,l
|
||||
double precision :: integral
|
||||
double precision, allocatable :: A(:,:), B(:,:)
|
||||
double precision :: tmp_re, tmp_im
|
||||
|
||||
integer :: n_integrals_1, n_integrals_2
|
||||
integer(key_kind), allocatable :: buffer_i_1(:), buffer_i_2(:)
|
||||
real(integral_kind), allocatable :: buffer_values_1(:), buffer_values_2(:)
|
||||
logical :: use_map1
|
||||
integer(key_kind) :: idx_tmp
|
||||
double precision :: sign
|
||||
|
||||
|
||||
|
||||
provide ao_two_e_integrals_in_map
|
||||
allocate (A(ao_num,ao_num), B(ao_num,ao_num) )
|
||||
|
||||
A(1,1) = huge(1.d0)
|
||||
iunit = getunitandopen('E.qp','r')
|
||||
read (iunit,*,end=9) A(1,1)
|
||||
9 continue
|
||||
close(iunit)
|
||||
if (A(1,1) /= huge(1.d0)) then
|
||||
! call ezfio_set_nuclei_nuclear_repulsion(A(1,1))
|
||||
! call ezfio_set_nuclei_io_nuclear_repulsion("Read")
|
||||
print*, nuclear_repulsion,A(1,1)
|
||||
endif
|
||||
|
||||
A = 0.d0
|
||||
B = 0.d0
|
||||
! iunit = getunitandopen('T.qp','r')
|
||||
! do
|
||||
! read (iunit,*,end=10) i,j, tmp_re, tmp_im
|
||||
! A(i,j) = tmp_re
|
||||
! B(i,j) = tmp_im
|
||||
! print*,ao_kinetic_integrals(i,j),A(i,j)
|
||||
! print*,ao_kinetic_integrals_imag(i,j),B(i,j)
|
||||
! if (i.ne.j) then
|
||||
! A(j,i) = tmp_re
|
||||
! B(j,i) = -tmp_im
|
||||
! print*,ao_kinetic_integrals(j,i),A(j,i)
|
||||
! print*,ao_kinetic_integrals_imag(j,i),B(j,i)
|
||||
! endif
|
||||
! enddo
|
||||
! 10 continue
|
||||
! close(iunit)
|
||||
! call ezfio_set_ao_one_e_ints_ao_integrals_kinetic(A(1:ao_num, 1:ao_num))
|
||||
! call ezfio_set_ao_one_e_ints_ao_integrals_kinetic_imag(B(1:ao_num, 1:ao_num))
|
||||
! call ezfio_set_ao_one_e_ints_io_ao_integrals_kinetic("Read")
|
||||
|
||||
A = 0.d0
|
||||
B = 0.d0
|
||||
! iunit = getunitandopen('S.qp','r')
|
||||
! do
|
||||
! read (iunit,*,end=11) i,j, tmp_re, tmp_im
|
||||
! A(i,j) = tmp_re
|
||||
! B(i,j) = tmp_im
|
||||
! print*,real(ao_overlap_complex(i,j)),A(i,j)
|
||||
! print*,imag(ao_overlap_complex(i,j)),B(i,j)
|
||||
! print*,ao_overlap_imag(i,j),B(i,j)
|
||||
! if (i.ne.j) then
|
||||
! A(j,i) = tmp_re
|
||||
! B(j,i) = -tmp_im
|
||||
! print*,real(ao_overlap_complex(j,i)),A(j,i)
|
||||
! print*,imag(ao_overlap_complex(j,i)),B(j,i)
|
||||
! print*,ao_overlap_imag(j,i),B(j,i)
|
||||
! endif
|
||||
! enddo
|
||||
! 11 continue
|
||||
! close(iunit)
|
||||
! call ezfio_set_ao_one_e_ints_ao_integrals_overlap(A(1:ao_num, 1:ao_num))
|
||||
! call ezfio_set_ao_one_e_ints_ao_integrals_overlap_imag(B(1:ao_num, 1:ao_num))
|
||||
! call ezfio_set_ao_one_e_ints_io_ao_integrals_overlap("Read")
|
||||
|
||||
A = 0.d0
|
||||
B = 0.d0
|
||||
! iunit = getunitandopen('P.qp','r')
|
||||
! do
|
||||
! read (iunit,*,end=14) i,j, tmp_re, tmp_im
|
||||
! A(i,j) = tmp_re
|
||||
! B(i,j) = tmp_im
|
||||
! print*,ao_pseudo_integrals(i,j),A(i,j)
|
||||
! print*,ao_pseudo_integrals_imag(i,j),B(i,j)
|
||||
! ! print*,real(ao_integrals_pseudo(i,j)),A(i,j)
|
||||
! ! print*,imag(ao_integrals_pseudo(i,j)),B(i,j)
|
||||
! if (i.ne.j) then
|
||||
! A(j,i) = tmp_re
|
||||
! B(j,i) = -tmp_im
|
||||
! print*,ao_pseudo_integrals(j,i),A(j,i)
|
||||
! print*,ao_pseudo_integrals_imag(j,i),B(j,i)
|
||||
! ! print*,real(ao_integrals_pseudo(j,i)),A(j,i)
|
||||
! ! print*,imag(ao_integrals_pseudo(j,i)),B(j,i)
|
||||
! endif
|
||||
! enddo
|
||||
! 14 continue
|
||||
! close(iunit)
|
||||
! call ezfio_set_ao_one_e_ints_ao_integrals_pseudo(A(1:ao_num,1:ao_num))
|
||||
! call ezfio_set_ao_one_e_ints_ao_integrals_pseudo_imag(B(1:ao_num,1:ao_num))
|
||||
! call ezfio_set_ao_one_e_ints_io_ao_integrals_pseudo("Read")
|
||||
|
||||
A = 0.d0
|
||||
B = 0.d0
|
||||
! iunit = getunitandopen('V.qp','r')
|
||||
! do
|
||||
! read (iunit,*,end=12) i,j, tmp_re, tmp_im
|
||||
! A(i,j) = tmp_re
|
||||
! B(i,j) = tmp_im
|
||||
! print*,ao_integrals_n_e(i,j),A(i,j)
|
||||
! print*,ao_integrals_n_e_imag(i,j),B(i,j)
|
||||
! if (i.ne.j) then
|
||||
! A(j,i) = tmp_re
|
||||
! B(j,i) = -tmp_im
|
||||
! print*,ao_integrals_n_e(j,i),A(j,i)
|
||||
! print*,ao_integrals_n_e_imag(j,i),B(j,i)
|
||||
! endif
|
||||
! enddo
|
||||
! 12 continue
|
||||
! close(iunit)
|
||||
! call ezfio_set_ao_one_e_ints_ao_integrals_n_e(A(1:ao_num, 1:ao_num))
|
||||
! call ezfio_set_ao_one_e_ints_ao_integrals_n_e_imag(B(1:ao_num, 1:ao_num))
|
||||
! call ezfio_set_ao_one_e_ints_io_ao_integrals_n_e("Read")
|
||||
complex*16 :: int2e_tmp1,int2e_tmp2,get_ao_two_e_integral_periodic_simple,get_ao_two_e_integral_periodic
|
||||
double precision :: tmp3,tmp4,tmp5,tmp6
|
||||
allocate(buffer_i_1(ao_num**3), buffer_values_1(ao_num**3))
|
||||
allocate(buffer_i_2(ao_num**3), buffer_values_2(ao_num**3))
|
||||
iunit = getunitandopen('W.qp','r')
|
||||
n_integrals_1=0
|
||||
n_integrals_2=0
|
||||
buffer_values_1 = 0.d0
|
||||
buffer_values_2 = 0.d0
|
||||
do
|
||||
read (iunit,*,end=13) i,j,k,l, tmp_re, tmp_im
|
||||
int2e_tmp1 = get_ao_two_e_integral_periodic_simple(i,j,k,l,ao_integrals_map,ao_integrals_map_2)
|
||||
int2e_tmp2 = get_ao_two_e_integral_periodic(i,j,k,l,ao_integrals_map,ao_integrals_map_2)
|
||||
print'(4(I4),3(E15.7))',i,j,k,l,tmp_re,real(int2e_tmp1),real(int2e_tmp2)
|
||||
print'(4(I4),3(E15.7))',i,j,k,l,tmp_im,imag(int2e_tmp1),imag(int2e_tmp2)
|
||||
call ao_two_e_integral_periodic_map_idx_sign(i,j,k,l,use_map1,idx_tmp,sign)
|
||||
print*,use_map1,idx_tmp,sign
|
||||
call map_get(ao_integrals_map,idx_tmp,tmp3)
|
||||
call map_get(ao_integrals_map_2,idx_tmp,tmp4)
|
||||
call map_get(ao_integrals_map,idx_tmp+1,tmp5)
|
||||
call map_get(ao_integrals_map_2,idx_tmp+1,tmp6)
|
||||
print*,tmp3,tmp4
|
||||
print*,tmp5,tmp6
|
||||
! if (use_map1) then
|
||||
! n_integrals_1 += 1
|
||||
! buffer_i_1(n_integrals_1-1)=idx_tmp
|
||||
! buffer_values_1(n_integrals_1-1)=tmp_re
|
||||
! if (sign.ne.0.d0) then
|
||||
! n_integrals_1 += 1
|
||||
! buffer_i_1(n_integrals_2)=idx_tmp+1
|
||||
! buffer_values_1(n_integrals_1)=tmp_im*sign
|
||||
! endif
|
||||
! if (n_integrals_1 >= size(buffer_i_1)-1) then
|
||||
!! call insert_into_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1)
|
||||
! n_integrals_1 = 0
|
||||
! endif
|
||||
! else
|
||||
! n_integrals_2 += 1
|
||||
! buffer_i_2(n_integrals_2-1)=idx_tmp
|
||||
! buffer_values_2(n_integrals_2-1)=tmp_re
|
||||
! if (sign.ne.0.d0) then
|
||||
! n_integrals_2 += 1
|
||||
! buffer_i_2(n_integrals_2)=idx_tmp+1
|
||||
! buffer_values_2(n_integrals_2)=tmp_im*sign
|
||||
! endif
|
||||
! if (n_integrals_2 >= size(buffer_i_2)-1) then
|
||||
!! call insert_into_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2)
|
||||
! n_integrals_2 = 0
|
||||
! endif
|
||||
! endif
|
||||
enddo
|
||||
13 continue
|
||||
close(iunit)
|
||||
|
||||
! if (n_integrals_1 > 0) then
|
||||
!! call insert_into_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1)
|
||||
! endif
|
||||
! if (n_integrals_2 > 0) then
|
||||
!! call insert_into_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2)
|
||||
! endif
|
||||
!
|
||||
! call map_sort(ao_integrals_map)
|
||||
! call map_unique(ao_integrals_map)
|
||||
! call map_sort(ao_integrals_map_2)
|
||||
! call map_unique(ao_integrals_map_2)
|
||||
!
|
||||
! call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_periodic_1',ao_integrals_map)
|
||||
! call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_periodic_2',ao_integrals_map_2)
|
||||
! call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read'
|
||||
print*,'map1'
|
||||
do i=0,ao_integrals_map%map_size
|
||||
print*,i,ao_integrals_map%map(i)%value(:)
|
||||
print*,i,ao_integrals_map%map(i)%key(:)
|
||||
enddo
|
||||
print*,'map2'
|
||||
do i=0,ao_integrals_map_2%map_size
|
||||
print*,i,ao_integrals_map_2%map(i)%value(:)
|
||||
print*,i,ao_integrals_map_2%map(i)%key(:)
|
||||
enddo
|
||||
end
|
||||
@ -17,9 +17,13 @@ subroutine run
|
||||
double precision, allocatable :: A(:,:), B(:,:)
|
||||
double precision :: tmp_re, tmp_im
|
||||
|
||||
integer :: n_integrals
|
||||
integer(key_kind), allocatable :: buffer_i(:)
|
||||
real(integral_kind), allocatable :: buffer_values(:)
|
||||
integer :: n_integrals_1, n_integrals_2
|
||||
integer(key_kind), allocatable :: buffer_i_1(:), buffer_i_2(:)
|
||||
real(integral_kind), allocatable :: buffer_values_1(:), buffer_values_2(:)
|
||||
logical :: use_map1
|
||||
integer(key_kind) :: idx_tmp
|
||||
double precision :: sign
|
||||
|
||||
|
||||
call ezfio_set_ao_basis_ao_num(ao_num)
|
||||
|
||||
@ -107,31 +111,66 @@ subroutine run
|
||||
call ezfio_set_ao_one_e_ints_ao_integrals_n_e_imag(B(1:ao_num, 1:ao_num))
|
||||
call ezfio_set_ao_one_e_ints_io_ao_integrals_n_e("Read")
|
||||
|
||||
! allocate(buffer_i(ao_num**3), buffer_values(ao_num**3))
|
||||
! iunit = getunitandopen('W.qp','r')
|
||||
! n_integrals=0
|
||||
! buffer_values = 0.d0
|
||||
! do
|
||||
! read (iunit,*,end=13) i,j,k,l, integral
|
||||
! n_integrals += 1
|
||||
! call two_e_integrals_index(i, j, k, l, buffer_i(n_integrals) )
|
||||
! buffer_values(n_integrals) = integral
|
||||
! if (n_integrals == size(buffer_i)) then
|
||||
! call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_values)
|
||||
! n_integrals = 0
|
||||
! endif
|
||||
! enddo
|
||||
! 13 continue
|
||||
! close(iunit)
|
||||
!
|
||||
! if (n_integrals > 0) then
|
||||
! call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_values)
|
||||
! endif
|
||||
!
|
||||
! call map_sort(ao_integrals_map)
|
||||
! call map_unique(ao_integrals_map)
|
||||
!
|
||||
! call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map)
|
||||
! call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read')
|
||||
allocate(buffer_i_1(ao_num**3), buffer_values_1(ao_num**3))
|
||||
allocate(buffer_i_2(ao_num**3), buffer_values_2(ao_num**3))
|
||||
iunit = getunitandopen('W.qp','r')
|
||||
n_integrals_1=0
|
||||
n_integrals_2=0
|
||||
buffer_values_1 = 0.d0
|
||||
buffer_values_2 = 0.d0
|
||||
do
|
||||
read (iunit,*,end=13) i,j,k,l, tmp_re, tmp_im
|
||||
call ao_two_e_integral_periodic_map_idx_sign(i,j,k,l,use_map1,idx_tmp,sign)
|
||||
print'(4(I4),(L3),(I6),(F7.1))',i,j,k,l,use_map1,idx_tmp,sign
|
||||
if (use_map1) then
|
||||
n_integrals_1 += 1
|
||||
buffer_i_1(n_integrals_1)=idx_tmp
|
||||
buffer_values_1(n_integrals_1)=tmp_re
|
||||
print'(A,4(I4),(I6),(E15.7))','map1',i,j,k,l,idx_tmp,tmp_re
|
||||
if (sign.ne.0.d0) then
|
||||
n_integrals_1 += 1
|
||||
buffer_i_1(n_integrals_1)=idx_tmp+1
|
||||
buffer_values_1(n_integrals_1)=tmp_im*sign
|
||||
print'(A,4(I4),(I6),(E15.7))','map1',i,j,k,l,idx_tmp+1,tmp_im*sign
|
||||
endif
|
||||
if (n_integrals_1 >= size(buffer_i_1)-1) then
|
||||
call insert_into_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1)
|
||||
n_integrals_1 = 0
|
||||
endif
|
||||
else
|
||||
n_integrals_2 += 1
|
||||
buffer_i_2(n_integrals_2)=idx_tmp
|
||||
buffer_values_2(n_integrals_2)=tmp_re
|
||||
print'(A,4(I4),(I6),(E15.7))','map2',i,j,k,l,idx_tmp,tmp_re
|
||||
if (sign.ne.0.d0) then
|
||||
n_integrals_2 += 1
|
||||
buffer_i_2(n_integrals_2)=idx_tmp+1
|
||||
buffer_values_2(n_integrals_2)=tmp_im*sign
|
||||
print'(A,4(I4),(I6),(E15.7))','map2',i,j,k,l,idx_tmp+1,tmp_im*sign
|
||||
endif
|
||||
if (n_integrals_2 >= size(buffer_i_2)-1) then
|
||||
call insert_into_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2)
|
||||
n_integrals_2 = 0
|
||||
endif
|
||||
endif
|
||||
enddo
|
||||
13 continue
|
||||
close(iunit)
|
||||
|
||||
if (n_integrals_1 > 0) then
|
||||
call insert_into_ao_integrals_map(n_integrals_1,buffer_i_1,buffer_values_1)
|
||||
endif
|
||||
if (n_integrals_2 > 0) then
|
||||
call insert_into_ao_integrals_map_2(n_integrals_2,buffer_i_2,buffer_values_2)
|
||||
endif
|
||||
|
||||
call map_sort(ao_integrals_map)
|
||||
call map_unique(ao_integrals_map)
|
||||
call map_sort(ao_integrals_map_2)
|
||||
call map_unique(ao_integrals_map_2)
|
||||
|
||||
call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_periodic_1',ao_integrals_map)
|
||||
call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_periodic_2',ao_integrals_map_2)
|
||||
call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read')
|
||||
|
||||
end
|
||||
|
||||
Loading…
Reference in New Issue
Block a user