diff --git a/src/ao_one_e_ints/screening.irp.f b/src/ao_one_e_ints/screening.irp.f index 1bbe3c73..bc95ea86 100644 --- a/src/ao_one_e_ints/screening.irp.f +++ b/src/ao_one_e_ints/screening.irp.f @@ -3,7 +3,7 @@ logical function ao_one_e_integral_zero(i,k) integer, intent(in) :: i,k ao_one_e_integral_zero = .False. - if (.not.((io_ao_integrals_overlap/='None').or.is_periodic)) then + if (.not.((io_ao_integrals_overlap/='None').or.is_complex)) then if (ao_overlap_abs(i,k) < ao_integrals_threshold) then ao_one_e_integral_zero = .True. return diff --git a/src/ao_two_e_ints/map_integrals.irp.f b/src/ao_two_e_ints/map_integrals.irp.f index 3cf3eadb..3a0a2659 100644 --- a/src/ao_two_e_ints/map_integrals.irp.f +++ b/src/ao_two_e_ints/map_integrals.irp.f @@ -217,111 +217,111 @@ double precision function get_ao_two_e_integral(i,j,k,l,map) result(result) result = tmp end -BEGIN_PROVIDER [ complex*16, ao_integrals_cache_periodic, (0:64*64*64*64) ] - implicit none - BEGIN_DOC - ! Cache of AO integrals for fast access - END_DOC - PROVIDE ao_two_e_integrals_in_map - integer :: i,j,k,l,ii - integer(key_kind) :: idx1, idx2 - real(integral_kind) :: tmp_re, tmp_im - integer(key_kind) :: idx_re,idx_im - complex(integral_kind) :: integral +!BEGIN_PROVIDER [ complex*16, ao_integrals_cache_periodic, (0:64*64*64*64) ] +! implicit none +! BEGIN_DOC +! ! Cache of AO integrals for fast access +! END_DOC +! PROVIDE ao_two_e_integrals_in_map +! integer :: i,j,k,l,ii +! integer(key_kind) :: idx1, idx2 +! real(integral_kind) :: tmp_re, tmp_im +! integer(key_kind) :: idx_re,idx_im +! complex(integral_kind) :: integral +! +! +! !$OMP PARALLEL DO PRIVATE (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 +! 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 +! +! 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_periodic(ii) = integral +! enddo +! enddo +! enddo +! enddo +! !$OMP END PARALLEL DO +! +!END_PROVIDER - !$OMP PARALLEL DO PRIVATE (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 - 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 - - 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_periodic(ii) = integral - enddo - enddo - enddo - enddo - !$OMP END PARALLEL DO - -END_PROVIDER - - -complex*16 function get_ao_two_e_integral_periodic(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) :: idx1,idx2 - real(integral_kind) :: tmp_re, tmp_im - integer(key_kind) :: idx_re,idx_im - type(map_type), intent(inout) :: map - integer :: ii - complex(integral_kind) :: tmp - PROVIDE ao_two_e_integrals_in_map ao_integrals_cache_periodic ao_integrals_cache_min - !DIR$ FORCEINLINE - logical, external :: ao_two_e_integral_zero - if (ao_two_e_integral_zero(i,j,k,l)) then - tmp = (0.d0,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 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 - 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_periodic(ii) - endif - result = tmp - endif -end +!complex*16 function get_ao_two_e_integral_periodic(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) :: idx1,idx2 +! real(integral_kind) :: tmp_re, tmp_im +! integer(key_kind) :: idx_re,idx_im +! type(map_type), intent(inout) :: map +! integer :: ii +! complex(integral_kind) :: tmp +! PROVIDE ao_two_e_integrals_in_map ao_integrals_cache_periodic ao_integrals_cache_min +! !DIR$ FORCEINLINE +! logical, external :: ao_two_e_integral_zero +! if (ao_two_e_integral_zero(i,j,k,l)) then +! tmp = (0.d0,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 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 +! 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_periodic(ii) +! endif +! result = tmp +! endif +!end subroutine get_ao_two_e_integrals(j,k,l,sze,out_val) @@ -353,33 +353,33 @@ subroutine get_ao_two_e_integrals(j,k,l,sze,out_val) end -subroutine get_ao_two_e_integrals_periodic(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. - ! physicist convention : - END_DOC - implicit none - integer, intent(in) :: j,k,l, sze - complex(integral_kind), intent(out) :: out_val(sze) - - integer :: i - integer(key_kind) :: hash - logical, external :: ao_one_e_integral_zero - PROVIDE ao_two_e_integrals_in_map ao_integrals_map - - if (ao_one_e_integral_zero(j,l)) then - out_val = 0.d0 - return - endif - - double precision :: get_ao_two_e_integral - do i=1,sze - out_val(i) = get_ao_two_e_integral(i,j,k,l,ao_integrals_map) - enddo - -end +!subroutine get_ao_two_e_integrals_periodic(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. +! ! physicist convention : +! END_DOC +! implicit none +! integer, intent(in) :: j,k,l, sze +! complex(integral_kind), intent(out) :: out_val(sze) +! +! integer :: i +! integer(key_kind) :: hash +! logical, external :: ao_one_e_integral_zero +! PROVIDE ao_two_e_integrals_in_map ao_integrals_map +! +! if (ao_one_e_integral_zero(j,l)) then +! out_val = 0.d0 +! return +! endif +! +! double precision :: get_ao_two_e_integral +! do i=1,sze +! out_val(i) = get_ao_two_e_integral(i,j,k,l,ao_integrals_map) +! enddo +! +!end subroutine get_ao_two_e_integrals_non_zero(j,k,l,sze,out_val,out_val_index,non_zero_int) use map_module diff --git a/src/ao_two_e_ints/map_integrals_cplx.irp.f b/src/ao_two_e_ints/map_integrals_cplx.irp.f index 449bca02..a5926e59 100644 --- a/src/ao_two_e_ints/map_integrals_cplx.irp.f +++ b/src/ao_two_e_ints/map_integrals_cplx.irp.f @@ -343,11 +343,10 @@ complex*16 function get_ao_two_e_integral_complex(i,j,k,l,map,map2) result(resul ! a.le.c, b.le.d, tri(a,c).le.tri(b,d) PROVIDE ao_two_e_integrals_in_map ao_integrals_cache_complex 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) -! else if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < ao_integrals_threshold) then -! tmp = (0.d0,0.d0) -! else + !logical, external :: ao_two_e_integral_zero + !if (ao_two_e_integral_zero(i,j,k,l)) then + ! tmp = (0.d0,0.d0) + !else if (.True.) then ii = l-ao_integrals_cache_min ii = ior(ii, k-ao_integrals_cache_min) @@ -362,8 +361,8 @@ complex*16 function get_ao_two_e_integral_complex(i,j,k,l,map,map2) result(resul ii = ior( shiftl(ii,6), i-ao_integrals_cache_min) tmp = ao_integrals_cache_complex(ii) endif - result = tmp endif + result = tmp end @@ -380,14 +379,13 @@ subroutine get_ao_two_e_integrals_complex(j,k,l,sze,out_val) integer :: i integer(key_kind) :: hash - double precision :: thresh + !logical, external :: ao_one_e_integral_zero PROVIDE ao_two_e_integrals_in_map ao_integrals_map - thresh = ao_integrals_threshold - if (ao_overlap_abs(j,l) < thresh) then - out_val = (0.d0,0.d0) - return - endif + !if (ao_one_e_integral_zero(j,l)) then + ! out_val = (0.d0,0.d0) + ! return + !endif complex*16 :: get_ao_two_e_integral_complex do i=1,sze @@ -397,17 +395,17 @@ subroutine get_ao_two_e_integrals_complex(j,k,l,sze,out_val) end subroutine get_ao_two_e_integrals_non_zero_complex(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 + complex(integral_kind), intent(out) :: out_val(sze) + integer, intent(out) :: out_val_index(sze),non_zero_int print*,'not implemented for periodic',irp_here stop -1 -! 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 @@ -445,18 +443,18 @@ end subroutine get_ao_two_e_integrals_non_zero_jl_complex(j,l,thresh,sze_max,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 + double precision, intent(in) :: thresh + integer, intent(in) :: j,l, sze,sze_max + complex(integral_kind), intent(out) :: out_val(sze_max) + integer, intent(out) :: out_val_index(2,sze_max),non_zero_int print*,'not implemented for periodic',irp_here stop -1 -! 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 -! double precision, intent(in) :: thresh -! integer, intent(in) :: j,l, sze,sze_max -! real(integral_kind), intent(out) :: out_val(sze_max) -! integer, intent(out) :: out_val_index(2,sze_max),non_zero_int ! ! integer :: i,k ! integer(key_kind) :: hash @@ -496,19 +494,19 @@ end subroutine get_ao_two_e_integrals_non_zero_jl_from_list_complex(j,l,thresh,list,n_list,sze_max,out_val,out_val_index,non_zero_int) + use map_module + implicit none + BEGIN_DOC + ! Gets multiple AO two-electron integrals from the AO map . + ! All non-zero i are retrieved for j,k,l fixed. + END_DOC + double precision, intent(in) :: thresh + integer, intent(in) :: sze_max + integer, intent(in) :: j,l, n_list,list(2,sze_max) + complex(integral_kind), intent(out) :: out_val(sze_max) + integer, intent(out) :: out_val_index(2,sze_max),non_zero_int print*,'not implemented for periodic',irp_here stop -1 -! use map_module -! implicit none -! BEGIN_DOC -! ! Gets multiple AO two-electron integrals from the AO map . -! ! All non-zero i are retrieved for j,k,l fixed. -! END_DOC -! double precision, intent(in) :: thresh -! integer, intent(in) :: sze_max -! integer, intent(in) :: j,l, n_list,list(2,sze_max) -! real(integral_kind), intent(out) :: out_val(sze_max) -! integer, intent(out) :: out_val_index(2,sze_max),non_zero_int ! ! integer :: i,k ! integer(key_kind) :: hash diff --git a/src/ao_two_e_ints/screening.irp.f b/src/ao_two_e_ints/screening.irp.f index d3230370..eebe0043 100644 --- a/src/ao_two_e_ints/screening.irp.f +++ b/src/ao_two_e_ints/screening.irp.f @@ -3,7 +3,7 @@ logical function ao_two_e_integral_zero(i,j,k,l) integer, intent(in) :: i,j,k,l ao_two_e_integral_zero = .False. - if (.not.(read_ao_two_e_integrals.or.is_periodic)) then + if (.not.(read_ao_two_e_integrals.or.is_complex)) then if (ao_overlap_abs(j,l)*ao_overlap_abs(i,k) < ao_integrals_threshold) then ao_two_e_integral_zero = .True. return diff --git a/src/mo_one_e_ints/pot_mo_ints_cplx.irp.f b/src/mo_one_e_ints/pot_mo_ints_cplx.irp.f index a9f793d9..f472a8ff 100644 --- a/src/mo_one_e_ints/pot_mo_ints_cplx.irp.f +++ b/src/mo_one_e_ints/pot_mo_ints_cplx.irp.f @@ -6,8 +6,8 @@ BEGIN_PROVIDER [complex*16, mo_integrals_n_e_complex, (mo_num,mo_num)] integer :: i,j print *, 'Providing MO N-e integrals' - if (read_mo_integrals_e_n) then - call ezfio_get_mo_one_e_ints_mo_integrals_e_n_complex(mo_integrals_n_e_complex) + if (read_mo_integrals_n_e) then + call ezfio_get_mo_one_e_ints_mo_integrals_n_e_complex(mo_integrals_n_e_complex) print *, 'MO N-e integrals read from disk' else print *, 'Providing MO N-e integrals from AO N-e integrals' @@ -18,8 +18,8 @@ BEGIN_PROVIDER [complex*16, mo_integrals_n_e_complex, (mo_num,mo_num)] size(mo_integrals_n_e_complex,1) & ) endif - if (write_mo_integrals_e_n) then - call ezfio_set_mo_one_e_ints_mo_integrals_e_n_complex(mo_integrals_n_e_complex) + if (write_mo_integrals_n_e) then + call ezfio_set_mo_one_e_ints_mo_integrals_n_e_complex(mo_integrals_n_e_complex) print *, 'MO N-e integrals written to disk' endif @@ -39,8 +39,8 @@ BEGIN_PROVIDER [complex*16, mo_integrals_n_e_kpts, (mo_num_per_kpt,mo_num_per_kp integer :: i,j print *, 'Providing MO N-e integrals' - if (read_mo_integrals_e_n) then - call ezfio_get_mo_one_e_ints_mo_integrals_e_n_kpts(mo_integrals_n_e_kpts) + if (read_mo_integrals_n_e) then + call ezfio_get_mo_one_e_ints_mo_integrals_n_e_kpts(mo_integrals_n_e_kpts) print *, 'MO N-e integrals read from disk' else print *, 'Providing MO N-e integrals from AO N-e integrals' @@ -51,8 +51,8 @@ BEGIN_PROVIDER [complex*16, mo_integrals_n_e_kpts, (mo_num_per_kpt,mo_num_per_kp size(mo_integrals_n_e_kpts,1) & ) endif - if (write_mo_integrals_e_n) then - call ezfio_set_mo_one_e_ints_mo_integrals_e_n_kpts(mo_integrals_n_e_kpts) + if (write_mo_integrals_n_e) then + call ezfio_set_mo_one_e_ints_mo_integrals_n_e_kpts(mo_integrals_n_e_kpts) print *, 'MO N-e integrals written to disk' endif diff --git a/src/scf_utils/huckel_cplx.irp.f b/src/scf_utils/huckel_cplx.irp.f index ec504d14..f5dee3a4 100644 --- a/src/scf_utils/huckel_cplx.irp.f +++ b/src/scf_utils/huckel_cplx.irp.f @@ -19,7 +19,7 @@ subroutine huckel_guess_complex enddo A(j,j) = ao_one_e_integrals_diag_complex(j) + dble(ao_two_e_integral_alpha_complex(j,j)) if (dabs(dimag(ao_two_e_integral_alpha_complex(j,j))) .gt. 1.0d-10) then - stop 'diagonal elements of ao_bi_elec_integral_alpha should be real' + stop 'diagonal elements of ao_two_e_integral_alpha should be real' endif enddo @@ -67,7 +67,7 @@ subroutine huckel_guess_kpts enddo A(j,j) = ao_one_e_integrals_diag_kpts(j,k) + dble(ao_two_e_integral_alpha_kpts(j,j,k)) if (dabs(dimag(ao_two_e_integral_alpha_kpts(j,j,k))) .gt. 1.0d-10) then - stop 'diagonal elements of ao_bi_elec_integral_alpha should be real' + stop 'diagonal elements of ao_two_e_integral_alpha should be real' endif enddo diff --git a/src/utils_complex/create_ezfio_complex_3idx.py b/src/utils_complex/create_ezfio_complex_3idx.py index 0360cfe8..3fef73d2 100755 --- a/src/utils_complex/create_ezfio_complex_3idx.py +++ b/src/utils_complex/create_ezfio_complex_3idx.py @@ -145,12 +145,12 @@ def convert_kpts(filename,qph5path): ezfio.set_mo_one_e_ints_mo_integrals_kinetic_kpts(kin_mo_reim) ezfio.set_mo_one_e_ints_mo_integrals_overlap_kpts(ovlp_mo_reim) #ezfio.set_mo_one_e_ints_mo_integrals_n_e_complex(ne_mo_reim) - ezfio.set_mo_one_e_ints_mo_integrals_e_n_kpts(ne_mo_reim) + ezfio.set_mo_one_e_ints_mo_integrals_n_e_kpts(ne_mo_reim) ezfio.set_mo_one_e_ints_io_mo_integrals_kinetic('Read') ezfio.set_mo_one_e_ints_io_mo_integrals_overlap('Read') #ezfio.set_mo_one_e_ints_io_mo_integrals_n_e('Read') - ezfio.set_mo_one_e_ints_io_mo_integrals_e_n('Read') + ezfio.set_mo_one_e_ints_io_mo_integrals_n_e('Read') ########################################## # # @@ -331,12 +331,12 @@ def convert_cplx(filename,qph5path): ezfio.set_mo_one_e_ints_mo_integrals_kinetic_complex(kin_mo_reim) #ezfio.set_mo_one_e_ints_mo_integrals_overlap_complex(ovlp_mo_reim) #ezfio.set_mo_one_e_ints_mo_integrals_n_e_complex(ne_mo_reim) - ezfio.set_mo_one_e_ints_mo_integrals_e_n_complex(ne_mo_reim) + ezfio.set_mo_one_e_ints_mo_integrals_n_e_complex(ne_mo_reim) ezfio.set_mo_one_e_ints_io_mo_integrals_kinetic('Read') #ezfio.set_mo_one_e_ints_io_mo_integrals_overlap('Read') #ezfio.set_mo_one_e_ints_io_mo_integrals_n_e('Read') - ezfio.set_mo_one_e_ints_io_mo_integrals_e_n('Read') + ezfio.set_mo_one_e_ints_io_mo_integrals_n_e('Read') ########################################## # #