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

Merge pull request #114 from scemama/master

Browse Source 
Integrals
This commit is contained in:
Anthony Scemama 2015-11-17 00:30:01 +01:00
commit d9e73e500a
13 changed files with 332 additions and 104 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
plugins/Hartree_Fock/EZFIO.cfg
View File

@ -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 to improve SCF convergence
interface: ezfio,provider,ocaml
default: 0.0
[mo_guess_type]
type: MO_guess
doc: Initial MO guess. Can be [ Huckel | HCore ]

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

@ -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
@ -108,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
@ -121,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), &
@ -152,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
@ -168,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
@ -315,7 +333,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

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

2
plugins/QmcChem/target_pt2_qmc.irp.f
View File

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

5
scripts/compilation/qp_create_ninja.py
View File

@ -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("")

6
src/Determinants/EZFIO.cfg
View File

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

16
src/Determinants/davidson.irp.f
View File

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

12
src/Determinants/s2.irp.f
View File

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

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

76
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
@ -378,8 +386,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 +427,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
@ -665,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

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

@ -291,19 +291,42 @@ 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
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
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
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

136
src/Integrals_Bielec/test_integrals.irp.f Normal file
View File

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