LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-26 15:12:14 +02:00
Code Issues Releases Wiki Activity

Introduced Schwartz screening in map_get

Browse Source
This commit is contained in:
Anthony Scemama 2015-11-16 22:08:36 +01:00
parent 1d2631f7df
commit 62a8253244
6 changed files with 147 additions and 77 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

41
plugins/Hartree_Fock/Fock_matrix.irp.f
View File

@ -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

2
plugins/Hartree_Fock/SCF.irp.f
View File

@ -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

72
src/Determinants/slater_rules.irp.f
View File

@ -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

68
src/Integrals_Bielec/ao_bi_integrals.irp.f
View File

@ -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

25
src/Integrals_Bielec/map_integrals.irp.f
View File

@ -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 <ij|kl> 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 <ij|kl> 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

16
src/Integrals_Bielec/mo_bi_integrals.irp.f
View File

@ -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