From 8842c1ff76f35d73b0a3fad789d6b1d5e3714bdf Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 22 Sep 2014 21:02:10 +0200 Subject: [PATCH] Direct integral driven HF --- src/BiInts/ao_bi_integrals.irp.f | 4 +- src/BiInts/map_integrals.irp.f | 1 + src/Hartree_Fock/Fock_matrix.irp.f | 88 ++++++++++++++++++++---------- 3 files changed, 61 insertions(+), 32 deletions(-) diff --git a/src/BiInts/ao_bi_integrals.irp.f b/src/BiInts/ao_bi_integrals.irp.f index 272f5ad8..3093243b 100644 --- a/src/BiInts/ao_bi_integrals.irp.f +++ b/src/BiInts/ao_bi_integrals.irp.f @@ -215,7 +215,7 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ] integer :: jl_pairs(2,ao_num*(ao_num+1)/2), kk, m, j1, i1, lmax PROVIDE gauleg_t2 ao_integrals_map all_utils - PROVIDE ao_bielec_integral_schwartz +! PROVIDE ao_bielec_integral_schwartz integral = ao_bielec_integral(1,1,1,1) real :: map_mb @@ -357,7 +357,7 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, ao_bielec_integral_schwartz,(ao_num,ao_num) ] implicit none BEGIN_DOC - ! Needed to compuet Schwartz inequalities + ! Needed to compute Schwartz inequalities END_DOC integer :: i,k diff --git a/src/BiInts/map_integrals.irp.f b/src/BiInts/map_integrals.irp.f index d82c88df..dc650294 100644 --- a/src/BiInts/map_integrals.irp.f +++ b/src/BiInts/map_integrals.irp.f @@ -9,6 +9,7 @@ BEGIN_PROVIDER [ type(map_type), ao_integrals_map ] ! AO integrals END_DOC integer*8 :: sze + PROVIDE gauleg_t2 call bielec_integrals_index(ao_num,ao_num,ao_num,ao_num,sze) call map_init(ao_integrals_map,sze) write(output_BiInts,*) 'AO map initialized' diff --git a/src/Hartree_Fock/Fock_matrix.irp.f b/src/Hartree_Fock/Fock_matrix.irp.f index 6802d0fb..a1553ebd 100644 --- a/src/Hartree_Fock/Fock_matrix.irp.f +++ b/src/Hartree_Fock/Fock_matrix.irp.f @@ -107,44 +107,73 @@ END_PROVIDER ! Alpha Fock matrix in AO basis set END_DOC - integer :: i,j,k,l,k1 + integer :: i,j,k,l,k1,r,s + integer*8 :: p,q double precision :: integral double precision :: ao_bielec_integral if (do_direct_integrals) then - PROVIDE all_utils ao_overlap_abs ao_integrals_threshold + PROVIDE all_utils ao_overlap_abs ao_integrals_threshold gauleg_t2 PROVIDE HF_density_matrix_ao_alpha HF_density_matrix_ao_beta PROVIDE ao_bi_elec_integral_alpha - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(i,j,l,k1,k,integral) & - !$OMP SHARED(ao_num,ao_bi_elec_integral_alpha,ao_mono_elec_integral,& - !$OMP ao_num_align,ao_bi_elec_integral_beta,HF_density_matrix_ao_alpha, & - !$OMP HF_density_matrix_ao_beta) - !$OMP DO SCHEDULE(GUIDED) - do j=1,ao_num - do i=1,j - ao_bi_elec_integral_alpha(i,j) = 0.d0 - ao_bi_elec_integral_beta (i,j) = 0.d0 - do l=1,ao_num - do k=1,ao_num - if ((abs(HF_density_matrix_ao_alpha(k,l)) > 1.d-9).or. & - (abs(HF_density_matrix_ao_beta (k,l)) > 1.d-9)) then - integral = (HF_density_matrix_ao_alpha(k,l)+HF_density_matrix_ao_beta (k,l)) * ao_bielec_integral(k,l,i,j) + + ao_bi_elec_integral_alpha = 0.d0 + ao_bi_elec_integral_beta = 0.d0 + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,p,q,r,s)& + !$OMP SHARED(ao_num,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) & + !$OMP REDUCTION(+:ao_bi_elec_integral_alpha,ao_bi_elec_integral_beta) + + allocate(keys(1), values(1)) + + q = ao_num*ao_num*ao_num*ao_num + !$OMP DO SCHEDULE(dynamic) + do p=1_8,q + call bielec_integrals_index_reverse(kk,ii,ll,jj,p) + if ( (kk(1)>ao_num).or. & + (ii(1)>ao_num).or. & + (jj(1)>ao_num).or. & + (ll(1)>ao_num) ) then + cycle + endif + k = kk(1) + i = ii(1) + l = ll(1) + j = jj(1) + + if (ao_overlap_abs(k,l)*ao_overlap_abs(i,j) & + < 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) +! values(1) = ao_bielec_integral(k,i,l,j) + if (abs(values(1)) < ao_integrals_threshold) then + cycle + endif + do k2=1,8 + if (kk(k2)==0) then + cycle + endif + i = ii(k2) + 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) ao_bi_elec_integral_alpha(i,j) += integral ao_bi_elec_integral_beta (i,j) += integral - - integral = ao_bielec_integral(k,j,i,l) - ao_bi_elec_integral_alpha(i,j) -= HF_density_matrix_ao_alpha(k,l)*integral - ao_bi_elec_integral_beta (i,j) -= HF_density_matrix_ao_beta (k,l)*integral - endif - enddo - enddo - ao_bi_elec_integral_alpha(j,i) = ao_bi_elec_integral_alpha(i,j) - ao_bi_elec_integral_beta (j,i) = ao_bi_elec_integral_beta (i,j) - enddo + integral = values(1) + ao_bi_elec_integral_alpha(l,j) -= HF_density_matrix_ao_alpha(k,i) * integral + ao_bi_elec_integral_beta (l,j) -= HF_density_matrix_ao_beta (k,i) * integral + enddo enddo - !$OMP END DO NOWAIT + !$OMP END DO + deallocate(keys,values) !$OMP END PARALLEL - else PROVIDE ao_bielec_integrals_in_map @@ -166,7 +195,6 @@ END_PROVIDER call get_cache_map_n_elements_max(ao_integrals_map,n_elements_max) allocate(keys(n_elements_max), values(n_elements_max)) - !$OMP BARRIER !$OMP DO SCHEDULE(dynamic) do i8=0_8,ao_integrals_map%map_size n_elements = n_elements_max