LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2025-04-25 09:44:43 +02:00
Code Releases Activity

added two_e_ints_keywords and substitute ao_cart in ao_two_e_ints

Browse Source
This commit is contained in:
eginer 2025-04-18 16:21:02 +02:00
parent 25e7b33bae
commit afe0baa60a
17 changed files with 5872 additions and 292 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

275
plugins/local/tc_keywords/EZFIO.cfg
View File

@ -1,275 +0,0 @@
[read_rl_eigv]
type: logical
doc: If |true|, read the right/left eigenvectors from ezfio
interface: ezfio,provider,ocaml
default: False
[comp_left_eigv]
type: logical
doc: If |true|, computes also the left-eigenvector
interface: ezfio,provider,ocaml
default: False
[three_body_h_tc]
type: logical
doc: If |true|, three-body terms are included
interface: ezfio,provider,ocaml
default: False
[three_e_3_idx_term]
type: logical
doc: If |true|, the diagonal 3-idx terms of the 3-e interaction are taken
interface: ezfio,provider,ocaml
default: True
[three_e_4_idx_term]
type: logical
doc: If |true|, the off-diagonal 4-idx terms of the 3-e interaction are taken
interface: ezfio,provider,ocaml
default: True
[three_e_5_idx_term]
type: logical
doc: If |true|, the off-diagonal 5-idx terms of the 3-e interaction are taken
interface: ezfio,provider,ocaml
default: True
[pure_three_body_h_tc]
type: logical
doc: If |true|, pure triple excitation three-body terms are included
interface: ezfio,provider,ocaml
default: False
[double_normal_ord]
type: logical
doc: If |true|, contracted double excitation three-body terms are included
interface: ezfio,provider,ocaml
default: False
[noL_standard]
type: logical
doc: If |true|, standard normal-ordering for L (to be used with three_body_h_tc |false|)
interface: ezfio,provider,ocaml
default: True
[core_tc_op]
type: logical
doc: If |true|, takes the usual Hamiltonian for core orbitals (assumed to be doubly occupied)
interface: ezfio,provider,ocaml
default: False
[full_tc_h_solver]
type: logical
doc: If |true|, you diagonalize the full TC H matrix
interface: ezfio,provider,ocaml
default: False
[thresh_it_dav]
type: Threshold
doc: Thresholds on the energy for iterative Davidson used in TC
interface: ezfio,provider,ocaml
default: 1.e-5
[thrsh_cycle_tc]
type: Threshold
doc: Thresholds to cycle the integrals with the envelop
interface: ezfio,provider,ocaml
default: 1.e-10
[max_it_dav]
type: integer
doc: nb max of iteration in Davidson used in TC
interface: ezfio,provider,ocaml
default: 1000
[thresh_psi_r]
type: Threshold
doc: Thresholds on the coefficients of the right-eigenvector. Used for PT2 computation.
interface: ezfio,provider,ocaml
default: 0.000005
[thresh_psi_r_norm]
type: logical
doc: If |true|, you prune the WF to compute the PT1 coef based on the norm. If False, the pruning is done through the amplitude on the right-coefficient.
interface: ezfio,provider,ocaml
default: False
[state_following_tc]
type: logical
doc: If |true|, the states are re-ordered to match the input states
default: False
interface: ezfio,provider,ocaml
[symmetric_fock_tc]
type: logical
doc: If |true|, using F+F^t as Fock TC
interface: ezfio,provider,ocaml
default: False
[selection_tc]
type: integer
doc: if +1: only positive is selected, -1: only negative is selected, :0 both positive and negative
interface: ezfio,provider,ocaml
default: 0
[mu_r_ct]
type: double precision
doc: a parameter used to define mu(r)
interface: ezfio, provider, ocaml
default: 1.5
[beta_rho_power]
type: double precision
doc: a parameter used to define mu(r)
interface: ezfio, provider, ocaml
default: 0.33333
[zeta_erf_mu_of_r]
type: double precision
doc: a parameter used to define mu(r)
interface: ezfio, provider, ocaml
default: 1.
[thr_degen_tc]
type: Threshold
doc: Threshold to determine if two orbitals are degenerate in TCSCF in order to avoid random quasi orthogonality between the right- and left-eigenvector for the same eigenvalue
interface: ezfio,provider,ocaml
default: 1.e-6
[maxovl_tc]
type: logical
doc: If |true|, maximize the overlap between orthogonalized left- and right eigenvectors
interface: ezfio,provider,ocaml
default: False
[test_cycle_tc]
type: logical
doc: If |true|, the integrals of the three-body jastrow are computed with cycles
interface: ezfio,provider,ocaml
default: False
[thresh_biorthog_diag]
type: Threshold
doc: Threshold to determine if diagonal elements of the bi-orthogonal condition L.T x R are close enouph to 1
interface: ezfio,provider,ocaml
default: 1.e-6
[thresh_lr_angle]
type: double precision
doc: Maximum value of the angle between the couple of left and right orbital for the rotations
interface: ezfio,provider,ocaml
default: 20.0
[thresh_biorthog_nondiag]
type: Threshold
doc: Threshold to determine if non-diagonal elements of L.T x R are close enouph to 0
interface: ezfio,provider,ocaml
default: 1.e-6
[var_tc]
type: logical
doc: If |true|, use VAR-TC
interface: ezfio,provider,ocaml
default: False
[io_tc_integ]
type: Disk_access
doc: Read/Write integrals int2_grad1_u12_ao, tc_grad_square_ao and tc_grad_and_lapl_ao from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml
default: None
[io_tc_norm_ord]
type: Disk_access
doc: Read/Write normal_two_body_bi_orth from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml
default: None
[only_spin_tc_right]
type: logical
doc: If |true|, only the right part of WF is used to compute spin dens
interface: ezfio,provider,ocaml
default: False
[save_sorted_tc_wf]
type: logical
doc: If |true|, save the bi-ortho wave functions in a sorted way
interface: ezfio,provider,ocaml
default: True
[use_ipp]
type: logical
doc: If |true|, use Manu IPP
interface: ezfio,provider,ocaml
default: True
[tc_grid1_a]
type: integer
doc: size of angular grid over r1: [ 6 | 14 | 26 | 38 | 50 | 74 | 86 | 110 | 146 | 170 | 194 | 230 | 266 | 302 | 350 | 434 | 590 | 770 | 974 | 1202 | 1454 | 1730 | 2030 | 2354 | 2702 | 3074 | 3470 | 3890 | 4334 | 4802 | 5294 | 5810 ]
interface: ezfio,provider,ocaml
default: 50
[tc_grid1_r]
type: integer
doc: size of radial grid over r1
interface: ezfio,provider,ocaml
default: 30
[tc_grid2_a]
type: integer
doc: size of angular grid over r2: [ 6 | 14 | 26 | 38 | 50 | 74 | 86 | 110 | 146 | 170 | 194 | 230 | 266 | 302 | 350 | 434 | 590 | 770 | 974 | 1202 | 1454 | 1730 | 2030 | 2354 | 2702 | 3074 | 3470 | 3890 | 4334 | 4802 | 5294 | 5810 ]
interface: ezfio,provider,ocaml
default: 266
[tc_grid2_r]
type: integer
doc: size of radial grid over r2
interface: ezfio,provider,ocaml
default: 70
[tc_integ_type]
type: character*(32)
doc: approach used to evaluate TC integrals [ analytic | numeric | semi-analytic ]
interface: ezfio,ocaml,provider
default: numeric
[minimize_lr_angles]
type: logical
doc: If |true|, you minimize the angle between the left and right vectors associated to degenerate orbitals
interface: ezfio,provider,ocaml
default: False
[thresh_de_tc_angles]
type: Threshold
doc: Thresholds on delta E for changing angles between orbitals
interface: ezfio,provider,ocaml
default: 1.e-6
[ao_to_mo_tc_n3]
type: logical
doc: If |true|, memory scale of TC ao -> mo: O(N3)
interface: ezfio,provider,ocaml
default: False
[tc_save_mem_loops]
type: logical
doc: If |true|, use loops to save memory TC
interface: ezfio,provider,ocaml
default: False
[tc_save_mem]
type: logical
doc: If |true|, more calc but less mem
interface: ezfio,provider,ocaml
default: False
[im_thresh_tc]
type: Threshold
doc: Thresholds on the Imag part of TC energy
interface: ezfio,provider,ocaml
default: 1.e-7
[transpose_two_e_int]
type: logical
doc: If |true|, you duplicate the two-electron TC integrals with the transpose matrix. Acceleates the PT2.
interface: ezfio,provider,ocaml
default: False

2
plugins/local/tc_keywords/NEED
View File

@ -1,2 +0,0 @@
ezfio_files
nuclei

6
src/ao_cart_two_e_ints/NEED Normal file
View File

@ -0,0 +1,6 @@
hamiltonian
ao_one_e_ints
pseudo
bitmask
ao_basis
two_e_ints_keywords

17
src/ao_cart_two_e_ints/README.rst Normal file
View File

@ -0,0 +1,17 @@
==================
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) = <ik|jl> = <12|12>

499
src/ao_cart_two_e_ints/cholesky.irp.f Normal file
View File

@ -0,0 +1,499 @@
double precision function get_ao_cart_integ_chol(i,j,k,l)
implicit none
BEGIN_DOC
! CHOLESKY representation of the integral of the AO basis <ik|jl> 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)
end
BEGIN_PROVIDER [ double precision, cholesky_ao_cart_transp, (cholesky_ao_cart_num, ao_cart_num, ao_cart_num) ]
implicit none
BEGIN_DOC
! Transposed of the Cholesky vectors in AO basis set
END_DOC
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)
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) ]
use mmap_module
implicit none
BEGIN_DOC
! Cholesky vectors in AO basis: (ik|a):
! <ij|kl> = (ik|jl) = sum_a (ik|a).(a|jl)
!
! Last dimension of cholesky_ao is cholesky_ao_cart_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
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_cart_two_e_integral
logical, external :: ao_cart_two_e_integral_zero
double precision, external :: ao_cart_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_cart_num*ao_cart_num*1_8+1
double precision :: wall0,wall1
type(mmap_type) :: map
PROVIDE nproc ao_cart_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_cart_cholesky) then
print *, 'Reading Cholesky AO vectors from disk...'
iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao', '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
else
call set_multiple_levels_omp(.False.)
if (do_direct_integrals) then
if (ao_cart_two_e_integral(1,1,1,1) < huge(1.d0)) then
! Trigger providers inside ao_cart_two_e_integral
continue
endif
else
PROVIDE ao_cart_two_e_integrals_in_map
endif
tau = ao_cart_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_cart_num
do i=1,ao_cart_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_cart_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_cart_two_e_integral(addr1(i8), addr1(i8), &
addr2(i8), addr2(i8), ao_cart_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_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)),&
addr1(Dset(m)), addr2(Dset(m)))
endif
enddo
!$OMP END PARALLEL DO
else
PROVIDE ao_cart_integrals_map
!$OMP PARALLEL DO PRIVATE(k) SCHEDULE(dynamic,21)
do k=1,np
Delta_col(k) = 0.d0
if (.not.ao_cart_two_e_integral_zero( addr1(Lset(k)), addr1(Dset(m)),&
addr2(Lset(k)), addr2(Dset(m)) ) ) then
Delta_col(k) = &
get_ao_cart_two_e_integral( addr1(Lset(k)), addr1(Dset(m)),&
addr2(Lset(k)), addr2(Dset(m)), ao_cart_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_cart_num,ao_cart_num,rank), stat=ierr)
if (ierr /= 0) then
call print_memory_usage()
print *, irp_here, ': Allocation failed'
stop -1
endif
! 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_cart_num
cholesky_ao(1:ao_cart_num,j,k) = L((j-1_8)*ao_cart_num+1_8:1_8*j*ao_cart_num,rank-k+1)
enddo
enddo
!$OMP END PARALLEL DO
call mmap_destroy(map)
cholesky_ao_cart_num = rank
if (write_ao_cart_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
close(iunit)
call ezfio_set_ao_cart_two_e_ints_io_ao_cart_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 *, ''
call wall_time(wall1)
print*,'Time to provide AO cholesky vectors = ',(wall1-wall0)/60.d0, ' min'
END_PROVIDER

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

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

703
src/ao_cart_two_e_ints/map_integrals.irp.f Normal file
View File

@ -0,0 +1,703 @@
use map_module
!! AO Map
!! ======
BEGIN_PROVIDER [ type(map_type), ao_cart_integrals_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_map,sze)
print*, 'AO map initialized : ', sze
END_PROVIDER
subroutine two_e_integrals_index(i,j,k,l,i1)
use map_module
implicit none
BEGIN_DOC
! Gives a unique index for i,j,k,l using permtuation symmetry.
! i <-> k, j <-> l, and (i,k) <-> (j,l) for non-periodic systems
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
subroutine two_e_integrals_index_reverse(i,j,k,l,i1)
use map_module
implicit none
BEGIN_DOC
! Computes the 4 indices $i,j,k,l$ from a unique index $i_1$.
! For 2 indices $i,j$ and $i \le j$, we have
! $p = i(i-1)/2 + j$.
! The key point is that because $j < i$,
! $i(i-1)/2 < p \le i(i+1)/2$. So $i$ can be found by solving
! $i^2 - i - 2p=0$. One obtains $i=1 + \sqrt{1+8p}/2$
! and $j = p - i(i-1)/2$.
! This rule is applied 3 times. First for the symmetry of the
! pairs (i,k) and (j,l), and then for the symmetry within each pair.
END_DOC
integer, intent(out) :: i(8),j(8),k(8),l(8)
integer(key_kind), intent(in) :: i1
integer(key_kind) :: i2,i3
i = 0
i2 = ceiling(0.5d0*(dsqrt(dble(shiftl(i1,3)+1))-1.d0))
l(1) = ceiling(0.5d0*(dsqrt(dble(shiftl(i2,3)+1))-1.d0))
i3 = i1 - shiftr(i2*i2-i2,1)
k(1) = ceiling(0.5d0*(dsqrt(dble(shiftl(i3,3)+1))-1.d0))
j(1) = int(i2 - shiftr(l(1)*l(1)-l(1),1),4)
i(1) = int(i3 - shiftr(k(1)*k(1)-k(1),1),4)
!ijkl
i(2) = i(1) !ilkj
j(2) = l(1)
k(2) = k(1)
l(2) = j(1)
i(3) = k(1) !kjil
j(3) = j(1)
k(3) = i(1)
l(3) = l(1)
i(4) = k(1) !klij
j(4) = l(1)
k(4) = i(1)
l(4) = j(1)
i(5) = j(1) !jilk
j(5) = i(1)
k(5) = l(1)
l(5) = k(1)
i(6) = j(1) !jkli
j(6) = k(1)
k(6) = l(1)
l(6) = i(1)
i(7) = l(1) !lijk
j(7) = i(1)
k(7) = j(1)
l(7) = k(1)
i(8) = l(1) !lkji
j(8) = k(1)
k(8) = j(1)
l(8) = i(1)
integer :: ii, jj
do ii=2,8
do jj=1,ii-1
if ( (i(ii) == i(jj)).and. &
(j(ii) == j(jj)).and. &
(k(ii) == k(jj)).and. &
(l(ii) == l(jj)) ) then
i(ii) = 0
exit
endif
enddo
enddo
! This has been tested with up to 1000 AOs, and all the reverse indices are
! correct ! We can remove the test
! do ii=1,8
! if (i(ii) /= 0) then
! call two_e_integrals_index(i(ii),j(ii),k(ii),l(ii),i2)
! if (i1 /= i2) then
! print *, i1, i2
! print *, i(ii), j(ii), k(ii), l(ii)
! stop 'two_e_integrals_index_reverse failed'
! endif
! endif
! enddo
end
subroutine ao_cart_idx2_sq(i,j,ij)
implicit none
integer, intent(in) :: i,j
integer, intent(out) :: ij
if (i<j) then
ij=(j-1)*(j-1)+2*i-mod(j+1,2)
else if (i>j) then
ij=(i-1)*(i-1)+2*j-mod(i,2)
else
ij=i*i
endif
end
subroutine idx2_tri_int(i,j,ij)
implicit none
integer, intent(in) :: i,j
integer, intent(out) :: ij
integer :: p,q
p = max(i,j)
q = min(i,j)
ij = q+ishft(p*p-p,-1)
end
subroutine ao_cart_idx2_tri_key(i,j,ij)
use map_module
implicit none
integer, intent(in) :: i,j
integer(key_kind), intent(out) :: ij
integer(key_kind) :: p,q
p = max(i,j)
q = min(i,j)
ij = q+ishft(p*p-p,-1)
end
subroutine two_e_integrals_index_2fold(i,j,k,l,i1)
use map_module
implicit none
integer, intent(in) :: i,j,k,l
integer(key_kind), intent(out) :: i1
integer :: ik,jl
call ao_cart_idx2_sq(i,k,ik)
call ao_cart_idx2_sq(j,l,jl)
call ao_cart_idx2_tri_key(ik,jl,i1)
end
subroutine ao_cart_idx2_sq_rev(i,k,ik)
BEGIN_DOC
! reverse square compound index
END_DOC
! p = ceiling(dsqrt(dble(ik)))
! q = ceiling(0.5d0*(dble(ik)-dble((p-1)*(p-1))))
! if (mod(ik,2)==0) then
! k=p
! i=q
! else
! i=p
! k=q
! endif
integer, intent(in) :: ik
integer, intent(out) :: i,k
integer :: pq(0:1),i1,i2
pq(0) = ceiling(dsqrt(dble(ik)))
pq(1) = ceiling(0.5d0*(dble(ik)-dble((pq(0)-1)*(pq(0)-1))))
i1=mod(ik,2)
i2=mod(ik+1,2)
k=pq(i1)
i=pq(i2)
end
subroutine ao_cart_idx2_tri_rev_key(i,k,ik)
use map_module
BEGIN_DOC
!return i<=k
END_DOC
integer(key_kind), intent(in) :: ik
integer, intent(out) :: i,k
integer(key_kind) :: tmp_k
k = ceiling(0.5d0*(dsqrt(8.d0*dble(ik)+1.d0)-1.d0))
tmp_k = k
i = int(ik - ishft(tmp_k*tmp_k-tmp_k,-1))
end
subroutine idx2_tri_rev_int(i,k,ik)
BEGIN_DOC
!return i<=k
END_DOC
integer, intent(in) :: ik
integer, intent(out) :: i,k
k = ceiling(0.5d0*(dsqrt(8.d0*dble(ik)+1.d0)-1.d0))
i = int(ik - ishft(k*k-k,-1))
end
subroutine two_e_integrals_index_reverse_2fold(i,j,k,l,i1)
use map_module
implicit none
integer, intent(out) :: i(2),j(2),k(2),l(2)
integer(key_kind), intent(in) :: i1
integer(key_kind) :: i0
integer :: i2,i3
i = 0
call ao_cart_idx2_tri_rev_key(i3,i2,i1)
call ao_cart_idx2_sq_rev(j(1),l(1),i2)
call ao_cart_idx2_sq_rev(i(1),k(1),i3)
!ijkl
i(2) = j(1) !jilk
j(2) = i(1)
k(2) = l(1)
l(2) = k(1)
! i(3) = k(1) !klij complex conjugate
! j(3) = l(1)
! k(3) = i(1)
! l(3) = j(1)
!
! i(4) = l(1) !lkji complex conjugate
! j(4) = k(1)
! k(4) = j(1)
! l(4) = i(1)
integer :: ii
if ( (i(1)==i(2)).and. &
(j(1)==j(2)).and. &
(k(1)==k(2)).and. &
(l(1)==l(2)) ) then
i(2) = 0
endif
! This has been tested with up to 1000 AOs, and all the reverse indices are
! correct ! We can remove the test
! do ii=1,2
! if (i(ii) /= 0) then
! call two_e_integrals_index_2fold(i(ii),j(ii),k(ii),l(ii),i0)
! if (i1 /= i0) then
! print *, i1, i0
! print *, i(ii), j(ii), k(ii), l(ii)
! stop 'two_e_integrals_index_reverse_2fold failed'
! endif
! endif
! enddo
end
BEGIN_PROVIDER [ integer, ao_cart_integrals_cache_min ]
&BEGIN_PROVIDER [ integer, ao_cart_integrals_cache_max ]
implicit none
BEGIN_DOC
! Min and max values of the AOs for which the integrals are in the cache
END_DOC
ao_cart_integrals_cache_min = max(1,ao_cart_num - 63)
ao_cart_integrals_cache_max = ao_cart_num
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_cart_integrals_cache, (0:64*64*64*64) ]
implicit none
BEGIN_DOC
! Cache of AO integrals for fast access
END_DOC
PROVIDE ao_cart_two_e_integrals_in_map
integer :: i,j,k,l,ii
integer(key_kind) :: idx, idx2
real(integral_kind) :: integral
real(integral_kind) :: tmp_re, tmp_im
integer(key_kind) :: idx_re,idx_im
!$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral)
do l=ao_cart_integrals_cache_min,ao_cart_integrals_cache_max
do k=ao_cart_integrals_cache_min,ao_cart_integrals_cache_max
do j=ao_cart_integrals_cache_min,ao_cart_integrals_cache_max
do i=ao_cart_integrals_cache_min,ao_cart_integrals_cache_max
!DIR$ FORCEINLINE
call two_e_integrals_index(i,j,k,l,idx)
!DIR$ FORCEINLINE
call map_get(ao_cart_integrals_map,idx,integral)
ii = l-ao_cart_integrals_cache_min
ii = ior( shiftl(ii,6), k-ao_cart_integrals_cache_min)
ii = ior( shiftl(ii,6), j-ao_cart_integrals_cache_min)
ii = ior( shiftl(ii,6), i-ao_cart_integrals_cache_min)
ao_cart_integrals_cache(ii) = integral
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
! ---
double precision function get_ao_cart_two_e_integral(i, j, k, l, map) result(result)
use map_module
implicit none
BEGIN_DOC
! Gets one AO bi-electronic integral from the AO map in PHYSICIST NOTATION
!
! <1:k, 2:l |1:i, 2:j>
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_in_map ao_cart_integrals_cache ao_cart_integrals_cache_min
!DIR$ FORCEINLINE
if (ao_cart_two_e_integral_zero(i,j,k,l)) then
tmp = 0.d0
else
ii = l-ao_cart_integrals_cache_min
ii = ior(ii, k-ao_cart_integrals_cache_min)
ii = ior(ii, j-ao_cart_integrals_cache_min)
ii = ior(ii, i-ao_cart_integrals_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)
else
ii = l-ao_cart_integrals_cache_min
ii = ior( shiftl(ii,6), k-ao_cart_integrals_cache_min)
ii = ior( shiftl(ii,6), j-ao_cart_integrals_cache_min)
ii = ior( shiftl(ii,6), i-ao_cart_integrals_cache_min)
tmp = ao_cart_integrals_cache(ii)
endif
endif
result = tmp
end
BEGIN_PROVIDER [ complex*16, ao_cart_integrals_cache_periodic, (0:64*64*64*64) ]
implicit none
BEGIN_DOC
! Cache of AO integrals for fast access
END_DOC
PROVIDE ao_cart_two_e_integrals_in_map
integer :: i,j,k,l,ii
integer(key_kind) :: idx1, idx2
real(integral_kind) :: tmp_re, tmp_im
integer(key_kind) :: idx_re,idx_im
complex(integral_kind) :: integral
!$OMP PARALLEL DO PRIVATE (i,j,k,l,idx1,idx2,tmp_re,tmp_im,idx_re,idx_im,ii,integral)
do l=ao_cart_integrals_cache_min,ao_cart_integrals_cache_max
do k=ao_cart_integrals_cache_min,ao_cart_integrals_cache_max
do j=ao_cart_integrals_cache_min,ao_cart_integrals_cache_max
do i=ao_cart_integrals_cache_min,ao_cart_integrals_cache_max
!DIR$ FORCEINLINE
call two_e_integrals_index_2fold(i,j,k,l,idx1)
!DIR$ FORCEINLINE
call two_e_integrals_index_2fold(k,l,i,j,idx2)
idx_re = min(idx1,idx2)
idx_im = max(idx1,idx2)
!DIR$ FORCEINLINE
call map_get(ao_cart_integrals_map,idx_re,tmp_re)
if (idx_re /= idx_im) then
call map_get(ao_cart_integrals_map,idx_im,tmp_im)
if (idx1 < idx2) then
integral = dcmplx(tmp_re,tmp_im)
else
integral = dcmplx(tmp_re,-tmp_im)
endif
else
tmp_im = 0.d0
integral = dcmplx(tmp_re,tmp_im)
endif
ii = l-ao_cart_integrals_cache_min
ii = ior( shiftl(ii,6), k-ao_cart_integrals_cache_min)
ii = ior( shiftl(ii,6), j-ao_cart_integrals_cache_min)
ii = ior( shiftl(ii,6), i-ao_cart_integrals_cache_min)
ao_cart_integrals_cache_periodic(ii) = integral
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
complex*16 function get_ao_cart_two_e_integral_periodic(i,j,k,l,map) result(result)
use map_module
implicit none
BEGIN_DOC
! Gets one AO bi-electronic integral from the AO map
END_DOC
integer, intent(in) :: i,j,k,l
integer(key_kind) :: idx1,idx2
real(integral_kind) :: tmp_re, tmp_im
integer(key_kind) :: idx_re,idx_im
type(map_type), intent(inout) :: map
integer :: ii
complex(integral_kind) :: tmp
PROVIDE ao_cart_two_e_integrals_in_map ao_cart_integrals_cache_periodic ao_cart_integrals_cache_min
!DIR$ FORCEINLINE
logical, external :: ao_cart_two_e_integral_zero
if (ao_cart_two_e_integral_zero(i,j,k,l)) then
tmp = (0.d0,0.d0)
else
ii = l-ao_cart_integrals_cache_min
ii = ior(ii, k-ao_cart_integrals_cache_min)
ii = ior(ii, j-ao_cart_integrals_cache_min)
ii = ior(ii, i-ao_cart_integrals_cache_min)
if (iand(ii, -64) /= 0) then
!DIR$ FORCEINLINE
call two_e_integrals_index_2fold(i,j,k,l,idx1)
!DIR$ FORCEINLINE
call two_e_integrals_index_2fold(k,l,i,j,idx2)
idx_re = min(idx1,idx2)
idx_im = max(idx1,idx2)
!DIR$ FORCEINLINE
call map_get(ao_cart_integrals_map,idx_re,tmp_re)
if (idx_re /= idx_im) then
call map_get(ao_cart_integrals_map,idx_im,tmp_im)
if (idx1 < idx2) then
tmp = dcmplx(tmp_re,tmp_im)
else
tmp = dcmplx(tmp_re,-tmp_im)
endif
else
tmp_im = 0.d0
tmp = dcmplx(tmp_re,tmp_im)
endif
else
ii = l-ao_cart_integrals_cache_min
ii = ior( shiftl(ii,6), k-ao_cart_integrals_cache_min)
ii = ior( shiftl(ii,6), j-ao_cart_integrals_cache_min)
ii = ior( shiftl(ii,6), i-ao_cart_integrals_cache_min)
tmp = ao_cart_integrals_cache_periodic(ii)
endif
result = tmp
endif
end
subroutine get_ao_cart_two_e_integrals(j,k,l,sze,out_val)
use map_module
BEGIN_DOC
! Gets multiple AO bi-electronic integral from the AO map .
! All i are retrieved for j,k,l fixed.
! physicist convention : <ij|kl>
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
logical, external :: ao_cart_one_e_integral_zero
PROVIDE ao_cart_two_e_integrals_in_map ao_cart_integrals_map
if (ao_cart_one_e_integral_zero(j,l)) then
out_val(1:sze) = 0.d0
return
endif
double precision :: get_ao_cart_two_e_integral
do i=1,sze
out_val(i) = get_ao_cart_two_e_integral(i,j,k,l,ao_cart_integrals_map)
enddo
end
subroutine get_ao_cart_two_e_integrals_periodic(j,k,l,sze,out_val)
use map_module
BEGIN_DOC
! Gets multiple AO bi-electronic integral from the AO map .
! All i are retrieved for j,k,l fixed.
! physicist convention : <ij|kl>
END_DOC
implicit none
integer, intent(in) :: j,k,l, sze
complex(integral_kind), intent(out) :: out_val(sze)
integer :: i
integer(key_kind) :: hash
logical, external :: ao_cart_one_e_integral_zero
PROVIDE ao_cart_two_e_integrals_in_map ao_cart_integrals_map
if (ao_cart_one_e_integral_zero(j,l)) then
out_val = 0.d0
return
endif
double precision :: get_ao_cart_two_e_integral
do i=1,sze
out_val(i) = get_ao_cart_two_e_integral(i,j,k,l,ao_cart_integrals_map)
enddo
end
subroutine get_ao_cart_two_e_integrals_non_zero(j,k,l,sze,out_val,out_val_index,non_zero_int)
use map_module
implicit none
BEGIN_DOC
! Gets multiple AO bi-electronic integral from the AO map .
! All non-zero i are retrieved for j,k,l fixed.
END_DOC
integer, intent(in) :: j,k,l, sze
real(integral_kind), intent(out) :: out_val(sze)
integer, intent(out) :: out_val_index(sze),non_zero_int
integer :: i
integer(key_kind) :: hash
double precision :: tmp
logical, external :: ao_cart_one_e_integral_zero
logical, external :: ao_cart_two_e_integral_zero
PROVIDE ao_cart_two_e_integrals_in_map
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
!DIR$ FORCEINLINE
if (ao_cart_two_e_integral_zero(i,j,k,l)) then
cycle
endif
call two_e_integrals_index(i,j,k,l,hash)
call map_get(ao_cart_integrals_map, hash,tmp)
if (dabs(tmp) < ao_cart_integrals_threshold) cycle
non_zero_int = non_zero_int+1
out_val_index(non_zero_int) = i
out_val(non_zero_int) = tmp
enddo
end
subroutine get_ao_cart_two_e_integrals_non_zero_jl(j,l,thresh,sze_max,sze,out_val,out_val_index,non_zero_int)
use map_module
implicit none
BEGIN_DOC
! Gets multiple AO bi-electronic integral from the AO map .
! All non-zero i are retrieved for j,k,l fixed.
END_DOC
double precision, intent(in) :: thresh
integer, intent(in) :: j,l, sze,sze_max
real(integral_kind), intent(out) :: out_val(sze_max)
integer, intent(out) :: out_val_index(2,sze_max),non_zero_int
integer :: i,k
integer(key_kind) :: hash
double precision :: tmp
logical, external :: ao_cart_one_e_integral_zero
logical, external :: ao_cart_two_e_integral_zero
PROVIDE ao_cart_two_e_integrals_in_map
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 k = 1, sze
do i = 1, sze
integer, external :: ao_cart_l4
double precision, external :: ao_cart_two_e_integral
!DIR$ FORCEINLINE
if (ao_cart_two_e_integral_zero(i,j,k,l)) then
cycle
endif
call two_e_integrals_index(i,j,k,l,hash)
call map_get(ao_cart_integrals_map, hash,tmp)
if (dabs(tmp) < thresh ) cycle
non_zero_int = non_zero_int+1
out_val_index(1,non_zero_int) = i
out_val_index(2,non_zero_int) = k
out_val(non_zero_int) = tmp
enddo
enddo
end
subroutine get_ao_cart_two_e_integrals_non_zero_jl_from_list(j,l,thresh,list,n_list,sze_max,out_val,out_val_index,non_zero_int)
use map_module
implicit none
BEGIN_DOC
! Gets multiple AO two-electron integrals from the AO map .
! All non-zero i are retrieved for j,k,l fixed.
END_DOC
double precision, intent(in) :: thresh
integer, intent(in) :: sze_max
integer, intent(in) :: j,l, n_list,list(2,sze_max)
real(integral_kind), intent(out) :: out_val(sze_max)
integer, intent(out) :: out_val_index(2,sze_max),non_zero_int
integer :: i,k
integer(key_kind) :: hash
double precision :: tmp
logical, external :: ao_cart_one_e_integral_zero
logical, external :: ao_cart_two_e_integral_zero
PROVIDE ao_cart_two_e_integrals_in_map
non_zero_int = 0
if (ao_cart_one_e_integral_zero(j,l)) then
out_val = 0.d0
return
endif
non_zero_int = 0
integer :: kk
do kk = 1, n_list
k = list(1,kk)
i = list(2,kk)
integer, external :: ao_cart_l4
double precision, external :: ao_cart_two_e_integral
!DIR$ FORCEINLINE
if (ao_cart_two_e_integral_zero(i,j,k,l)) then
cycle
endif
call two_e_integrals_index(i,j,k,l,hash)
call map_get(ao_cart_integrals_map, hash,tmp)
if (dabs(tmp) < thresh ) cycle
non_zero_int = non_zero_int+1
out_val_index(1,non_zero_int) = i
out_val_index(2,non_zero_int) = k
out_val(non_zero_int) = tmp
enddo
end
function get_ao_cart_map_size()
implicit none
integer (map_size_kind) :: get_ao_cart_map_size
BEGIN_DOC
! Returns the number of elements in the AO map
END_DOC
get_ao_cart_map_size = ao_cart_integrals_map % n_elements
end
subroutine clear_ao_cart_map
implicit none
BEGIN_DOC
! Frees the memory of the AO map
END_DOC
call map_deinit(ao_cart_integrals_map)
FREE ao_cart_integrals_map
end
subroutine insert_into_ao_cart_integrals_map(n_integrals,buffer_i, buffer_values)
use map_module
implicit none
BEGIN_DOC
! Create new entry into AO map
END_DOC
integer, intent(in) :: n_integrals
integer(key_kind), intent(inout) :: buffer_i(n_integrals)
real(integral_kind), intent(inout) :: buffer_values(n_integrals)
call map_append(ao_cart_integrals_map, buffer_i, buffer_values, n_integrals)
end

288
src/ao_cart_two_e_ints/map_integrals_erf.irp.f Normal file
View File

@ -0,0 +1,288 @@
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

112
src/ao_cart_two_e_ints/providers_ao_erf.irp.f Normal file
View File

@ -0,0 +1,112 @@
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

15
src/ao_cart_two_e_ints/screening.irp.f Normal file
View File

@ -0,0 +1,15 @@
logical function ao_cart_two_e_integral_zero(i,j,k,l)
implicit none
integer, intent(in) :: i,j,k,l
ao_cart_two_e_integral_zero = .False.
if (.not.(read_ao_cart_two_e_integrals.or.is_periodic.or.use_cgtos)) then
if (ao_cart_overlap_abs(j,l)*ao_cart_overlap_abs(i,k) < ao_cart_integrals_threshold) then
ao_cart_two_e_integral_zero = .True.
return
endif
if (ao_cart_two_e_integral_schwartz(j,l)*ao_cart_two_e_integral_schwartz(i,k) < ao_cart_integrals_threshold) then
ao_cart_two_e_integral_zero = .True.
endif
endif
end

17
src/ao_cart_two_e_ints/thresh.irp.f Normal file
View File

@ -0,0 +1,17 @@
BEGIN_PROVIDER [ double precision, ao_cart_integrals_threshold]
implicit none
ao_cart_integrals_threshold = ao_integrals_threshold
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_cart_cholesky_threshold]
implicit none
ao_cart_cholesky_threshold = ao_cholesky_threshold
END_PROVIDER
BEGIN_PROVIDER [ logical, do_ao_cart_cholesky]
implicit none
do_ao_cart_cholesky = do_ao_cholesky
END_PROVIDER

1674
src/ao_cart_two_e_ints/two_e_coul_integrals_cgtos.irp.f Normal file
View File

File diff suppressed because it is too large Load Diff

1790
src/ao_cart_two_e_ints/two_e_integrals.irp.f Normal file
View File

File diff suppressed because it is too large Load Diff

658
src/ao_cart_two_e_ints/two_e_integrals_erf.irp.f Normal file
View File

@ -0,0 +1,658 @@
double precision function ao_cart_two_e_integral_erf(i,j,k,l)
implicit none
BEGIN_DOC
! integral of the AO basis <ik|jl> 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_cart_two_e_integral_schwartz_accel_erf
provide mu_erf
if (ao_cart_prim_num(i) * ao_cart_prim_num(j) * ao_cart_prim_num(k) * ao_cart_prim_num(l) > 1024 ) then
ao_cart_two_e_integral_erf = ao_cart_two_e_integral_schwartz_accel_erf(i,j,k,l)
return
endif
dim1 = n_pt_max_integrals
num_i = ao_cart_nucl(i)
num_j = ao_cart_nucl(j)
num_k = ao_cart_nucl(k)
num_l = ao_cart_nucl(l)
ao_cart_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_cart_power(i,p)
J_power(p) = ao_cart_power(j,p)
K_power(p) = ao_cart_power(k,p)
L_power(p) = ao_cart_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_cart_prim_num(i)
coef1 = ao_cart_coef_normalized_ordered_transp(p,i)
do q = 1, ao_cart_prim_num(j)
coef2 = coef1*ao_cart_coef_normalized_ordered_transp(q,j)
call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,&
ao_cart_expo_ordered_transp(p,i),ao_cart_expo_ordered_transp(q,j), &
I_power,J_power,I_center,J_center,dim1)
p_inv = 1.d0/pp
do r = 1, ao_cart_prim_num(k)
coef3 = coef2*ao_cart_coef_normalized_ordered_transp(r,k)
do s = 1, ao_cart_prim_num(l)
coef4 = coef3*ao_cart_coef_normalized_ordered_transp(s,l)
call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,&
ao_cart_expo_ordered_transp(r,k),ao_cart_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_cart_two_e_integral_erf = ao_cart_two_e_integral_erf + coef4 * integral
enddo ! s
enddo ! r
enddo ! q
enddo ! p
else
do p = 1, 3
I_power(p) = ao_cart_power(i,p)
J_power(p) = ao_cart_power(j,p)
K_power(p) = ao_cart_power(k,p)
L_power(p) = ao_cart_power(l,p)
enddo
double precision :: ERI_erf
do p = 1, ao_cart_prim_num(i)
coef1 = ao_cart_coef_normalized_ordered_transp(p,i)
do q = 1, ao_cart_prim_num(j)
coef2 = coef1*ao_cart_coef_normalized_ordered_transp(q,j)
do r = 1, ao_cart_prim_num(k)
coef3 = coef2*ao_cart_coef_normalized_ordered_transp(r,k)
do s = 1, ao_cart_prim_num(l)
coef4 = coef3*ao_cart_coef_normalized_ordered_transp(s,l)
integral = ERI_erf( &
ao_cart_expo_ordered_transp(p,i),ao_cart_expo_ordered_transp(q,j),ao_cart_expo_ordered_transp(r,k),ao_cart_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_cart_two_e_integral_erf = ao_cart_two_e_integral_erf + coef4 * integral
enddo ! s
enddo ! r
enddo ! q
enddo ! p
endif
end
double precision function ao_cart_two_e_integral_schwartz_accel_erf(i,j,k,l)
implicit none
BEGIN_DOC
! integral of the AO basis <ik|jl> 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_cart_nucl(i)
num_j = ao_cart_nucl(j)
num_k = ao_cart_nucl(k)
num_l = ao_cart_nucl(l)
ao_cart_two_e_integral_schwartz_accel_erf = 0.d0
double precision :: thr
thr = ao_cart_integrals_threshold*ao_cart_integrals_threshold
allocate(schwartz_kl(0:ao_cart_prim_num(l),0:ao_cart_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_cart_power(i,p)
J_power(p) = ao_cart_power(j,p)
K_power(p) = ao_cart_power(k,p)
L_power(p) = ao_cart_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_cart_prim_num(k)
coef1 = ao_cart_coef_normalized_ordered_transp(r,k)*ao_cart_coef_normalized_ordered_transp(r,k)
schwartz_kl(0,r) = 0.d0
do s = 1, ao_cart_prim_num(l)
coef2 = coef1 * ao_cart_coef_normalized_ordered_transp(s,l) * ao_cart_coef_normalized_ordered_transp(s,l)
call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,&
ao_cart_expo_ordered_transp(r,k),ao_cart_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_cart_prim_num(i)
coef1 = ao_cart_coef_normalized_ordered_transp(p,i)
do q = 1, ao_cart_prim_num(j)
coef2 = coef1*ao_cart_coef_normalized_ordered_transp(q,j)
call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,&
ao_cart_expo_ordered_transp(p,i),ao_cart_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_cart_prim_num(k)
if (schwartz_kl(0,r)*schwartz_ij < thr) then
cycle
endif
coef3 = coef2*ao_cart_coef_normalized_ordered_transp(r,k)
do s = 1, ao_cart_prim_num(l)
if (schwartz_kl(s,r)*schwartz_ij < thr) then
cycle
endif
coef4 = coef3*ao_cart_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_cart_expo_ordered_transp(r,k),ao_cart_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_cart_two_e_integral_schwartz_accel_erf = ao_cart_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_cart_power(i,p)
J_power(p) = ao_cart_power(j,p)
K_power(p) = ao_cart_power(k,p)
L_power(p) = ao_cart_power(l,p)
enddo
double precision :: ERI_erf
schwartz_kl(0,0) = 0.d0
do r = 1, ao_cart_prim_num(k)
coef1 = ao_cart_coef_normalized_ordered_transp(r,k)*ao_cart_coef_normalized_ordered_transp(r,k)
schwartz_kl(0,r) = 0.d0
do s = 1, ao_cart_prim_num(l)
coef2 = coef1*ao_cart_coef_normalized_ordered_transp(s,l)*ao_cart_coef_normalized_ordered_transp(s,l)
schwartz_kl(s,r) = ERI_erf( &
ao_cart_expo_ordered_transp(r,k),ao_cart_expo_ordered_transp(s,l),ao_cart_expo_ordered_transp(r,k),ao_cart_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_cart_prim_num(i)
coef1 = ao_cart_coef_normalized_ordered_transp(p,i)
do q = 1, ao_cart_prim_num(j)
coef2 = coef1*ao_cart_coef_normalized_ordered_transp(q,j)
schwartz_ij = ERI_erf( &
ao_cart_expo_ordered_transp(p,i),ao_cart_expo_ordered_transp(q,j),ao_cart_expo_ordered_transp(p,i),ao_cart_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_cart_prim_num(k)
if (schwartz_kl(0,r)*schwartz_ij < thr) then
cycle
endif
coef3 = coef2*ao_cart_coef_normalized_ordered_transp(r,k)
do s = 1, ao_cart_prim_num(l)
if (schwartz_kl(s,r)*schwartz_ij < thr) then
cycle
endif
coef4 = coef3*ao_cart_coef_normalized_ordered_transp(s,l)
integral = ERI_erf( &
ao_cart_expo_ordered_transp(p,i),ao_cart_expo_ordered_transp(q,j),ao_cart_expo_ordered_transp(r,k),ao_cart_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_cart_two_e_integral_schwartz_accel_erf = ao_cart_two_e_integral_schwartz_accel_erf + coef4 * integral
enddo ! s
enddo ! r
enddo ! q
enddo ! p
endif
deallocate (schwartz_kl)
end
subroutine compute_ao_cart_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_cart_two_e_integral_erf
integer :: i
logical, external :: ao_cart_one_e_integral_zero
logical, external :: ao_cart_two_e_integral_zero
if (ao_cart_one_e_integral_zero(j,l)) then
buffer_value = 0._integral_kind
return
endif
if (ao_cart_two_e_integral_erf_schwartz(j,l) < thresh ) then
buffer_value = 0._integral_kind
return
endif
do i = 1, ao_cart_num
if (ao_cart_two_e_integral_zero(i,j,k,l)) then
buffer_value(i) = 0._integral_kind
cycle
endif
if (ao_cart_two_e_integral_erf_schwartz(i,k)*ao_cart_two_e_integral_erf_schwartz(j,l) < thresh ) then
buffer_value(i) = 0._integral_kind
cycle
endif
!DIR$ FORCEINLINE
buffer_value(i) = ao_cart_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 <pq|rs> 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_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

15
src/ao_two_e_ints/screening.irp.f
View File

@ -1,15 +0,0 @@
logical function ao_two_e_integral_zero(i,j,k,l)
implicit none
integer, intent(in) :: i,j,k,l
ao_two_e_integral_zero = .False.
if (.not.(read_ao_two_e_integrals.or.is_periodic.or.use_cgtos)) then
if (ao_overlap_abs(j,l)*ao_overlap_abs(i,k) < ao_integrals_threshold) then
ao_two_e_integral_zero = .True.
return
endif
if (ao_two_e_integral_schwartz(j,l)*ao_two_e_integral_schwartz(i,k) < ao_integrals_threshold) then
ao_two_e_integral_zero = .True.
endif
endif
end

25
src/two_e_ints_keywords/EZFIO.cfg Normal file
View File

@ -0,0 +1,25 @@
[ao_integrals_threshold]
type: Threshold
doc: If | (pq|rs) | < `ao_integrals_threshold` then (pq|rs) is zero
interface: ezfio,provider,ocaml
default: 1.e-15
ezfio_name: threshold_ao
[ao_cholesky_threshold]
type: Threshold
doc: If | (ii|jj) | < `ao_cholesky_threshold` then (ii|jj) is zero
interface: ezfio,provider,ocaml
default: 1.e-12
[do_ao_cholesky]
type: logical
doc: Perform Cholesky decomposition of AO integrals
interface: ezfio,provider,ocaml
default: True
[use_only_lr]
type: logical
doc: If true, use only the long range part of the two-electron integrals instead of 1/r12
interface: ezfio, provider, ocaml
default: False

11
src/two_e_ints_keywords/direct.irp.f Normal file
View File

@ -0,0 +1,11 @@
BEGIN_PROVIDER [ logical, do_direct_integrals ]
implicit none
BEGIN_DOC
! Compute integrals on the fly
END_DOC
do_direct_integrals = do_ao_cholesky
END_PROVIDER