From 0c84015982f554740b0597a0a8bb2a37e104ec84 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 16 Nov 2015 19:00:49 +0100 Subject: [PATCH 1/6] Added level shift in Hartree-Fock --- plugins/Hartree_Fock/EZFIO.cfg | 6 ++++++ plugins/Hartree_Fock/Fock_matrix.irp.f | 11 +++++++++-- 2 files changed, 15 insertions(+), 2 deletions(-) diff --git a/plugins/Hartree_Fock/EZFIO.cfg b/plugins/Hartree_Fock/EZFIO.cfg index c39c3483..e3cb3673 100644 --- a/plugins/Hartree_Fock/EZFIO.cfg +++ b/plugins/Hartree_Fock/EZFIO.cfg @@ -10,6 +10,12 @@ doc: Maximum number of SCF iterations interface: ezfio,provider,ocaml default: 200 +[level_shift] +type: Positive_float +doc: Energy shift on the virtual MOs +interface: ezfio,provider,ocaml +default: 0.0 + [mo_guess_type] type: MO_guess doc: Initial MO guess. Can be [ Huckel | HCore ] diff --git a/plugins/Hartree_Fock/Fock_matrix.irp.f b/plugins/Hartree_Fock/Fock_matrix.irp.f index 2561ad03..60413ab9 100644 --- a/plugins/Hartree_Fock/Fock_matrix.irp.f +++ b/plugins/Hartree_Fock/Fock_matrix.irp.f @@ -73,8 +73,13 @@ enddo endif + ! Introduce level shift here + do i = elec_alpha_num+1, mo_tot_num + Fock_matrix_mo(i,i) += level_shift + enddo + do i = 1, mo_tot_num - Fock_matrix_diag_mo(i) = Fock_matrix_mo(i,i) + Fock_matrix_diag_mo(i) = Fock_matrix_mo(i,i) enddo END_PROVIDER @@ -315,7 +320,9 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num_align, ao_num) ] ! Fock matrix in AO basis set END_DOC - if (elec_alpha_num == elec_beta_num) then + if ( (elec_alpha_num == elec_beta_num).and. & + (level_shift == 0.) ) & + then integer :: i,j do j=1,ao_num !DIR$ VECTOR ALIGNED From 1e336aeca1201ab3b5ddf0a623c5d7d91bdf8346 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 16 Nov 2015 20:13:06 +0100 Subject: [PATCH 2/6] Added Schwartz in AOs --- plugins/Hartree_Fock/EZFIO.cfg | 2 +- src/Integrals_Bielec/ao_bi_integrals.irp.f | 8 ++++++-- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/plugins/Hartree_Fock/EZFIO.cfg b/plugins/Hartree_Fock/EZFIO.cfg index e3cb3673..1b8486ee 100644 --- a/plugins/Hartree_Fock/EZFIO.cfg +++ b/plugins/Hartree_Fock/EZFIO.cfg @@ -12,7 +12,7 @@ default: 200 [level_shift] type: Positive_float -doc: Energy shift on the virtual MOs +doc: Energy shift on the virtual MOs to improve SCF convergence interface: ezfio,provider,ocaml default: 0.0 diff --git a/src/Integrals_Bielec/ao_bi_integrals.irp.f b/src/Integrals_Bielec/ao_bi_integrals.irp.f index aa186454..4d1b2f95 100644 --- a/src/Integrals_Bielec/ao_bi_integrals.irp.f +++ b/src/Integrals_Bielec/ao_bi_integrals.irp.f @@ -378,8 +378,9 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ] !$OMP DEFAULT(NONE) & !$OMP SHARED (ao_num, jl_pairs, ao_integrals_map,thresh, & !$OMP cpu_1,wall_1,lock, lmax,n_centers,ao_nucl, & - !$OMP ao_overlap_abs,ao_overlap,abort_here, & - !$OMP wall_0,progress_bar,progress_value) + !$OMP ao_overlap_abs,ao_overlap,abort_here, & + !$OMP wall_0,progress_bar,progress_value, & + !$OMP ao_bielec_integral_schwartz) allocate(buffer_i(size_buffer)) allocate(buffer_value(size_buffer)) @@ -418,6 +419,9 @@ IRP_ENDIF if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then cycle endif + if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < thresh ) then + cycle + endif !DIR$ FORCEINLINE integral = ao_bielec_integral(i,k,j,l) if (abs(integral) < thresh) then From 1d2631f7dfe56c731f3f1185aaaabf87e4f329e3 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 16 Nov 2015 20:14:25 +0100 Subject: [PATCH 3/6] Print in QMC --- plugins/QmcChem/target_pt2_qmc.irp.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/plugins/QmcChem/target_pt2_qmc.irp.f b/plugins/QmcChem/target_pt2_qmc.irp.f index 6e6d2a60..228bcb5e 100644 --- a/plugins/QmcChem/target_pt2_qmc.irp.f +++ b/plugins/QmcChem/target_pt2_qmc.irp.f @@ -62,7 +62,7 @@ program e_curve endif enddo call compute_energy(psi_bilinear_matrix_values_save,E,m,norm) - print '(E9.1,2X,I8,2X,F10.2,2X,F10.8,2X,F12.6)', norm_sort(n), m, & + print '(E9.1,2X,I8,2X,F10.2,2X,F10.6,2X,F12.6)', norm_sort(n), m, & dble( elec_alpha_num**3 + elec_alpha_num**2 * m ) / & dble( elec_alpha_num**3 + elec_alpha_num**2 * n ), norm, E if (E < target_energy) then From 62a8253244a9cc3c083c195bfc7c61b938180dca Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 16 Nov 2015 22:08:36 +0100 Subject: [PATCH 4/6] Introduced Schwartz screening in map_get --- plugins/Hartree_Fock/Fock_matrix.irp.f | 41 +++++++----- plugins/Hartree_Fock/SCF.irp.f | 2 +- src/Determinants/slater_rules.irp.f | 72 +++++++++++----------- src/Integrals_Bielec/ao_bi_integrals.irp.f | 68 ++++++++++++-------- src/Integrals_Bielec/map_integrals.irp.f | 25 +++++++- src/Integrals_Bielec/mo_bi_integrals.irp.f | 16 +++++ 6 files changed, 147 insertions(+), 77 deletions(-) diff --git a/plugins/Hartree_Fock/Fock_matrix.irp.f b/plugins/Hartree_Fock/Fock_matrix.irp.f index 60413ab9..82540de3 100644 --- a/plugins/Hartree_Fock/Fock_matrix.irp.f +++ b/plugins/Hartree_Fock/Fock_matrix.irp.f @@ -113,9 +113,10 @@ END_PROVIDER END_DOC integer :: i,j,k,l,k1,r,s + integer :: i0,j0,k0,l0 integer*8 :: p,q - double precision :: integral - double precision :: ao_bielec_integral + double precision :: integral, c0, c1, c2 + double precision :: ao_bielec_integral, local_threshold double precision, allocatable :: ao_bi_elec_integral_alpha_tmp(:,:) double precision, allocatable :: ao_bi_elec_integral_beta_tmp(:,:) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: ao_bi_elec_integral_beta_tmp @@ -126,11 +127,12 @@ END_PROVIDER if (do_direct_integrals) then !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,p,q,r,s, & - !$OMP ao_bi_elec_integral_alpha_tmp,ao_bi_elec_integral_beta_tmp)& + !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,p,q,r,s,i0,j0,k0,l0, & + !$OMP ao_bi_elec_integral_alpha_tmp,ao_bi_elec_integral_beta_tmp, c0, c1, c2, & + !$OMP local_threshold)& !$OMP SHARED(ao_num,ao_num_align,HF_density_matrix_ao_alpha,HF_density_matrix_ao_beta,& !$OMP ao_integrals_map,ao_integrals_threshold, ao_bielec_integral_schwartz, & - !$OMP ao_overlap_abs, ao_bi_elec_integral_alpha, ao_bi_elec_integral_beta) + !$OMP ao_overlap_abs, ao_bi_elec_integral_alpha, ao_bi_elec_integral_beta) allocate(keys(1), values(1)) allocate(ao_bi_elec_integral_alpha_tmp(ao_num_align,ao_num), & @@ -157,14 +159,16 @@ END_PROVIDER < ao_integrals_threshold) then cycle endif - if (ao_bielec_integral_schwartz(k,l)*ao_bielec_integral_schwartz(i,j) & - < ao_integrals_threshold) then - cycle - endif - values(1) = ao_bielec_integral(k,l,i,j) - if (abs(values(1)) < ao_integrals_threshold) then + local_threshold = ao_bielec_integral_schwartz(k,l)*ao_bielec_integral_schwartz(i,j) + if (local_threshold < ao_integrals_threshold) then cycle endif + i0 = i + j0 = j + k0 = k + l0 = l + values(1) = 0.d0 + local_threshold = ao_integrals_threshold/local_threshold do k2=1,8 if (kk(k2)==0) then cycle @@ -173,12 +177,21 @@ END_PROVIDER j = jj(k2) k = kk(k2) l = ll(k2) - integral = (HF_density_matrix_ao_alpha(k,l)+HF_density_matrix_ao_beta(k,l)) * values(1) + c0 = HF_density_matrix_ao_alpha(k,l)+HF_density_matrix_ao_beta(k,l) + c1 = HF_density_matrix_ao_alpha(k,i) + c2 = HF_density_matrix_ao_beta(k,i) + if ( dabs(c0)+dabs(c1)+dabs(c2) < local_threshold) then + cycle + endif + if (values(1) == 0.d0) then + values(1) = ao_bielec_integral(k0,l0,i0,j0) + endif + integral = c0 * values(1) ao_bi_elec_integral_alpha_tmp(i,j) += integral ao_bi_elec_integral_beta_tmp (i,j) += integral integral = values(1) - ao_bi_elec_integral_alpha_tmp(l,j) -= HF_density_matrix_ao_alpha(k,i) * integral - ao_bi_elec_integral_beta_tmp (l,j) -= HF_density_matrix_ao_beta (k,i) * integral + ao_bi_elec_integral_alpha_tmp(l,j) -= c1 * integral + ao_bi_elec_integral_beta_tmp (l,j) -= c2 * integral enddo enddo !$OMP END DO NOWAIT diff --git a/plugins/Hartree_Fock/SCF.irp.f b/plugins/Hartree_Fock/SCF.irp.f index f21c7399..dead61ee 100644 --- a/plugins/Hartree_Fock/SCF.irp.f +++ b/plugins/Hartree_Fock/SCF.irp.f @@ -42,7 +42,7 @@ subroutine run BEGIN_DOC ! Run SCF calculation END_DOC - double precision :: SCF_energy_before,SCF_energy_after,diag_H_mat_elem,get_mo_bielec_integral + double precision :: SCF_energy_before,SCF_energy_after,diag_H_mat_elem double precision :: E0 integer :: i_it, i, j, k diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index 4cb1afff..7a122793 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -365,7 +365,7 @@ subroutine i_H_j(key_i,key_j,Nint,hij) integer :: exc(0:2,2,2) integer :: degree - double precision :: get_mo_bielec_integral + double precision :: get_mo_bielec_integral_schwartz integer :: m,n,p,q integer :: i,j,k integer :: occ(Nint*bit_kind_size,2) @@ -390,31 +390,31 @@ subroutine i_H_j(key_i,key_j,Nint,hij) call get_double_excitation(key_i,key_j,exc,phase,Nint) if (exc(0,1,1) == 1) then ! Mono alpha, mono beta - hij = phase*get_mo_bielec_integral( & + hij = phase*get_mo_bielec_integral_schwartz( & exc(1,1,1), & exc(1,1,2), & exc(1,2,1), & exc(1,2,2) ,mo_integrals_map) else if (exc(0,1,1) == 2) then ! Double alpha - hij = phase*(get_mo_bielec_integral( & + hij = phase*(get_mo_bielec_integral_schwartz( & exc(1,1,1), & exc(2,1,1), & exc(1,2,1), & exc(2,2,1) ,mo_integrals_map) - & - get_mo_bielec_integral( & + get_mo_bielec_integral_schwartz( & exc(1,1,1), & exc(2,1,1), & exc(2,2,1), & exc(1,2,1) ,mo_integrals_map) ) else if (exc(0,1,2) == 2) then ! Double beta - hij = phase*(get_mo_bielec_integral( & + hij = phase*(get_mo_bielec_integral_schwartz( & exc(1,1,2), & exc(2,1,2), & exc(1,2,2), & exc(2,2,2) ,mo_integrals_map) - & - get_mo_bielec_integral( & + get_mo_bielec_integral_schwartz( & exc(1,1,2), & exc(2,1,2), & exc(2,2,2), & @@ -432,15 +432,15 @@ subroutine i_H_j(key_i,key_j,Nint,hij) do k = 1, elec_alpha_num i = occ(k,1) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) has_mipi(i) = .True. endif enddo do k = 1, elec_beta_num i = occ(k,2) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) has_mipi(i) = .True. endif enddo @@ -459,15 +459,15 @@ subroutine i_H_j(key_i,key_j,Nint,hij) do k = 1, elec_beta_num i = occ(k,2) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) has_mipi(i) = .True. endif enddo do k = 1, elec_alpha_num i = occ(k,1) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) has_mipi(i) = .True. endif enddo @@ -501,7 +501,7 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree) integer,intent(out) :: exc(0:2,2,2) integer,intent(out) :: degree - double precision :: get_mo_bielec_integral + double precision :: get_mo_bielec_integral_schwartz integer :: m,n,p,q integer :: i,j,k integer :: occ(Nint*bit_kind_size,2) @@ -526,31 +526,31 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree) call get_double_excitation(key_i,key_j,exc,phase,Nint) if (exc(0,1,1) == 1) then ! Mono alpha, mono beta - hij = phase*get_mo_bielec_integral( & + hij = phase*get_mo_bielec_integral_schwartz( & exc(1,1,1), & exc(1,1,2), & exc(1,2,1), & exc(1,2,2) ,mo_integrals_map) else if (exc(0,1,1) == 2) then ! Double alpha - hij = phase*(get_mo_bielec_integral( & + hij = phase*(get_mo_bielec_integral_schwartz( & exc(1,1,1), & exc(2,1,1), & exc(1,2,1), & exc(2,2,1) ,mo_integrals_map) - & - get_mo_bielec_integral( & + get_mo_bielec_integral_schwartz( & exc(1,1,1), & exc(2,1,1), & exc(2,2,1), & exc(1,2,1) ,mo_integrals_map) ) else if (exc(0,1,2) == 2) then ! Double beta - hij = phase*(get_mo_bielec_integral( & + hij = phase*(get_mo_bielec_integral_schwartz( & exc(1,1,2), & exc(2,1,2), & exc(1,2,2), & exc(2,2,2) ,mo_integrals_map) - & - get_mo_bielec_integral( & + get_mo_bielec_integral_schwartz( & exc(1,1,2), & exc(2,1,2), & exc(2,2,2), & @@ -568,15 +568,15 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree) do k = 1, elec_alpha_num i = occ(k,1) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) has_mipi(i) = .True. endif enddo do k = 1, elec_beta_num i = occ(k,2) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) has_mipi(i) = .True. endif enddo @@ -595,15 +595,15 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree) do k = 1, elec_beta_num i = occ(k,2) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) has_mipi(i) = .True. endif enddo do k = 1, elec_alpha_num i = occ(k,1) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) has_mipi(i) = .True. endif enddo @@ -637,7 +637,7 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble) integer :: exc(0:2,2,2) integer :: degree - double precision :: get_mo_bielec_integral + double precision :: get_mo_bielec_integral_schwartz integer :: m,n,p,q integer :: i,j,k integer :: occ(Nint*bit_kind_size,2) @@ -664,31 +664,31 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble) call get_double_excitation(key_i,key_j,exc,phase,Nint) if (exc(0,1,1) == 1) then ! Mono alpha, mono beta - hij = phase*get_mo_bielec_integral( & + hij = phase*get_mo_bielec_integral_schwartz( & exc(1,1,1), & exc(1,1,2), & exc(1,2,1), & exc(1,2,2) ,mo_integrals_map) else if (exc(0,1,1) == 2) then ! Double alpha - hij = phase*(get_mo_bielec_integral( & + hij = phase*(get_mo_bielec_integral_schwartz( & exc(1,1,1), & exc(2,1,1), & exc(1,2,1), & exc(2,2,1) ,mo_integrals_map) - & - get_mo_bielec_integral( & + get_mo_bielec_integral_schwartz( & exc(1,1,1), & exc(2,1,1), & exc(2,2,1), & exc(1,2,1) ,mo_integrals_map) ) else if (exc(0,1,2) == 2) then ! Double beta - hij = phase*(get_mo_bielec_integral( & + hij = phase*(get_mo_bielec_integral_schwartz( & exc(1,1,2), & exc(2,1,2), & exc(1,2,2), & exc(2,2,2) ,mo_integrals_map) - & - get_mo_bielec_integral( & + get_mo_bielec_integral_schwartz( & exc(1,1,2), & exc(2,1,2), & exc(2,2,2), & @@ -706,15 +706,15 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble) do k = 1, elec_alpha_num i = occ(k,1) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) has_mipi(i) = .True. endif enddo do k = 1, elec_beta_num i = occ(k,2) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) has_mipi(i) = .True. endif enddo @@ -733,15 +733,15 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble) do k = 1, elec_beta_num i = occ(k,2) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) has_mipi(i) = .True. endif enddo do k = 1, elec_alpha_num i = occ(k,1) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) has_mipi(i) = .True. endif enddo diff --git a/src/Integrals_Bielec/ao_bi_integrals.irp.f b/src/Integrals_Bielec/ao_bi_integrals.irp.f index 4d1b2f95..19f4d2a7 100644 --- a/src/Integrals_Bielec/ao_bi_integrals.irp.f +++ b/src/Integrals_Bielec/ao_bi_integrals.irp.f @@ -204,7 +204,7 @@ double precision function ao_bielec_integral_schwartz_accel(i,j,k,l) 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_bielec_integral_schwartz_accel += coef4 * integral + ao_bielec_integral_schwartz_accel = ao_bielec_integral_schwartz_accel + coef4 * integral enddo ! s enddo ! r enddo ! q @@ -264,7 +264,7 @@ double precision function ao_bielec_integral_schwartz_accel(i,j,k,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_bielec_integral_schwartz_accel += coef4 * integral + ao_bielec_integral_schwartz_accel = ao_bielec_integral_schwartz_accel + coef4 * integral enddo ! s enddo ! r enddo ! q @@ -307,12 +307,20 @@ subroutine compute_ao_bielec_integrals(j,k,l,sze,buffer_value) buffer_value = 0._integral_kind return endif + if (ao_bielec_integral_schwartz(j,l) < thresh ) then + buffer_value = 0._integral_kind + return + endif do i = 1, ao_num if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then buffer_value(i) = 0._integral_kind cycle endif + if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < thresh ) then + buffer_value(i) = 0._integral_kind + cycle + endif !DIR$ FORCEINLINE buffer_value(i) = ao_bielec_integral(i,k,j,l) enddo @@ -669,32 +677,44 @@ double precision function ERI(alpha,beta,delta,gama,a_x,b_x,c_x,d_x,a_y,b_y,c_y, integer :: n_pt_sup double precision :: p,q,denom,coeff double precision :: I_f + integer :: nx,ny,nz include 'Utils/constants.include.F' - if(iand(a_x+b_x+c_x+d_x,1).eq.1.or.iand(a_y+b_y+c_y+d_y,1).eq.1.or.iand(a_z+b_z+c_z+d_z,1).eq.1)then + nx = a_x+b_x+c_x+d_x + if(iand(nx,1) == 1) then ERI = 0.d0 return - else - - 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) - coeff = pi_5_2 / (p * q * dsqrt(p+q)) - !DIR$ FORCEINLINE - n_pt = 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) - - 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 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 = ishft( 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 diff --git a/src/Integrals_Bielec/map_integrals.irp.f b/src/Integrals_Bielec/map_integrals.irp.f index 83950c5c..2a0618c9 100644 --- a/src/Integrals_Bielec/map_integrals.irp.f +++ b/src/Integrals_Bielec/map_integrals.irp.f @@ -295,15 +295,36 @@ double precision function get_mo_bielec_integral(i,j,k,l,map) get_mo_bielec_integral = dble(tmp) end +double precision function get_mo_bielec_integral_schwartz(i,j,k,l,map) + use map_module + implicit none + BEGIN_DOC + ! Returns one integral in the MO basis + END_DOC + integer, intent(in) :: i,j,k,l + integer(key_kind) :: idx + type(map_type), intent(inout) :: map + real(integral_kind) :: tmp + PROVIDE mo_bielec_integrals_in_map + !DIR$ FORCEINLINE + if (mo_bielec_integral_schwartz(i,k)*mo_bielec_integral_schwartz(j,l) < mo_integrals_threshold) then + tmp = 0.d0 + else + call bielec_integrals_index(i,j,k,l,idx) + call map_get(map,idx,tmp) + endif + get_mo_bielec_integral_schwartz = dble(tmp) +end + double precision function mo_bielec_integral(i,j,k,l) implicit none BEGIN_DOC ! Returns one integral in the MO basis END_DOC integer, intent(in) :: i,j,k,l - double precision :: get_mo_bielec_integral + double precision :: get_mo_bielec_integral_schwartz PROVIDE mo_bielec_integrals_in_map - mo_bielec_integral = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + mo_bielec_integral = get_mo_bielec_integral_schwartz(i,j,k,l,mo_integrals_map) return end diff --git a/src/Integrals_Bielec/mo_bi_integrals.irp.f b/src/Integrals_Bielec/mo_bi_integrals.irp.f index 7d772ac3..98737e2a 100644 --- a/src/Integrals_Bielec/mo_bi_integrals.irp.f +++ b/src/Integrals_Bielec/mo_bi_integrals.irp.f @@ -488,3 +488,19 @@ END_PROVIDER enddo END_PROVIDER + +BEGIN_PROVIDER [ double precision, mo_bielec_integral_schwartz,(mo_tot_num,mo_tot_num) ] + implicit none + BEGIN_DOC + ! Needed to compute Schwartz inequalities + END_DOC + + integer :: i,k + + do i=1,mo_tot_num + do k=1,mo_tot_num + mo_bielec_integral_schwartz(k,i) = dsqrt(mo_bielec_integral_jj(k,i)) + enddo + enddo + +END_PROVIDER From b3d0c0620933f8f45490b74131d94065f7a09c79 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 17 Nov 2015 00:16:00 +0100 Subject: [PATCH 5/6] Added integrals access benchmark --- scripts/compilation/qp_create_ninja.py | 5 +- src/Integrals_Bielec/map_integrals.irp.f | 10 +- src/Integrals_Bielec/test_integrals.irp.f | 136 ++++++++++++++++++++++ 3 files changed, 145 insertions(+), 6 deletions(-) create mode 100644 src/Integrals_Bielec/test_integrals.irp.f diff --git a/scripts/compilation/qp_create_ninja.py b/scripts/compilation/qp_create_ninja.py index 094ec7a6..b0fe2f8c 100755 --- a/scripts/compilation/qp_create_ninja.py +++ b/scripts/compilation/qp_create_ninja.py @@ -35,7 +35,8 @@ except ImportError: from qp_path import QP_ROOT, QP_SRC, QP_EZFIO -EZFIO_LIB = join(QP_ROOT, "lib", "libezfio.a") +LIB = "" # join(QP_ROOT, "lib", "rdtsc.o") +EZFIO_LIB = join(QP_ROOT, "lib", "libezfio.a") ROOT_BUILD_NINJA = join(QP_ROOT, "config", "build.ninja") header = r"""# @@ -94,7 +95,7 @@ def ninja_create_env_variable(pwd_config_file): l_string.append(str_) lib_lapack = get_compilation_option(pwd_config_file, "LAPACK_LIB") - l_string.append("{0} = {1} {2}".format("LIB", lib_lapack, EZFIO_LIB)) + l_string.append("LIB = {0} {1} {2}".format(LIB, lib_lapack, EZFIO_LIB)) l_string.append("") diff --git a/src/Integrals_Bielec/map_integrals.irp.f b/src/Integrals_Bielec/map_integrals.irp.f index 2a0618c9..181075bb 100644 --- a/src/Integrals_Bielec/map_integrals.irp.f +++ b/src/Integrals_Bielec/map_integrals.irp.f @@ -291,6 +291,7 @@ double precision function get_mo_bielec_integral(i,j,k,l,map) PROVIDE mo_bielec_integrals_in_map !DIR$ FORCEINLINE call bielec_integrals_index(i,j,k,l,idx) + !DIR$ FORCEINLINE call map_get(map,idx,tmp) get_mo_bielec_integral = dble(tmp) end @@ -306,12 +307,13 @@ double precision function get_mo_bielec_integral_schwartz(i,j,k,l,map) type(map_type), intent(inout) :: map real(integral_kind) :: tmp PROVIDE mo_bielec_integrals_in_map - !DIR$ FORCEINLINE - if (mo_bielec_integral_schwartz(i,k)*mo_bielec_integral_schwartz(j,l) < mo_integrals_threshold) then - tmp = 0.d0 - else + if (mo_bielec_integral_schwartz(i,k)*mo_bielec_integral_schwartz(j,l) > mo_integrals_threshold) then + !DIR$ FORCEINLINE call bielec_integrals_index(i,j,k,l,idx) + !DIR$ FORCEINLINE call map_get(map,idx,tmp) + else + tmp = 0.d0 endif get_mo_bielec_integral_schwartz = dble(tmp) end diff --git a/src/Integrals_Bielec/test_integrals.irp.f b/src/Integrals_Bielec/test_integrals.irp.f new file mode 100644 index 00000000..8d590f27 --- /dev/null +++ b/src/Integrals_Bielec/test_integrals.irp.f @@ -0,0 +1,136 @@ +program bench_maps + implicit none + BEGIN_DOC +! Performs timing benchmarks on integral access + END_DOC + integer :: i,j,k,l + integer*8 :: ii,jj + double precision :: r, cpu + integer*8 :: cpu0, cpu1, count_rate, count_max + + PROVIDE mo_bielec_integrals_in_map + print *, mo_tot_num, 'MOs' + + ! Time random function + call system_clock(cpu0, count_rate, count_max) + do ii=1,100000000_8 + call random_number(r) + i = int(r*mo_tot_num)+1 + call random_number(r) + j = int(r*mo_tot_num)+1 + call random_number(r) + k = int(r*mo_tot_num)+1 + call random_number(r) + l = int(r*mo_tot_num)+1 + enddo + call system_clock(cpu1, count_rate, count_max) + cpu = (cpu1-cpu0)/count_rate + print *, 'Random function : ', cpu/dble(ii) + + call system_clock(cpu0, count_rate, count_max) + do ii=1,100000000_8 + call random_number(r) + i = int(r*mo_tot_num)+1 + call random_number(r) + j = int(r*mo_tot_num)+1 + call random_number(r) + k = int(r*mo_tot_num)+1 + call random_number(r) + l = int(r*mo_tot_num)+1 + call get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + enddo + call system_clock(cpu1, count_rate, count_max) + cpu = -cpu + (cpu1 - cpu0)/count_rate + print *, 'Random access : ', cpu/dble(ii) + + ii=0_8 + call system_clock(cpu0, count_rate, count_max) + do jj=1,10 + do l=1,mo_tot_num + do k=1,mo_tot_num + do j=1,mo_tot_num + do i=1,mo_tot_num + ii += 1 + call get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + enddo + enddo + enddo + enddo + enddo + call system_clock(cpu1, count_rate, count_max) + cpu = (cpu1 - cpu0)/count_rate + print *, 'loop ijkl : ', cpu/dble(ii) + + ii=0 + call system_clock(cpu0, count_rate, count_max) + do jj=1,10 + do l=1,mo_tot_num + do j=1,mo_tot_num + do k=1,mo_tot_num + do i=1,mo_tot_num + ii += 1 + call get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + enddo + enddo + enddo + enddo + enddo + call system_clock(cpu1, count_rate, count_max) + cpu = (cpu1 - cpu0)/count_rate + print *, 'loop ikjl : ', cpu/dble(ii) + + ii=0 + call system_clock(cpu0, count_rate, count_max) + do jj=1,10 + do k=1,mo_tot_num + do l=1,mo_tot_num + do j=1,mo_tot_num + do i=1,mo_tot_num + ii += 1 + call get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + enddo + enddo + enddo + enddo + enddo + call system_clock(cpu1, count_rate, count_max) + cpu = (cpu1 - cpu0)/count_rate + print *, 'loop ijlk : ', cpu/dble(ii) + + ii=0 + call system_clock(cpu0, count_rate, count_max) + do jj=1,10 + do i=1,mo_tot_num + do j=1,mo_tot_num + do k=1,mo_tot_num + do l=1,mo_tot_num + ii += 1 + call get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + enddo + enddo + enddo + enddo + enddo + call system_clock(cpu1, count_rate, count_max) + cpu = (cpu1 - cpu0)/count_rate + print *, 'loop lkji : ', cpu/dble(ii) + + ii=0 + call system_clock(cpu0, count_rate, count_max) + do jj=1,10 + do j=1,mo_tot_num + do i=1,mo_tot_num + do k=1,mo_tot_num + do l=1,mo_tot_num + ii += 1 + call get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + enddo + enddo + enddo + enddo + enddo + call system_clock(cpu1, count_rate, count_max) + cpu = (cpu1 - cpu0)/count_rate + print *, 'loop lkij : ', cpu/dble(ii) + +end From 4679023df94fcf070cd6e0beaf23d9e414ac8553 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 17 Nov 2015 00:16:15 +0100 Subject: [PATCH 6/6] Added Davidson threshold in EZFIO --- src/Determinants/EZFIO.cfg | 6 ++++++ src/Determinants/davidson.irp.f | 16 +++++++--------- src/Determinants/s2.irp.f | 12 +++++------- src/Determinants/slater_rules.irp.f | 8 ++++---- 4 files changed, 22 insertions(+), 20 deletions(-) diff --git a/src/Determinants/EZFIO.cfg b/src/Determinants/EZFIO.cfg index 9fe2f23c..9613c6c1 100644 --- a/src/Determinants/EZFIO.cfg +++ b/src/Determinants/EZFIO.cfg @@ -40,6 +40,12 @@ doc: Force the wave function to be an eigenfunction of S^2 interface: ezfio,provider,ocaml default: False +[threshold_davidson] +type: Threshold +doc: Thresholds of Davidson's algorithm +interface: ezfio,provider,ocaml +default: 1.e-8 + [threshold_generators] type: Threshold doc: Thresholds on generators (fraction of the norm) diff --git a/src/Determinants/davidson.irp.f b/src/Determinants/davidson.irp.f index fd2efcf9..5de5d379 100644 --- a/src/Determinants/davidson.irp.f +++ b/src/Determinants/davidson.irp.f @@ -591,14 +591,12 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun abort_here = abort_all end - BEGIN_PROVIDER [ character(64), davidson_criterion ] -&BEGIN_PROVIDER [ double precision, davidson_threshold ] +BEGIN_PROVIDER [ character(64), davidson_criterion ] implicit none BEGIN_DOC ! Can be : [ energy | residual | both | wall_time | cpu_time | iterations ] END_DOC davidson_criterion = 'residual' - davidson_threshold = 1.d-10 END_PROVIDER subroutine davidson_converged(energy,residual,wall,iterations,cpu,N_st,converged) @@ -621,20 +619,20 @@ subroutine davidson_converged(energy,residual,wall,iterations,cpu,N_st,converged E = energy - energy_old energy_old = energy if (davidson_criterion == 'energy') then - converged = dabs(maxval(E(1:N_st))) < davidson_threshold + converged = dabs(maxval(E(1:N_st))) < threshold_davidson else if (davidson_criterion == 'residual') then - converged = dabs(maxval(residual(1:N_st))) < davidson_threshold + converged = dabs(maxval(residual(1:N_st))) < threshold_davidson else if (davidson_criterion == 'both') then converged = dabs(maxval(residual(1:N_st))) + dabs(maxval(E(1:N_st)) ) & - < davidson_threshold + < threshold_davidson else if (davidson_criterion == 'wall_time') then call wall_time(time) - converged = time - wall > davidson_threshold + converged = time - wall > threshold_davidson else if (davidson_criterion == 'cpu_time') then call cpu_time(time) - converged = time - cpu > davidson_threshold + converged = time - cpu > threshold_davidson else if (davidson_criterion == 'iterations') then - converged = iterations >= int(davidson_threshold) + converged = iterations >= int(threshold_davidson) endif converged = converged.or.abort_here end diff --git a/src/Determinants/s2.irp.f b/src/Determinants/s2.irp.f index 573f094f..577e50f2 100644 --- a/src/Determinants/s2.irp.f +++ b/src/Determinants/s2.irp.f @@ -122,15 +122,13 @@ subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2) integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi, pass double precision :: davidson_threshold_bis - !PROVIDE davidson_threshold - s2 = 0.d0 - davidson_threshold_bis = davidson_threshold + davidson_threshold_bis = threshold_davidson call sort_dets_ab_v(psi_keys_tmp, sorted, sort_idx, shortcut, version, n, N_int) !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(i,j,s2_tmp,sh, sh2, ni, exa, ext, org_i, org_j, endi, pass)& - !$OMP SHARED(n,psi_coefs_tmp,psi_keys_tmp,N_int,davidson_threshold,shortcut,sorted,sort_idx,version)& + !$OMP SHARED(n,psi_coefs_tmp,psi_keys_tmp,N_int,threshold_davidson,shortcut,sorted,sort_idx,version)& !$OMP REDUCTION(+:s2) !$OMP DO SCHEDULE(dynamic) @@ -162,7 +160,7 @@ subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2) org_j = sort_idx(j) if ( dabs(psi_coefs_tmp(org_j)) + dabs(psi_coefs_tmp(org_i))& - > davidson_threshold ) then + > threshold_davidson ) then call get_s2(psi_keys_tmp(1,1,org_i),psi_keys_tmp(1,1,org_j),s2_tmp,N_int) s2 = s2 + psi_coefs_tmp(org_i)*psi_coefs_tmp(org_j)*s2_tmp endif @@ -179,7 +177,7 @@ subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2) !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(i,j,s2_tmp,sh, sh2, ni, exa, ext, org_i, org_j, endi, pass)& - !$OMP SHARED(n,psi_coefs_tmp,psi_keys_tmp,N_int,davidson_threshold,shortcut,sorted,sort_idx,version)& + !$OMP SHARED(n,psi_coefs_tmp,psi_keys_tmp,N_int,threshold_davidson,shortcut,sorted,sort_idx,version)& !$OMP REDUCTION(+:s2) !$OMP DO SCHEDULE(dynamic) @@ -195,7 +193,7 @@ subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2) org_j = sort_idx(j) if ( dabs(psi_coefs_tmp(org_j)) + dabs(psi_coefs_tmp(org_i))& - > davidson_threshold ) then + > threshold_davidson ) then call get_s2(psi_keys_tmp(1,1,org_i),psi_keys_tmp(1,1,org_j),s2_tmp,N_int) s2 = s2 + psi_coefs_tmp(org_i)*psi_coefs_tmp(org_j)*s2_tmp endif diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index 7a122793..4eaf2d31 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -1256,7 +1256,7 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint) !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(i,hij,j,k,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,local_threshold,sorted_i)& - !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,davidson_threshold,sorted,shortcut,sort_idx,version) + !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,threshold_davidson,sorted,shortcut,sort_idx,version) allocate(vt(n)) Vt = 0.d0 @@ -1273,7 +1273,7 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint) do i=shortcut(sh),shortcut(sh+1)-1 org_i = sort_idx(i) - local_threshold = davidson_threshold - dabs(u_0(org_i)) + local_threshold = threshold_davidson - dabs(u_0(org_i)) if(sh==sh2) then endi = i-1 else @@ -1315,14 +1315,14 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint) !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(i,hij,j,k,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,local_threshold)& - !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,davidson_threshold,sorted,shortcut,sort_idx,version) + !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,threshold_davidson,sorted,shortcut,sort_idx,version) allocate(vt(n)) Vt = 0.d0 !$OMP DO SCHEDULE(dynamic) do sh=1,shortcut(0) do i=shortcut(sh),shortcut(sh+1)-1 - local_threshold = davidson_threshold - dabs(u_0(org_i)) + local_threshold = threshold_davidson - dabs(u_0(org_i)) org_i = sort_idx(i) do j=shortcut(sh),i-1 org_j = sort_idx(j)