diff --git a/plugins/local/tc_keywords/EZFIO.cfg b/plugins/local/tc_keywords/EZFIO.cfg deleted file mode 100644 index b858fa5b..00000000 --- a/plugins/local/tc_keywords/EZFIO.cfg +++ /dev/null @@ -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 diff --git a/plugins/local/tc_keywords/NEED b/plugins/local/tc_keywords/NEED deleted file mode 100644 index f1c051ff..00000000 --- a/plugins/local/tc_keywords/NEED +++ /dev/null @@ -1,2 +0,0 @@ -ezfio_files -nuclei diff --git a/src/ao_cart_two_e_ints/NEED b/src/ao_cart_two_e_ints/NEED new file mode 100644 index 00000000..dba58c3d --- /dev/null +++ b/src/ao_cart_two_e_ints/NEED @@ -0,0 +1,6 @@ +hamiltonian +ao_one_e_ints +pseudo +bitmask +ao_basis +two_e_ints_keywords diff --git a/src/ao_cart_two_e_ints/README.rst b/src/ao_cart_two_e_ints/README.rst new file mode 100644 index 00000000..f2dd5206 --- /dev/null +++ b/src/ao_cart_two_e_ints/README.rst @@ -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) = = <12|12> + + + diff --git a/src/ao_cart_two_e_ints/cholesky.irp.f b/src/ao_cart_two_e_ints/cholesky.irp.f new file mode 100644 index 00000000..8de035b6 --- /dev/null +++ b/src/ao_cart_two_e_ints/cholesky.irp.f @@ -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 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): + ! = (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 + diff --git a/src/ao_cart_two_e_ints/gauss_legendre.irp.f b/src/ao_cart_two_e_ints/gauss_legendre.irp.f new file mode 100644 index 00000000..4bdadb6e --- /dev/null +++ b/src/ao_cart_two_e_ints/gauss_legendre.irp.f @@ -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 + diff --git a/src/ao_cart_two_e_ints/map_integrals.irp.f b/src/ao_cart_two_e_ints/map_integrals.irp.f new file mode 100644 index 00000000..bf0791f1 --- /dev/null +++ b/src/ao_cart_two_e_ints/map_integrals.irp.f @@ -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 (ij) 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 : + 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 : + 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 + + 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 new file mode 100644 index 00000000..f24f8080 --- /dev/null +++ b/src/ao_cart_two_e_ints/map_integrals_erf.irp.f @@ -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 + + + + 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 new file mode 100644 index 00000000..a01334da --- /dev/null +++ b/src/ao_cart_two_e_ints/providers_ao_erf.irp.f @@ -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 + + diff --git a/src/ao_cart_two_e_ints/screening.irp.f b/src/ao_cart_two_e_ints/screening.irp.f new file mode 100644 index 00000000..805ef279 --- /dev/null +++ b/src/ao_cart_two_e_ints/screening.irp.f @@ -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 diff --git a/src/ao_cart_two_e_ints/thresh.irp.f b/src/ao_cart_two_e_ints/thresh.irp.f new file mode 100644 index 00000000..5dd7c251 --- /dev/null +++ b/src/ao_cart_two_e_ints/thresh.irp.f @@ -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 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_coul_integrals_cgtos.irp.f new file mode 100644 index 00000000..193573c3 --- /dev/null +++ b/src/ao_cart_two_e_ints/two_e_coul_integrals_cgtos.irp.f @@ -0,0 +1,1674 @@ + +! --- + +double precision function ao_cart_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_cart_2e_cgtos_schwartz_accel + complex*16, external :: ERI_cgtos + complex*16, external :: general_primitive_integral_cgtos + + + + 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_cgtos = ao_cart_2e_cgtos_schwartz_accel(i, j, k, l) + + return + endif + + dim1 = n_pt_max_integrals + + ii = ao_cart_nucl(i) + jj = ao_cart_nucl(j) + kk = ao_cart_nucl(k) + ll = ao_cart_nucl(l) + + do m = 1, 3 + I_power(m) = ao_cart_power(i,m) + J_power(m) = ao_cart_power(j,m) + K_power(m) = ao_cart_power(k,m) + L_power(m) = ao_cart_power(l,m) + enddo + + + ao_cart_two_e_integral_cgtos = 0.d0 + + if(use_pw .or. ii /= jj .or. kk /= ll .or. jj /= kk) then + + do p = 1, ao_cart_prim_num(i) + + coef1 = ao_cart_coef_cgtos_norm_ord_transp(p,i) + expo1 = ao_cart_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_cart_expo_pw_ord_transp(m,p,i) + enddo + phiI = ao_cart_expo_phase_ord_transp(4,p,i) + KI2 = ao_cart_expo_pw_ord_transp(4,p,i) + + do q = 1, ao_cart_prim_num(j) + + coef2 = coef1 * ao_cart_coef_cgtos_norm_ord_transp(q,j) + expo2 = ao_cart_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_cart_expo_pw_ord_transp(m,q,j) + enddo + phiJ = ao_cart_expo_phase_ord_transp(4,q,j) + KJ2 = ao_cart_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_cart_prim_num(k) + + coef3 = coef2 * ao_cart_coef_cgtos_norm_ord_transp(r,k) + expo3 = ao_cart_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_cart_expo_pw_ord_transp(m,r,k) + enddo + phiK = ao_cart_expo_phase_ord_transp(4,r,k) + KK2 = ao_cart_expo_pw_ord_transp(4,r,k) + + do s = 1, ao_cart_prim_num(l) + + coef4 = coef3 * ao_cart_coef_cgtos_norm_ord_transp(s,l) + expo4 = ao_cart_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_cart_expo_pw_ord_transp(m,s,l) + enddo + phiL = ao_cart_expo_phase_ord_transp(4,s,l) + KL2 = ao_cart_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_cart_two_e_integral_cgtos = ao_cart_two_e_integral_cgtos + coef4 * 2.d0 * real(int_tot) + enddo ! s + enddo ! r + enddo ! q + enddo ! p + + else + + do p = 1, ao_cart_prim_num(i) + + coef1 = ao_cart_coef_cgtos_norm_ord_transp(p,i) + expo1 = ao_cart_expo_cgtos_ord_transp(p,i) + phiI = ao_cart_expo_phase_ord_transp(4,p,i) + + do q = 1, ao_cart_prim_num(j) + + coef2 = coef1 * ao_cart_coef_cgtos_norm_ord_transp(q,j) + expo2 = ao_cart_expo_cgtos_ord_transp(q,j) + phiJ = ao_cart_expo_phase_ord_transp(4,q,j) + + do r = 1, ao_cart_prim_num(k) + + coef3 = coef2 * ao_cart_coef_cgtos_norm_ord_transp(r,k) + expo3 = ao_cart_expo_cgtos_ord_transp(r,k) + phiK = ao_cart_expo_phase_ord_transp(4,r,k) + + do s = 1, ao_cart_prim_num(l) + + coef4 = coef3 * ao_cart_coef_cgtos_norm_ord_transp(s,l) + expo4 = ao_cart_expo_cgtos_ord_transp(s,l) + phiL = ao_cart_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_cart_two_e_integral_cgtos = ao_cart_two_e_integral_cgtos + coef4 * 2.d0 * real(int_tot) + enddo ! s + enddo ! r + enddo ! q + enddo ! p + + endif + +end + +! --- + +double precision function ao_cart_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_cart_2e_cgtos_schwartz_accel = 0.d0 + + dim1 = n_pt_max_integrals + + ii = ao_cart_nucl(i) + jj = ao_cart_nucl(j) + kk = ao_cart_nucl(k) + ll = ao_cart_nucl(l) + + do m = 1, 3 + I_power(m) = ao_cart_power(i,m) + J_power(m) = ao_cart_power(j,m) + K_power(m) = ao_cart_power(k,m) + L_power(m) = ao_cart_power(l,m) + enddo + + + thr = ao_cart_integrals_threshold*ao_cart_integrals_threshold + + allocate(schwartz_kl(0:ao_cart_prim_num(l),0:ao_cart_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_cart_prim_num(k) + + coef1 = ao_cart_coef_cgtos_norm_ord_transp(r,k) * ao_cart_coef_cgtos_norm_ord_transp(r,k) + expo1 = ao_cart_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_cart_expo_pw_ord_transp(m,r,k) + enddo + phiK = ao_cart_expo_phase_ord_transp(4,r,k) + KK2 = ao_cart_expo_pw_ord_transp(4,r,k) + + schwartz_kl(0,r) = 0.d0 + do s = 1, ao_cart_prim_num(l) + + coef2 = coef1 * ao_cart_coef_cgtos_norm_ord_transp(s,l) * ao_cart_coef_cgtos_norm_ord_transp(s,l) + expo2 = ao_cart_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_cart_expo_pw_ord_transp(m,s,l) + enddo + phiL = ao_cart_expo_phase_ord_transp(4,s,l) + KL2 = ao_cart_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_cart_prim_num(i) + + coef1 = ao_cart_coef_cgtos_norm_ord_transp(p,i) + expo1 = ao_cart_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_cart_expo_pw_ord_transp(m,p,i) + enddo + phiI = ao_cart_expo_phase_ord_transp(4,p,i) + KI2 = ao_cart_expo_pw_ord_transp(4,p,i) + + do q = 1, ao_cart_prim_num(j) + + coef2 = coef1 * ao_cart_coef_cgtos_norm_ord_transp(q,j) + expo2 = ao_cart_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_cart_expo_pw_ord_transp(m,q,j) + enddo + phiJ = ao_cart_expo_phase_ord_transp(4,q,j) + KJ2 = ao_cart_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_cart_prim_num(k) + if(schwartz_kl(0,r)*schwartz_ij < thr) cycle + + coef3 = coef2 * ao_cart_coef_cgtos_norm_ord_transp(r,k) + expo3 = ao_cart_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_cart_expo_pw_ord_transp(m,r,k) + enddo + phiK = ao_cart_expo_phase_ord_transp(4,r,k) + KK2 = ao_cart_expo_pw_ord_transp(4,r,k) + + do s = 1, ao_cart_prim_num(l) + if(schwartz_kl(s,r)*schwartz_ij < thr) cycle + + coef4 = coef3 * ao_cart_coef_cgtos_norm_ord_transp(s,l) + expo4 = ao_cart_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_cart_expo_pw_ord_transp(m,s,l) + enddo + phiL = ao_cart_expo_phase_ord_transp(4,s,l) + KL2 = ao_cart_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_cart_2e_cgtos_schwartz_accel = ao_cart_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_cart_prim_num(k) + + coef1 = ao_cart_coef_cgtos_norm_ord_transp(r,k) * ao_cart_coef_cgtos_norm_ord_transp(r,k) + expo1 = ao_cart_expo_cgtos_ord_transp(r,k) + phiK = ao_cart_expo_phase_ord_transp(4,r,k) + + schwartz_kl(0,r) = 0.d0 + do s = 1, ao_cart_prim_num(l) + + coef2 = coef1 * ao_cart_coef_cgtos_norm_ord_transp(s,l) * ao_cart_coef_cgtos_norm_ord_transp(s,l) + expo2 = ao_cart_expo_cgtos_ord_transp(s,l) + phiL = ao_cart_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_cart_prim_num(i) + + coef1 = ao_cart_coef_cgtos_norm_ord_transp(p,i) + expo1 = ao_cart_expo_cgtos_ord_transp(p,i) + phiI = ao_cart_expo_phase_ord_transp(4,p,i) + + do q = 1, ao_cart_prim_num(j) + + coef2 = coef1 * ao_cart_coef_cgtos_norm_ord_transp(q,j) + expo2 = ao_cart_expo_cgtos_ord_transp(q,j) + phiJ = ao_cart_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_cart_prim_num(k) + if(schwartz_kl(0,r)*schwartz_ij < thr) cycle + + coef3 = coef2 * ao_cart_coef_cgtos_norm_ord_transp(r,k) + expo3 = ao_cart_expo_cgtos_ord_transp(r,k) + phiK = ao_cart_expo_phase_ord_transp(4,r,k) + + do s = 1, ao_cart_prim_num(l) + if(schwartz_kl(s,r)*schwartz_ij < thr) cycle + + coef4 = coef3 * ao_cart_coef_cgtos_norm_ord_transp(s,l) + expo4 = ao_cart_expo_cgtos_ord_transp(s,l) + phiL = ao_cart_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_cart_2e_cgtos_schwartz_accel = ao_cart_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_cart_2e_cgtos_schwartz, (ao_cart_num, ao_cart_num)] + + BEGIN_DOC + ! Needed to compute Schwartz inequalities + END_DOC + + implicit none + integer :: i, k + double precision :: ao_cart_two_e_integral_cgtos + + ao_cart_2e_cgtos_schwartz(1,1) = ao_cart_two_e_integral_cgtos(1, 1, 1, 1) + + !$OMP PARALLEL DO PRIVATE(i,k) & + !$OMP DEFAULT(NONE) & + !$OMP SHARED(ao_cart_num, ao_cart_2e_cgtos_schwartz) & + !$OMP SCHEDULE(dynamic) + do i = 1, ao_cart_num + do k = 1, i + ao_cart_2e_cgtos_schwartz(i,k) = dsqrt(ao_cart_two_e_integral_cgtos(i, i, k, k)) + ao_cart_2e_cgtos_schwartz(k,i) = ao_cart_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_cart_two_e_ints/two_e_integrals.irp.f b/src/ao_cart_two_e_ints/two_e_integrals.irp.f new file mode 100644 index 00000000..3df0d015 --- /dev/null +++ b/src/ao_cart_two_e_ints/two_e_integrals.irp.f @@ -0,0 +1,1790 @@ + +! --- + +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_cart_cholesky) then + do_schwartz_accel = .False. + else + do_schwartz_accel = (ao_cart_prim_num(i) * ao_cart_prim_num(j) * & + ao_cart_prim_num(k) * ao_cart_prim_num(l) > 1024 ) + endif + +end function + + +double precision function ao_cart_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_cart_two_e_integral_erf + double precision, external :: ao_cart_two_e_integral_cgtos + double precision, external :: ao_cart_two_e_integral_schwartz_accel + double precision, external :: ao_cart_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_cart_two_e_integral = ao_cart_two_e_integral_cgtos(i, j, k, l) + return + + else if (use_only_lr) then + + ao_cart_two_e_integral = ao_cart_two_e_integral_erf(i, j, k, l) + return + + else if (do_schwartz_accel(i,j,k,l)) then + + ao_cart_two_e_integral = ao_cart_two_e_integral_schwartz_accel(i,j,k,l) + + else + + 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 = 0.d0 + + if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k) then + + ao_cart_two_e_integral = ao_cart_two_e_integral_general(i,j,k,l,general_primitive_integral) + + 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 + + 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( & + 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 = ao_cart_two_e_integral + coef4 * integral + enddo ! s + enddo ! r + enddo ! q + enddo ! p + + endif + + endif + +end + +! --- + +double precision function ao_cart_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_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_general = 0.d0 + + 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 + + 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 = 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_cart_two_e_integral_general = ao_cart_two_e_integral_general + coef4 * integral + enddo ! s + enddo ! r + enddo ! q + enddo ! p + +end + +double precision function ao_cart_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_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 = 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))) + + + 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(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) + double precision :: coef1 + coef1 = ao_cart_coef_normalized_ordered_transp(p,i) + do q = 1, ao_cart_prim_num(j) + double precision :: coef2 + coef2 = coef1*ao_cart_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_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(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 + double precision :: coef3 + coef3 = coef2*ao_cart_coef_normalized_ordered_transp(r,k) + do s = 1, ao_cart_prim_num(l) + double precision :: coef4 + 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 + 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(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 = ao_cart_two_e_integral_schwartz_accel + 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 + + 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( & + 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( & + 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( & + 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 = ao_cart_two_e_integral_schwartz_accel + coef4 * integral + enddo ! s + enddo ! r + enddo ! q + enddo ! p + + endif + deallocate (schwartz_kl) + +end + + +integer function ao_cart_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_cart_l4 = ao_cart_l(i)*ao_cart_l(j)*ao_cart_l(k)*ao_cart_l(l) +end + + + +subroutine compute_ao_cart_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_cart_two_e_integral + + 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 + + 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 + !DIR$ FORCEINLINE + buffer_value(i) = ao_cart_two_e_integral(i,k,j,l) + enddo + +end + +BEGIN_PROVIDER [ logical, ao_cart_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_cart_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_cart_two_e_integrals io_ao_cart_two_e_integrals ao_cart_integrals_map + + if (read_ao_cart_two_e_integrals) then + print*,'Reading the AO integrals' + call map_load_from_disk(trim(ezfio_filename)//'/work/ao_cart_ints',ao_cart_integrals_map) + print*, 'AO integrals provided' + ao_cart_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_cart_two_e_integral(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_jl(j,l,n_integrals,buffer_i,buffer_value) + call insert_into_ao_cart_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_cart_integrals_map) + call cpu_time(cpu_2) + call wall_time(wall_2) + integer(map_size_kind) :: get_ao_cart_map_size, ao_cart_map_size + ao_cart_map_size = get_ao_cart_map_size() + + print*, 'AO integrals provided:' + print*, ' Size of AO map : ', map_mb(ao_cart_integrals_map) ,'MB' + print*, ' Number of AO integrals :', ao_cart_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_in_map = .True. + + if (write_ao_cart_two_e_integrals.and.mpi_master) then + call ezfio_set_work_empty(.False.) + call map_save_to_disk(trim(ezfio_filename)//'/work/ao_cart_ints',ao_cart_integrals_map) + call ezfio_set_ao_cart_two_e_ints_io_ao_cart_two_e_integrals('Read') + endif + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, ao_cart_two_e_integral_schwartz, (ao_cart_num, ao_cart_num) ] + + BEGIN_DOC + ! Needed to compute Schwartz inequalities + END_DOC + + implicit none + integer :: i, k + double precision :: ao_cart_two_e_integral,cpu_1,cpu_2, wall_1, wall_2 + + ao_cart_two_e_integral_schwartz(1,1) = ao_cart_two_e_integral(1,1,1,1) + !$OMP PARALLEL DO PRIVATE(i,k) & + !$OMP DEFAULT(NONE) & + !$OMP SHARED (ao_cart_num,ao_cart_two_e_integral_schwartz) & + !$OMP SCHEDULE(dynamic) + do i=ao_cart_num,1,-1 + do k=1,i + ao_cart_two_e_integral_schwartz(i,k) = dsqrt(ao_cart_two_e_integral(i,i,k,k)) + ao_cart_two_e_integral_schwartz(k,i) = ao_cart_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_cart_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_cart_num*ao_cart_num) + real(integral_kind),intent(out) :: buffer_value(ao_cart_num*ao_cart_num) + logical, external :: ao_cart_two_e_integral_zero + + integer :: i,k + double precision, external :: ao_cart_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_cart_integrals_threshold + + n_integrals = 0 + + j1 = j+shiftr(l*l-l,1) + do k = 1, ao_cart_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_cart_two_e_integral_zero(i,j,k,l)) then + cycle + endif + !DIR$ FORCEINLINE + integral = ao_cart_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_cart_two_e_ints/two_e_integrals_erf.irp.f b/src/ao_cart_two_e_ints/two_e_integrals_erf.irp.f new file mode 100644 index 00000000..9a06475c --- /dev/null +++ b/src/ao_cart_two_e_ints/two_e_integrals_erf.irp.f @@ -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 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 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 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 diff --git a/src/ao_two_e_ints/screening.irp.f b/src/ao_two_e_ints/screening.irp.f deleted file mode 100644 index 8c80af6c..00000000 --- a/src/ao_two_e_ints/screening.irp.f +++ /dev/null @@ -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 diff --git a/src/two_e_ints_keywords/EZFIO.cfg b/src/two_e_ints_keywords/EZFIO.cfg new file mode 100644 index 00000000..1be271ac --- /dev/null +++ b/src/two_e_ints_keywords/EZFIO.cfg @@ -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 diff --git a/src/two_e_ints_keywords/direct.irp.f b/src/two_e_ints_keywords/direct.irp.f new file mode 100644 index 00000000..ce093a57 --- /dev/null +++ b/src/two_e_ints_keywords/direct.irp.f @@ -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 + +