diff --git a/src/ao_basis/ao_transfo_matrix.irp.f b/src/ao_basis/ao_transfo_matrix.irp.f new file mode 100644 index 00000000..2405ea94 --- /dev/null +++ b/src/ao_basis/ao_transfo_matrix.irp.f @@ -0,0 +1,107 @@ +BEGIN_PROVIDER [ integer, ao_num ] + implicit none + BEGIN_DOC +! Number of |AOs| + END_DOC + + logical :: has + PROVIDE ezfio_filename + if (mpi_master) then + + call ezfio_has_ao_basis_ao_num(has) + if (has) then +! write(6,'(A)') '.. >>>>> [ IO READ: ao_num ] <<<<< ..' + call ezfio_get_ao_basis_ao_num(ao_num) + else + print *, 'Using the transformation matrix' + if(.True.)then + ao_num = ao_cart_num + endif + endif + endif + IRP_IF MPI_DEBUG + print *, irp_here, mpi_rank + call MPI_BARRIER(MPI_COMM_WORLD, ierr) + IRP_ENDIF + IRP_IF MPI + include 'mpif.h' + integer :: ierr + call MPI_BCAST( ao_num, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read ao_num with MPI' + endif + IRP_ENDIF + +! call write_time(6) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_cart_to_ao_basis_mat, (ao_num, ao_cart_num)] + implicit none + if(.True.)then + ao_cart_to_ao_basis_mat=0.d0 + integer :: i + do i = 1, ao_num + ao_cart_to_ao_basis_mat(i,i) = 1.d0 + enddo + endif + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_cart_to_ao_basis_mat_transp, (ao_cart_num, ao_num)] + implicit none + integer :: i,j + do i = 1, ao_num + do j = 1, ao_cart_num + ao_cart_to_ao_basis_mat_transp(j,i) = ao_cart_to_ao_basis_mat(i,j) + enddo + enddo +END_PROVIDER + +subroutine ao_cart_to_ao_basis(A_cart, LDA_cart, A_ao_basis, LDA_ao_basis) + implicit none + BEGIN_DOC + ! transform a matrix "A_cart" expressed in the cartesian AO basis + ! + ! to the "AO" basis with the transformation matrix "ao_cart_to_ao_basis_mat" + ! + ! A_cart = ao_cart_to_ao_basis_mat A_cart ao_cart_to_ao_basis_mat^T + END_DOC + double precision, intent(in) :: A_cart(LDA_cart, ao_cart_num) + double precision, intent(out :: A_ao_basis(LDA_ao_basis, ao_num) + integer, intent(in) :: LDA_cart, LDA_ao_basis + double precision, allocatable :: tmp(:,:) + allocate (tmp(ao_num,ao_cart_num)) + + call dgemm('N','N',ao_num,ao_cart_num,ao_cart_num, 1.d0, & + ao_cart_to_ao_basis_mat,size(ao_cart_to_ao_basis_mat,1), & + A_cart,size(A_cart,1), 0.d0, & + tmp, size(tmp,1)) + + call dgemm('N','T',ao_num,ao_num,ao_cart_num, 1.d0, & + tmp, size(tmp,1), & + ao_cart_to_ao_basis_mat,size(ao_cart_to_ao_basis_mat,1), 0.d0, & + A_ao_basis,size(A_ao_basis,1)) + + deallocate(tmp) + +end + +subroutine ao_cart_to_ao_basis_vec(V_cart, V_ao_basis) + implicit none + BEGIN_DOC + ! transform a VECTOR "V_cart" expressed in the cartesian AO basis + ! + ! to the "AO" basis with the transformation matrix "ao_cart_to_ao_basis_mat" + ! + ! V_cart = ao_cart_to_ao_basis_mat V_cart + END_DOC + double precision, intent(in) :: V_cart(ao_cart_num) + double precision, intent(out :: A_ao_basis(ao_num) + + call dgemv('N',ao_num,ao_cart_num, 1.d0, & + ao_cart_to_ao_basis_mat,size(ao_cart_to_ao_basis_mat,1), & + V_cart,1, 0.d0, V_ao_basis, 1) + +end + diff --git a/src/ao_cart_two_e_ints/README.rst b/src/ao_cart_two_e_ints/README.rst index f2dd5206..5bd6813e 100644 --- a/src/ao_cart_two_e_ints/README.rst +++ b/src/ao_cart_two_e_ints/README.rst @@ -3,12 +3,6 @@ ao_cart_two_e_ints ================== Here, all two-electron integrals (:math:`1/r_{12}`) are computed. -As they have 4 indices and many are zero, they are stored in a map, as defined -in :file:`utils/map_module.f90`. - -To fetch an |AO| integral, use the -`get_ao_cart_two_e_integral(i,j,k,l,ao_cart_integrals_map)` function. - The conventions are: * For |AO| integrals : (ij|kl) = (11|22) = = <12|12> diff --git a/src/ao_cart_two_e_ints/cholesky.irp.f b/src/ao_cart_two_e_ints/cholesky.irp.f index 2f68d11a..d565b810 100644 --- a/src/ao_cart_two_e_ints/cholesky.irp.f +++ b/src/ao_cart_two_e_ints/cholesky.irp.f @@ -1,16 +1,18 @@ -double precision function get_ao_cart_integ_chol(i,j,k,l) +BEGIN_TEMPLATE +double precision function get_ao_cart$_erf_integ_chol(i,j,k,l) implicit none BEGIN_DOC ! CHOLESKY representation of the integral of the AO basis or (ij|kl) ! i(r1) j(r1) 1/r12 k(r2) l(r2) + ! END_DOC integer, intent(in) :: i,j,k,l double precision, external :: ddot - get_ao_cart_integ_chol = ddot(cholesky_ao_cart_num, cholesky_ao_cart_transp(1,i,j), 1, cholesky_ao_cart_transp(1,k,l), 1) + get_ao_cart$_erf_integ_chol = ddot(cholesky_ao_cart$_erf_num, cholesky_ao_cart$_erf_transp(1,i,j), 1, cholesky_ao_cart$_erf_transp(1,k,l), 1) end -BEGIN_PROVIDER [ double precision, cholesky_ao_cart_transp, (cholesky_ao_cart_num, ao_cart_num, ao_cart_num) ] +BEGIN_PROVIDER [ double precision, cholesky_ao_cart_transp, (cholesky_ao$_erf_cart_num, ao_cart_num, ao_cart_num) ] implicit none BEGIN_DOC ! Transposed of the Cholesky vectors in AO basis set @@ -18,23 +20,23 @@ BEGIN_PROVIDER [ double precision, cholesky_ao_cart_transp, (cholesky_ao_cart_nu integer :: i,j,k do j=1,ao_cart_num do i=1,ao_cart_num - do k=1,cholesky_ao_cart_num - cholesky_ao_cart_transp(k,i,j) = cholesky_ao(i,j,k) + do k=1,cholesky_ao_cart$_erf_num + cholesky_ao_cart$_erf_transp(k,i,j) = cholesky_ao(i,j,k) enddo enddo enddo END_PROVIDER - BEGIN_PROVIDER [ integer, cholesky_ao_cart_num ] -&BEGIN_PROVIDER [ double precision, cholesky_ao, (ao_cart_num, ao_cart_num, 1) ] + BEGIN_PROVIDER [ integer, cholesky_ao_cart$_erf_num ] +&BEGIN_PROVIDER [ double precision, cholesky_ao_cart$_erf, (ao_cart_num, ao_cart_num, 1) ] use mmap_module implicit none BEGIN_DOC ! Cholesky vectors in AO basis: (ik|a): ! = (ik|jl) = sum_a (ik|a).(a|jl) ! - ! Last dimension of cholesky_ao is cholesky_ao_cart_num + ! Last dimension of cholesky_ao is cholesky_ao_cart$_erf_num ! ! https://mogp-emulator.readthedocs.io/en/latest/methods/proc/ProcPivotedCholesky.html ! @@ -60,10 +62,9 @@ END_PROVIDER integer :: N, np, nq double precision :: Dmax, Dmin, Qmax, f - double precision, external :: get_ao_cart_two_e_integral logical, external :: ao_cart_two_e_integral_zero - double precision, external :: ao_cart_two_e_integral + double precision, external :: ao_cart_two_e_integral$_erf integer :: block_size, iblock double precision :: mem, mem0 @@ -78,7 +79,7 @@ END_PROVIDER type(mmap_type) :: map - PROVIDE nproc ao_cart_cholesky_threshold do_direct_integrals qp_max_mem + PROVIDE nproc ao_cholesky_threshold do_direct_integrals qp_max_mem PROVIDE nucl_coord call set_multiple_levels_omp(.False.) @@ -87,25 +88,25 @@ END_PROVIDER ! Will be reallocated at the end deallocate(cholesky_ao) - if (read_ao_cart_cholesky) then - print *, 'Reading Cholesky AO vectors from disk...' - iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao', 'R') + if (read_ao_cart$_erf_cholesky) then + print *, 'Reading Cholesky AO$_erf vectors from disk...' + iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao_cart$_erf', 'R') read(iunit) rank allocate(cholesky_ao(ao_cart_num,ao_cart_num,rank), stat=ierr) read(iunit) cholesky_ao close(iunit) - cholesky_ao_cart_num = rank + cholesky_ao_cart$_erf_num = rank else call set_multiple_levels_omp(.False.) - if (ao_cart_two_e_integral(1,1,1,1) < huge(1.d0)) then + if (ao_cart_two_e_integral$_erf(1,1,1,1) < huge(1.d0)) then ! Trigger providers inside ao_cart_two_e_integral continue endif - tau = ao_cart_cholesky_threshold + tau = ao_cholesky_threshold tau2 = tau*tau rank = 0 @@ -118,7 +119,7 @@ END_PROVIDER call print_memory_usage() print *, '' - print *, 'Cholesky decomposition of AO integrals' + print *, 'Cholesky decomposition of AO$_erf integrals' print *, '======================================' print *, '' print *, '============ =============' @@ -138,7 +139,7 @@ END_PROVIDER !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i8) SCHEDULE(dynamic,21) do i8=ndim8-1,1,-1 - D(i8) = ao_cart_two_e_integral(addr1(i8), addr2(i8), & + D(i8) = ao_cart_two_e_integral$_erf(addr1(i8), addr2(i8), & addr1(i8), addr2(i8)) enddo !$OMP END PARALLEL DO @@ -345,7 +346,7 @@ END_PROVIDER if (.not.ao_cart_two_e_integral_zero( addr1(Lset(k)), addr1(Dset(m)),& addr2(Lset(k)), addr2(Dset(m)) ) ) then Delta_col(k) = & - ao_cart_two_e_integral(addr1(Lset(k)), addr2(Lset(k)),& + ao_cart_two_e_integral$_erf(addr1(Lset(k)), addr2(Lset(k)),& addr1(Dset(m)), addr2(Dset(m))) endif enddo @@ -448,24 +449,31 @@ END_PROVIDER call mmap_destroy(map) - cholesky_ao_cart_num = rank + cholesky_ao_cart$_erf_num = rank - if (write_ao_cart_cholesky) then - print *, 'Writing Cholesky AO vectors to disk...' - iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao', 'W') + if (write_ao_cart$_erf_cholesky) then + print *, 'Writing Cholesky AO$_erf vectors to disk...' + iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao_cart', 'W') write(iunit) rank - write(iunit) cholesky_ao + write(iunit) cholesky_ao_cart close(iunit) - call ezfio_set_ao_cart_two_e_ints_io_ao_cart_cholesky('Read') + call ezfio_set_ao_cart_two_e_ints_io_ao_cart$_erf_cholesky('Read') endif endif - print *, 'Rank : ', cholesky_ao_cart_num, '(', 100.d0*dble(cholesky_ao_cart_num)/dble(ao_cart_num*ao_cart_num), ' %)' + print *, 'Rank : ', cholesky_ao_cart$_erf_num, '(', 100.d0*dble(cholesky_ao_cart$_erf_num)/dble(ao_cart_num*ao_cart_num), ' %)' print *, '' call wall_time(wall1) - print*,'Time to provide AO cholesky vectors = ',(wall1-wall0)/60.d0, ' min' + print*,'Time to provide AO$_erf cholesky vectors = ',(wall1-wall0)/60.d0, ' min' END_PROVIDER +SUBST [ _erf ] + +;; +_erf;; +_cgtos;; + +END_TEMPLATE diff --git a/src/ao_cart_two_e_ints/map_integrals_erf.irp.f b/src/ao_cart_two_e_ints/map_integrals_erf.irp.f deleted file mode 100644 index f24f8080..00000000 --- a/src/ao_cart_two_e_ints/map_integrals_erf.irp.f +++ /dev/null @@ -1,288 +0,0 @@ -use map_module - -!! AO Map -!! ====== - -BEGIN_PROVIDER [ type(map_type), ao_cart_integrals_erf_map ] - implicit none - BEGIN_DOC - ! |AO| integrals - END_DOC - integer(key_kind) :: key_max - integer(map_size_kind) :: sze - call two_e_integrals_index(ao_cart_num,ao_cart_num,ao_cart_num,ao_cart_num,key_max) - sze = key_max - call map_init(ao_cart_integrals_erf_map,sze) - print*, 'AO map initialized : ', sze -END_PROVIDER - - BEGIN_PROVIDER [ integer, ao_cart_integrals_erf_cache_min ] -&BEGIN_PROVIDER [ integer, ao_cart_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_cart_integrals_erf_cache_min = max(1,ao_cart_num - 63) - ao_cart_integrals_erf_cache_max = ao_cart_num - -END_PROVIDER - -BEGIN_PROVIDER [ double precision, ao_cart_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_cart_two_e_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_cart_integrals_erf_cache_min,ao_cart_integrals_erf_cache_max - do k=ao_cart_integrals_erf_cache_min,ao_cart_integrals_erf_cache_max - do j=ao_cart_integrals_erf_cache_min,ao_cart_integrals_erf_cache_max - do i=ao_cart_integrals_erf_cache_min,ao_cart_integrals_erf_cache_max - !DIR$ FORCEINLINE - call two_e_integrals_index(i,j,k,l,idx) - !DIR$ FORCEINLINE - call map_get(ao_cart_integrals_erf_map,idx,integral) - ii = l-ao_cart_integrals_erf_cache_min - ii = ior( ishft(ii,6), k-ao_cart_integrals_erf_cache_min) - ii = ior( ishft(ii,6), j-ao_cart_integrals_erf_cache_min) - ii = ior( ishft(ii,6), i-ao_cart_integrals_erf_cache_min) - ao_cart_integrals_erf_cache(ii) = integral - enddo - enddo - enddo - enddo - !$OMP END PARALLEL DO - -END_PROVIDER - - -subroutine insert_into_ao_cart_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_cart_integrals_erf_map, buffer_i, buffer_values, n_integrals) -end - -double precision function get_ao_cart_two_e_integral_erf(i,j,k,l,map) result(result) - use map_module - implicit none - BEGIN_DOC - ! Gets one |AO| two-electron 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 - logical, external :: ao_cart_two_e_integral_zero - PROVIDE ao_cart_two_e_integrals_erf_in_map ao_cart_integrals_erf_cache ao_cart_integrals_erf_cache_min - !DIR$ FORCEINLINE - if (ao_cart_two_e_integral_zero(i,j,k,l)) then - tmp = 0.d0 - else if (ao_cart_two_e_integral_erf_schwartz(i,k)*ao_cart_two_e_integral_erf_schwartz(j,l) < ao_cart_integrals_threshold) then - tmp = 0.d0 - else - ii = l-ao_cart_integrals_erf_cache_min - ii = ior(ii, k-ao_cart_integrals_erf_cache_min) - ii = ior(ii, j-ao_cart_integrals_erf_cache_min) - ii = ior(ii, i-ao_cart_integrals_erf_cache_min) - if (iand(ii, -64) /= 0) then - !DIR$ FORCEINLINE - call two_e_integrals_index(i,j,k,l,idx) - !DIR$ FORCEINLINE - call map_get(map,idx,tmp) - tmp = tmp - else - ii = l-ao_cart_integrals_erf_cache_min - ii = ior( ishft(ii,6), k-ao_cart_integrals_erf_cache_min) - ii = ior( ishft(ii,6), j-ao_cart_integrals_erf_cache_min) - ii = ior( ishft(ii,6), i-ao_cart_integrals_erf_cache_min) - tmp = ao_cart_integrals_erf_cache(ii) - endif - endif - result = tmp -end - - -subroutine get_ao_cart_two_e_integrals_erf(j,k,l,sze,out_val) - use map_module - BEGIN_DOC - ! Gets multiple |AO| two-electron 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 - logical, external :: ao_cart_one_e_integral_zero - PROVIDE ao_cart_two_e_integrals_erf_in_map ao_cart_integrals_erf_map - thresh = ao_cart_integrals_threshold - - if (ao_cart_one_e_integral_zero(j,l)) then - out_val = 0.d0 - return - endif - - double precision :: get_ao_cart_two_e_integral_erf - do i=1,sze - out_val(i) = get_ao_cart_two_e_integral_erf(i,j,k,l,ao_cart_integrals_erf_map) - enddo - -end - -subroutine get_ao_cart_two_e_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| two-electron integrals 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 - logical, external :: ao_cart_one_e_integral_zero - PROVIDE ao_cart_two_e_integrals_erf_in_map - thresh = ao_cart_integrals_threshold - - non_zero_int = 0 - if (ao_cart_one_e_integral_zero(j,l)) then - out_val = 0.d0 - return - endif - - non_zero_int = 0 - do i=1,sze - integer, external :: ao_cart_l4 - double precision, external :: ao_cart_two_e_integral_erf - !DIR$ FORCEINLINE - if (ao_cart_two_e_integral_erf_schwartz(i,k)*ao_cart_two_e_integral_erf_schwartz(j,l) < thresh) then - cycle - endif - call two_e_integrals_index(i,j,k,l,hash) - call map_get(ao_cart_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_cart_erf_map_size() - implicit none - integer (map_size_kind) :: get_ao_cart_erf_map_size - BEGIN_DOC - ! Returns the number of elements in the |AO| map - END_DOC - get_ao_cart_erf_map_size = ao_cart_integrals_erf_map % n_elements -end - -subroutine clear_ao_cart_erf_map - implicit none - BEGIN_DOC - ! Frees the memory of the |AO| map - END_DOC - call map_deinit(ao_cart_integrals_erf_map) - FREE ao_cart_integrals_erf_map -end - - - -subroutine dump_ao_cart_integrals_erf(filename) - use map_module - implicit none - BEGIN_DOC - ! Save to disk the |AO| erf 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_cart_integrals_erf_map%sorted, ao_cart_integrals_erf_map%map_size, & - ao_cart_integrals_erf_map%n_elements - do i=0_8,ao_cart_integrals_erf_map%map_size - write(66) ao_cart_integrals_erf_map%map(i)%sorted, ao_cart_integrals_erf_map%map(i)%map_size,& - ao_cart_integrals_erf_map%map(i)%n_elements - enddo - do i=0_8,ao_cart_integrals_erf_map%map_size - key => ao_cart_integrals_erf_map%map(i)%key - val => ao_cart_integrals_erf_map%map(i)%value - n = ao_cart_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_cart_integrals_erf(filename) - implicit none - BEGIN_DOC - ! Read from disk the |AO| erf 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_cart_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_cart_integrals_erf_map%sorted, ao_cart_integrals_erf_map%map_size,& - ao_cart_integrals_erf_map%n_elements - do i=0_8, ao_cart_integrals_erf_map%map_size - read(66,err=99,end=99) ao_cart_integrals_erf_map%map(i)%sorted, & - ao_cart_integrals_erf_map%map(i)%map_size, ao_cart_integrals_erf_map%map(i)%n_elements - call cache_map_reallocate(ao_cart_integrals_erf_map%map(i),ao_cart_integrals_erf_map%map(i)%map_size) - enddo - do i=0_8, ao_cart_integrals_erf_map%map_size - key => ao_cart_integrals_erf_map%map(i)%key - val => ao_cart_integrals_erf_map%map(i)%value - n = ao_cart_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_cart_integrals_erf_map) - load_ao_cart_integrals_erf = 0 - return - 99 continue - call map_deinit(ao_cart_integrals_erf_map) - 98 continue - stop 'Problem reading ao_cart_integrals_erf_map file in work/' - -end - - - - diff --git a/src/ao_cart_two_e_ints/providers_ao_erf.irp.f b/src/ao_cart_two_e_ints/providers_ao_erf.irp.f deleted file mode 100644 index a01334da..00000000 --- a/src/ao_cart_two_e_ints/providers_ao_erf.irp.f +++ /dev/null @@ -1,112 +0,0 @@ -BEGIN_PROVIDER [ logical, ao_cart_two_e_integrals_erf_in_map ] - implicit none - use map_module - BEGIN_DOC - ! Map of Atomic integrals - ! i(r1) j(r2) 1/r12 k(r1) l(r2) - END_DOC - - integer :: i,j,k,l - double precision :: ao_cart_two_e_integral_erf,cpu_1,cpu_2, wall_1, wall_2 - double precision :: integral, wall_0 - include 'utils/constants.include.F' - - ! For integrals file - integer(key_kind),allocatable :: buffer_i(:) - integer :: size_buffer - real(integral_kind),allocatable :: buffer_value(:) - - integer :: n_integrals, rc - integer :: kk, m, j1, i1, lmax - character*(64) :: fmt - - double precision :: map_mb - PROVIDE read_ao_cart_two_e_integrals_erf io_ao_cart_two_e_integrals_erf ao_cart_integrals_erf_map - - if (read_ao_cart_two_e_integrals_erf) then - print*,'Reading the AO ERF integrals' - call map_load_from_disk(trim(ezfio_filename)//'/work/ao_cart_ints_erf',ao_cart_integrals_erf_map) - print*, 'AO ERF integrals provided' - ao_cart_two_e_integrals_erf_in_map = .True. - return - endif - - print*, 'Providing the AO ERF integrals' - call wall_time(wall_0) - call wall_time(wall_1) - call cpu_time(cpu_1) - - if (.True.) then - ! Avoid openMP - integral = ao_cart_two_e_integral_erf(1,1,1,1) - endif - - size_buffer = ao_cart_num*ao_cart_num - !$OMP PARALLEL DEFAULT(shared) private(j,l) & - !$OMP PRIVATE(buffer_i, buffer_value, n_integrals) - allocate(buffer_i(size_buffer), buffer_value(size_buffer)) - n_integrals = 0 - !$OMP DO COLLAPSE(1) SCHEDULE(dynamic) - do l=1,ao_cart_num - do j=1,l - call compute_ao_cart_integrals_erf_jl(j,l,n_integrals,buffer_i,buffer_value) - call insert_into_ao_cart_integrals_erf_map(n_integrals,buffer_i,buffer_value) - enddo - enddo - !$OMP END DO - deallocate(buffer_i, buffer_value) - !$OMP END PARALLEL - - - - print*, 'Sorting the map' - call map_sort(ao_cart_integrals_erf_map) - call cpu_time(cpu_2) - call wall_time(wall_2) - integer(map_size_kind) :: get_ao_cart_erf_map_size, ao_cart_erf_map_size - ao_cart_erf_map_size = get_ao_cart_erf_map_size() - - print*, 'AO ERF integrals provided:' - print*, ' Size of AO ERF map : ', map_mb(ao_cart_integrals_erf_map) ,'MB' - print*, ' Number of AO ERF integrals :', ao_cart_erf_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_cart_two_e_integrals_erf_in_map = .True. - - if (write_ao_cart_two_e_integrals_erf) then - call ezfio_set_work_empty(.False.) - call map_save_to_disk(trim(ezfio_filename)//'/work/ao_cart_ints_erf',ao_cart_integrals_erf_map) - call ezfio_set_ao_cart_two_e_ints_io_ao_cart_two_e_integrals_erf('Read') - endif - -END_PROVIDER - - - - -BEGIN_PROVIDER [ double precision, ao_cart_two_e_integral_erf_schwartz,(ao_cart_num,ao_cart_num) ] - implicit none - BEGIN_DOC - ! Needed to compute Schwartz inequalities - END_DOC - - integer :: i,k - double precision :: ao_cart_two_e_integral_erf,cpu_1,cpu_2, wall_1, wall_2 - - ao_cart_two_e_integral_erf_schwartz(1,1) = ao_cart_two_e_integral_erf(1,1,1,1) - !$OMP PARALLEL DO PRIVATE(i,k) & - !$OMP DEFAULT(NONE) & - !$OMP SHARED (ao_cart_num,ao_cart_two_e_integral_erf_schwartz) & - !$OMP SCHEDULE(dynamic) - do i=1,ao_cart_num - do k=1,i - ao_cart_two_e_integral_erf_schwartz(i,k) = dsqrt(ao_cart_two_e_integral_erf(i,k,i,k)) - ao_cart_two_e_integral_erf_schwartz(k,i) = ao_cart_two_e_integral_erf_schwartz(i,k) - enddo - enddo - !$OMP END PARALLEL DO - -END_PROVIDER - - diff --git a/src/ao_cart_two_e_ints/two_e_coul_integrals_cgtos.irp.f b/src/ao_cart_two_e_ints/two_e_integrals_cgtos.irp.f similarity index 100% rename from src/ao_cart_two_e_ints/two_e_coul_integrals_cgtos.irp.f rename to src/ao_cart_two_e_ints/two_e_integrals_cgtos.irp.f diff --git a/src/ao_cart_two_e_ints/two_e_integrals_erf.irp.f b/src/ao_cart_two_e_ints/two_e_integrals_erf.irp.f index 9a06475c..2f19fb41 100644 --- a/src/ao_cart_two_e_ints/two_e_integrals_erf.irp.f +++ b/src/ao_cart_two_e_ints/two_e_integrals_erf.irp.f @@ -601,58 +601,3 @@ 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 end - - -subroutine compute_ao_cart_integrals_erf_jl(j,l,n_integrals,buffer_i,buffer_value) - implicit none - use map_module - BEGIN_DOC - ! Parallel client for AO integrals - END_DOC - - integer, intent(in) :: j,l - integer,intent(out) :: n_integrals - integer(key_kind),intent(out) :: buffer_i(ao_cart_num*ao_cart_num) - real(integral_kind),intent(out) :: buffer_value(ao_cart_num*ao_cart_num) - - integer :: i,k - double precision :: ao_cart_two_e_integral_erf,cpu_1,cpu_2, wall_1, wall_2 - double precision :: integral, wall_0 - double precision :: thr - integer :: kk, m, j1, i1 - logical, external :: ao_cart_two_e_integral_zero - - thr = ao_cart_integrals_threshold - - n_integrals = 0 - - j1 = j+ishft(l*l-l,-1) - do k = 1, ao_cart_num ! r1 - i1 = ishft(k*k-k,-1) - if (i1 > j1) then - exit - endif - do i = 1, k - i1 += 1 - if (i1 > j1) then - exit - endif - if (ao_cart_two_e_integral_zero(i,j,k,l)) then - cycle - endif - if (ao_cart_two_e_integral_erf_schwartz(i,k)*ao_cart_two_e_integral_erf_schwartz(j,l) < thr ) then - cycle - endif - !DIR$ FORCEINLINE - integral = ao_cart_two_e_integral_erf(i,k,j,l) ! i,k : r1 j,l : r2 - if (abs(integral) < thr) then - cycle - endif - n_integrals += 1 - !DIR$ FORCEINLINE - call two_e_integrals_index(i,j,k,l,buffer_i(n_integrals)) - buffer_value(n_integrals) = integral - enddo - enddo - -end diff --git a/src/ao_two_e_ints/ao_bi_integrals.irp.f b/src/ao_two_e_ints/ao_bi_integrals.irp.f new file mode 100644 index 00000000..f0ad316b --- /dev/null +++ b/src/ao_two_e_ints/ao_bi_integrals.irp.f @@ -0,0 +1,150 @@ +subroutine ao_two_e_integrals_index(i,j,k,l,i1) + use map_module + implicit none + BEGIN_DOC + ! Computes an unique index for i,j,k,l integrals + END_DOC + 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 + +BEGIN_TEMPLATE + +BEGIN_PROVIDER [ logical, ao_two_e_integrals$_erf_in_map ] + use map_module + implicit none + BEGIN_DOC + ! If True, the map of AO two-electron integrals$_erf is provided + END_DOC + integer(bit_kind) :: mask_ijkl(N_int,4) + integer(bit_kind) :: mask_ijk(N_int,3) + double precision :: cpu_1, cpu_2, wall_1, wall_2 + + + ao_two_e_integrals$_erf_in_map = .True. + if (read_ao_two_e_integrals$_erf) then + print*,'Reading the AO integrals$_erf' + call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints$_erf',ao_integrals$_erf_map) + print*, 'AO integrals$_erf provided' + return + endif + + call wall_time(wall_1) + call cpu_time(cpu_1) + + if (do_ao_cholesky) then + PROVIDE cholesky_ao_transp + else + call add_integrals$_erf_to_map_cholesky_ao + integer*8 :: get_ao_map_size, ao_map_size + ao_map_size = get_ao_map_size() + double precision, external :: map_mb + print*,'Atomic integrals$_erf provided:' + print*,' Size of AO map ', map_mb(ao_integrals$_erf_map) ,'MB' + print*,' Number of AO integrals$_erf: ', ao_map_size + endif + + call wall_time(wall_2) + call cpu_time(cpu_2) + + 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), ')' + + if (write_ao_two_e_integrals$_erf.and.mpi_master) 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_ao_two_e_ints_io_ao_two_e_integrals$_erf('Read') + endif + +END_PROVIDER + +subroutine add_integrals$_erf_to_map_cholesky_ao + use bitmasks + implicit none + + BEGIN_DOC + ! Adds integrals$_erf to the AO map using Cholesky vectors + END_DOC + + integer :: i,j,k,l,m + integer :: size_buffer, n_integrals$_erf + size_buffer = min(ao_num*ao_num*ao_num,16000000) + + double precision, allocatable :: Vtmp(:,:,:) + integer(key_kind) , allocatable :: buffer_i(:) + real(integral_kind), allocatable :: buffer_value(:) + + PROVIDE cholesky_ao_transp + call set_multiple_levels_omp(.False.) + + !$OMP PARALLEL DEFAULT(SHARED) & + !$OMP PRIVATE(i,j,k,l,n_integrals$_erf,buffer_value, buffer_i, Vtmp) + allocate (buffer_i(size_buffer), buffer_value(size_buffer)) + allocate (Vtmp(ao_num,ao_num,ao_num)) + n_integrals$_erf = 0 + + !$OMP DO SCHEDULE(dynamic) + do l=1,ao_num + call dgemm('T','N',ao_num*ao_num,ao_num,cholesky_ao_num,1.d0, & + cholesky_ao_transp, cholesky_ao_num, & + cholesky_ao_transp(1,1,l), cholesky_ao_num, 0.d0, & + Vtmp, ao_num*ao_num) + + do k=1,l + do j=1,ao_num + do i=1,j + if (dabs(Vtmp(i,j,k)) > ao_integrals_threshold) then + n_integrals$_erf = n_integrals$_erf + 1 + buffer_value(n_integrals$_erf) = Vtmp(i,j,k) + !DIR$ FORCEINLINE + call ao_two_e_integrals$_erf_index(i,k,j,l,buffer_i(n_integrals$_erf)) + if (n_integrals$_erf == size_buffer) then + call map_append(ao_integrals$_erf_map, buffer_i, buffer_value, n_integrals$_erf) + n_integrals$_erf = 0 + endif + endif + enddo + enddo + enddo + enddo + !$OMP END DO NOWAIT + + if (n_integrals$_erf > 0) then + call map_append(ao_integrals$_erf_map, buffer_i, buffer_value, n_integrals$_erf) + endif + deallocate(buffer_i, buffer_value, Vtmp) + !$OMP BARRIER + !$OMP END PARALLEL + + call map_sort(ao_integrals$_erf_map) + call map_unique(ao_integrals$_erf_map) + +end + +subroutine clear_ao_map + implicit none + BEGIN_DOC + ! Frees the meaory of the AO map + END_DOC + call map_deinit(ao_integrals$_erf_map) + FREE ao_integrals$_erf_map + FREE ao_two_e_integrals$_erf_in_map +end + +SUBST [ _erf ] + +;; +_erf;; +_cgtos;; + +END_TEMPLATE diff --git a/src/ao_two_e_ints/cholesky.irp.f b/src/ao_two_e_ints/cholesky.irp.f index 5c83cec6..172e33d1 100644 --- a/src/ao_two_e_ints/cholesky.irp.f +++ b/src/ao_two_e_ints/cholesky.irp.f @@ -1,499 +1,163 @@ -double precision function get_ao_integ_chol(i,j,k,l) +BEGIN_TEMPLATE + +double precision function get_ao$_erf_integ_chol(i,j,k,l) implicit none BEGIN_DOC - ! CHOLESKY representation of the integral of the AO basis or (ij|kl) + ! CHOLESKY representation of the integral of the AO$_erf basis or (ij|kl) ! i(r1) j(r1) 1/r12 k(r2) l(r2) END_DOC integer, intent(in) :: i,j,k,l double precision, external :: ddot - get_ao_integ_chol = ddot(cholesky_ao_num, cholesky_ao_transp(1,i,j), 1, cholesky_ao_transp(1,k,l), 1) + get_ao$_erf_integ_chol = ddot(cholesky_ao$_erf_num, cholesky_ao$_erf_transp(1,i,j), 1, cholesky_ao$_erf_transp(1,k,l), 1) end -BEGIN_PROVIDER [ double precision, cholesky_ao_transp, (cholesky_ao_num, ao_num, ao_num) ] + BEGIN_PROVIDER [ integer, cholesky_ao$_erf_num ] implicit none BEGIN_DOC -! Transposed of the Cholesky vectors in AO basis set + ! Number of Cholesky vectors in MO basis END_DOC - integer :: i,j,k - do j=1,ao_num - do i=1,ao_num - do k=1,cholesky_ao_num - cholesky_ao_transp(k,i,j) = cholesky_ao(i,j,k) - enddo + integer, external :: getUnitAndOpen + integer :: iunit + if (read_ao$_erf_cholesky) then + iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao$_erf', 'R') + read(iunit) cholesky_ao$_erf_num + close(iunit) + else + cholesky_ao$_erf_num = cholesky_ao$_erf_cart_num + endif +END_PROVIDER + +BEGIN_PROVIDER [ double precision, cholesky_ao$_erf, (ao_num, ao_num, cholesky_ao$_erf_num) ] + implicit none + BEGIN_DOC + ! Cholesky vectors in ao basis + END_DOC + + integer :: k, i, j + + call set_multiple_levels_omp(.False.) + !$OMP PARALLEL DO PRIVATE(k) + do k=1,cholesky_ao$_erf_num + do j=1,ao_num + do i=1,ao_num + cholesky_ao$_erf(i,j,k) = cholesky_ao$_erf_transp(k,i,j) + enddo enddo enddo + !$OMP END PARALLEL DO + free cholesky_ao$_erf_transp + END_PROVIDER +BEGIN_PROVIDER [ double precision, cholesky_ao$_erf_transp, (cholesky_ao$_erf_num, ao_num, ao_num) ] + implicit none + BEGIN_DOC + ! Cholesky vectors in ao basis. Warning: it is transposed wrt cholesky_ao$_erf_cart: + ! + ! - cholesky_ao$_erf_cart is (ao_cart_num^2 x cholesky_ao$_erf_cart_num) + ! + ! - cholesky_ao$_erf_transp is (cholesky_ao$_erf_num x ao_num^2) + END_DOC - BEGIN_PROVIDER [ integer, cholesky_ao_num ] -&BEGIN_PROVIDER [ double precision, cholesky_ao, (ao_num, ao_num, 1) ] - use mmap_module - implicit none - BEGIN_DOC - ! Cholesky vectors in AO basis: (ik|a): - ! = (ik|jl) = sum_a (ik|a).(a|jl) - ! - ! Last dimension of cholesky_ao is cholesky_ao_num - ! - ! https://mogp-emulator.readthedocs.io/en/latest/methods/proc/ProcPivotedCholesky.html - ! - ! https://doi.org/10.1016/j.apnum.2011.10.001 : Page 4, Algorithm 1 - ! - ! https://www.diva-portal.org/smash/get/diva2:396223/FULLTEXT01.pdf - END_DOC + double precision, allocatable :: X(:,:,:) + double precision :: wall0, wall1 + integer, external :: getUnitAndOpen + integer :: iunit, ierr, rank - integer*8 :: ndim8 - integer :: rank - double precision :: tau, tau2 - double precision, pointer :: L(:,:) - - double precision :: s - - double precision, allocatable :: D(:), Ltmp_p(:,:), Ltmp_q(:,:), D_sorted(:), Delta_col(:), Delta(:,:) - integer, allocatable :: addr1(:), addr2(:) - integer*8, allocatable :: Lset(:), Dset(:) - logical, allocatable :: computed(:) - - integer :: i,j,k,m,p,q, dj, p2, q2, ii, jj - integer*8 :: i8, j8, p8, qj8, rank_max, np8 - integer :: N, np, nq - - double precision :: Dmax, Dmin, Qmax, f - double precision, external :: get_ao_two_e_integral - logical, external :: ao_two_e_integral_zero - - double precision, external :: ao_two_e_integral - integer :: block_size, iblock - - double precision :: mem, mem0 - double precision, external :: memory_of_double, memory_of_int - double precision, external :: memory_of_double8, memory_of_int8 - - integer, external :: getUnitAndOpen - integer :: iunit, ierr - - ndim8 = ao_num*ao_num*1_8+1 - double precision :: wall0,wall1 - - type(mmap_type) :: map - - PROVIDE nproc ao_cholesky_threshold do_direct_integrals qp_max_mem - PROVIDE nucl_coord - call set_multiple_levels_omp(.False.) - - call wall_time(wall0) - - ! Will be reallocated at the end - deallocate(cholesky_ao) - - if (read_ao_cholesky) then - print *, 'Reading Cholesky AO vectors from disk...' - iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao', 'R') + if (read_ao$_erf_cholesky) then + print *, 'Reading Cholesky ao$_erf vectors from disk...' + iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao$_erf_transp', 'R') read(iunit) rank - allocate(cholesky_ao(ao_num,ao_num,rank), stat=ierr) - read(iunit) cholesky_ao + if (cholesky_ao$_erf_num /= rank) then + stop 'inconsistent rank' + endif + read(iunit) cholesky_ao$_erf_transp close(iunit) - cholesky_ao_num = rank + else + print *, '' + print *, 'ao$_erf_cart->ao$_erf Transformation of Cholesky vectors' + print *, '-----------------------------------------' + print *, '' - else - - call set_multiple_levels_omp(.False.) - - if (do_direct_integrals) then - if (ao_two_e_integral(1,1,1,1) < huge(1.d0)) then - ! Trigger providers inside ao_two_e_integral - continue - endif - else - PROVIDE ao_two_e_integrals_in_map - endif - - tau = ao_cholesky_threshold - tau2 = tau*tau - - rank = 0 - - allocate( D(ndim8), Lset(ndim8), Dset(ndim8), D_sorted(ndim8)) - allocate( addr1(ndim8), addr2(ndim8), Delta_col(ndim8), computed(ndim8) ) - - call resident_memory(mem0) - - call print_memory_usage() - - print *, '' - print *, 'Cholesky decomposition of AO integrals' - print *, '======================================' - print *, '' - print *, '============ =============' - print *, ' Rank Threshold' - print *, '============ =============' - - - ! 1. - i8=0 - do j=1,ao_num - do i=1,ao_num - i8 = i8+1 - addr1(i8) = i - addr2(i8) = j - enddo - enddo - - if (do_direct_integrals) then - !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i8) SCHEDULE(dynamic,21) - do i8=ndim8-1,1,-1 - D(i8) = ao_two_e_integral(addr1(i8), addr2(i8), & - addr1(i8), addr2(i8)) - enddo - !$OMP END PARALLEL DO - else - !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i8) SCHEDULE(dynamic,21) - do i8=ndim8-1,1,-1 - D(i8) = get_ao_two_e_integral(addr1(i8), addr1(i8), & - addr2(i8), addr2(i8), ao_integrals_map) - enddo - !$OMP END PARALLEL DO - endif - ! Just to guarentee termination - D(ndim8) = 0.d0 - - D_sorted(:) = -D(:) - call dsort_noidx_big(D_sorted,ndim8) - D_sorted(:) = -D_sorted(:) - Dmax = D_sorted(1) - - ! 2. - np8=0_8 - do p8=1,ndim8 - if ( Dmax*D(p8) >= tau2 ) then - np8 = np8+1_8 - Lset(np8) = p8 - endif - enddo - if (np8 > ndim8) stop 'np>ndim8' - np = int(np8,4) - if (np <= 0) stop 'np<=0' - - rank_max = np - ! Avoid too large arrays when there are many electrons - if (elec_num > 10) then - rank_max = min(np,25*elec_num*elec_num) - endif - - call mmap_create_d('', (/ ndim8, rank_max /), .False., .True., map) - L => map%d2 - - ! 3. - N = 0 - - ! 4. - i = 0 - - mem = memory_of_double(np) & ! Delta(np,nq) - + (np+1)*memory_of_double(block_size) ! Ltmp_p(np,block_size) + Ltmp_q(nq,block_size) - -! call check_mem(mem) - ! 5. - do while ( (Dmax > tau).and.(np > 0) ) - ! a. - i = i+1 - - block_size = max(N,24) - - ! Determine nq so that Delta fits in memory - - s = 0.1d0 - Dmin = max(s*Dmax,tau) - do nq=2,np-1 - if (D_sorted(nq) < Dmin) exit - enddo - - do while (.True.) - - mem = mem0 & - + np*memory_of_double(nq) & ! Delta(np,nq) - + (np+nq)*memory_of_double(block_size) ! Ltmp_p(np,block_size) + Ltmp_q(nq,block_size) - - if (mem > qp_max_mem*0.5d0) then - Dmin = D_sorted(nq/2) - do ii=nq/2,np-1 - if (D_sorted(ii) < Dmin) then - nq = ii - exit - endif - enddo - else - exit - endif - - enddo -!call print_memory_usage -!print *, 'np, nq, Predicted memory: ', np, nq, mem - - if (nq <= 0) then - print *, nq - stop 'bug in cholesky: nq <= 0' - endif - - Dmin = D_sorted(nq) - nq=0 - do p=1,np - if ( D(Lset(p)) >= Dmin ) then - nq = nq+1 - Dset(nq) = Lset(p) - endif - enddo - - - allocate(Delta(np,nq)) - allocate(Ltmp_p(np,block_size), stat=ierr) - - if (ierr /= 0) then - call print_memory_usage() - print *, irp_here, ': allocation failed : (Ltmp_p(np,block_size))' - stop -1 - endif - - allocate(Ltmp_q(nq,block_size), stat=ierr) - - if (ierr /= 0) then - call print_memory_usage() - print *, irp_here, ': allocation failed : (Ltmp_q(nq,block_size))' - stop -1 - endif - - - computed(1:nq) = .False. - - - !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(k,p,q) - do k=1,N - !$OMP DO - do p=1,np - Ltmp_p(p,k) = L(Lset(p),k) - enddo - !$OMP END DO NOWAIT - - !$OMP DO - do q=1,nq - Ltmp_q(q,k) = L(Dset(q),k) - enddo - !$OMP END DO NOWAIT - enddo - !$OMP BARRIER - !$OMP END PARALLEL - - if (N>0) then - - call dgemm('N', 'T', np, nq, N, -1.d0, & - Ltmp_p(1,1), np, Ltmp_q(1,1), nq, 0.d0, Delta, np) - - else - - !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(q,j) - do q=1,nq - Delta(:,q) = 0.d0 - enddo - !$OMP END PARALLEL DO - - endif - - ! f. - Qmax = D(Dset(1)) - do q=1,nq - Qmax = max(Qmax, D(Dset(q))) - enddo - - if (Qmax < Dmin) exit - - ! g. - - iblock = 0 - - do j=1,nq - - if ( (Qmax < Dmin).or.(N+j*1_8 > ndim8) ) exit - - ! i. - rank = N+j - if (rank == rank_max) then - print *, 'cholesky: rank_max reached' - exit - endif - - if (iblock == block_size) then - - call dgemm('N','T',np,nq,block_size,-1.d0, & - Ltmp_p, np, Ltmp_q, nq, 1.d0, Delta, np) - - iblock = 0 - - endif - - ! ii. - do dj=1,nq - qj8 = Dset(dj) - if (D(qj8) == Qmax) then - exit - endif - enddo - - do i8=1,ndim8 - L(i8, rank) = 0.d0 - enddo - - iblock = iblock+1 - !$OMP PARALLEL DO PRIVATE(p) - do p=1,np - Ltmp_p(p,iblock) = Delta(p,dj) - enddo - !$OMP END PARALLEL DO - - if (.not.computed(dj)) then - m = dj - if (do_direct_integrals) then - !$OMP PARALLEL DO PRIVATE(k) SCHEDULE(dynamic,21) - do k=1,np - Delta_col(k) = 0.d0 - if (.not.ao_two_e_integral_zero( addr1(Lset(k)), addr1(Dset(m)),& - addr2(Lset(k)), addr2(Dset(m)) ) ) then - Delta_col(k) = & - ao_two_e_integral(addr1(Lset(k)), addr2(Lset(k)),& - addr1(Dset(m)), addr2(Dset(m))) - endif - enddo - !$OMP END PARALLEL DO - else - PROVIDE ao_integrals_map - !$OMP PARALLEL DO PRIVATE(k) SCHEDULE(dynamic,21) - do k=1,np - Delta_col(k) = 0.d0 - if (.not.ao_two_e_integral_zero( addr1(Lset(k)), addr1(Dset(m)),& - addr2(Lset(k)), addr2(Dset(m)) ) ) then - Delta_col(k) = & - get_ao_two_e_integral( addr1(Lset(k)), addr1(Dset(m)),& - addr2(Lset(k)), addr2(Dset(m)), ao_integrals_map) - endif - enddo - !$OMP END PARALLEL DO - endif - - !$OMP PARALLEL DO PRIVATE(p) - do p=1,np - Ltmp_p(p,iblock) = Ltmp_p(p,iblock) + Delta_col(p) - Delta(p,dj) = Ltmp_p(p,iblock) - enddo - !$OMP END PARALLEL DO - - computed(dj) = .True. - endif - - ! iv. - if (iblock > 1) then - call dgemv('N', np, iblock-1, -1.d0, Ltmp_p, np, Ltmp_q(dj,1), nq, 1.d0,& - Ltmp_p(1,iblock), 1) - endif - - ! iii. - f = 1.d0/dsqrt(Qmax) - - !$OMP PARALLEL PRIVATE(p,q) DEFAULT(shared) - !$OMP DO - do p=1,np - Ltmp_p(p,iblock) = Ltmp_p(p,iblock) * f - L(Lset(p), rank) = Ltmp_p(p,iblock) - D(Lset(p)) = D(Lset(p)) - Ltmp_p(p,iblock) * Ltmp_p(p,iblock) - enddo - !$OMP END DO - - !$OMP DO - do q=1,nq - Ltmp_q(q,iblock) = L(Dset(q), rank) - enddo - !$OMP END DO - !$OMP END PARALLEL - - Qmax = D(Dset(1)) - do q=1,nq - Qmax = max(Qmax, D(Dset(q))) - enddo - - enddo - - print '(I10, 4X, ES12.3)', rank, Qmax - - deallocate(Ltmp_p) - deallocate(Ltmp_q) - deallocate(Delta) - - ! i. - N = rank - - ! j. - D_sorted(:) = -D(:) - call dsort_noidx_big(D_sorted,ndim8) - D_sorted(:) = -D_sorted(:) - - Dmax = D_sorted(1) - - np8=0_8 - do p8=1,ndim8 - if ( Dmax*D(p8) >= tau2 ) then - np8 = np8+1_8 - Lset(np8) = p8 - endif - enddo - np = int(np8,4) - - enddo - - - print *, '============ =============' - print *, '' - - deallocate( D, Lset, Dset, D_sorted ) - deallocate( addr1, addr2, Delta_col, computed ) - - - allocate(cholesky_ao(ao_num,ao_num,rank), stat=ierr) + call wall_time(wall0) + allocate(X(ao_num,cholesky_ao$_erf_num,ao_cart_num), stat=ierr) if (ierr /= 0) then - call print_memory_usage() - print *, irp_here, ': Allocation failed' - stop -1 + print *, irp_here, ': Allocation failed' endif + call dgemm('T','N', ao_cart_num*cholesky_ao$_erf_num, ao_num, ao_cart_num, 1.d0, & + cholesky_ao$_erf_cart, ao_cart_num, ao_cart_to_ao_basis_mat_transp, ao_cart_num, & + 0.d0, X, ao_cart_num*cholesky_ao$_erf_num) + call dgemm('T','N', cholesky_ao$_erf_num*ao_num, ao_num, ao_cart_num, 1.d0, & + X, ao_cart_num, ao_cart_to_ao_basis_mat_transp, ao_cart_num, 0.d0, & + cholesky_ao$_erf_transp, cholesky_ao$_erf_num*ao_num) + deallocate(X) + call wall_time(wall1) + print*,'Time to provide ao$_erf cholesky vectors = ',(wall1-wall0)/60.d0, ' min' - ! Reverse order of Cholesky vectors to increase precision in dot products - !$OMP PARALLEL DO PRIVATE(k,j) - do k=1,rank - do j=1,ao_num - cholesky_ao(1:ao_num,j,k) = L((j-1_8)*ao_num+1_8:1_8*j*ao_num,rank-k+1) - enddo - enddo - !$OMP END PARALLEL DO - - call mmap_destroy(map) - - cholesky_ao_num = rank - - if (write_ao_cholesky) then - print *, 'Writing Cholesky AO vectors to disk...' - iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao', 'W') - write(iunit) rank - write(iunit) cholesky_ao + if (write_ao$_erf_cholesky) then + print *, 'Writing Cholesky ao$_erf vectors to disk...' + iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao$_erf_transp', 'W') + write(iunit) cholesky_ao$_erf_num + write(iunit) cholesky_ao$_erf_transp close(iunit) - call ezfio_set_ao_two_e_ints_io_ao_cholesky('Read') + call ezfio_set_ao_two_e_ints_io_ao$_erf_cholesky('Read') endif - - endif - - print *, 'Rank : ', cholesky_ao_num, '(', 100.d0*dble(cholesky_ao_num)/dble(ao_num*ao_num), ' %)' - print *, '' - call wall_time(wall1) - print*,'Time to provide AO cholesky vectors = ',(wall1-wall0)/60.d0, ' min' - + endif END_PROVIDER + +BEGIN_PROVIDER [ real, cholesky_ao$_erf_sp, (ao_num, ao_num, cholesky_ao$_erf_num) ] + implicit none + BEGIN_DOC + ! Cholesky vectors in ao$_erf basis in stored in single precision + END_DOC + + integer :: k, i, j + + call set_multiple_levels_omp(.False.) + !$OMP PARALLEL DO PRIVATE(k) + do k=1,cholesky_ao$_erf_num + do j=1,ao_num + do i=1,ao_num + cholesky_ao$_erf_sp(i,j,k) = cholesky_ao$_erf_transp_sp(k,i,j) + enddo + enddo + enddo + !$OMP END PARALLEL DO + +END_PROVIDER + +BEGIN_PROVIDER [ real, cholesky_ao$_erf_transp_sp, (cholesky_ao$_erf_num, ao_num, ao_num) ] + implicit none + BEGIN_DOC + ! Cholesky vectors in ao$_erf basis in s. Warning: it is transposed wrt cholesky_ao$_erf_cart: + ! + ! - cholesky_ao$_erf_cart is (ao_cart_num^2 x cholesky_ao$_erf_cart_num) + ! + ! - cholesky_ao$_erf_transp is (cholesky_ao$_erf_num x ao_num^2) + END_DOC + + integer :: i,j,k + !$OMP PARALLEL DO PRIVATE(k) + do j=1,ao_num + do i=1,ao_num + do k=1,cholesky_ao$_erf_num + cholesky_ao$_erf_transp_sp(k,i,j) = cholesky_ao$_erf_transp(k,i,j) + enddo + enddo + enddo + !$OMP END PARALLEL DO + +END_PROVIDER + +SUBST [ _erf ] + +;; +_erf;; +_cgtos;; + +END_TEMPLATE diff --git a/src/ao_two_e_ints/providers_ao_erf.irp.f b/src/ao_two_e_ints/providers_ao_erf.irp.f deleted file mode 100644 index 8468fb0a..00000000 --- a/src/ao_two_e_ints/providers_ao_erf.irp.f +++ /dev/null @@ -1,112 +0,0 @@ -BEGIN_PROVIDER [ logical, ao_two_e_integrals_erf_in_map ] - implicit none - use map_module - BEGIN_DOC - ! Map of Atomic integrals - ! i(r1) j(r2) 1/r12 k(r1) l(r2) - END_DOC - - integer :: i,j,k,l - double precision :: ao_two_e_integral_erf,cpu_1,cpu_2, wall_1, wall_2 - double precision :: integral, wall_0 - include 'utils/constants.include.F' - - ! For integrals file - integer(key_kind),allocatable :: buffer_i(:) - integer :: size_buffer - real(integral_kind),allocatable :: buffer_value(:) - - integer :: n_integrals, rc - integer :: kk, m, j1, i1, lmax - character*(64) :: fmt - - double precision :: map_mb - PROVIDE read_ao_two_e_integrals_erf io_ao_two_e_integrals_erf ao_integrals_erf_map - - if (read_ao_two_e_integrals_erf) then - print*,'Reading the AO ERF integrals' - call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints_erf',ao_integrals_erf_map) - print*, 'AO ERF integrals provided' - ao_two_e_integrals_erf_in_map = .True. - return - endif - - print*, 'Providing the AO ERF integrals' - call wall_time(wall_0) - call wall_time(wall_1) - call cpu_time(cpu_1) - - if (.True.) then - ! Avoid openMP - integral = ao_two_e_integral_erf(1,1,1,1) - endif - - size_buffer = ao_num*ao_num - !$OMP PARALLEL DEFAULT(shared) private(j,l) & - !$OMP PRIVATE(buffer_i, buffer_value, n_integrals) - allocate(buffer_i(size_buffer), buffer_value(size_buffer)) - n_integrals = 0 - !$OMP DO COLLAPSE(1) SCHEDULE(dynamic) - do l=1,ao_num - do j=1,l - call compute_ao_integrals_erf_jl(j,l,n_integrals,buffer_i,buffer_value) - call insert_into_ao_integrals_erf_map(n_integrals,buffer_i,buffer_value) - enddo - enddo - !$OMP END DO - deallocate(buffer_i, buffer_value) - !$OMP END PARALLEL - - - - print*, 'Sorting the map' - call map_sort(ao_integrals_erf_map) - call cpu_time(cpu_2) - call wall_time(wall_2) - integer(map_size_kind) :: get_ao_erf_map_size, ao_erf_map_size - ao_erf_map_size = get_ao_erf_map_size() - - print*, 'AO ERF integrals provided:' - print*, ' Size of AO ERF map : ', map_mb(ao_integrals_erf_map) ,'MB' - print*, ' Number of AO ERF integrals :', ao_erf_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_erf_in_map = .True. - - if (write_ao_two_e_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_ao_two_e_ints_io_ao_two_e_integrals_erf('Read') - endif - -END_PROVIDER - - - - -BEGIN_PROVIDER [ double precision, ao_two_e_integral_erf_schwartz,(ao_num,ao_num) ] - implicit none - BEGIN_DOC - ! Needed to compute Schwartz inequalities - END_DOC - - integer :: i,k - double precision :: ao_two_e_integral_erf,cpu_1,cpu_2, wall_1, wall_2 - - ao_two_e_integral_erf_schwartz(1,1) = ao_two_e_integral_erf(1,1,1,1) - !$OMP PARALLEL DO PRIVATE(i,k) & - !$OMP DEFAULT(NONE) & - !$OMP SHARED (ao_num,ao_two_e_integral_erf_schwartz) & - !$OMP SCHEDULE(dynamic) - do i=1,ao_num - do k=1,i - ao_two_e_integral_erf_schwartz(i,k) = dsqrt(ao_two_e_integral_erf(i,k,i,k)) - ao_two_e_integral_erf_schwartz(k,i) = ao_two_e_integral_erf_schwartz(i,k) - enddo - enddo - !$OMP END PARALLEL DO - -END_PROVIDER - - diff --git a/src/ao_two_e_ints/routines_save_integrals_erf.irp.f b/src/ao_two_e_ints/routines_save_integrals_erf.irp.f index d98d96ce..2bf20a78 100644 --- a/src/ao_two_e_ints/routines_save_integrals_erf.irp.f +++ b/src/ao_two_e_ints/routines_save_integrals_erf.irp.f @@ -12,6 +12,7 @@ subroutine save_erf_two_e_ints_ao_into_ints_ao integer :: i,j,k,l PROVIDE ao_two_e_integrals_erf_in_map call ezfio_set_work_empty(.False.) + ! intentionally ao_ints in order to trick the reading call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_erf_map) call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read') call ezfio_set_ao_two_e_ints_do_ao_cholesky('Read') diff --git a/src/ao_two_e_ints/two_e_coul_integrals_cgtos.irp.f b/src/ao_two_e_ints/two_e_coul_integrals_cgtos.irp.f deleted file mode 100644 index 30e4824a..00000000 --- a/src/ao_two_e_ints/two_e_coul_integrals_cgtos.irp.f +++ /dev/null @@ -1,1674 +0,0 @@ - -! --- - -double precision function ao_two_e_integral_cgtos(i, j, k, l) - - BEGIN_DOC - ! integral of the AO basis or (ij|kl) - ! i(r1) j(r1) 1/r12 k(r2) l(r2) - END_DOC - - implicit none - include 'utils/constants.include.F' - - integer, intent(in) :: i, j, k, l - - integer :: p, q, r, s, m - integer :: ii, jj, kk, ll, dim1, I_power(3), J_power(3), K_power(3), L_power(3) - integer :: iorder_p1(3), iorder_p2(3), iorder_q1(3), iorder_q2(3) - double precision :: coef1, coef2, coef3, coef4 - double precision :: KI2, phiI - double precision :: KJ2, phiJ - double precision :: KK2, phiK - double precision :: KL2, phiL - complex*16 :: expo1, expo1_inv, Ie_center(3), Ip_center(3) - complex*16 :: expo2, expo2_inv, Je_center(3), Jp_center(3) - complex*16 :: expo3, expo3_inv, Ke_center(3), Kp_center(3) - complex*16 :: expo4, expo4_inv, Le_center(3), Lp_center(3) - complex*16 :: P1_new(0:max_dim,3), P1_center(3), fact_p1, pp1, p1_inv - complex*16 :: P2_new(0:max_dim,3), P2_center(3), fact_p2, pp2, p2_inv - complex*16 :: Q1_new(0:max_dim,3), Q1_center(3), fact_q1, qq1, q1_inv - complex*16 :: Q2_new(0:max_dim,3), Q2_center(3), fact_q2, qq2, q2_inv - complex*16 :: int1, int2, int3, int4 - complex*16 :: int5, int6, int7, int8 - complex*16 :: C1, C2, C3, C4, C5, C6, C7, C8 - complex*16 :: int_tot - - double precision, external :: ao_2e_cgtos_schwartz_accel - complex*16, external :: ERI_cgtos - complex*16, external :: general_primitive_integral_cgtos - - - - if(ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024) then - - ao_two_e_integral_cgtos = ao_2e_cgtos_schwartz_accel(i, j, k, l) - - return - endif - - dim1 = n_pt_max_integrals - - ii = ao_nucl(i) - jj = ao_nucl(j) - kk = ao_nucl(k) - ll = ao_nucl(l) - - do m = 1, 3 - I_power(m) = ao_power(i,m) - J_power(m) = ao_power(j,m) - K_power(m) = ao_power(k,m) - L_power(m) = ao_power(l,m) - enddo - - - ao_two_e_integral_cgtos = 0.d0 - - if(use_pw .or. ii /= jj .or. kk /= ll .or. jj /= kk) then - - do p = 1, ao_prim_num(i) - - coef1 = ao_coef_cgtos_norm_ord_transp(p,i) - expo1 = ao_expo_cgtos_ord_transp(p,i) - expo1_inv = (1.d0, 0.d0) / expo1 - do m = 1, 3 - Ip_center(m) = nucl_coord(ii,m) - Ie_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,p,i) - enddo - phiI = ao_expo_phase_ord_transp(4,p,i) - KI2 = ao_expo_pw_ord_transp(4,p,i) - - do q = 1, ao_prim_num(j) - - coef2 = coef1 * ao_coef_cgtos_norm_ord_transp(q,j) - expo2 = ao_expo_cgtos_ord_transp(q,j) - expo2_inv = (1.d0, 0.d0) / expo2 - do m = 1, 3 - Jp_center(m) = nucl_coord(jj,m) - Je_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,q,j) - enddo - phiJ = ao_expo_phase_ord_transp(4,q,j) - KJ2 = ao_expo_pw_ord_transp(4,q,j) - - call give_explicit_cpoly_and_cgaussian(P1_new, P1_center, pp1, fact_p1, iorder_p1, & - expo1, expo2, I_power, J_power, Ie_center, Je_center, Ip_center, Jp_center, dim1) - - p1_inv = (1.d0, 0.d0) / pp1 - - call give_explicit_cpoly_and_cgaussian(P2_new, P2_center, pp2, fact_p2, iorder_p2, & - conjg(expo1), expo2, I_power, J_power, conjg(Ie_center), Je_center, conjg(Ip_center), Jp_center, dim1) - - p2_inv = (1.d0, 0.d0) / pp2 - - do r = 1, ao_prim_num(k) - - coef3 = coef2 * ao_coef_cgtos_norm_ord_transp(r,k) - expo3 = ao_expo_cgtos_ord_transp(r,k) - expo3_inv = (1.d0, 0.d0) / expo3 - do m = 1, 3 - Kp_center(m) = nucl_coord(kk,m) - Ke_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo3_inv * ao_expo_pw_ord_transp(m,r,k) - enddo - phiK = ao_expo_phase_ord_transp(4,r,k) - KK2 = ao_expo_pw_ord_transp(4,r,k) - - do s = 1, ao_prim_num(l) - - coef4 = coef3 * ao_coef_cgtos_norm_ord_transp(s,l) - expo4 = ao_expo_cgtos_ord_transp(s,l) - expo4_inv = (1.d0, 0.d0) / expo4 - do m = 1, 3 - Lp_center(m) = nucl_coord(ll,m) - Le_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo4_inv * ao_expo_pw_ord_transp(m,s,l) - enddo - phiL = ao_expo_phase_ord_transp(4,s,l) - KL2 = ao_expo_pw_ord_transp(4,s,l) - - call give_explicit_cpoly_and_cgaussian(Q1_new, Q1_center, qq1, fact_q1, iorder_q1, & - expo3, expo4, K_power, L_power, Ke_center, Le_center, Kp_center, Lp_center, dim1) - - q1_inv = (1.d0, 0.d0) / qq1 - - call give_explicit_cpoly_and_cgaussian(Q2_new, Q2_center, qq2, fact_q2, iorder_q2, & - conjg(expo3), expo4, K_power, L_power, conjg(Ke_center), Le_center, conjg(Kp_center), Lp_center, dim1) - - q2_inv = (1.d0, 0.d0) / qq2 - - C1 = zexp((0.d0, 1.d0) * (-phiI - phiJ - phiK - phiL) & - - 0.25d0 * (expo1_inv * KI2 + expo2_inv * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) - C2 = zexp((0.d0, 1.d0) * (-phiI - phiJ + phiK - phiL) & - - 0.25d0 * (expo1_inv * KI2 + expo2_inv * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) - C3 = zexp((0.d0, 1.d0) * ( phiI - phiJ - phiK - phiL) & - - 0.25d0 * (conjg(expo1_inv) * KI2 + expo2_inv * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) - C4 = zexp((0.d0, 1.d0) * ( phiI - phiJ + phiK - phiL) & - - 0.25d0 * (conjg(expo1_inv) * KI2 + expo2_inv * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) - C5 = zexp((0.d0, 1.d0) * (-phiI + phiJ - phiK - phiL) & - - 0.25d0 * (expo1_inv * KI2 + conjg(expo2_inv) * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) - C6 = zexp((0.d0, 1.d0) * (-phiI + phiJ + phiK - phiL) & - - 0.25d0 * (expo1_inv * KI2 + conjg(expo2_inv) * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) - C7 = zexp((0.d0, 1.d0) * ( phiI + phiJ - phiK - phiL) & - - 0.25d0 * (conjg(expo1_inv) * KI2 + conjg(expo2_inv) * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) - C8 = zexp((0.d0, 1.d0) * ( phiI + phiJ + phiK - phiL) & - - 0.25d0 * (conjg(expo1_inv) * KI2 + conjg(expo2_inv) * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) - - int1 = general_primitive_integral_cgtos(dim1, & - P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & - Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) - - int2 = general_primitive_integral_cgtos(dim1, & - P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & - Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) - - int3 = general_primitive_integral_cgtos(dim1, & - P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & - Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) - - int4 = general_primitive_integral_cgtos(dim1, & - P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & - Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) - - int5 = general_primitive_integral_cgtos(dim1, & - conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, & - Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) - - int6 = general_primitive_integral_cgtos(dim1, & - conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, & - Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) - - int7 = general_primitive_integral_cgtos(dim1, & - conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, & - Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) - - int8 = general_primitive_integral_cgtos(dim1, & - conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, & - Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) - - int_tot = C1 * int1 + C2 * int2 + C3 * int3 + C4 * int4 + C5 * int5 + C6 * int6 + C7 * int7 + C8 * int8 - - ao_two_e_integral_cgtos = ao_two_e_integral_cgtos + coef4 * 2.d0 * real(int_tot) - enddo ! s - enddo ! r - enddo ! q - enddo ! p - - else - - do p = 1, ao_prim_num(i) - - coef1 = ao_coef_cgtos_norm_ord_transp(p,i) - expo1 = ao_expo_cgtos_ord_transp(p,i) - phiI = ao_expo_phase_ord_transp(4,p,i) - - do q = 1, ao_prim_num(j) - - coef2 = coef1 * ao_coef_cgtos_norm_ord_transp(q,j) - expo2 = ao_expo_cgtos_ord_transp(q,j) - phiJ = ao_expo_phase_ord_transp(4,q,j) - - do r = 1, ao_prim_num(k) - - coef3 = coef2 * ao_coef_cgtos_norm_ord_transp(r,k) - expo3 = ao_expo_cgtos_ord_transp(r,k) - phiK = ao_expo_phase_ord_transp(4,r,k) - - do s = 1, ao_prim_num(l) - - coef4 = coef3 * ao_coef_cgtos_norm_ord_transp(s,l) - expo4 = ao_expo_cgtos_ord_transp(s,l) - phiL = ao_expo_phase_ord_transp(4,s,l) - - C1 = zexp((0.d0, 1.d0) * (-phiI - phiJ - phiK - phiL)) - C2 = zexp((0.d0, 1.d0) * (-phiI - phiJ + phiK - phiL)) - C3 = zexp((0.d0, 1.d0) * ( phiI - phiJ - phiK - phiL)) - C4 = zexp((0.d0, 1.d0) * ( phiI - phiJ + phiK - phiL)) - C5 = zexp((0.d0, 1.d0) * (-phiI + phiJ - phiK - phiL)) - C6 = zexp((0.d0, 1.d0) * (-phiI + phiJ + phiK - phiL)) - C7 = zexp((0.d0, 1.d0) * ( phiI + phiJ - phiK - phiL)) - C8 = zexp((0.d0, 1.d0) * ( phiI + phiJ + phiK - phiL)) - - int1 = ERI_cgtos(expo1, expo2, expo3, expo4, & - I_power(1), J_power(1), K_power(1), L_power(1), & - I_power(2), J_power(2), K_power(2), L_power(2), & - I_power(3), J_power(3), K_power(3), L_power(3)) - - int2 = ERI_cgtos(expo1, expo2, conjg(expo3), expo4, & - I_power(1), J_power(1), K_power(1), L_power(1), & - I_power(2), J_power(2), K_power(2), L_power(2), & - I_power(3), J_power(3), K_power(3), L_power(3)) - - int3 = ERI_cgtos(conjg(expo1), expo2, expo3, expo4, & - I_power(1), J_power(1), K_power(1), L_power(1), & - I_power(2), J_power(2), K_power(2), L_power(2), & - I_power(3), J_power(3), K_power(3), L_power(3)) - - int4 = ERI_cgtos(conjg(expo1), expo2, conjg(expo3), expo4, & - I_power(1), J_power(1), K_power(1), L_power(1), & - I_power(2), J_power(2), K_power(2), L_power(2), & - I_power(3), J_power(3), K_power(3), L_power(3)) - - int5 = ERI_cgtos(expo1, conjg(expo2), expo3, expo4, & - I_power(1), J_power(1), K_power(1), L_power(1), & - I_power(2), J_power(2), K_power(2), L_power(2), & - I_power(3), J_power(3), K_power(3), L_power(3)) - - int6 = ERI_cgtos(expo1, conjg(expo2), conjg(expo3), expo4, & - I_power(1), J_power(1), K_power(1), L_power(1), & - I_power(2), J_power(2), K_power(2), L_power(2), & - I_power(3), J_power(3), K_power(3), L_power(3)) - - int7 = ERI_cgtos(conjg(expo1), conjg(expo2), expo3, expo4, & - I_power(1), J_power(1), K_power(1), L_power(1), & - I_power(2), J_power(2), K_power(2), L_power(2), & - I_power(3), J_power(3), K_power(3), L_power(3)) - - int8 = ERI_cgtos(conjg(expo1), conjg(expo2), conjg(expo3), expo4, & - I_power(1), J_power(1), K_power(1), L_power(1), & - I_power(2), J_power(2), K_power(2), L_power(2), & - I_power(3), J_power(3), K_power(3), L_power(3)) - - int_tot = C1 * int1 + C2 * int2 + C3 * int3 + C4 * int4 & - + C5 * int5 + C6 * int6 + C7 * int7 + C8 * int8 - - ao_two_e_integral_cgtos = ao_two_e_integral_cgtos + coef4 * 2.d0 * real(int_tot) - enddo ! s - enddo ! r - enddo ! q - enddo ! p - - endif - -end - -! --- - -double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) - - BEGIN_DOC - ! integral of the AO basis or (ij|kl) - ! i(r1) j(r1) 1/r12 k(r2) l(r2) - END_DOC - - implicit none - include 'utils/constants.include.F' - - integer, intent(in) :: i, j, k, l - - integer :: p, q, r, s, m - integer :: ii, jj, kk, ll, dim1, I_power(3), J_power(3), K_power(3), L_power(3) - integer :: iorder_p1(3), iorder_p2(3), iorder_q1(3), iorder_q2(3) - double precision :: coef1, coef2, coef3, coef4 - double precision :: KI2, phiI - double precision :: KJ2, phiJ - double precision :: KK2, phiK - double precision :: KL2, phiL - complex*16 :: expo1, expo1_inv, Ie_center(3), Ip_center(3) - complex*16 :: expo2, expo2_inv, Je_center(3), Jp_center(3) - complex*16 :: expo3, expo3_inv, Ke_center(3), Kp_center(3) - complex*16 :: expo4, expo4_inv, Le_center(3), Lp_center(3) - complex*16 :: P1_new(0:max_dim,3), P1_center(3), fact_p1, pp1, p1_inv - complex*16 :: P2_new(0:max_dim,3), P2_center(3), fact_p2, pp2, p2_inv - complex*16 :: Q1_new(0:max_dim,3), Q1_center(3), fact_q1, qq1, q1_inv - complex*16 :: Q2_new(0:max_dim,3), Q2_center(3), fact_q2, qq2, q2_inv - complex*16 :: int1, int2, int3, int4 - complex*16 :: int5, int6, int7, int8 - complex*16 :: C1, C2, C3, C4, C5, C6, C7, C8 - complex*16 :: int_tot - - double precision, allocatable :: schwartz_kl(:,:) - double precision :: thr - double precision :: schwartz_ij - - complex*16, external :: ERI_cgtos - complex*16, external :: general_primitive_integral_cgtos - - ao_2e_cgtos_schwartz_accel = 0.d0 - - dim1 = n_pt_max_integrals - - ii = ao_nucl(i) - jj = ao_nucl(j) - kk = ao_nucl(k) - ll = ao_nucl(l) - - do m = 1, 3 - I_power(m) = ao_power(i,m) - J_power(m) = ao_power(j,m) - K_power(m) = ao_power(k,m) - L_power(m) = ao_power(l,m) - enddo - - - thr = ao_integrals_threshold*ao_integrals_threshold - - allocate(schwartz_kl(0:ao_prim_num(l),0:ao_prim_num(k))) - - if(use_pw .or. ii /= jj .or. kk /= ll .or. jj /= kk) then - - schwartz_kl(0,0) = 0.d0 - do r = 1, ao_prim_num(k) - - coef1 = ao_coef_cgtos_norm_ord_transp(r,k) * ao_coef_cgtos_norm_ord_transp(r,k) - expo1 = ao_expo_cgtos_ord_transp(r,k) - expo1_inv = (1.d0, 0.d0) / expo1 - do m = 1, 3 - Kp_center(m) = nucl_coord(kk,m) - Ke_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,r,k) - enddo - phiK = ao_expo_phase_ord_transp(4,r,k) - KK2 = ao_expo_pw_ord_transp(4,r,k) - - schwartz_kl(0,r) = 0.d0 - do s = 1, ao_prim_num(l) - - coef2 = coef1 * ao_coef_cgtos_norm_ord_transp(s,l) * ao_coef_cgtos_norm_ord_transp(s,l) - expo2 = ao_expo_cgtos_ord_transp(s,l) - expo2_inv = (1.d0, 0.d0) / expo2 - do m = 1, 3 - Lp_center(m) = nucl_coord(ll,m) - Le_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,s,l) - enddo - phiL = ao_expo_phase_ord_transp(4,s,l) - KL2 = ao_expo_pw_ord_transp(4,s,l) - - call give_explicit_cpoly_and_cgaussian(P1_new, P1_center, pp1, fact_p1, iorder_p1, & - expo1, expo2, K_power, L_power, Ke_center, Le_center, Kp_center, Lp_center, dim1) - - p1_inv = (1.d0, 0.d0) / pp1 - - call give_explicit_cpoly_and_cgaussian(P2_new, P2_center, pp2, fact_p2, iorder_p2, & - conjg(expo1), expo2, K_power, L_power, conjg(Ke_center), Le_center, conjg(Kp_center), Lp_center, dim1) - - p2_inv = (1.d0, 0.d0) / pp2 - - C1 = zexp(-(0.d0, 2.d0) * (phiK + phiL) - 0.5d0 * (expo1_inv * KK2 + expo2_inv * KL2)) - C2 = zexp(-(0.d0, 2.d0) * phiL - 0.5d0 * (real(expo1_inv) * KK2 + expo2_inv * KL2)) - !C3 = C2 - C4 = zexp((0.d0, 2.d0) * (phiK - phiL) - 0.5d0 * (conjg(expo1_inv) * KK2 + expo2_inv * KL2)) - C5 = zexp(-(0.d0, 2.d0) * phiK - 0.5d0 * (expo1_inv * KK2 + real(expo2_inv) * KL2)) - C6 = zexp(-(0.5d0, 0.d0) * (real(expo1_inv) * KK2 + real(expo2_inv) * KL2)) - !C7 = C6 - !C8 = conjg(C5) - - int1 = general_primitive_integral_cgtos(dim1, & - P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & - P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) - - int2 = general_primitive_integral_cgtos(dim1, & - P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & - P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) - - ! int3 = int2 - !int3 = general_primitive_integral_cgtos(dim1, & - ! P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & - ! P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) - - int4 = general_primitive_integral_cgtos(dim1, & - P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & - P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) - - int5 = general_primitive_integral_cgtos(dim1, & - conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, & - P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) - - int6 = general_primitive_integral_cgtos(dim1, & - conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, & - P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) - - int7 = general_primitive_integral_cgtos(dim1, & - conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, & - P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) - - !int8 = conjg(int5) - !int8 = general_primitive_integral_cgtos(dim1, & - ! conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, & - ! P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) - - !int_tot = C1 * int1 + C2 * int2 + C3 * int3 + C4 * int4 + C5 * int5 + C6 * int6 + C7 * int7 + C8 * int8 - int_tot = C1 * int1 + 2.d0 * C2 * int2 + C4 * int4 + 2.d0 * real(C5 * int5) + C6 * (int6 + int7) - - schwartz_kl(s,r) = coef2 * 2.d0 * real(int_tot) - - schwartz_kl(0,r) = max(schwartz_kl(0,r), schwartz_kl(s,r)) - enddo - - schwartz_kl(0,0) = max(schwartz_kl(0,r), schwartz_kl(0,0)) - enddo - - - do p = 1, ao_prim_num(i) - - coef1 = ao_coef_cgtos_norm_ord_transp(p,i) - expo1 = ao_expo_cgtos_ord_transp(p,i) - expo1_inv = (1.d0, 0.d0) / expo1 - do m = 1, 3 - Ip_center(m) = nucl_coord(ii,m) - Ie_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,p,i) - enddo - phiI = ao_expo_phase_ord_transp(4,p,i) - KI2 = ao_expo_pw_ord_transp(4,p,i) - - do q = 1, ao_prim_num(j) - - coef2 = coef1 * ao_coef_cgtos_norm_ord_transp(q,j) - expo2 = ao_expo_cgtos_ord_transp(q,j) - expo2_inv = (1.d0, 0.d0) / expo2 - do m = 1, 3 - Jp_center(m) = nucl_coord(jj,m) - Je_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,q,j) - enddo - phiJ = ao_expo_phase_ord_transp(4,q,j) - KJ2 = ao_expo_pw_ord_transp(4,q,j) - - call give_explicit_cpoly_and_cgaussian(P1_new, P1_center, pp1, fact_p1, iorder_p1, & - expo1, expo2, I_power, J_power, Ie_center, Je_center, Ip_center, Jp_center, dim1) - - p1_inv = (1.d0, 0.d0) / pp1 - - call give_explicit_cpoly_and_cgaussian(P2_new, P2_center, pp2, fact_p2, iorder_p2, & - conjg(expo1), expo2, I_power, J_power, conjg(Ie_center), Je_center, conjg(Ip_center), Jp_center, dim1) - p2_inv = (1.d0, 0.d0) / pp2 - - C1 = zexp(-(0.d0, 2.d0) * (phiI + phiJ) - 0.5d0 * (expo1_inv * KI2 + expo2_inv * KJ2)) - C2 = zexp(-(0.d0, 2.d0) * phiJ - 0.5d0 * (real(expo1_inv) * KI2 + expo2_inv * KJ2)) - !C3 = C2 - C4 = zexp((0.d0, 2.d0) * (phiI - phiJ) - 0.5d0 * (conjg(expo1_inv) * KI2 + expo2_inv * KJ2)) - C5 = zexp(-(0.d0, 2.d0) * phiI - 0.5d0 * (expo1_inv * KI2 + real(expo2_inv) * KJ2)) - C6 = zexp(-(0.5d0, 0.d0) * (real(expo1_inv) * KI2 + real(expo2_inv) * KJ2)) - !C7 = C6 - !C8 = conjg(C5) - - int1 = general_primitive_integral_cgtos(dim1, & - P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & - P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) - - int2 = general_primitive_integral_cgtos(dim1, & - P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & - P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) - - !int3 = int1 - !int3 = general_primitive_integral_cgtos(dim1, & - ! P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & - ! P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) - - int4 = general_primitive_integral_cgtos(dim1, & - P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & - P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) - - int5 = general_primitive_integral_cgtos(dim1, & - conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, & - P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) - - int6 = general_primitive_integral_cgtos(dim1, & - conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, & - P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) - - int7 = general_primitive_integral_cgtos(dim1, & - conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, & - P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) - - !int8 = conjg(I5) - !int8 = general_primitive_integral_cgtos(dim1, & - ! conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, & - ! P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) - - !int_tot = C1 * int1 + C2 * int2 + C3 * int3 + C4 * int4 + C5 * int5 + C6 * int6 + C7 * int7 + C8 * int8 - int_tot = C1 * int1 + 2.d0 * C2 * int2 + C4 * int4 + 2.d0 * real(C5 * int5) + C6 * (int6 + int7) - - schwartz_ij = coef2 * coef2 * 2.d0 * real(int_tot) - - if(schwartz_kl(0,0)*schwartz_ij < thr) cycle - - do r = 1, ao_prim_num(k) - if(schwartz_kl(0,r)*schwartz_ij < thr) cycle - - coef3 = coef2 * ao_coef_cgtos_norm_ord_transp(r,k) - expo3 = ao_expo_cgtos_ord_transp(r,k) - expo3_inv = (1.d0, 0.d0) / expo3 - do m = 1, 3 - Kp_center(m) = nucl_coord(kk,m) - Ke_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo3_inv * ao_expo_pw_ord_transp(m,r,k) - enddo - phiK = ao_expo_phase_ord_transp(4,r,k) - KK2 = ao_expo_pw_ord_transp(4,r,k) - - do s = 1, ao_prim_num(l) - if(schwartz_kl(s,r)*schwartz_ij < thr) cycle - - coef4 = coef3 * ao_coef_cgtos_norm_ord_transp(s,l) - expo4 = ao_expo_cgtos_ord_transp(s,l) - expo4_inv = (1.d0, 0.d0) / expo4 - do m = 1, 3 - Lp_center(m) = nucl_coord(ll,m) - Le_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo4_inv * ao_expo_pw_ord_transp(m,s,l) - enddo - phiL = ao_expo_phase_ord_transp(4,s,l) - KL2 = ao_expo_pw_ord_transp(4,s,l) - - call give_explicit_cpoly_and_cgaussian(Q1_new, Q1_center, qq1, fact_q1, iorder_q1, & - expo3, expo4, K_power, L_power, Ke_center, Le_center, Kp_center, Lp_center, dim1) - - q1_inv = (1.d0, 0.d0) / qq1 - - call give_explicit_cpoly_and_cgaussian(Q2_new, Q2_center, qq2, fact_q2, iorder_q2, & - conjg(expo3), expo4, K_power, L_power, conjg(Ke_center), Le_center, conjg(Kp_center), Lp_center, dim1) - - q2_inv = (1.d0, 0.d0) / qq2 - - C1 = zexp((0.d0, 1.d0) * (-phiI - phiJ - phiK - phiL) & - - 0.25d0 * (expo1_inv * KI2 + expo2_inv * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) - C2 = zexp((0.d0, 1.d0) * (-phiI - phiJ + phiK - phiL) & - - 0.25d0 * (expo1_inv * KI2 + expo2_inv * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) - C3 = zexp((0.d0, 1.d0) * ( phiI - phiJ - phiK - phiL) & - - 0.25d0 * (conjg(expo1_inv) * KI2 + expo2_inv * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) - C4 = zexp((0.d0, 1.d0) * ( phiI - phiJ + phiK - phiL) & - - 0.25d0 * (conjg(expo1_inv) * KI2 + expo2_inv * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) - C5 = zexp((0.d0, 1.d0) * (-phiI + phiJ - phiK - phiL) & - - 0.25d0 * (expo1_inv * KI2 + conjg(expo2_inv) * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) - C6 = zexp((0.d0, 1.d0) * (-phiI + phiJ + phiK - phiL) & - - 0.25d0 * (expo1_inv * KI2 + conjg(expo2_inv) * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) - C7 = zexp((0.d0, 1.d0) * ( phiI + phiJ - phiK - phiL) & - - 0.25d0 * (conjg(expo1_inv) * KI2 + conjg(expo2_inv) * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) - C8 = zexp((0.d0, 1.d0) * ( phiI + phiJ + phiK - phiL) & - - 0.25d0 * (conjg(expo1_inv) * KI2 + conjg(expo2_inv) * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) - - int1 = general_primitive_integral_cgtos(dim1, & - P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & - Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) - - int2 = general_primitive_integral_cgtos(dim1, & - P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & - Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) - - int3 = general_primitive_integral_cgtos(dim1, & - P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & - Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) - - int4 = general_primitive_integral_cgtos(dim1, & - P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & - Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) - - int5 = general_primitive_integral_cgtos(dim1, & - conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, & - Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) - - int6 = general_primitive_integral_cgtos(dim1, & - conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, & - Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) - - int7 = general_primitive_integral_cgtos(dim1, & - conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, & - Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) - - int8 = general_primitive_integral_cgtos(dim1, & - conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, & - Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) - - int_tot = C1 * int1 + C2 * int2 + C3 * int3 + C4 * int4 + C5 * int5 + C6 * int6 + C7 * int7 + C8 * int8 - - ao_2e_cgtos_schwartz_accel = ao_2e_cgtos_schwartz_accel + coef4 * 2.d0 * real(int_tot) - enddo ! s - enddo ! r - enddo ! q - enddo ! p - - else - - schwartz_kl(0,0) = 0.d0 - do r = 1, ao_prim_num(k) - - coef1 = ao_coef_cgtos_norm_ord_transp(r,k) * ao_coef_cgtos_norm_ord_transp(r,k) - expo1 = ao_expo_cgtos_ord_transp(r,k) - phiK = ao_expo_phase_ord_transp(4,r,k) - - schwartz_kl(0,r) = 0.d0 - do s = 1, ao_prim_num(l) - - coef2 = coef1 * ao_coef_cgtos_norm_ord_transp(s,l) * ao_coef_cgtos_norm_ord_transp(s,l) - expo2 = ao_expo_cgtos_ord_transp(s,l) - phiL = ao_expo_phase_ord_transp(4,s,l) - - C1 = zexp(-(0.d0, 2.d0) * (phiK + phiL)) - C2 = zexp(-(0.d0, 2.d0) * phiL) - !C3 = C2 - C4 = zexp((0.d0, 2.d0) * (phiK - phiL)) - C5 = zexp(-(0.d0, 2.d0) * phiK) - C6 = (1.d0, 0.d0) - !C7 = C6 - !C8 = conjg(C5) - - int1 = ERI_cgtos(expo1, expo2, expo1, expo2, & - K_power(1), L_power(1), K_power(1), L_power(1), & - K_power(2), L_power(2), K_power(2), L_power(2), & - K_power(3), L_power(3), K_power(3), L_power(3)) - - int2 = ERI_cgtos(expo1, expo2, conjg(expo1), expo2, & - K_power(1), L_power(1), K_power(1), L_power(1), & - K_power(2), L_power(2), K_power(2), L_power(2), & - K_power(3), L_power(3), K_power(3), L_power(3)) - - !int3 = int2 - !int3 = ERI_cgtos(conjg(expo1), expo2, expo1, expo2, & - ! K_power(1), L_power(1), K_power(1), L_power(1), & - ! K_power(2), L_power(2), K_power(2), L_power(2), & - ! K_power(3), L_power(3), K_power(3), L_power(3)) - - int4 = ERI_cgtos(conjg(expo1), expo2, conjg(expo1), expo2, & - K_power(1), L_power(1), K_power(1), L_power(1), & - K_power(2), L_power(2), K_power(2), L_power(2), & - K_power(3), L_power(3), K_power(3), L_power(3)) - - int5 = ERI_cgtos(expo1, conjg(expo2), expo1, expo2, & - K_power(1), L_power(1), K_power(1), L_power(1), & - K_power(2), L_power(2), K_power(2), L_power(2), & - K_power(3), L_power(3), K_power(3), L_power(3)) - - int6 = ERI_cgtos(expo1, conjg(expo2), conjg(expo1), expo2, & - K_power(1), L_power(1), K_power(1), L_power(1), & - K_power(2), L_power(2), K_power(2), L_power(2), & - K_power(3), L_power(3), K_power(3), L_power(3)) - - int7 = ERI_cgtos(conjg(expo1), conjg(expo2), expo1, expo2, & - K_power(1), L_power(1), K_power(1), L_power(1), & - K_power(2), L_power(2), K_power(2), L_power(2), & - K_power(3), L_power(3), K_power(3), L_power(3)) - - int8 = conjg(int5) - !int8 = ERI_cgtos(conjg(expo1), conjg(expo2), conjg(expo1), expo2, & - ! K_power(1), L_power(1), K_power(1), L_power(1), & - ! K_power(2), L_power(2), K_power(2), L_power(2), & - ! K_power(3), L_power(3), K_power(3), L_power(3)) - - !int_tot = C1 * int1 + C2 * int2 + C3 * int3 + C4 * int4 + C5 * int5 + C6 * int6 + C7 * int7 + C8 * int8 - int_tot = C1 * int1 + 2.d0 * C2 * int2 + C4 * int4 + 2.d0 * real(C5 * int5) + C6 * (int6 + int7) - - schwartz_kl(s,r) = coef2 * 2.d0 * real(int_tot) - - schwartz_kl(0,r) = max(schwartz_kl(0,r), schwartz_kl(s,r)) - enddo - schwartz_kl(0,0) = max(schwartz_kl(0,r), schwartz_kl(0,0)) - enddo - - do p = 1, ao_prim_num(i) - - coef1 = ao_coef_cgtos_norm_ord_transp(p,i) - expo1 = ao_expo_cgtos_ord_transp(p,i) - phiI = ao_expo_phase_ord_transp(4,p,i) - - do q = 1, ao_prim_num(j) - - coef2 = coef1 * ao_coef_cgtos_norm_ord_transp(q,j) - expo2 = ao_expo_cgtos_ord_transp(q,j) - phiJ = ao_expo_phase_ord_transp(4,q,j) - - C1 = zexp(-(0.d0, 2.d0) * (phiI + phiJ)) - C2 = zexp(-(0.d0, 2.d0) * phiJ) - !C3 = C2 - C4 = zexp((0.d0, 2.d0) * (phiI - phiJ)) - C5 = zexp(-(0.d0, 2.d0) * phiI) - C6 = (1.d0, 0.d0) - !C7 = C6 - !C8 = conjg(C5) - - int1 = ERI_cgtos(expo1, expo2, expo1, expo2, & - I_power(1), J_power(1), I_power(1), J_power(1), & - I_power(2), J_power(2), I_power(2), J_power(2), & - I_power(3), J_power(3), I_power(3), J_power(3)) - - int2 = ERI_cgtos(expo1, expo2, conjg(expo1), expo2, & - I_power(1), J_power(1), I_power(1), J_power(1), & - I_power(2), J_power(2), I_power(2), J_power(2), & - I_power(3), J_power(3), I_power(3), J_power(3)) - - !int2 = int2 - !int3 = ERI_cgtos(conjg(expo1), expo2, expo1, expo2, & - ! I_power(1), J_power(1), I_power(1), J_power(1), & - ! I_power(2), J_power(2), I_power(2), J_power(2), & - ! I_power(3), J_power(3), I_power(3), J_power(3)) - - int4 = ERI_cgtos(conjg(expo1), expo2, conjg(expo1), expo2, & - I_power(1), J_power(1), I_power(1), J_power(1), & - I_power(2), J_power(2), I_power(2), J_power(2), & - I_power(3), J_power(3), I_power(3), J_power(3)) - - int5 = ERI_cgtos(expo1, conjg(expo2), expo1, expo2, & - I_power(1), J_power(1), I_power(1), J_power(1), & - I_power(2), J_power(2), I_power(2), J_power(2), & - I_power(3), J_power(3), I_power(3), J_power(3)) - - int6 = ERI_cgtos(expo1, conjg(expo2), conjg(expo1), expo2, & - I_power(1), J_power(1), I_power(1), J_power(1), & - I_power(2), J_power(2), I_power(2), J_power(2), & - I_power(3), J_power(3), I_power(3), J_power(3)) - - int7 = ERI_cgtos(conjg(expo1), conjg(expo2), expo1, expo2, & - I_power(1), J_power(1), I_power(1), J_power(1), & - I_power(2), J_power(2), I_power(2), J_power(2), & - I_power(3), J_power(3), I_power(3), J_power(3)) - - !int8 = conjg(int5) - !int8 = ERI_cgtos(conjg(expo1), conjg(expo2), conjg(expo1), expo2, & - ! I_power(1), J_power(1), I_power(1), J_power(1), & - ! I_power(2), J_power(2), I_power(2), J_power(2), & - ! I_power(3), J_power(3), I_power(3), J_power(3)) - - !int_tot = C1 * int1 + C2 * int2 + C3 * int3 + C4 * int4 + C5 * int5 + C6 * int6 + C7 * int7 + C8 * int8 - int_tot = C1 * int1 + 2.d0 * C2 * int2 + C4 * int4 + 2.d0 * real(C5 * int5) + C6 * (int6 + int7) - - schwartz_ij = coef2 * coef2 * 2.d0 * real(int_tot) - - if(schwartz_kl(0,0)*schwartz_ij < thr) cycle - do r = 1, ao_prim_num(k) - if(schwartz_kl(0,r)*schwartz_ij < thr) cycle - - coef3 = coef2 * ao_coef_cgtos_norm_ord_transp(r,k) - expo3 = ao_expo_cgtos_ord_transp(r,k) - phiK = ao_expo_phase_ord_transp(4,r,k) - - do s = 1, ao_prim_num(l) - if(schwartz_kl(s,r)*schwartz_ij < thr) cycle - - coef4 = coef3 * ao_coef_cgtos_norm_ord_transp(s,l) - expo4 = ao_expo_cgtos_ord_transp(s,l) - phiL = ao_expo_phase_ord_transp(4,s,l) - - C1 = zexp((0.d0, 1.d0) * (-phiI - phiJ - phiK - phiL)) - C2 = zexp((0.d0, 1.d0) * (-phiI - phiJ + phiK - phiL)) - C3 = zexp((0.d0, 1.d0) * ( phiI - phiJ - phiK - phiL)) - C4 = zexp((0.d0, 1.d0) * ( phiI - phiJ + phiK - phiL)) - C5 = zexp((0.d0, 1.d0) * (-phiI + phiJ - phiK - phiL)) - C6 = zexp((0.d0, 1.d0) * (-phiI + phiJ + phiK - phiL)) - C7 = zexp((0.d0, 1.d0) * ( phiI + phiJ - phiK - phiL)) - C8 = zexp((0.d0, 1.d0) * ( phiI + phiJ + phiK - phiL)) - - int1 = ERI_cgtos(expo1, expo2, expo3, expo4, & - I_power(1), J_power(1), K_power(1), L_power(1), & - I_power(2), J_power(2), K_power(2), L_power(2), & - I_power(3), J_power(3), K_power(3), L_power(3)) - - int2 = ERI_cgtos(expo1, expo2, conjg(expo3), expo4, & - I_power(1), J_power(1), K_power(1), L_power(1), & - I_power(2), J_power(2), K_power(2), L_power(2), & - I_power(3), J_power(3), K_power(3), L_power(3)) - - int3 = ERI_cgtos(conjg(expo1), expo2, expo3, expo4, & - I_power(1), J_power(1), K_power(1), L_power(1), & - I_power(2), J_power(2), K_power(2), L_power(2), & - I_power(3), J_power(3), K_power(3), L_power(3)) - - int4 = ERI_cgtos(conjg(expo1), expo2, conjg(expo3), expo4, & - I_power(1), J_power(1), K_power(1), L_power(1), & - I_power(2), J_power(2), K_power(2), L_power(2), & - I_power(3), J_power(3), K_power(3), L_power(3)) - - int5 = ERI_cgtos(expo1, conjg(expo2), expo3, expo4, & - I_power(1), J_power(1), K_power(1), L_power(1), & - I_power(2), J_power(2), K_power(2), L_power(2), & - I_power(3), J_power(3), K_power(3), L_power(3)) - - int6 = ERI_cgtos(expo1, conjg(expo2), conjg(expo3), expo4, & - I_power(1), J_power(1), K_power(1), L_power(1), & - I_power(2), J_power(2), K_power(2), L_power(2), & - I_power(3), J_power(3), K_power(3), L_power(3)) - - int7 = ERI_cgtos(conjg(expo1), conjg(expo2), expo3, expo4, & - I_power(1), J_power(1), K_power(1), L_power(1), & - I_power(2), J_power(2), K_power(2), L_power(2), & - I_power(3), J_power(3), K_power(3), L_power(3)) - - int8 = ERI_cgtos(conjg(expo1), conjg(expo2), conjg(expo3), expo4, & - I_power(1), J_power(1), K_power(1), L_power(1), & - I_power(2), J_power(2), K_power(2), L_power(2), & - I_power(3), J_power(3), K_power(3), L_power(3)) - - int_tot = C1 * int1 + C2 * int2 + C3 * int3 + C4 * int4 + C5 * int5 + C6 * int6 + C7 * int7 + C8 * int8 - - ao_2e_cgtos_schwartz_accel = ao_2e_cgtos_schwartz_accel + coef4 * 2.d0 * real(int_tot) - enddo ! s - enddo ! r - enddo ! q - enddo ! p - - endif - - deallocate(schwartz_kl) - -end - -! --- - -BEGIN_PROVIDER [double precision, ao_2e_cgtos_schwartz, (ao_num, ao_num)] - - BEGIN_DOC - ! Needed to compute Schwartz inequalities - END_DOC - - implicit none - integer :: i, k - double precision :: ao_two_e_integral_cgtos - - ao_2e_cgtos_schwartz(1,1) = ao_two_e_integral_cgtos(1, 1, 1, 1) - - !$OMP PARALLEL DO PRIVATE(i,k) & - !$OMP DEFAULT(NONE) & - !$OMP SHARED(ao_num, ao_2e_cgtos_schwartz) & - !$OMP SCHEDULE(dynamic) - do i = 1, ao_num - do k = 1, i - ao_2e_cgtos_schwartz(i,k) = dsqrt(ao_two_e_integral_cgtos(i, i, k, k)) - ao_2e_cgtos_schwartz(k,i) = ao_2e_cgtos_schwartz(i,k) - enddo - enddo - !$OMP END PARALLEL DO - -END_PROVIDER - -! --- - -complex*16 function general_primitive_integral_cgtos(dim, P_new, P_center, fact_p, p, p_inv, iorder_p, & - Q_new, Q_center, fact_q, q, q_inv, iorder_q) - - BEGIN_DOC - ! - ! Computes the integral where p,q,r,s are cos-cGTOS primitives - ! - END_DOC - - implicit none - include 'utils/constants.include.F' - - integer, intent(in) :: dim - integer, intent(in) :: iorder_p(3), iorder_q(3) - complex*16, intent(in) :: P_new(0:max_dim,3), P_center(3), fact_p, p, p_inv - complex*16, intent(in) :: Q_new(0:max_dim,3), Q_center(3), fact_q, q, q_inv - - integer :: i, j, nx, ny, nz, n_Ix, n_Iy, n_Iz, iorder, n_pt_tmp, n_pt_out - complex*16 :: pq, pq_inv, pq_inv_2, p01_1, p01_2, p10_1, p10_2, ppq, sq_ppq - complex*16 :: rho, dist, const - complex*16 :: accu, tmp_p, tmp_q - complex*16 :: dx(0:max_dim), Ix_pol(0:max_dim), dy(0:max_dim), Iy_pol(0:max_dim), dz(0:max_dim), Iz_pol(0:max_dim) - complex*16 :: d1(0:max_dim), d_poly(0:max_dim) - - complex*16, external :: crint_sum - - - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: dx, Ix_pol, dy, Iy_pol, dz, Iz_pol - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: d1, d_poly - - general_primitive_integral_cgtos = (0.d0, 0.d0) - - pq = (0.5d0, 0.d0) * p_inv * q_inv - pq_inv = (0.5d0, 0.d0) / (p + q) - pq_inv_2 = pq_inv + pq_inv - p10_1 = q * pq ! 1/(2p) - p01_1 = p * pq ! 1/(2q) - p10_2 = pq_inv_2 * p10_1 * q ! 0.5d0*q/(pq + p*p) - p01_2 = pq_inv_2 * p01_1 * p ! 0.5d0*p/(q*q + pq) - - sq_ppq = zsqrt(p + q) - - ! --- - - iorder = iorder_p(1) + iorder_q(1) + iorder_p(1) + iorder_q(1) - - do i = 0, iorder - Ix_pol(i) = (0.d0, 0.d0) - enddo - - n_Ix = 0 - do i = 0, iorder_p(1) - - tmp_p = P_new(i,1) - if(zabs(tmp_p) < thresh) cycle - - do j = 0, iorder_q(1) - - tmp_q = tmp_p * Q_new(j,1) - if(zabs(tmp_q) < thresh) cycle - - !DIR$ FORCEINLINE - call give_cpolynom_mult_center_x(P_center(1), Q_center(1), i, j, p, q, iorder, pq_inv, pq_inv_2, p10_1, p01_1, p10_2, p01_2, dx, nx) - !DIR$ FORCEINLINE - call add_cpoly_multiply(dx, nx, tmp_q, Ix_pol, n_Ix) - enddo - enddo - if(n_Ix == -1) then - return - endif - - ! --- - - iorder = iorder_p(2) + iorder_q(2) + iorder_p(2) + iorder_q(2) - - do i = 0, iorder - Iy_pol(i) = (0.d0, 0.d0) - enddo - - n_Iy = 0 - do i = 0, iorder_p(2) - - tmp_p = P_new(i,2) - if(zabs(tmp_p) < thresh) cycle - - do j = 0, iorder_q(2) - - tmp_q = tmp_p * Q_new(j,2) - if(zabs(tmp_q) < thresh) cycle - - !DIR$ FORCEINLINE - call give_cpolynom_mult_center_x(P_center(2), Q_center(2), i, j, p, q, iorder, pq_inv, pq_inv_2, p10_1, p01_1, p10_2, p01_2, dy, ny) - !DIR$ FORCEINLINE - call add_cpoly_multiply(dy, ny, tmp_q, Iy_pol, n_Iy) - enddo - enddo - - if(n_Iy == -1) then - return - endif - - ! --- - - iorder = iorder_p(3) + iorder_q(3) + iorder_p(3) + iorder_q(3) - - do i = 0, iorder - Iz_pol(i) = (0.d0, 0.d0) - enddo - - n_Iz = 0 - do i = 0, iorder_p(3) - - tmp_p = P_new(i,3) - if(zabs(tmp_p) < thresh) cycle - - do j = 0, iorder_q(3) - - tmp_q = tmp_p * Q_new(j,3) - if(zabs(tmp_q) < thresh) cycle - - !DIR$ FORCEINLINE - call give_cpolynom_mult_center_x(P_center(3), Q_center(3), i, j, p, q, iorder, pq_inv, pq_inv_2, p10_1, p01_1, p10_2, p01_2, dz, nz) - !DIR$ FORCEINLINE - call add_cpoly_multiply(dz, nz, tmp_q, Iz_pol, n_Iz) - enddo - enddo - - if(n_Iz == -1) then - return - endif - - ! --- - - rho = p * q * pq_inv_2 - dist = (P_center(1) - Q_center(1)) * (P_center(1) - Q_center(1)) & - + (P_center(2) - Q_center(2)) * (P_center(2) - Q_center(2)) & - + (P_center(3) - Q_center(3)) * (P_center(3) - Q_center(3)) - const = dist * rho - - n_pt_tmp = n_Ix + n_Iy - do i = 0, n_pt_tmp - d_poly(i) = (0.d0, 0.d0) - enddo - - !DIR$ FORCEINLINE - call multiply_cpoly(Ix_pol, n_Ix, Iy_pol, n_Iy, d_poly, n_pt_tmp) - if(n_pt_tmp == -1) then - return - endif - n_pt_out = n_pt_tmp + n_Iz - do i = 0, n_pt_out - d1(i) = (0.d0, 0.d0) - enddo - - !DIR$ FORCEINLINE - call multiply_cpoly(d_poly, n_pt_tmp, Iz_pol, n_Iz, d1, n_pt_out) - if(n_pt_out == -1) then - return - endif - - accu = crint_sum(n_pt_out, const, d1) - - general_primitive_integral_cgtos = fact_p * fact_q * accu * pi_5_2 * p_inv * q_inv / sq_ppq - -end - -! --- - -complex*16 function ERI_cgtos(alpha, beta, delta, gama, a_x, b_x, c_x, d_x, a_y, b_y, c_y, d_y, a_z, b_z, c_z, d_z) - - BEGIN_DOC - ! ATOMIC PRIMTIVE two-electron integral between the 4 primitives :: - ! primitive_1 = x1**(a_x) y1**(a_y) z1**(a_z) exp(-alpha * r1**2) - ! primitive_2 = x1**(b_x) y1**(b_y) z1**(b_z) exp(- beta * r1**2) - ! primitive_3 = x2**(c_x) y2**(c_y) z2**(c_z) exp(-delta * r2**2) - ! primitive_4 = x2**(d_x) y2**(d_y) z2**(d_z) exp(- gama * r2**2) - END_DOC - - implicit none - include 'utils/constants.include.F' - - integer, intent(in) :: a_x, b_x, c_x, d_x, a_y, b_y, c_y, d_y, a_z, b_z, c_z, d_z - complex*16, intent(in) :: delta, gama, alpha, beta - - integer :: a_x_2, b_x_2, c_x_2, d_x_2, a_y_2, b_y_2, c_y_2, d_y_2, a_z_2, b_z_2, c_z_2, d_z_2 - integer :: i, j, k, l, n_pt - integer :: nx, ny, nz - complex*16 :: p, q, ppq, sq_ppq, coeff, I_f - - ERI_cgtos = (0.d0, 0.d0) - - ASSERT (real(alpha) >= 0.d0) - ASSERT (real(beta ) >= 0.d0) - ASSERT (real(delta) >= 0.d0) - ASSERT (real(gama ) >= 0.d0) - - nx = a_x + b_x + c_x + d_x - if(iand(nx, 1) == 1) then - ERI_cgtos = (0.d0, 0.d0) - return - endif - - ny = a_y + b_y + c_y + d_y - if(iand(ny, 1) == 1) then - ERI_cgtos = (0.d0, 0.d0) - return - endif - - nz = a_z + b_z + c_z + d_z - if(iand(nz, 1) == 1) then - ERI_cgtos = (0.d0, 0.d0) - return - endif - - n_pt = shiftl(nx+ny+nz, 1) - - p = alpha + beta - q = delta + gama - - sq_ppq = zsqrt(p + q) - - coeff = pi_5_2 / (p * q * sq_ppq) - if(n_pt == 0) then - ERI_cgtos = coeff - return - endif - - call integrale_new_cgtos(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, p, q, n_pt) - - ERI_cgtos = I_f * coeff - -end - -! --- - -subroutine integrale_new_cgtos(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, p, q, n_pt) - - BEGIN_DOC - ! Calculates the integral of the polynomial : - ! - ! $I_{x_1}(a_x+b_x, c_x+d_x, p, q) \, I_{x_1}(a_y+b_y, c_y+d_y, p, q) \, I_{x_1}(a_z+b_z, c_z+d_z, p, q)$ - ! in $( 0 ; 1)$ - END_DOC - - implicit none - include 'utils/constants.include.F' - - integer, intent(in) :: n_pt - integer, intent(in) :: a_x, b_x, c_x, d_x, a_y, b_y, c_y, d_y, a_z, b_z, c_z, d_z - complex*16, intent(out) :: I_f - - integer :: i, j, ix, iy, iz, jx, jy, jz, sx, sy, sz - complex*16 :: p, q - complex*16 :: pq_inv, p10_1, p10_2, p01_1, p01_2, pq_inv_2 - complex*16 :: B00(n_pt_max_integrals), B10(n_pt_max_integrals), B01(n_pt_max_integrals) - complex*16 :: t1(n_pt_max_integrals), t2(n_pt_max_integrals) - - - ASSERT (n_pt > 1) - - j = shiftr(n_pt, 1) - - pq_inv = (0.5d0, 0.d0) / (p + q) - p10_1 = (0.5d0, 0.d0) / p - p01_1 = (0.5d0, 0.d0) / q - p10_2 = (0.5d0, 0.d0) * q /(p * q + p * p) - p01_2 = (0.5d0, 0.d0) * p /(q * q + q * p) - pq_inv_2 = pq_inv + pq_inv - - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: t1, t2, B10, B01, B00 - ix = a_x + b_x - jx = c_x + d_x - iy = a_y + b_y - jy = c_y + d_y - iz = a_z + b_z - jz = c_z + d_z - sx = ix + jx - sy = iy + jy - sz = iz + jz - - do i = 1, n_pt - B10(i) = p10_1 - gauleg_t2(i, j) * p10_2 - B01(i) = p01_1 - gauleg_t2(i, j) * p01_2 - B00(i) = gauleg_t2(i, j) * pq_inv - enddo - - if(sx > 0) then - call I_x1_new_cgtos(ix, jx, B10, B01, B00, t1, n_pt) - else - do i = 1, n_pt - t1(i) = (1.d0, 0.d0) - enddo - endif - - if(sy > 0) then - call I_x1_new_cgtos(iy, jy, B10, B01, B00, t2, n_pt) - do i = 1, n_pt - t1(i) = t1(i) * t2(i) - enddo - endif - - if(sz > 0) then - call I_x1_new_cgtos(iz, jz, B10, B01, B00, t2, n_pt) - do i = 1, n_pt - t1(i) = t1(i) * t2(i) - enddo - endif - - I_f = (0.d0, 0.d0) - do i = 1, n_pt - I_f += gauleg_w(i, j) * t1(i) - enddo - -end - -! --- - -recursive subroutine I_x1_new_cgtos(a, c, B_10, B_01, B_00, res, n_pt) - - BEGIN_DOC - ! recursive function involved in the two-electron integral - END_DOC - - implicit none - include 'utils/constants.include.F' - - integer, intent(in) :: a, c, n_pt - complex*16, intent(in) :: B_10(n_pt_max_integrals), B_01(n_pt_max_integrals), B_00(n_pt_max_integrals) - complex*16, intent(out) :: res(n_pt_max_integrals) - - integer :: i - complex*16 :: res2(n_pt_max_integrals) - - if(c < 0) then - - do i = 1, n_pt - res(i) = (0.d0, 0.d0) - enddo - - else if (a == 0) then - - call I_x2_new_cgtos(c, B_10, B_01, B_00, res, n_pt) - - else if (a == 1) then - - call I_x2_new_cgtos(c-1, B_10, B_01, B_00, res, n_pt) - do i = 1, n_pt - res(i) = dble(c) * B_00(i) * res(i) - enddo - - else - - call I_x1_new_cgtos(a-2, c , B_10, B_01, B_00, res , n_pt) - call I_x1_new_cgtos(a-1, c-1, B_10, B_01, B_00, res2, n_pt) - do i = 1, n_pt - res(i) = dble(a-1) * B_10(i) * res(i) + dble(c) * B_00(i) * res2(i) - enddo - - endif - -end - -! --- - -recursive subroutine I_x2_new_cgtos(c, B_10, B_01, B_00, res, n_pt) - - BEGIN_DOC - ! recursive function involved in the two-electron integral - END_DOC - - implicit none - include 'utils/constants.include.F' - - integer, intent(in) :: c, n_pt - complex*16, intent(in) :: B_10(n_pt_max_integrals), B_01(n_pt_max_integrals), B_00(n_pt_max_integrals) - complex*16, intent(out) :: res(n_pt_max_integrals) - - integer :: i - - if(c == 1) then - - do i = 1, n_pt - res(i) = (0.d0, 0.d0) - enddo - - elseif(c == 0) then - - do i = 1, n_pt - res(i) = (1.d0, 0.d0) - enddo - - else - - call I_x1_new_cgtos(0, c-2, B_10, B_01, B_00, res, n_pt) - do i = 1, n_pt - res(i) = dble(c-1) * B_01(i) * res(i) - enddo - - endif - -end - -! --- - -subroutine give_cpolynom_mult_center_x(P_center, Q_center, a_x, d_x, p, q, n_pt_in, & - pq_inv, pq_inv_2, p10_1, p01_1, p10_2, p01_2, d, n_pt_out) - - BEGIN_DOC - ! subroutine that returns the explicit polynom in term of the "t" - ! variable of the following polynoms : - ! - ! $I_{x_1}(a_x,d_x,p,q) \, I_{x_1}(a_y,d_y,p,q) \ I_{x_1}(a_z,d_z,p,q)$ - END_DOC - - implicit none - include 'utils/constants.include.F' - - integer, intent(in) :: n_pt_in, a_x, d_x - complex*16, intent(in) :: P_center, Q_center, p, q, pq_inv, p10_1, p01_1, p10_2, p01_2, pq_inv_2 - integer, intent(out) :: n_pt_out - complex*16, intent(out) :: d(0:max_dim) - - integer :: n_pt1, i - complex*16 :: B10(0:2), B01(0:2), B00(0:2), C00(0:2), D00(0:2) - - ASSERT (n_pt_in >= 0) - - B10(0) = p10_1 - B10(1) = (0.d0, 0.d0) - B10(2) = -p10_2 - - B01(0) = p01_1 - B01(1) = (0.d0, 0.d0) - B01(2) = -p01_2 - - B00(0) = (0.d0, 0.d0) - B00(1) = (0.d0, 0.d0) - B00(2) = pq_inv - - C00(0) = (0.d0, 0.d0) - C00(1) = (0.d0, 0.d0) - C00(2) = -q * (P_center - Q_center) * pq_inv_2 - - D00(0) = (0.d0, 0.d0) - D00(1) = (0.d0, 0.d0) - D00(2) = -p * (Q_center - P_center) * pq_inv_2 - - do i = 0, n_pt_in - d(i) = (0.d0, 0.d0) - enddo - - n_pt1 = n_pt_in - - !DIR$ FORCEINLINE - call I_x1_pol_mult_cgtos(a_x, d_x, B10, B01, B00, C00, D00, d, n_pt1, n_pt_in) - n_pt_out = n_pt1 - - if(n_pt1 < 0) then - n_pt_out = -1 - do i = 0, n_pt_in - d(i) = (0.d0, 0.d0) - enddo - return - endif - -end - -! --- - -subroutine I_x1_pol_mult_cgtos(a, c, B_10, B_01, B_00, C_00, D_00, d, nd, n_pt_in) - - BEGIN_DOC - ! Recursive function involved in the two-electron integral - END_DOC - - implicit none - include 'utils/constants.include.F' - - integer, intent(in) :: n_pt_in, a, c - complex*16, intent(in) :: B_10(0:2), B_01(0:2), B_00(0:2), C_00(0:2), D_00(0:2) - integer, intent(inout) :: nd - complex*16, intent(inout) :: d(0:max_dim) - - if( (c >= 0) .and. (nd >= 0) ) then - - if(a == 1) then - call I_x1_pol_mult_a1_cgtos(c, B_10, B_01, B_00, C_00, D_00, d, nd, n_pt_in) - else if(a == 2) then - call I_x1_pol_mult_a2_cgtos(c, B_10, B_01, B_00, C_00, D_00, d, nd, n_pt_in) - else if(a > 2) then - call I_x1_pol_mult_recurs_cgtos(a, c, B_10, B_01, B_00, C_00, D_00, d, nd, n_pt_in) - else ! a == 0 - - if(c == 0)then - nd = 0 - d(0) = (1.d0, 0.d0) - return - endif - - call I_x2_pol_mult_cgtos(c, B_10, B_01, B_00, C_00, D_00, d, nd, n_pt_in) - endif - - else - - nd = -1 - - endif - -end - -! --- - -recursive subroutine I_x1_pol_mult_recurs_cgtos(a, c, B_10, B_01, B_00, C_00, D_00, d, nd, n_pt_in) - - BEGIN_DOC - ! Recursive function involved in the two-electron integral - END_DOC - - implicit none - include 'utils/constants.include.F' - - integer, intent(in) :: n_pt_in, a, c - complex*16, intent(in) :: B_10(0:2), B_01(0:2), B_00(0:2), C_00(0:2), D_00(0:2) - integer, intent(inout) :: nd - complex*16, intent(inout) :: d(0:max_dim) - - integer :: nx, ix, iy, ny - complex*16 :: X(0:max_dim) - complex*16 :: Y(0:max_dim) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X,Y - - ASSERT (a > 2) - - !DIR$ LOOP COUNT(8) - do ix = 0, n_pt_in - X(ix) = (0.d0, 0.d0) - enddo - - nx = 0 - if(a == 3) then - call I_x1_pol_mult_a1_cgtos(c, B_10, B_01, B_00, C_00, D_00, X, nx, n_pt_in) - elseif(a == 4) then - call I_x1_pol_mult_a2_cgtos(c, B_10, B_01, B_00, C_00, D_00, X, nx, n_pt_in) - else - ASSERT (a >= 5) - call I_x1_pol_mult_recurs_cgtos(a-2, c, B_10, B_01, B_00, C_00, D_00, X, nx, n_pt_in) - endif - - !DIR$ LOOP COUNT(8) - do ix = 0, nx - X(ix) *= dble(a-1) - enddo - - !DIR$ FORCEINLINE - call multiply_cpoly(X, nx, B_10, 2, d, nd) - nx = nd - - !DIR$ LOOP COUNT(8) - do ix = 0, n_pt_in - X(ix) = (0.d0, 0.d0) - enddo - - if(c > 0) then - - if(a == 3) then - call I_x1_pol_mult_a2_cgtos(c-1, B_10, B_01, B_00, C_00, D_00, X, nx, n_pt_in) - else - ASSERT(a >= 4) - call I_x1_pol_mult_recurs_cgtos(a-1, c-1, B_10, B_01, B_00, C_00, D_00, X, nx, n_pt_in) - endif - - if(c > 1) then - !DIR$ LOOP COUNT(8) - do ix = 0, nx - X(ix) *= dble(c) - enddo - endif - !DIR$ FORCEINLINE - call multiply_cpoly(X, nx, B_00, 2, d, nd) - - endif - - ny = 0 - - !DIR$ LOOP COUNT(8) - do ix = 0, n_pt_in - Y(ix) = (0.d0, 0.d0) - enddo - - ASSERT (a > 2) - - if(a == 3) then - call I_x1_pol_mult_a2_cgtos(c, B_10, B_01, B_00, C_00, D_00, Y, ny, n_pt_in) - else - ASSERT(a >= 4) - call I_x1_pol_mult_recurs_cgtos(a-1, c, B_10, B_01, B_00, C_00, D_00, Y, ny, n_pt_in) - endif - - !DIR$ FORCEINLINE - call multiply_cpoly(Y, ny, C_00, 2, d, nd) - -end - -! --- - -recursive subroutine I_x1_pol_mult_a1_cgtos(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) - - BEGIN_DOC - ! Recursive function involved in the two-electron integral - END_DOC - - implicit none - include 'utils/constants.include.F' - - integer, intent(in) :: n_pt_in, c - complex*16, intent(in) :: B_10(0:2), B_01(0:2), B_00(0:2), C_00(0:2), D_00(0:2) - integer, intent(inout) :: nd - complex*16, intent(inout) :: d(0:max_dim) - - integer :: nx, ix, iy, ny - complex*16 :: X(0:max_dim) - complex*16 :: Y(0:max_dim) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X,Y - - if((c < 0) .or. (nd < 0)) then - nd = -1 - return - endif - - nx = nd - !DIR$ LOOP COUNT(8) - do ix = 0, n_pt_in - X(ix) = (0.d0, 0.d0) - enddo - call I_x2_pol_mult_cgtos(c-1, B_10, B_01, B_00, C_00, D_00, X, nx, n_pt_in) - - if(c > 1) then - !DIR$ LOOP COUNT(8) - do ix = 0, nx - X(ix) *= dble(c) - enddo - endif - - !DIR$ FORCEINLINE - call multiply_cpoly(X, nx, B_00, 2, d, nd) - - ny = 0 - - !DIR$ LOOP COUNT(8) - do ix = 0, n_pt_in - Y(ix) = (0.d0, 0.d0) - enddo - call I_x2_pol_mult_cgtos(c, B_10, B_01, B_00, C_00, D_00, Y, ny, n_pt_in) - - !DIR$ FORCEINLINE - call multiply_cpoly(Y, ny, C_00, 2, d, nd) - -end - -! --- - -recursive subroutine I_x1_pol_mult_a2_cgtos(c, B_10, B_01, B_00, C_00, D_00, d, nd, n_pt_in) - - BEGIN_DOC - ! Recursive function involved in the two-electron integral - END_DOC - - implicit none - include 'utils/constants.include.F' - - integer, intent(in) :: n_pt_in, c - complex*16, intent(in) :: B_10(0:2), B_01(0:2), B_00(0:2), C_00(0:2), D_00(0:2) - integer, intent(inout) :: nd - complex*16, intent(inout) :: d(0:max_dim) - - integer :: nx, ix, iy, ny - complex*16 :: X(0:max_dim) - complex*16 :: Y(0:max_dim) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X,Y - - !DIR$ LOOP COUNT(8) - do ix = 0, n_pt_in - X(ix) = (0.d0, 0.d0) - enddo - - nx = 0 - call I_x2_pol_mult_cgtos(c, B_10, B_01, B_00, C_00, D_00, X, nx, n_pt_in) - - !DIR$ FORCEINLINE - call multiply_cpoly(X, nx, B_10, 2, d, nd) - - nx = nd - !DIR$ LOOP COUNT(8) - do ix = 0, n_pt_in - X(ix) = (0.d0, 0.d0) - enddo - - !DIR$ FORCEINLINE - call I_x1_pol_mult_a1_cgtos(c-1, B_10, B_01, B_00, C_00, D_00, X, nx, n_pt_in) - - if (c>1) then - !DIR$ LOOP COUNT(8) - do ix = 0, nx - X(ix) *= dble(c) - enddo - endif - - !DIR$ FORCEINLINE - call multiply_cpoly(X, nx, B_00, 2, d, nd) - - ny = 0 - !DIR$ LOOP COUNT(8) - do ix = 0, n_pt_in - Y(ix) = 0.d0 - enddo - !DIR$ FORCEINLINE - call I_x1_pol_mult_a1_cgtos(c, B_10, B_01, B_00, C_00, D_00, Y, ny, n_pt_in) - - !DIR$ FORCEINLINE - call multiply_cpoly(Y, ny, C_00, 2, d, nd) - -end - -! --- - -recursive subroutine I_x2_pol_mult_cgtos(c, B_10, B_01, B_00, C_00, D_00, d, nd, dim) - - BEGIN_DOC - ! Recursive function involved in the two-electron integral - END_DOC - - implicit none - include 'utils/constants.include.F' - - integer, intent(in) :: dim, c - complex*16, intent(in) :: B_10(0:2), B_01(0:2), B_00(0:2), C_00(0:2), D_00(0:2) - integer, intent(inout) :: nd - complex*16, intent(inout) :: d(0:max_dim) - - integer :: i - integer :: nx, ix, ny - complex*16 :: X(0:max_dim), Y(0:max_dim) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X, Y - - select case (c) - - case (0) - nd = 0 - d(0) = (1.d0, 0.d0) - return - - case (:-1) - nd = -1 - return - - case (1) - nd = 2 - d(0) = D_00(0) - d(1) = D_00(1) - d(2) = D_00(2) - return - - case (2) - nd = 2 - d(0) = B_01(0) - d(1) = B_01(1) - d(2) = B_01(2) - - ny = 2 - Y(0) = D_00(0) - Y(1) = D_00(1) - Y(2) = D_00(2) - - !DIR$ FORCEINLINE - call multiply_cpoly(Y, ny, D_00, 2, d, nd) - return - - case default - - !DIR$ LOOP COUNT(6) - do ix = 0, c+c - X(ix) = (0.d0, 0.d0) - enddo - nx = 0 - call I_x2_pol_mult_cgtos(c-2, B_10, B_01, B_00, C_00, D_00, X, nx, dim) - - !DIR$ LOOP COUNT(6) - do ix = 0, nx - X(ix) *= dble(c-1) - enddo - - !DIR$ FORCEINLINE - call multiply_cpoly(X, nx, B_01, 2, d, nd) - - ny = 0 - !DIR$ LOOP COUNT(6) - do ix = 0, c+c - Y(ix) = 0.d0 - enddo - call I_x2_pol_mult_cgtos(c-1, B_10, B_01, B_00, C_00, D_00, Y, ny, dim) - - !DIR$ FORCEINLINE - call multiply_cpoly(Y, ny, D_00, 2, d, nd) - - end select - -end - diff --git a/src/ao_two_e_ints/two_e_integrals.irp.f b/src/ao_two_e_ints/two_e_integrals.irp.f deleted file mode 100644 index 1cb7617e..00000000 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ /dev/null @@ -1,1790 +0,0 @@ - -! --- - -logical function do_schwartz_accel(i,j,k,l) - implicit none - BEGIN_DOC - ! If true, use Schwatrz to accelerate direct integral calculation - END_DOC - integer, intent(in) :: i, j, k, l - if (do_ao_cholesky) then - do_schwartz_accel = .False. - else - do_schwartz_accel = (ao_prim_num(i) * ao_prim_num(j) * & - ao_prim_num(k) * ao_prim_num(l) > 1024 ) - endif - -end function - - -double precision function ao_two_e_integral(i, j, k, l) - - BEGIN_DOC - ! integral of the AO basis or (ij|kl) - ! i(r1) j(r1) 1/r12 k(r2) l(r2) - END_DOC - - implicit none - include 'utils/constants.include.F' - - integer, intent(in) :: i, j, k, l - - integer :: p, q, r, s - integer :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3) - integer :: iorder_p(3), iorder_q(3) - double precision :: I_center(3), J_center(3), K_center(3), L_center(3) - double precision :: integral - double precision :: P_new(0:max_dim,3),P_center(3),fact_p,pp - double precision :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq - - double precision, external :: ao_two_e_integral_erf - double precision, external :: ao_two_e_integral_cgtos - double precision, external :: ao_two_e_integral_schwartz_accel - double precision, external :: ao_two_e_integral_general - double precision, external :: general_primitive_integral - - logical, external :: do_schwartz_accel - double precision :: coef1, coef2, coef3, coef4 - - if(use_cgtos) then - - ao_two_e_integral = ao_two_e_integral_cgtos(i, j, k, l) - return - - else if (use_only_lr) then - - ao_two_e_integral = ao_two_e_integral_erf(i, j, k, l) - return - - else if (do_schwartz_accel(i,j,k,l)) then - - ao_two_e_integral = ao_two_e_integral_schwartz_accel(i,j,k,l) - - else - - num_i = ao_nucl(i) - num_j = ao_nucl(j) - num_k = ao_nucl(k) - num_l = ao_nucl(l) - ao_two_e_integral = 0.d0 - - if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k) then - - ao_two_e_integral = ao_two_e_integral_general(i,j,k,l,general_primitive_integral) - - else - - do p = 1, 3 - I_power(p) = ao_power(i,p) - J_power(p) = ao_power(j,p) - K_power(p) = ao_power(k,p) - L_power(p) = ao_power(l,p) - enddo - double precision :: ERI - - do p = 1, ao_prim_num(i) - coef1 = ao_coef_normalized_ordered_transp(p,i) - do q = 1, ao_prim_num(j) - coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) - do r = 1, ao_prim_num(k) - coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) - do s = 1, ao_prim_num(l) - coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) - integral = ERI( & - ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),& - I_power(1),J_power(1),K_power(1),L_power(1), & - I_power(2),J_power(2),K_power(2),L_power(2), & - I_power(3),J_power(3),K_power(3),L_power(3)) - ao_two_e_integral = ao_two_e_integral + coef4 * integral - enddo ! s - enddo ! r - enddo ! q - enddo ! p - - endif - - endif - -end - -! --- - -double precision function ao_two_e_integral_general(i, j, k, l, op) - - BEGIN_DOC - ! integral of the AO basis or (ij|kl) - ! i(r1) j(r1) 1/r12 k(r2) l(r2) - END_DOC - - implicit none - include 'utils/constants.include.F' - - integer, intent(in) :: i, j, k, l - double precision, external :: op ! Operator function - - integer :: p, q, r, s - integer :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3) - integer :: iorder_p(3), iorder_q(3) - double precision :: I_center(3), J_center(3), K_center(3), L_center(3) - double precision :: integral - double precision :: P_new(0:max_dim,3),P_center(3),fact_p,pp - double precision :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq - - dim1 = n_pt_max_integrals - - num_i = ao_nucl(i) - num_j = ao_nucl(j) - num_k = ao_nucl(k) - num_l = ao_nucl(l) - ao_two_e_integral_general = 0.d0 - - do p = 1, 3 - I_power(p) = ao_power(i,p) - J_power(p) = ao_power(j,p) - K_power(p) = ao_power(k,p) - L_power(p) = ao_power(l,p) - I_center(p) = nucl_coord(num_i,p) - J_center(p) = nucl_coord(num_j,p) - K_center(p) = nucl_coord(num_k,p) - L_center(p) = nucl_coord(num_l,p) - enddo - - double precision :: coef1, coef2, coef3, coef4 - double precision :: p_inv,q_inv - - do p = 1, ao_prim_num(i) - coef1 = ao_coef_normalized_ordered_transp(p,i) - do q = 1, ao_prim_num(j) - coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) - call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,& - ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j), & - I_power,J_power,I_center,J_center,dim1) - p_inv = 1.d0/pp - do r = 1, ao_prim_num(k) - coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) - do s = 1, ao_prim_num(l) - coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) - call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,& - ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), & - K_power,L_power,K_center,L_center,dim1) - q_inv = 1.d0/qq - integral = op(dim1, & - P_new,P_center,fact_p,pp,p_inv,iorder_p, & - Q_new,Q_center,fact_q,qq,q_inv,iorder_q) - ao_two_e_integral_general = ao_two_e_integral_general + coef4 * integral - enddo ! s - enddo ! r - enddo ! q - enddo ! p - -end - -double precision function ao_two_e_integral_schwartz_accel(i,j,k,l) - implicit none - BEGIN_DOC - ! integral of the AO basis or (ij|kl) - ! i(r1) j(r1) 1/r12 k(r2) l(r2) - END_DOC - integer,intent(in) :: i,j,k,l - integer :: p,q,r,s - double precision :: I_center(3),J_center(3),K_center(3),L_center(3) - integer :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3) - double precision :: integral - include 'utils/constants.include.F' - double precision :: P_new(0:max_dim,3),P_center(3),fact_p,pp - double precision :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq - integer :: iorder_p(3), iorder_q(3) - double precision, allocatable :: schwartz_kl(:,:) - double precision :: schwartz_ij - - dim1 = n_pt_max_integrals - - num_i = ao_nucl(i) - num_j = ao_nucl(j) - num_k = ao_nucl(k) - num_l = ao_nucl(l) - ao_two_e_integral_schwartz_accel = 0.d0 - double precision :: thr - thr = ao_integrals_threshold*ao_integrals_threshold - - allocate(schwartz_kl(0:ao_prim_num(l),0:ao_prim_num(k))) - - - if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k)then - do p = 1, 3 - I_power(p) = ao_power(i,p) - J_power(p) = ao_power(j,p) - K_power(p) = ao_power(k,p) - L_power(p) = ao_power(l,p) - I_center(p) = nucl_coord(num_i,p) - J_center(p) = nucl_coord(num_j,p) - K_center(p) = nucl_coord(num_k,p) - L_center(p) = nucl_coord(num_l,p) - enddo - - schwartz_kl(0,0) = 0.d0 - do r = 1, ao_prim_num(k) - coef1 = ao_coef_normalized_ordered_transp(r,k)*ao_coef_normalized_ordered_transp(r,k) - schwartz_kl(0,r) = 0.d0 - do s = 1, ao_prim_num(l) - coef2 = coef1 * ao_coef_normalized_ordered_transp(s,l) * ao_coef_normalized_ordered_transp(s,l) - call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,& - ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), & - K_power,L_power,K_center,L_center,dim1) - q_inv = 1.d0/qq - schwartz_kl(s,r) = general_primitive_integral(dim1, & - Q_new,Q_center,fact_q,qq,q_inv,iorder_q, & - Q_new,Q_center,fact_q,qq,q_inv,iorder_q) & - * coef2 - schwartz_kl(0,r) = max(schwartz_kl(0,r),schwartz_kl(s,r)) - enddo - schwartz_kl(0,0) = max(schwartz_kl(0,r),schwartz_kl(0,0)) - enddo - - do p = 1, ao_prim_num(i) - double precision :: coef1 - coef1 = ao_coef_normalized_ordered_transp(p,i) - do q = 1, ao_prim_num(j) - double precision :: coef2 - coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) - double precision :: p_inv,q_inv - call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,& - ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j), & - I_power,J_power,I_center,J_center,dim1) - p_inv = 1.d0/pp - schwartz_ij = general_primitive_integral(dim1, & - P_new,P_center,fact_p,pp,p_inv,iorder_p, & - P_new,P_center,fact_p,pp,p_inv,iorder_p) * & - coef2*coef2 - if (schwartz_kl(0,0)*schwartz_ij < thr) then - cycle - endif - do r = 1, ao_prim_num(k) - if (schwartz_kl(0,r)*schwartz_ij < thr) then - cycle - endif - double precision :: coef3 - coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) - do s = 1, ao_prim_num(l) - double precision :: coef4 - if (schwartz_kl(s,r)*schwartz_ij < thr) then - cycle - endif - coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) - double precision :: general_primitive_integral - call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,& - ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), & - K_power,L_power,K_center,L_center,dim1) - q_inv = 1.d0/qq - integral = general_primitive_integral(dim1, & - P_new,P_center,fact_p,pp,p_inv,iorder_p, & - Q_new,Q_center,fact_q,qq,q_inv,iorder_q) - ao_two_e_integral_schwartz_accel = ao_two_e_integral_schwartz_accel + coef4 * integral - enddo ! s - enddo ! r - enddo ! q - enddo ! p - - else - - do p = 1, 3 - I_power(p) = ao_power(i,p) - J_power(p) = ao_power(j,p) - K_power(p) = ao_power(k,p) - L_power(p) = ao_power(l,p) - enddo - double precision :: ERI - - schwartz_kl(0,0) = 0.d0 - do r = 1, ao_prim_num(k) - coef1 = ao_coef_normalized_ordered_transp(r,k)*ao_coef_normalized_ordered_transp(r,k) - schwartz_kl(0,r) = 0.d0 - do s = 1, ao_prim_num(l) - coef2 = coef1*ao_coef_normalized_ordered_transp(s,l)*ao_coef_normalized_ordered_transp(s,l) - schwartz_kl(s,r) = ERI( & - ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),& - K_power(1),L_power(1),K_power(1),L_power(1), & - K_power(2),L_power(2),K_power(2),L_power(2), & - K_power(3),L_power(3),K_power(3),L_power(3)) * & - coef2 - schwartz_kl(0,r) = max(schwartz_kl(0,r),schwartz_kl(s,r)) - enddo - schwartz_kl(0,0) = max(schwartz_kl(0,r),schwartz_kl(0,0)) - enddo - - do p = 1, ao_prim_num(i) - coef1 = ao_coef_normalized_ordered_transp(p,i) - do q = 1, ao_prim_num(j) - coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) - schwartz_ij = ERI( & - ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),& - I_power(1),J_power(1),I_power(1),J_power(1), & - I_power(2),J_power(2),I_power(2),J_power(2), & - I_power(3),J_power(3),I_power(3),J_power(3))*coef2*coef2 - if (schwartz_kl(0,0)*schwartz_ij < thr) then - cycle - endif - do r = 1, ao_prim_num(k) - if (schwartz_kl(0,r)*schwartz_ij < thr) then - cycle - endif - coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) - do s = 1, ao_prim_num(l) - if (schwartz_kl(s,r)*schwartz_ij < thr) then - cycle - endif - coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) - integral = ERI( & - ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),& - I_power(1),J_power(1),K_power(1),L_power(1), & - I_power(2),J_power(2),K_power(2),L_power(2), & - I_power(3),J_power(3),K_power(3),L_power(3)) - ao_two_e_integral_schwartz_accel = ao_two_e_integral_schwartz_accel + coef4 * integral - enddo ! s - enddo ! r - enddo ! q - enddo ! p - - endif - deallocate (schwartz_kl) - -end - - -integer function ao_l4(i,j,k,l) - implicit none - BEGIN_DOC -! Computes the product of l values of i,j,k,and l - END_DOC - integer, intent(in) :: i,j,k,l - ao_l4 = ao_l(i)*ao_l(j)*ao_l(k)*ao_l(l) -end - - - -subroutine compute_ao_two_e_integrals(j,k,l,sze,buffer_value) - implicit none - use map_module - - BEGIN_DOC - ! Compute AO 1/r12 integrals for all i and fixed j,k,l - END_DOC - - include 'utils/constants.include.F' - integer, intent(in) :: j,k,l,sze - real(integral_kind), intent(out) :: buffer_value(sze) - double precision :: ao_two_e_integral - - integer :: i - logical, external :: ao_one_e_integral_zero - logical, external :: ao_two_e_integral_zero - - - if (ao_one_e_integral_zero(j,l)) then - buffer_value = 0._integral_kind - return - endif - - do i = 1, ao_num - if (ao_two_e_integral_zero(i,j,k,l)) then - buffer_value(i) = 0._integral_kind - cycle - endif - !DIR$ FORCEINLINE - buffer_value(i) = ao_two_e_integral(i,k,j,l) - enddo - -end - -BEGIN_PROVIDER [ logical, ao_two_e_integrals_in_map ] - implicit none - use map_module - BEGIN_DOC - ! Map of Atomic integrals - ! i(r1) j(r2) 1/r12 k(r1) l(r2) - END_DOC - - integer :: i,j,k,l - double precision :: ao_two_e_integral,cpu_1,cpu_2, wall_1, wall_2 - double precision :: integral, wall_0 - include 'utils/constants.include.F' - - ! For integrals file - integer(key_kind),allocatable :: buffer_i(:) - integer :: size_buffer - real(integral_kind),allocatable :: buffer_value(:) - - integer :: n_integrals, rc - integer :: kk, m, j1, i1, lmax - character*(64) :: fmt - - double precision :: map_mb - PROVIDE read_ao_two_e_integrals io_ao_two_e_integrals ao_integrals_map - - 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) - - if (.True.) then - ! Avoid openMP - integral = ao_two_e_integral(1,1,1,1) - endif - - size_buffer = ao_num*ao_num - !$OMP PARALLEL DEFAULT(shared) private(j,l) & - !$OMP PRIVATE(buffer_i, buffer_value, n_integrals) - allocate(buffer_i(size_buffer), buffer_value(size_buffer)) - n_integrals = 0 - !$OMP DO COLLAPSE(1) SCHEDULE(dynamic) - do l=1,ao_num - do j=1,l - call compute_ao_integrals_jl(j,l,n_integrals,buffer_i,buffer_value) - call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value) - enddo - enddo - !$OMP END DO - deallocate(buffer_i, buffer_value) - !$OMP END PARALLEL - - 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)), ' )' - - 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') - endif - -END_PROVIDER - -! --- - -BEGIN_PROVIDER [ double precision, ao_two_e_integral_schwartz, (ao_num, ao_num) ] - - BEGIN_DOC - ! Needed to compute Schwartz inequalities - END_DOC - - implicit none - integer :: i, k - double precision :: ao_two_e_integral,cpu_1,cpu_2, wall_1, wall_2 - - ao_two_e_integral_schwartz(1,1) = ao_two_e_integral(1,1,1,1) - !$OMP PARALLEL DO PRIVATE(i,k) & - !$OMP DEFAULT(NONE) & - !$OMP SHARED (ao_num,ao_two_e_integral_schwartz) & - !$OMP SCHEDULE(dynamic) - do i=ao_num,1,-1 - do k=1,i - ao_two_e_integral_schwartz(i,k) = dsqrt(ao_two_e_integral(i,i,k,k)) - ao_two_e_integral_schwartz(k,i) = ao_two_e_integral_schwartz(i,k) - enddo - enddo - !$OMP END PARALLEL DO - -END_PROVIDER - -! --- - -double precision function general_primitive_integral(dim, & - P_new,P_center,fact_p,p,p_inv,iorder_p, & - Q_new,Q_center,fact_q,q,q_inv,iorder_q) - implicit none - BEGIN_DOC - ! Computes the integral where p,q,r,s are Gaussian primitives - END_DOC - integer,intent(in) :: dim - include 'utils/constants.include.F' - double precision, intent(in) :: P_new(0:max_dim,3),P_center(3),fact_p,p,p_inv - double precision, intent(in) :: Q_new(0:max_dim,3),Q_center(3),fact_q,q,q_inv - integer, intent(in) :: iorder_p(3) - integer, intent(in) :: iorder_q(3) - - double precision :: r_cut,gama_r_cut,rho,dist - double precision :: dx(0:max_dim),Ix_pol(0:max_dim),dy(0:max_dim),Iy_pol(0:max_dim),dz(0:max_dim),Iz_pol(0:max_dim) - integer :: n_Ix,n_Iy,n_Iz,nx,ny,nz - double precision :: bla - integer :: ix,iy,iz,jx,jy,jz,i - double precision :: a,b,c,d,e,f,accu,pq,const - double precision :: pq_inv, p10_1, p10_2, p01_1, p01_2,pq_inv_2 - integer :: n_pt_tmp,n_pt_out, iorder - double precision :: d1(0:max_dim),d_poly(0:max_dim),d1_screened(0:max_dim) - - general_primitive_integral = 0.d0 - - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: dx,Ix_pol,dy,Iy_pol,dz,Iz_pol - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: d1, d_poly - - ! Gaussian Product - ! ---------------- - - pq = p_inv*0.5d0*q_inv - pq_inv = 0.5d0/(p+q) - p10_1 = q*pq ! 1/(2p) - p01_1 = p*pq ! 1/(2q) - pq_inv_2 = pq_inv+pq_inv - p10_2 = pq_inv_2 * p10_1*q !0.5d0*q/(pq + p*p) - p01_2 = pq_inv_2 * p01_1*p !0.5d0*p/(q*q + pq) - - - accu = 0.d0 - iorder = iorder_p(1)+iorder_q(1)+iorder_p(1)+iorder_q(1) - do ix=0,iorder - Ix_pol(ix) = 0.d0 - enddo - n_Ix = 0 - do ix = 0, iorder_p(1) - if (abs(P_new(ix,1)) < thresh) cycle - a = P_new(ix,1) - do jx = 0, iorder_q(1) - d = a*Q_new(jx,1) - if (abs(d) < thresh) cycle - !DIR$ FORCEINLINE - call give_polynom_mult_center_x(P_center(1),Q_center(1),ix,jx,p,q,iorder,pq_inv,pq_inv_2,p10_1,p01_1,p10_2,p01_2,dx,nx) - !DIR$ FORCEINLINE - call add_poly_multiply(dx,nx,d,Ix_pol,n_Ix) - enddo - enddo - if (n_Ix == -1) then - return - endif - iorder = iorder_p(2)+iorder_q(2)+iorder_p(2)+iorder_q(2) - do ix=0, iorder - Iy_pol(ix) = 0.d0 - enddo - n_Iy = 0 - do iy = 0, iorder_p(2) - if (abs(P_new(iy,2)) > thresh) then - b = P_new(iy,2) - do jy = 0, iorder_q(2) - e = b*Q_new(jy,2) - if (abs(e) < thresh) cycle - !DIR$ FORCEINLINE - call give_polynom_mult_center_x(P_center(2),Q_center(2),iy,jy,p,q,iorder,pq_inv,pq_inv_2,p10_1,p01_1,p10_2,p01_2,dy,ny) - !DIR$ FORCEINLINE - call add_poly_multiply(dy,ny,e,Iy_pol,n_Iy) - enddo - endif - enddo - if (n_Iy == -1) then - return - endif - - iorder = iorder_p(3)+iorder_q(3)+iorder_p(3)+iorder_q(3) - do ix=0,iorder - Iz_pol(ix) = 0.d0 - enddo - n_Iz = 0 - do iz = 0, iorder_p(3) - if (abs(P_new(iz,3)) > thresh) then - c = P_new(iz,3) - do jz = 0, iorder_q(3) - f = c*Q_new(jz,3) - if (abs(f) < thresh) cycle - !DIR$ FORCEINLINE - call give_polynom_mult_center_x(P_center(3),Q_center(3),iz,jz,p,q,iorder,pq_inv,pq_inv_2,p10_1,p01_1,p10_2,p01_2,dz,nz) - !DIR$ FORCEINLINE - call add_poly_multiply(dz,nz,f,Iz_pol,n_Iz) - enddo - endif - enddo - if (n_Iz == -1) then - return - endif - - rho = p*q *pq_inv_2 - dist = (P_center(1) - Q_center(1))*(P_center(1) - Q_center(1)) + & - (P_center(2) - Q_center(2))*(P_center(2) - Q_center(2)) + & - (P_center(3) - Q_center(3))*(P_center(3) - Q_center(3)) - const = dist*rho - - n_pt_tmp = n_Ix+n_Iy - do i=0,n_pt_tmp - d_poly(i)=0.d0 - enddo - -! call multiply_poly(Ix_pol,n_Ix,Iy_pol,n_Iy,d_poly,n_pt_tmp) - integer :: ib, ic - if (ior(n_Ix,n_Iy) >= 0) then - do ib=0,n_Ix - do ic = 0,n_Iy - d_poly(ib+ic) = d_poly(ib+ic) + Iy_pol(ic) * Ix_pol(ib) - enddo - enddo - - do n_pt_tmp = n_Ix+n_Iy, 0, -1 - if (d_poly(n_pt_tmp) /= 0.d0) exit - enddo - endif - - if (n_pt_tmp == -1) then - return - endif - n_pt_out = n_pt_tmp+n_Iz - do i=0,n_pt_out - d1(i)=0.d0 - enddo - -! call multiply_poly(d_poly ,n_pt_tmp ,Iz_pol,n_Iz,d1,n_pt_out) - if (ior(n_pt_tmp,n_Iz) >= 0) then - ! Bottleneck here - do ib=0,n_pt_tmp - do ic = 0,n_Iz - d1(ib+ic) = d1(ib+ic) + Iz_pol(ic) * d_poly(ib) - enddo - enddo - - do n_pt_out = n_pt_tmp+n_Iz, 0, -1 - if (d1(n_pt_out) /= 0.d0) exit - enddo - endif - - - double precision :: rint_sum - accu = accu + rint_sum(n_pt_out,const,d1) - - general_primitive_integral = fact_p * fact_q * accu *pi_5_2*p_inv*q_inv/dsqrt(p+q) -end - - -double precision function ERI(alpha,beta,delta,gama,a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z) - implicit none - BEGIN_DOC - ! ATOMIC PRIMTIVE two-electron integral between the 4 primitives :: - ! primitive_1 = x1**(a_x) y1**(a_y) z1**(a_z) exp(-alpha * r1**2) - ! primitive_2 = x1**(b_x) y1**(b_y) z1**(b_z) exp(- beta * r1**2) - ! primitive_3 = x2**(c_x) y2**(c_y) z2**(c_z) exp(-delta * r2**2) - ! primitive_4 = x2**(d_x) y2**(d_y) z2**(d_z) exp(- gama * r2**2) - END_DOC - double precision, intent(in) :: delta,gama,alpha,beta - integer, intent(in) :: 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 :: a_x_2,b_x_2,c_x_2,d_x_2,a_y_2,b_y_2,c_y_2,d_y_2,a_z_2,b_z_2,c_z_2,d_z_2 - integer :: i,j,k,l,n_pt - integer :: n_pt_sup - double precision :: p,q,denom,coeff - double precision :: I_f - integer :: nx,ny,nz - include 'utils/constants.include.F' - nx = a_x+b_x+c_x+d_x - if(iand(nx,1) == 1) then - ERI = 0.d0 - return - endif - - ny = a_y+b_y+c_y+d_y - if(iand(ny,1) == 1) then - ERI = 0.d0 - return - endif - - nz = a_z+b_z+c_z+d_z - if(iand(nz,1) == 1) then - ERI = 0.d0 - return - endif - - ASSERT (alpha >= 0.d0) - ASSERT (beta >= 0.d0) - ASSERT (delta >= 0.d0) - ASSERT (gama >= 0.d0) - p = alpha + beta - q = delta + gama - ASSERT (p+q >= 0.d0) - n_pt = shiftl( nx+ny+nz,1 ) - - coeff = pi_5_2 / (p * q * dsqrt(p+q)) - if (n_pt == 0) then - ERI = coeff - return - endif - - call integrale_new(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,p,q,n_pt) - - ERI = I_f * coeff -end - - -subroutine integrale_new(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,p,q,n_pt) - BEGIN_DOC - ! Calculates the integral of the polynomial : - ! - ! $I_{x_1}(a_x+b_x,c_x+d_x,p,q) \, I_{x_1}(a_y+b_y,c_y+d_y,p,q) \, I_{x_1}(a_z+b_z,c_z+d_z,p,q)$ - ! in $( 0 ; 1)$ - END_DOC - - - implicit none - 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_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 - - j = shiftr(n_pt,1) - ASSERT (n_pt > 1) - pq_inv = 0.5d0/(p+q) - pq_inv_2 = pq_inv + pq_inv - p10_1 = 0.5d0/p - p01_1 = 0.5d0/q - p10_2 = 0.5d0 * q /(p * q + p * p) - p01_2 = 0.5d0 * p /(q * q + q * p) - double precision :: B00(n_pt_max_integrals) - double precision :: B10(n_pt_max_integrals), B01(n_pt_max_integrals) - double precision :: t1(n_pt_max_integrals), t2(n_pt_max_integrals) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: t1, t2, B10, B01, B00 - ix = a_x+b_x - jx = c_x+d_x - iy = a_y+b_y - jy = c_y+d_y - iz = a_z+b_z - jz = c_z+d_z - sx = ix+jx - sy = iy+jy - sz = iz+jz - - do i = 1,n_pt - B10(i) = p10_1 - gauleg_t2(i,j)* p10_2 - B01(i) = p01_1 - gauleg_t2(i,j)* p01_2 - B00(i) = gauleg_t2(i,j)*pq_inv - enddo - if (sx > 0) then - call I_x1_new(ix,jx,B10,B01,B00,t1,n_pt) - else - do i = 1,n_pt - t1(i) = 1.d0 - enddo - endif - if (sy > 0) then - call I_x1_new(iy,jy,B10,B01,B00,t2,n_pt) - do i = 1,n_pt - t1(i) = t1(i)*t2(i) - enddo - endif - if (sz > 0) then - call I_x1_new(iz,jz,B10,B01,B00,t2,n_pt) - do i = 1,n_pt - t1(i) = t1(i)*t2(i) - enddo - endif - I_f= 0.d0 - do i = 1,n_pt - I_f += gauleg_w(i,j)*t1(i) - enddo - - - -end - -recursive subroutine I_x1_new(a,c,B_10,B_01,B_00,res,n_pt) - BEGIN_DOC - ! recursive function involved in the two-electron integral - END_DOC - implicit none - include 'utils/constants.include.F' - integer, intent(in) :: a,c,n_pt - double precision, intent(in) :: B_10(n_pt_max_integrals),B_01(n_pt_max_integrals),B_00(n_pt_max_integrals) - double precision, intent(out) :: res(n_pt_max_integrals) - double precision :: res2(n_pt_max_integrals) - integer :: i - - if(c<0)then - do i=1,n_pt - res(i) = 0.d0 - enddo - else if (a==0) then - call I_x2_new(c,B_10,B_01,B_00,res,n_pt) - else if (a==1) then - call I_x2_new(c-1,B_10,B_01,B_00,res,n_pt) - do i=1,n_pt - res(i) = c * B_00(i) * res(i) - enddo - else - call I_x1_new(a-2,c,B_10,B_01,B_00,res,n_pt) - call I_x1_new(a-1,c-1,B_10,B_01,B_00,res2,n_pt) - do i=1,n_pt - res(i) = (a-1) * B_10(i) * res(i) & - + c * B_00(i) * res2(i) - enddo - endif -end - -recursive subroutine I_x2_new(c,B_10,B_01,B_00,res,n_pt) - implicit none - BEGIN_DOC - ! recursive function involved in the two-electron integral - END_DOC - include 'utils/constants.include.F' - integer, intent(in) :: c, n_pt - double precision, intent(in) :: B_10(n_pt_max_integrals),B_01(n_pt_max_integrals),B_00(n_pt_max_integrals) - double precision, intent(out) :: res(n_pt_max_integrals) - integer :: i - - if(c==1)then - do i=1,n_pt - res(i) = 0.d0 - enddo - elseif(c==0) then - do i=1,n_pt - res(i) = 1.d0 - enddo - else - call I_x1_new(0,c-2,B_10,B_01,B_00,res,n_pt) - do i=1,n_pt - res(i) = (c-1) * B_01(i) * res(i) - enddo - endif -end - - -integer function n_pt_sup(a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z) - implicit none - BEGIN_DOC - ! Returns the upper boundary of the degree of the polynomial involved in the - ! two-electron integral : - ! - ! $I_x(a_x,b_x,c_x,d_x) \, I_y(a_y,b_y,c_y,d_y) \, I_z(a_z,b_z,c_z,d_z)$ - END_DOC - 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 - n_pt_sup = shiftl( a_x+b_x+c_x+d_x + a_y+b_y+c_y+d_y + a_z+b_z+c_z+d_z,1 ) -end - - - - -subroutine give_polynom_mult_center_x(P_center,Q_center,a_x,d_x,p,q,n_pt_in,pq_inv,pq_inv_2,p10_1,p01_1,p10_2,p01_2,d,n_pt_out) - implicit none - BEGIN_DOC - ! subroutine that returns the explicit polynom in term of the "t" - ! variable of the following polynomw : - ! - ! $I_{x_1}(a_x,d_x,p,q) \, I_{x_1}(a_y,d_y,p,q) \ I_{x_1}(a_z,d_z,p,q)$ - END_DOC - integer, intent(in) :: n_pt_in - integer,intent(out) :: n_pt_out - integer, intent(in) :: a_x,d_x - double precision, intent(in) :: P_center, Q_center - double precision, intent(in) :: p,q,pq_inv,p10_1,p01_1,p10_2,p01_2,pq_inv_2 - include 'utils/constants.include.F' - double precision,intent(out) :: d(0:max_dim) - double precision :: accu - accu = 0.d0 - ASSERT (n_pt_in >= 0) - ! pq_inv = 0.5d0/(p+q) - ! pq_inv_2 = 1.d0/(p+q) - ! p10_1 = 0.5d0/p - ! p01_1 = 0.5d0/q - ! p10_2 = 0.5d0 * q /(p * q + p * p) - ! p01_2 = 0.5d0 * p /(q * q + q * p) - double precision :: B10(0:2), B01(0:2), B00(0:2),C00(0:2),D00(0:2) - B10(0) = p10_1 - B10(1) = 0.d0 - B10(2) = - p10_2 - ! B10 = p01_1 - t**2 * p10_2 - B01(0) = p01_1 - B01(1) = 0.d0 - B01(2) = - p01_2 - ! B01 = p01_1- t**2 * pq_inv - B00(0) = 0.d0 - B00(1) = 0.d0 - B00(2) = pq_inv - ! B00 = t**2 * pq_inv - do i = 0,n_pt_in - d(i) = 0.d0 - enddo - integer :: n_pt1,dim,i - n_pt1 = n_pt_in - ! C00 = -q/(p+q)*(Px-Qx) * t^2 - C00(0) = 0.d0 - C00(1) = 0.d0 - C00(2) = -q*(P_center-Q_center) * pq_inv_2 - ! D00 = -p/(p+q)*(Px-Qx) * t^2 - D00(0) = 0.d0 - D00(1) = 0.d0 - D00(2) = -p*(Q_center-P_center) * pq_inv_2 - !D00(2) = -p*(Q_center(1)-P_center(1)) /(p+q) - !DIR$ FORCEINLINE - call I_x1_pol_mult(a_x,d_x,B10,B01,B00,C00,D00,d,n_pt1,n_pt_in) - n_pt_out = n_pt1 - if(n_pt1<0)then - n_pt_out = -1 - do i = 0,n_pt_in - d(i) = 0.d0 - enddo - return - endif - -end - -subroutine I_x1_pol_mult(a,c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) - implicit none - BEGIN_DOC - ! Recursive function involved in the two-electron integral - END_DOC - integer , intent(in) :: n_pt_in - include 'utils/constants.include.F' - double precision,intent(inout) :: d(0:max_dim) - integer,intent(inout) :: nd - integer, intent(in) :: a,c - double precision, intent(in) :: B_10(0:2),B_01(0:2),B_00(0:2),C_00(0:2),D_00(0:2) - if( (c>=0).and.(nd>=0) )then - - if (a==1) then - call I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) - else if (a==2) then - call I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) - else if (a>2) then - call I_x1_pol_mult_recurs(a,c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) - else ! a == 0 - - if( c==0 )then - nd = 0 - d(0) = 1.d0 - return - endif - - call I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) - endif - else - nd = -1 - endif -end - -recursive subroutine I_x1_pol_mult_recurs(a,c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) - implicit none - BEGIN_DOC - ! Recursive function involved in the two-electron integral - END_DOC - integer , intent(in) :: n_pt_in - include 'utils/constants.include.F' - double precision,intent(inout) :: d(0:max_dim) - integer,intent(inout) :: nd - integer, intent(in) :: a,c - double precision, intent(in) :: B_10(0:2),B_01(0:2),B_00(0:2),C_00(0:2),D_00(0:2) - double precision :: X(0:max_dim) - double precision :: Y(0:max_dim) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X,Y - integer :: nx, ix,iy,ny,ib - - ASSERT (a>2) - !DIR$ LOOP COUNT(8) - do ix=0,n_pt_in - X(ix) = 0.d0 - enddo - nx = 0 - if (a==3) then - call I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,X,nx,n_pt_in) - else if (a==4) then - call I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,X,nx,n_pt_in) - else - ASSERT (a>=5) - call I_x1_pol_mult_recurs(a-2,c,B_10,B_01,B_00,C_00,D_00,X,nx,n_pt_in) - endif - - !DIR$ LOOP COUNT(8) - do ix=0,nx - X(ix) *= dble(a-1) - enddo - -! !DIR$ FORCEINLINE -! call multiply_poly_c2_inline_2e(X,nx,B_10,d,nd) - if (nx >= 0) then - select case (nx) - case (0) - d(0) = d(0) + B_10(0) * X(0) - d(1) = d(1) + B_10(1) * X(0) - d(2) = d(2) + B_10(2) * X(0) - - case (1) - d(0) = d(0) + B_10(0) * X(0) - d(1) = d(1) + B_10(0) * X(1) + B_10(1) * X(0) - d(2) = d(2) + B_10(1) * X(1) + B_10(2) * X(0) - d(3) = d(3) + B_10(2) * X(1) - - case (2) - d(0) = d(0) + B_10(0) * X(0) - d(1) = d(1) + B_10(0) * X(1) + B_10(1) * X(0) - d(2) = d(2) + B_10(0) * X(2) + B_10(1) * X(1) + B_10(2) * X(0) - d(3) = d(3) + B_10(1) * X(2) + B_10(2) * X(1) - d(4) = d(4) + B_10(2) * X(2) - - case default - d(0) = d(0) + B_10(0) * X(0) - d(1) = d(1) + B_10(0) * X(1) + B_10(1) * X(0) - do ib=2,nx - d(ib) = d(ib) + B_10(0) * X(ib) + B_10(1) * X(ib-1) + B_10(2) * X(ib-2) - enddo - d(nx+1) = d(nx+1) + B_10(1) * X(nx) + B_10(2) * X(nx-1) - d(nx+2) = d(nx+2) + B_10(2) * X(nx) - - end select - - do nd = nx+2,0,-1 - if (d(nd) /= 0.d0) exit - enddo - - endif - - nx = nd - !DIR$ LOOP COUNT(8) - do ix=0,n_pt_in - X(ix) = 0.d0 - enddo - - if (c>0) then - if (a==3) then - call I_x1_pol_mult_a2(c-1,B_10,B_01,B_00,C_00,D_00,X,nx,n_pt_in) - else - ASSERT(a >= 4) - call I_x1_pol_mult_recurs(a-1,c-1,B_10,B_01,B_00,C_00,D_00,X,nx,n_pt_in) - endif - if (c>1) then - !DIR$ LOOP COUNT(8) - do ix=0,nx - X(ix) *= c - enddo - endif - -! !DIR$ FORCEINLINE -! call multiply_poly_c2_inline_2e(X,nx,B_00,d,nd) - if(nx >= 0) then - - select case (nx) - case (0) - d(0) = d(0) + B_00(0) * X(0) - d(1) = d(1) + B_00(1) * X(0) - d(2) = d(2) + B_00(2) * X(0) - - case (1) - d(0) = d(0) + B_00(0) * X(0) - d(1) = d(1) + B_00(0) * X(1) + B_00(1) * X(0) - d(2) = d(2) + B_00(1) * X(1) + B_00(2) * X(0) - d(3) = d(3) + B_00(2) * X(1) - - case (2) - d(0) = d(0) + B_00(0) * X(0) - d(1) = d(1) + B_00(0) * X(1) + B_00(1) * X(0) - d(2) = d(2) + B_00(0) * X(2) + B_00(1) * X(1) + B_00(2) * X(0) - d(3) = d(3) + B_00(1) * X(2) + B_00(2) * X(1) - d(4) = d(4) + B_00(2) * X(2) - - case default - d(0) = d(0) + B_00(0) * X(0) - d(1) = d(1) + B_00(0) * X(1) + B_00(1) * X(0) - do ib=2,nx - d(ib) = d(ib) + B_00(0) * X(ib) + B_00(1) * X(ib-1) + B_00(2) * X(ib-2) - enddo - d(nx+1) = d(nx+1) + B_00(1) * X(nx) + B_00(2) * X(nx-1) - d(nx+2) = d(nx+2) + B_00(2) * X(nx) - - end select - - do nd = nx+2,0,-1 - if (d(nd) /= 0.d0) exit - enddo - - endif - - endif - - ny=0 - - !DIR$ LOOP COUNT(8) - do ix=0,n_pt_in - Y(ix) = 0.d0 - enddo - ASSERT(a > 2) - if (a==3) then - call I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,Y,ny,n_pt_in) - else - ASSERT(a >= 4) - call I_x1_pol_mult_recurs(a-1,c,B_10,B_01,B_00,C_00,D_00,Y,ny,n_pt_in) - endif - -! !DIR$ FORCEINLINE -! call multiply_poly_c2_inline_2e(Y,ny,C_00,d,nd) - if(ny >= 0) then - - select case (ny) - case (0) - d(0) = d(0) + C_00(0) * Y(0) - d(1) = d(1) + C_00(1) * Y(0) - d(2) = d(2) + C_00(2) * Y(0) - - case (1) - d(0) = d(0) + C_00(0) * Y(0) - d(1) = d(1) + C_00(0) * Y(1) + C_00(1) * Y(0) - d(2) = d(2) + C_00(1) * Y(1) + C_00(2) * Y(0) - d(3) = d(3) + C_00(2) * Y(1) - - case (2) - d(0) = d(0) + C_00(0) * Y(0) - d(1) = d(1) + C_00(0) * Y(1) + C_00(1) * Y(0) - d(2) = d(2) + C_00(0) * Y(2) + C_00(1) * Y(1) + C_00(2) * Y(0) - d(3) = d(3) + C_00(1) * Y(2) + C_00(2) * Y(1) - d(4) = d(4) + C_00(2) * Y(2) - - case default - d(0) = d(0) + C_00(0) * Y(0) - d(1) = d(1) + C_00(0) * Y(1) + C_00(1) * Y(0) - do ib=2,ny - d(ib) = d(ib) + C_00(0) * Y(ib) + C_00(1) * Y(ib-1) + C_00(2) * Y(ib-2) - enddo - d(ny+1) = d(ny+1) + C_00(1) * Y(ny) + C_00(2) * Y(ny-1) - d(ny+2) = d(ny+2) + C_00(2) * Y(ny) - - end select - - do nd = ny+2,0,-1 - if (d(nd) /= 0.d0) exit - enddo - - endif - -end - -recursive subroutine I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) - implicit none - BEGIN_DOC - ! Recursive function involved in the two-electron integral - END_DOC - integer , intent(in) :: n_pt_in - include 'utils/constants.include.F' - double precision,intent(inout) :: d(0:max_dim) - integer,intent(inout) :: nd - integer, intent(in) :: c - double precision, intent(in) :: B_10(0:2),B_01(0:2),B_00(0:2),C_00(0:2),D_00(0:2) - double precision :: X(0:max_dim) - double precision :: Y(0:max_dim) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X,Y - integer :: nx, ix,iy,ny,ib - - if( (c<0).or.(nd<0) )then - nd = -1 - return - endif - - nx = nd - !DIR$ LOOP COUNT(8) - do ix=0,n_pt_in - X(ix) = 0.d0 - enddo - call I_x2_pol_mult(c-1,B_10,B_01,B_00,C_00,D_00,X,nx,n_pt_in) - - if (c>1) then - !DIR$ LOOP COUNT(8) - do ix=0,nx - X(ix) *= dble(c) - enddo - endif - -! !DIR$ FORCEINLINE -! call multiply_poly_c2_inline_2e(X,nx,B_00,d,nd) - if(nx >= 0) then - - select case (nx) - case (0) - d(0) = d(0) + B_00(0) * X(0) - d(1) = d(1) + B_00(1) * X(0) - d(2) = d(2) + B_00(2) * X(0) - - case (1) - d(0) = d(0) + B_00(0) * X(0) - d(1) = d(1) + B_00(0) * X(1) + B_00(1) * X(0) - d(2) = d(2) + B_00(1) * X(1) + B_00(2) * X(0) - d(3) = d(3) + B_00(2) * X(1) - - case (2) - d(0) = d(0) + B_00(0) * X(0) - d(1) = d(1) + B_00(0) * X(1) + B_00(1) * X(0) - d(2) = d(2) + B_00(0) * X(2) + B_00(1) * X(1) + B_00(2) * X(0) - d(3) = d(3) + B_00(1) * X(2) + B_00(2) * X(1) - d(4) = d(4) + B_00(2) * X(2) - - case default - d(0) = d(0) + B_00(0) * X(0) - d(1) = d(1) + B_00(0) * X(1) + B_00(1) * X(0) - do ib=2,nx - d(ib) = d(ib) + B_00(0) * X(ib) + B_00(1) * X(ib-1) + B_00(2) * X(ib-2) - enddo - d(nx+1) = d(nx+1) + B_00(1) * X(nx) + B_00(2) * X(nx-1) - d(nx+2) = d(nx+2) + B_00(2) * X(nx) - - end select - - do nd = nx+2,0,-1 - if (d(nd) /= 0.d0) exit - enddo - - endif - - ny=0 - - !DIR$ LOOP COUNT(8) - do ix=0,n_pt_in - Y(ix) = 0.d0 - enddo - call I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,Y,ny,n_pt_in) - -! !DIR$ FORCEINLINE -! call multiply_poly_c2_inline_2e(Y,ny,C_00,d,nd) - if(ny >= 0) then - - select case (ny) - case (0) - d(0) = d(0) + C_00(0) * Y(0) - d(1) = d(1) + C_00(1) * Y(0) - d(2) = d(2) + C_00(2) * Y(0) - - case (1) - d(0) = d(0) + C_00(0) * Y(0) - d(1) = d(1) + C_00(0) * Y(1) + C_00(1) * Y(0) - d(2) = d(2) + C_00(1) * Y(1) + C_00(2) * Y(0) - d(3) = d(3) + C_00(2) * Y(1) - - case (2) - d(0) = d(0) + C_00(0) * Y(0) - d(1) = d(1) + C_00(0) * Y(1) + C_00(1) * Y(0) - d(2) = d(2) + C_00(0) * Y(2) + C_00(1) * Y(1) + C_00(2) * Y(0) - d(3) = d(3) + C_00(1) * Y(2) + C_00(2) * Y(1) - d(4) = d(4) + C_00(2) * Y(2) - - case default - d(0) = d(0) + C_00(0) * Y(0) - d(1) = d(1) + C_00(0) * Y(1) + C_00(1) * Y(0) - do ib=2,ny - d(ib) = d(ib) + C_00(0) * Y(ib) + C_00(1) * Y(ib-1) + C_00(2) * Y(ib-2) - enddo - d(ny+1) = d(ny+1) + C_00(1) * Y(ny) + C_00(2) * Y(ny-1) - d(ny+2) = d(ny+2) + C_00(2) * Y(ny) - - end select - - do nd = ny+2,0,-1 - if (d(nd) /= 0.d0) exit - enddo - - endif - -end - -recursive subroutine I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) - implicit none - BEGIN_DOC - ! Recursive function involved in the two-electron integral - END_DOC - integer , intent(in) :: n_pt_in - include 'utils/constants.include.F' - double precision,intent(inout) :: d(0:max_dim) - integer,intent(inout) :: nd - integer, intent(in) :: c - double precision, intent(in) :: B_10(0:2),B_01(0:2),B_00(0:2),C_00(0:2),D_00(0:2) - double precision :: X(0:max_dim) - double precision :: Y(0:max_dim) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X,Y - integer :: nx, ix,iy,ny,ib - - !DIR$ LOOP COUNT(8) - do ix=0,n_pt_in - X(ix) = 0.d0 - enddo - nx = 0 - call I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,X,nx,n_pt_in) - -! !DIR$ FORCEINLINE -! call multiply_poly_c2_inline_2e(X,nx,B_10,d,nd) - if(nx >= 0) then - - select case (nx) - case (0) - d(0) = d(0) + B_10(0) * X(0) - d(1) = d(1) + B_10(1) * X(0) - d(2) = d(2) + B_10(2) * X(0) - - case (1) - d(0) = d(0) + B_10(0) * X(0) - d(1) = d(1) + B_10(0) * X(1) + B_10(1) * X(0) - d(2) = d(2) + B_10(1) * X(1) + B_10(2) * X(0) - d(3) = d(3) + B_10(2) * X(1) - - case (2) - d(0) = d(0) + B_10(0) * X(0) - d(1) = d(1) + B_10(0) * X(1) + B_10(1) * X(0) - d(2) = d(2) + B_10(0) * X(2) + B_10(1) * X(1) + B_10(2) * X(0) - d(3) = d(3) + B_10(1) * X(2) + B_10(2) * X(1) - d(4) = d(4) + B_10(2) * X(2) - - case default - d(0) = d(0) + B_10(0) * X(0) - d(1) = d(1) + B_10(0) * X(1) + B_10(1) * X(0) - do ib=2,nx - d(ib) = d(ib) + B_10(0) * X(ib) + B_10(1) * X(ib-1) + B_10(2) * X(ib-2) - enddo - d(nx+1) = d(nx+1) + B_10(1) * X(nx) + B_10(2) * X(nx-1) - d(nx+2) = d(nx+2) + B_10(2) * X(nx) - - end select - - do nd = nx+2,0,-1 - if (d(nd) /= 0.d0) exit - enddo - - endif - - nx = nd - !DIR$ LOOP COUNT(8) - do ix=0,n_pt_in - X(ix) = 0.d0 - enddo - - !DIR$ FORCEINLINE - call I_x1_pol_mult_a1(c-1,B_10,B_01,B_00,C_00,D_00,X,nx,n_pt_in) - - if (c>1) then - !DIR$ LOOP COUNT(8) - do ix=0,nx - X(ix) *= dble(c) - enddo - endif - -! !DIR$ FORCEINLINE -! call multiply_poly_c2_inline_2e(X,nx,B_00,d,nd) - if(nx >= 0) then - - select case (nx) - case (0) - d(0) = d(0) + B_00(0) * X(0) - d(1) = d(1) + B_00(1) * X(0) - d(2) = d(2) + B_00(2) * X(0) - - case (1) - d(0) = d(0) + B_00(0) * X(0) - d(1) = d(1) + B_00(0) * X(1) + B_00(1) * X(0) - d(2) = d(2) + B_00(1) * X(1) + B_00(2) * X(0) - d(3) = d(3) + B_00(2) * X(1) - - case (2) - d(0) = d(0) + B_00(0) * X(0) - d(1) = d(1) + B_00(0) * X(1) + B_00(1) * X(0) - d(2) = d(2) + B_00(0) * X(2) + B_00(1) * X(1) + B_00(2) * X(0) - d(3) = d(3) + B_00(1) * X(2) + B_00(2) * X(1) - d(4) = d(4) + B_00(2) * X(2) - - case default - d(0) = d(0) + B_00(0) * X(0) - d(1) = d(1) + B_00(0) * X(1) + B_00(1) * X(0) - do ib=2,nx - d(ib) = d(ib) + B_00(0) * X(ib) + B_00(1) * X(ib-1) + B_00(2) * X(ib-2) - enddo - d(nx+1) = d(nx+1) + B_00(1) * X(nx) + B_00(2) * X(nx-1) - d(nx+2) = d(nx+2) + B_00(2) * X(nx) - - end select - - do nd = nx+2,0,-1 - if (d(nd) /= 0.d0) exit - enddo - - endif - - ny=0 - !DIR$ LOOP COUNT(8) - do ix=0,n_pt_in - Y(ix) = 0.d0 - enddo - !DIR$ FORCEINLINE - call I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,Y,ny,n_pt_in) - -! !DIR$ FORCEINLINE -! call multiply_poly_c2_inline_2e(Y,ny,C_00,d,nd) - if(ny >= 0) then - - select case (ny) - case (0) - d(0) = d(0) + C_00(0) * Y(0) - d(1) = d(1) + C_00(1) * Y(0) - d(2) = d(2) + C_00(2) * Y(0) - - case (1) - d(0) = d(0) + C_00(0) * Y(0) - d(1) = d(1) + C_00(0) * Y(1) + C_00(1) * Y(0) - d(2) = d(2) + C_00(1) * Y(1) + C_00(2) * Y(0) - d(3) = d(3) + C_00(2) * Y(1) - - case (2) - d(0) = d(0) + C_00(0) * Y(0) - d(1) = d(1) + C_00(0) * Y(1) + C_00(1) * Y(0) - d(2) = d(2) + C_00(0) * Y(2) + C_00(1) * Y(1) + C_00(2) * Y(0) - d(3) = d(3) + C_00(1) * Y(2) + C_00(2) * Y(1) - d(4) = d(4) + C_00(2) * Y(2) - - case default - d(0) = d(0) + C_00(0) * Y(0) - d(1) = d(1) + C_00(0) * Y(1) + C_00(1) * Y(0) - do ib=2,ny - d(ib) = d(ib) + C_00(0) * Y(ib) + C_00(1) * Y(ib-1) + C_00(2) * Y(ib-2) - enddo - d(ny+1) = d(ny+1) + C_00(1) * Y(ny) + C_00(2) * Y(ny-1) - d(ny+2) = d(ny+2) + C_00(2) * Y(ny) - - end select - - do nd = ny+2,0,-1 - if (d(nd) /= 0.d0) exit - enddo - - endif - -end - -recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim) - implicit none - BEGIN_DOC - ! Recursive function involved in the two-electron integral - END_DOC - integer , intent(in) :: dim - include 'utils/constants.include.F' - double precision :: d(0:max_dim) - integer,intent(inout) :: nd - integer, intent(in) :: c - double precision, intent(in) :: B_10(0:2),B_01(0:2),B_00(0:2),C_00(0:2),D_00(0:2) - integer :: nx, ix,ny - double precision :: X(0:max_dim),Y(0:max_dim) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X, Y - integer :: i, ib - - select case (c) - case (0) - nd = 0 - d(0) = 1.d0 - return - - case (:-1) - nd = -1 - return - - case (1) - nd = 2 - d(0) = D_00(0) - d(1) = D_00(1) - d(2) = D_00(2) - return - - case (2) - nd = 2 - d(0) = B_01(0) - d(1) = B_01(1) - d(2) = B_01(2) - - ny = 2 - Y(0) = D_00(0) - Y(1) = D_00(1) - Y(2) = D_00(2) - -! !DIR$ FORCEINLINE -! call multiply_poly_c2_inline_2e(Y,ny,D_00,d,nd) - if(ny >= 0) then - - select case (ny) - case (0) - d(0) = d(0) + D_00(0) * Y(0) - d(1) = d(1) + D_00(1) * Y(0) - d(2) = d(2) + D_00(2) * Y(0) - - case (1) - d(0) = d(0) + D_00(0) * Y(0) - d(1) = d(1) + D_00(0) * Y(1) + D_00(1) * Y(0) - d(2) = d(2) + D_00(1) * Y(1) + D_00(2) * Y(0) - d(3) = d(3) + D_00(2) * Y(1) - - case (2) - d(0) = d(0) + D_00(0) * Y(0) - d(1) = d(1) + D_00(0) * Y(1) + D_00(1) * Y(0) - d(2) = d(2) + D_00(0) * Y(2) + D_00(1) * Y(1) + D_00(2) * Y(0) - d(3) = d(3) + D_00(1) * Y(2) + D_00(2) * Y(1) - d(4) = d(4) + D_00(2) * Y(2) - - case default - d(0) = d(0) + D_00(0) * Y(0) - d(1) = d(1) + D_00(0) * Y(1) + D_00(1) * Y(0) - do ib=2,ny - d(ib) = d(ib) + D_00(0) * Y(ib) + D_00(1) * Y(ib-1) + D_00(2) * Y(ib-2) - enddo - d(ny+1) = d(ny+1) + D_00(1) * Y(ny) + D_00(2) * Y(ny-1) - d(ny+2) = d(ny+2) + D_00(2) * Y(ny) - - end select - - do nd = ny+2,0,-1 - if (d(nd) /= 0.d0) exit - enddo - - endif - - - return - - case default - - !DIR$ LOOP COUNT(6) - do ix=0,c+c - X(ix) = 0.d0 - enddo - nx = 0 - call I_x2_pol_mult(c-2,B_10,B_01,B_00,C_00,D_00,X,nx,dim) - - !DIR$ LOOP COUNT(6) - do ix=0,nx - X(ix) *= dble(c-1) - enddo - -! !DIR$ FORCEINLINE -! call multiply_poly_c2_inline_2e(X,nx,B_01,d,nd) - if(nx >= 0) then - - select case (nx) - case (0) - d(0) = d(0) + B_01(0) * X(0) - d(1) = d(1) + B_01(1) * X(0) - d(2) = d(2) + B_01(2) * X(0) - - case (1) - d(0) = d(0) + B_01(0) * X(0) - d(1) = d(1) + B_01(0) * X(1) + B_01(1) * X(0) - d(2) = d(2) + B_01(1) * X(1) + B_01(2) * X(0) - d(3) = d(3) + B_01(2) * X(1) - - case (2) - d(0) = d(0) + B_01(0) * X(0) - d(1) = d(1) + B_01(0) * X(1) + B_01(1) * X(0) - d(2) = d(2) + B_01(0) * X(2) + B_01(1) * X(1) + B_01(2) * X(0) - d(3) = d(3) + B_01(1) * X(2) + B_01(2) * X(1) - d(4) = d(4) + B_01(2) * X(2) - - case default - d(0) = d(0) + B_01(0) * X(0) - d(1) = d(1) + B_01(0) * X(1) + B_01(1) * X(0) - do ib=2,nx - d(ib) = d(ib) + B_01(0) * X(ib) + B_01(1) * X(ib-1) + B_01(2) * X(ib-2) - enddo - d(nx+1) = d(nx+1) + B_01(1) * X(nx) + B_01(2) * X(nx-1) - d(nx+2) = d(nx+2) + B_01(2) * X(nx) - - end select - - do nd = nx+2,0,-1 - if (d(nd) /= 0.d0) exit - enddo - - endif - - ny = 0 - !DIR$ LOOP COUNT(6) - do ix=0,c+c - Y(ix) = 0.d0 - enddo - call I_x2_pol_mult(c-1,B_10,B_01,B_00,C_00,D_00,Y,ny,dim) - -! !DIR$ FORCEINLINE -! call multiply_poly_c2_inline_2e(Y,ny,D_00,d,nd) - - if(ny >= 0) then - - select case (ny) - case (0) - d(0) = d(0) + D_00(0) * Y(0) - d(1) = d(1) + D_00(1) * Y(0) - d(2) = d(2) + D_00(2) * Y(0) - - case (1) - d(0) = d(0) + D_00(0) * Y(0) - d(1) = d(1) + D_00(0) * Y(1) + D_00(1) * Y(0) - d(2) = d(2) + D_00(1) * Y(1) + D_00(2) * Y(0) - d(3) = d(3) + D_00(2) * Y(1) - - case (2) - d(0) = d(0) + D_00(0) * Y(0) - d(1) = d(1) + D_00(0) * Y(1) + D_00(1) * Y(0) - d(2) = d(2) + D_00(0) * Y(2) + D_00(1) * Y(1) + D_00(2) * Y(0) - d(3) = d(3) + D_00(1) * Y(2) + D_00(2) * Y(1) - d(4) = d(4) + D_00(2) * Y(2) - - case default - d(0) = d(0) + D_00(0) * Y(0) - d(1) = d(1) + D_00(0) * Y(1) + D_00(1) * Y(0) - do ib=2,ny - d(ib) = d(ib) + D_00(0) * Y(ib) + D_00(1) * Y(ib-1) + D_00(2) * Y(ib-2) - enddo - d(ny+1) = d(ny+1) + D_00(1) * Y(ny) + D_00(2) * Y(ny-1) - d(ny+2) = d(ny+2) + D_00(2) * Y(ny) - - end select - - do nd = ny+2,0,-1 - if (d(nd) /= 0.d0) exit - enddo - - endif - - end select -end - - - - -subroutine compute_ao_integrals_jl(j,l,n_integrals,buffer_i,buffer_value) - implicit none - use map_module - BEGIN_DOC - ! Parallel client for AO integrals - END_DOC - - integer, intent(in) :: j,l - integer,intent(out) :: n_integrals - integer(key_kind),intent(out) :: buffer_i(ao_num*ao_num) - real(integral_kind),intent(out) :: buffer_value(ao_num*ao_num) - logical, external :: ao_two_e_integral_zero - - integer :: i,k - double precision, external :: ao_two_e_integral - double precision :: cpu_1,cpu_2, wall_1, wall_2 - double precision :: integral, wall_0 - double precision :: thr - integer :: kk, m, j1, i1 - - thr = ao_integrals_threshold - - n_integrals = 0 - - j1 = j+shiftr(l*l-l,1) - do k = 1, ao_num ! r1 - i1 = shiftr(k*k-k,1) - if (i1 > j1) then - exit - endif - do i = 1, k - i1 += 1 - if (i1 > j1) then - exit - endif - if (ao_two_e_integral_zero(i,j,k,l)) then - cycle - endif - !DIR$ FORCEINLINE - integral = ao_two_e_integral(i,k,j,l) ! i,k : r1 j,l : r2 - if (abs(integral) < thr) then - cycle - endif - n_integrals += 1 - !DIR$ FORCEINLINE - call two_e_integrals_index(i,j,k,l,buffer_i(n_integrals)) - buffer_value(n_integrals) = integral - enddo - enddo - -end - - -subroutine multiply_poly_local(b,nb,c,nc,d,nd) - implicit none - BEGIN_DOC - ! Multiply two polynomials - ! D(t) += B(t)*C(t) - END_DOC - - integer, intent(in) :: nb, nc - integer, intent(out) :: nd - double precision, intent(in) :: b(0:nb), c(0:nc) - double precision, intent(inout) :: d(0:nb+nc) - - integer :: ndtmp - integer :: ib, ic, id, k - if(ior(nc,nb) < 0) return !False if nc>=0 and nb>=0 - - do ib=0,nb - do ic = 0,nc - d(ib+ic) = d(ib+ic) + c(ic) * b(ib) - enddo - enddo - - do nd = nb+nc,0,-1 - if (d(nd) /= 0.d0) exit - enddo - -end - - -!DIR$ FORCEINLINE -subroutine multiply_poly_c2_inline_2e(b,nb,c,d,nd) - implicit none - BEGIN_DOC - ! Multiply two polynomials - ! D(t) += B(t)*C(t) - END_DOC - - integer, intent(in) :: nb - integer, intent(out) :: nd - double precision, intent(in) :: b(0:nb), c(0:2) - double precision, intent(inout) :: d(0:nb+2) - - integer :: ndtmp - integer :: ib, ic, id, k - if(nb < 0) return !False if nb>=0 - - select case (nb) - case (0) - d(0) = d(0) + c(0) * b(0) - d(1) = d(1) + c(1) * b(0) - d(2) = d(2) + c(2) * b(0) - - case (1) - d(0) = d(0) + c(0) * b(0) - d(1) = d(1) + c(0) * b(1) + c(1) * b(0) - d(2) = d(2) + c(1) * b(1) + c(2) * b(0) - d(3) = d(3) + c(2) * b(1) - - case (2) - d(0) = d(0) + c(0) * b(0) - d(1) = d(1) + c(0) * b(1) + c(1) * b(0) - d(2) = d(2) + c(0) * b(2) + c(1) * b(1) + c(2) * b(0) - d(3) = d(3) + c(1) * b(2) + c(2) * b(1) - d(4) = d(4) + c(2) * b(2) - - case default - d(0) = d(0) + c(0) * b(0) - d(1) = d(1) + c(0) * b(1) + c(1) * b(0) - do ib=2,nb - d(ib) = d(ib) + c(0) * b(ib) + c(1) * b(ib-1) + c(2) * b(ib-2) - enddo - d(nb+1) = d(nb+1) + c(1) * b(nb) + c(2) * b(nb-1) - d(nb+2) = d(nb+2) + c(2) * b(nb) - - end select - - do nd = nb+2,0,-1 - if (d(nd) /= 0.d0) exit - enddo - -end - diff --git a/src/ao_two_e_ints/two_e_integrals_erf.irp.f b/src/ao_two_e_ints/two_e_integrals_erf.irp.f deleted file mode 100644 index 88c74ac0..00000000 --- a/src/ao_two_e_ints/two_e_integrals_erf.irp.f +++ /dev/null @@ -1,658 +0,0 @@ -double precision function ao_two_e_integral_erf(i,j,k,l) - implicit none - BEGIN_DOC - ! integral of the AO basis or (ij|kl) - ! i(r1) j(r1) 1/r12 k(r2) l(r2) - END_DOC - - integer,intent(in) :: i,j,k,l - integer :: p,q,r,s - double precision :: I_center(3),J_center(3),K_center(3),L_center(3) - integer :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3) - double precision :: integral - include 'utils/constants.include.F' - double precision :: P_new(0:max_dim,3),P_center(3),fact_p,pp - double precision :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq - integer :: iorder_p(3), iorder_q(3) - double precision :: ao_two_e_integral_schwartz_accel_erf - - provide mu_erf - - if (ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then - ao_two_e_integral_erf = ao_two_e_integral_schwartz_accel_erf(i,j,k,l) - return - endif - - dim1 = n_pt_max_integrals - - num_i = ao_nucl(i) - num_j = ao_nucl(j) - num_k = ao_nucl(k) - num_l = ao_nucl(l) - ao_two_e_integral_erf = 0.d0 - - if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k)then - do p = 1, 3 - I_power(p) = ao_power(i,p) - J_power(p) = ao_power(j,p) - K_power(p) = ao_power(k,p) - L_power(p) = ao_power(l,p) - I_center(p) = nucl_coord(num_i,p) - J_center(p) = nucl_coord(num_j,p) - K_center(p) = nucl_coord(num_k,p) - L_center(p) = nucl_coord(num_l,p) - enddo - - double precision :: coef1, coef2, coef3, coef4 - double precision :: p_inv,q_inv - double precision :: general_primitive_integral_erf - - do p = 1, ao_prim_num(i) - coef1 = ao_coef_normalized_ordered_transp(p,i) - do q = 1, ao_prim_num(j) - coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) - call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,& - ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j), & - I_power,J_power,I_center,J_center,dim1) - p_inv = 1.d0/pp - do r = 1, ao_prim_num(k) - coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) - do s = 1, ao_prim_num(l) - coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) - call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,& - ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), & - K_power,L_power,K_center,L_center,dim1) - q_inv = 1.d0/qq - integral = general_primitive_integral_erf(dim1, & - P_new,P_center,fact_p,pp,p_inv,iorder_p, & - Q_new,Q_center,fact_q,qq,q_inv,iorder_q) - ao_two_e_integral_erf = ao_two_e_integral_erf + coef4 * integral - enddo ! s - enddo ! r - enddo ! q - enddo ! p - - else - - do p = 1, 3 - I_power(p) = ao_power(i,p) - J_power(p) = ao_power(j,p) - K_power(p) = ao_power(k,p) - L_power(p) = ao_power(l,p) - enddo - double precision :: ERI_erf - - do p = 1, ao_prim_num(i) - coef1 = ao_coef_normalized_ordered_transp(p,i) - do q = 1, ao_prim_num(j) - coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) - do r = 1, ao_prim_num(k) - coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) - do s = 1, ao_prim_num(l) - coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) - integral = ERI_erf( & - ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),& - I_power(1),J_power(1),K_power(1),L_power(1), & - I_power(2),J_power(2),K_power(2),L_power(2), & - I_power(3),J_power(3),K_power(3),L_power(3)) - ao_two_e_integral_erf = ao_two_e_integral_erf + coef4 * integral - enddo ! s - enddo ! r - enddo ! q - enddo ! p - - endif - -end - -double precision function ao_two_e_integral_schwartz_accel_erf(i,j,k,l) - implicit none - BEGIN_DOC - ! integral of the AO basis or (ij|kl) - ! i(r1) j(r1) 1/r12 k(r2) l(r2) - END_DOC - integer,intent(in) :: i,j,k,l - integer :: p,q,r,s - double precision :: I_center(3),J_center(3),K_center(3),L_center(3) - integer :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3) - double precision :: integral - include 'utils/constants.include.F' - double precision :: P_new(0:max_dim,3),P_center(3),fact_p,pp - double precision :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq - integer :: iorder_p(3), iorder_q(3) - double precision, allocatable :: schwartz_kl(:,:) - double precision :: schwartz_ij - - dim1 = n_pt_max_integrals - - num_i = ao_nucl(i) - num_j = ao_nucl(j) - num_k = ao_nucl(k) - num_l = ao_nucl(l) - ao_two_e_integral_schwartz_accel_erf = 0.d0 - double precision :: thr - thr = ao_integrals_threshold*ao_integrals_threshold - - allocate(schwartz_kl(0:ao_prim_num(l),0:ao_prim_num(k))) - - double precision :: coef3 - double precision :: coef2 - double precision :: p_inv,q_inv - double precision :: coef1 - double precision :: coef4 - - if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k)then - do p = 1, 3 - I_power(p) = ao_power(i,p) - J_power(p) = ao_power(j,p) - K_power(p) = ao_power(k,p) - L_power(p) = ao_power(l,p) - I_center(p) = nucl_coord(num_i,p) - J_center(p) = nucl_coord(num_j,p) - K_center(p) = nucl_coord(num_k,p) - L_center(p) = nucl_coord(num_l,p) - enddo - - schwartz_kl(0,0) = 0.d0 - do r = 1, ao_prim_num(k) - coef1 = ao_coef_normalized_ordered_transp(r,k)*ao_coef_normalized_ordered_transp(r,k) - schwartz_kl(0,r) = 0.d0 - do s = 1, ao_prim_num(l) - coef2 = coef1 * ao_coef_normalized_ordered_transp(s,l) * ao_coef_normalized_ordered_transp(s,l) - call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,& - ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), & - K_power,L_power,K_center,L_center,dim1) - q_inv = 1.d0/qq - schwartz_kl(s,r) = general_primitive_integral_erf(dim1, & - Q_new,Q_center,fact_q,qq,q_inv,iorder_q, & - Q_new,Q_center,fact_q,qq,q_inv,iorder_q) & - * coef2 - schwartz_kl(0,r) = max(schwartz_kl(0,r),schwartz_kl(s,r)) - enddo - schwartz_kl(0,0) = max(schwartz_kl(0,r),schwartz_kl(0,0)) - enddo - - do p = 1, ao_prim_num(i) - coef1 = ao_coef_normalized_ordered_transp(p,i) - do q = 1, ao_prim_num(j) - coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) - call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,& - ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j), & - I_power,J_power,I_center,J_center,dim1) - p_inv = 1.d0/pp - schwartz_ij = general_primitive_integral_erf(dim1, & - P_new,P_center,fact_p,pp,p_inv,iorder_p, & - P_new,P_center,fact_p,pp,p_inv,iorder_p) * & - coef2*coef2 - if (schwartz_kl(0,0)*schwartz_ij < thr) then - cycle - endif - do r = 1, ao_prim_num(k) - if (schwartz_kl(0,r)*schwartz_ij < thr) then - cycle - endif - coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) - do s = 1, ao_prim_num(l) - if (schwartz_kl(s,r)*schwartz_ij < thr) then - cycle - endif - coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) - double precision :: general_primitive_integral_erf - call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,& - ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), & - K_power,L_power,K_center,L_center,dim1) - q_inv = 1.d0/qq - integral = general_primitive_integral_erf(dim1, & - P_new,P_center,fact_p,pp,p_inv,iorder_p, & - Q_new,Q_center,fact_q,qq,q_inv,iorder_q) - ao_two_e_integral_schwartz_accel_erf = ao_two_e_integral_schwartz_accel_erf + coef4 * integral - enddo ! s - enddo ! r - enddo ! q - enddo ! p - - else - - do p = 1, 3 - I_power(p) = ao_power(i,p) - J_power(p) = ao_power(j,p) - K_power(p) = ao_power(k,p) - L_power(p) = ao_power(l,p) - enddo - double precision :: ERI_erf - - schwartz_kl(0,0) = 0.d0 - do r = 1, ao_prim_num(k) - coef1 = ao_coef_normalized_ordered_transp(r,k)*ao_coef_normalized_ordered_transp(r,k) - schwartz_kl(0,r) = 0.d0 - do s = 1, ao_prim_num(l) - coef2 = coef1*ao_coef_normalized_ordered_transp(s,l)*ao_coef_normalized_ordered_transp(s,l) - schwartz_kl(s,r) = ERI_erf( & - ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),& - K_power(1),L_power(1),K_power(1),L_power(1), & - K_power(2),L_power(2),K_power(2),L_power(2), & - K_power(3),L_power(3),K_power(3),L_power(3)) * & - coef2 - schwartz_kl(0,r) = max(schwartz_kl(0,r),schwartz_kl(s,r)) - enddo - schwartz_kl(0,0) = max(schwartz_kl(0,r),schwartz_kl(0,0)) - enddo - - do p = 1, ao_prim_num(i) - coef1 = ao_coef_normalized_ordered_transp(p,i) - do q = 1, ao_prim_num(j) - coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) - schwartz_ij = ERI_erf( & - ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),& - I_power(1),J_power(1),I_power(1),J_power(1), & - I_power(2),J_power(2),I_power(2),J_power(2), & - I_power(3),J_power(3),I_power(3),J_power(3))*coef2*coef2 - if (schwartz_kl(0,0)*schwartz_ij < thr) then - cycle - endif - do r = 1, ao_prim_num(k) - if (schwartz_kl(0,r)*schwartz_ij < thr) then - cycle - endif - coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) - do s = 1, ao_prim_num(l) - if (schwartz_kl(s,r)*schwartz_ij < thr) then - cycle - endif - coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) - integral = ERI_erf( & - ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),& - I_power(1),J_power(1),K_power(1),L_power(1), & - I_power(2),J_power(2),K_power(2),L_power(2), & - I_power(3),J_power(3),K_power(3),L_power(3)) - ao_two_e_integral_schwartz_accel_erf = ao_two_e_integral_schwartz_accel_erf + coef4 * integral - enddo ! s - enddo ! r - enddo ! q - enddo ! p - - endif - deallocate (schwartz_kl) - -end - - -subroutine compute_ao_two_e_integrals_erf(j,k,l,sze,buffer_value) - implicit none - use map_module - - BEGIN_DOC - ! Compute AO 1/r12 integrals for all i and fixed j,k,l - END_DOC - - include 'utils/constants.include.F' - integer, intent(in) :: j,k,l,sze - real(integral_kind), intent(out) :: buffer_value(sze) - double precision :: ao_two_e_integral_erf - - integer :: i - logical, external :: ao_one_e_integral_zero - logical, external :: ao_two_e_integral_zero - - if (ao_one_e_integral_zero(j,l)) then - buffer_value = 0._integral_kind - return - endif - if (ao_two_e_integral_erf_schwartz(j,l) < thresh ) then - buffer_value = 0._integral_kind - return - endif - - do i = 1, ao_num - if (ao_two_e_integral_zero(i,j,k,l)) then - buffer_value(i) = 0._integral_kind - cycle - endif - if (ao_two_e_integral_erf_schwartz(i,k)*ao_two_e_integral_erf_schwartz(j,l) < thresh ) then - buffer_value(i) = 0._integral_kind - cycle - endif - !DIR$ FORCEINLINE - buffer_value(i) = ao_two_e_integral_erf(i,k,j,l) - enddo - -end - -double precision function general_primitive_integral_erf(dim, & - P_new,P_center,fact_p,p,p_inv,iorder_p, & - Q_new,Q_center,fact_q,q,q_inv,iorder_q) - implicit none - BEGIN_DOC - ! Computes the integral where p,q,r,s are Gaussian primitives - END_DOC - integer,intent(in) :: dim - include 'utils/constants.include.F' - double precision, intent(in) :: P_new(0:max_dim,3),P_center(3),fact_p,p,p_inv - double precision, intent(in) :: Q_new(0:max_dim,3),Q_center(3),fact_q,q,q_inv - integer, intent(in) :: iorder_p(3) - integer, intent(in) :: iorder_q(3) - - double precision :: r_cut,gama_r_cut,rho,dist - double precision :: dx(0:max_dim),Ix_pol(0:max_dim),dy(0:max_dim),Iy_pol(0:max_dim),dz(0:max_dim),Iz_pol(0:max_dim) - integer :: n_Ix,n_Iy,n_Iz,nx,ny,nz - double precision :: bla - integer :: ix,iy,iz,jx,jy,jz,i - double precision :: a,b,c,d,e,f,accu,pq,const - double precision :: pq_inv, p10_1, p10_2, p01_1, p01_2,pq_inv_2 - integer :: n_pt_tmp,n_pt_out, iorder - double precision :: d1(0:max_dim),d_poly(0:max_dim),rint,d1_screened(0:max_dim) - - general_primitive_integral_erf = 0.d0 - - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: dx,Ix_pol,dy,Iy_pol,dz,Iz_pol - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: d1, d_poly - - ! Gaussian Product - ! ---------------- - double precision :: p_plus_q - p_plus_q = (p+q) * ((p*q)/(p+q) + mu_erf*mu_erf)/(mu_erf*mu_erf) - pq = p_inv*0.5d0*q_inv - - pq_inv = 0.5d0/p_plus_q - p10_1 = q*pq ! 1/(2p) - p01_1 = p*pq ! 1/(2q) - pq_inv_2 = pq_inv+pq_inv - p10_2 = pq_inv_2 * p10_1*q !0.5d0*q/(pq + p*p) - p01_2 = pq_inv_2 * p01_1*p !0.5d0*p/(q*q + pq) - - - accu = 0.d0 - iorder = iorder_p(1)+iorder_q(1)+iorder_p(1)+iorder_q(1) - !DIR$ VECTOR ALIGNED - do ix=0,iorder - Ix_pol(ix) = 0.d0 - enddo - n_Ix = 0 - do ix = 0, iorder_p(1) - if (abs(P_new(ix,1)) < thresh) cycle - a = P_new(ix,1) - do jx = 0, iorder_q(1) - d = a*Q_new(jx,1) - if (abs(d) < thresh) cycle - !DEC$ FORCEINLINE - call give_polynom_mult_center_x(P_center(1),Q_center(1),ix,jx,p,q,iorder,pq_inv,pq_inv_2,p10_1,p01_1,p10_2,p01_2,dx,nx) - !DEC$ FORCEINLINE - call add_poly_multiply(dx,nx,d,Ix_pol,n_Ix) - enddo - enddo - if (n_Ix == -1) then - return - endif - iorder = iorder_p(2)+iorder_q(2)+iorder_p(2)+iorder_q(2) - !DIR$ VECTOR ALIGNED - do ix=0, iorder - Iy_pol(ix) = 0.d0 - enddo - n_Iy = 0 - do iy = 0, iorder_p(2) - if (abs(P_new(iy,2)) > thresh) then - b = P_new(iy,2) - do jy = 0, iorder_q(2) - e = b*Q_new(jy,2) - if (abs(e) < thresh) cycle - !DEC$ FORCEINLINE - call give_polynom_mult_center_x(P_center(2),Q_center(2),iy,jy,p,q,iorder,pq_inv,pq_inv_2,p10_1,p01_1,p10_2,p01_2,dy,ny) - !DEC$ FORCEINLINE - call add_poly_multiply(dy,ny,e,Iy_pol,n_Iy) - enddo - endif - enddo - if (n_Iy == -1) then - return - endif - - iorder = iorder_p(3)+iorder_q(3)+iorder_p(3)+iorder_q(3) - do ix=0,iorder - Iz_pol(ix) = 0.d0 - enddo - n_Iz = 0 - do iz = 0, iorder_p(3) - if (abs(P_new(iz,3)) > thresh) then - c = P_new(iz,3) - do jz = 0, iorder_q(3) - f = c*Q_new(jz,3) - if (abs(f) < thresh) cycle - !DEC$ FORCEINLINE - call give_polynom_mult_center_x(P_center(3),Q_center(3),iz,jz,p,q,iorder,pq_inv,pq_inv_2,p10_1,p01_1,p10_2,p01_2,dz,nz) - !DEC$ FORCEINLINE - call add_poly_multiply(dz,nz,f,Iz_pol,n_Iz) - enddo - endif - enddo - if (n_Iz == -1) then - return - endif - - rho = p*q *pq_inv_2 ! le rho qui va bien - dist = (P_center(1) - Q_center(1))*(P_center(1) - Q_center(1)) + & - (P_center(2) - Q_center(2))*(P_center(2) - Q_center(2)) + & - (P_center(3) - Q_center(3))*(P_center(3) - Q_center(3)) - const = dist*rho - - n_pt_tmp = n_Ix+n_Iy - do i=0,n_pt_tmp - d_poly(i)=0.d0 - enddo - - !DEC$ FORCEINLINE - call multiply_poly(Ix_pol,n_Ix,Iy_pol,n_Iy,d_poly,n_pt_tmp) - if (n_pt_tmp == -1) then - return - endif - n_pt_out = n_pt_tmp+n_Iz - do i=0,n_pt_out - d1(i)=0.d0 - enddo - - !DEC$ FORCEINLINE - call multiply_poly(d_poly ,n_pt_tmp ,Iz_pol,n_Iz,d1,n_pt_out) - double precision :: rint_sum - accu = accu + rint_sum(n_pt_out,const,d1) - - ! change p+q in dsqrt - general_primitive_integral_erf = fact_p * fact_q * accu *pi_5_2*p_inv*q_inv/dsqrt(p_plus_q) -end - - -double precision function ERI_erf(alpha,beta,delta,gama,a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z) - implicit none - BEGIN_DOC - ! Atomic primtive two-electron integral between the 4 primitives : - ! - ! * primitive 1 : $x_1^{a_x} y_1^{a_y} z_1^{a_z} \exp(-\alpha * r1^2)$ - ! * primitive 2 : $x_1^{b_x} y_1^{b_y} z_1^{b_z} \exp(- \beta * r1^2)$ - ! * primitive 3 : $x_2^{c_x} y_2^{c_y} z_2^{c_z} \exp(-\delta * r2^2)$ - ! * primitive 4 : $x_2^{d_x} y_2^{d_y} z_2^{d_z} \exp(-\gamma * r2^2)$ - ! - END_DOC - double precision, intent(in) :: delta,gama,alpha,beta - integer, intent(in) :: 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 :: a_x_2,b_x_2,c_x_2,d_x_2,a_y_2,b_y_2,c_y_2,d_y_2,a_z_2,b_z_2,c_z_2,d_z_2 - integer :: i,j,k,l,n_pt - integer :: n_pt_sup - double precision :: p,q,denom,coeff - double precision :: I_f - integer :: nx,ny,nz - include 'utils/constants.include.F' - nx = a_x+b_x+c_x+d_x - if(iand(nx,1) == 1) then - ERI_erf = 0.d0 - return - endif - - ny = a_y+b_y+c_y+d_y - if(iand(ny,1) == 1) then - ERI_erf = 0.d0 - return - endif - - nz = a_z+b_z+c_z+d_z - if(iand(nz,1) == 1) then - ERI_erf = 0.d0 - return - endif - - ASSERT (alpha >= 0.d0) - ASSERT (beta >= 0.d0) - ASSERT (delta >= 0.d0) - ASSERT (gama >= 0.d0) - p = alpha + beta - q = delta + gama - double precision :: p_plus_q - p_plus_q = (p+q) * ((p*q)/(p+q) + mu_erf*mu_erf)/(mu_erf*mu_erf) - ASSERT (p+q >= 0.d0) - n_pt = ishft( nx+ny+nz,1 ) - - coeff = pi_5_2 / (p * q * dsqrt(p_plus_q)) - if (n_pt == 0) then - ERI_erf = coeff - return - endif - - call 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,p,q,n_pt) - - ERI_erf = I_f * coeff -end - - - -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,p,q,n_pt) - BEGIN_DOC - ! Calculate the integral of the polynomial : - ! - ! $I_x1(a_x+b_x, c_x+d_x,p,q) \, I_x1(a_y+b_y, c_y+d_y,p,q) \, I_x1(a_z+b_z, c_z+d_z,p,q)$ - ! - ! between $( 0 ; 1)$ - END_DOC - - - implicit none - 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_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 - - j = ishft(n_pt,-1) - ASSERT (n_pt > 1) - double precision :: p_plus_q - p_plus_q = (p+q) * ((p*q)/(p+q) + mu_erf*mu_erf)/(mu_erf*mu_erf) - - pq_inv = 0.5d0/(p_plus_q) - pq_inv_2 = pq_inv + pq_inv - p10_1 = 0.5d0/p - p01_1 = 0.5d0/q - p10_2 = 0.5d0 * q /(p * p_plus_q) - p01_2 = 0.5d0 * p /(q * p_plus_q) - double precision :: B00(n_pt_max_integrals) - double precision :: B10(n_pt_max_integrals), B01(n_pt_max_integrals) - double precision :: t1(n_pt_max_integrals), t2(n_pt_max_integrals) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: t1, t2, B10, B01, B00 - ix = a_x+b_x - jx = c_x+d_x - iy = a_y+b_y - jy = c_y+d_y - iz = a_z+b_z - jz = c_z+d_z - sx = ix+jx - sy = iy+jy - sz = iz+jz - - !DIR$ VECTOR ALIGNED - do i = 1,n_pt - B10(i) = p10_1 - gauleg_t2(i,j)* p10_2 - B01(i) = p01_1 - gauleg_t2(i,j)* p01_2 - B00(i) = gauleg_t2(i,j)*pq_inv - enddo - if (sx > 0) then - call I_x1_new(ix,jx,B10,B01,B00,t1,n_pt) - else - !DIR$ VECTOR ALIGNED - do i = 1,n_pt - t1(i) = 1.d0 - enddo - endif - if (sy > 0) then - call I_x1_new(iy,jy,B10,B01,B00,t2,n_pt) - !DIR$ VECTOR ALIGNED - do i = 1,n_pt - t1(i) = t1(i)*t2(i) - enddo - endif - if (sz > 0) then - call I_x1_new(iz,jz,B10,B01,B00,t2,n_pt) - !DIR$ VECTOR ALIGNED - do i = 1,n_pt - t1(i) = t1(i)*t2(i) - enddo - endif - I_f= 0.d0 - !DIR$ VECTOR ALIGNED - do i = 1,n_pt - I_f += gauleg_w(i,j)*t1(i) - enddo - - - -end - - -subroutine compute_ao_integrals_erf_jl(j,l,n_integrals,buffer_i,buffer_value) - implicit none - use map_module - BEGIN_DOC - ! Parallel client for AO integrals - END_DOC - - integer, intent(in) :: j,l - integer,intent(out) :: n_integrals - integer(key_kind),intent(out) :: buffer_i(ao_num*ao_num) - real(integral_kind),intent(out) :: buffer_value(ao_num*ao_num) - - integer :: i,k - double precision :: ao_two_e_integral_erf,cpu_1,cpu_2, wall_1, wall_2 - double precision :: integral, wall_0 - double precision :: thr - integer :: kk, m, j1, i1 - logical, external :: ao_two_e_integral_zero - - thr = ao_integrals_threshold - - n_integrals = 0 - - j1 = j+ishft(l*l-l,-1) - do k = 1, ao_num ! r1 - i1 = ishft(k*k-k,-1) - if (i1 > j1) then - exit - endif - do i = 1, k - i1 += 1 - if (i1 > j1) then - exit - endif - if (ao_two_e_integral_zero(i,j,k,l)) then - cycle - endif - if (ao_two_e_integral_erf_schwartz(i,k)*ao_two_e_integral_erf_schwartz(j,l) < thr ) then - cycle - endif - !DIR$ FORCEINLINE - integral = ao_two_e_integral_erf(i,k,j,l) ! i,k : r1 j,l : r2 - if (abs(integral) < thr) then - cycle - endif - n_integrals += 1 - !DIR$ FORCEINLINE - call two_e_integrals_index(i,j,k,l,buffer_i(n_integrals)) - buffer_value(n_integrals) = integral - enddo - enddo - -end