LCPQ
/
qp2
mirror of https://github.com/QuantumPackage/qp2.git
9
1
Fork 0
Code Releases Activity

Fixed schwartz screening when integrals are read

This commit is contained in:
Michel Caffarel 2020-05-12 22:46:39 +02:00
parent a62456b238
commit 329bcf805b
7 changed files with 137 additions and 153 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
src/ao_basis/aos.irp.f
View File

@ -3,7 +3,6 @@ BEGIN_PROVIDER [ integer, ao_prim_num_max ]
BEGIN_DOC
! Max number of primitives.
END_DOC
print *, 'XXXX', irp_here
ao_prim_num_max = maxval(ao_prim_num)
END_PROVIDER
@ -20,7 +19,6 @@ END_PROVIDER
C_A(1) = 0.d0
C_A(2) = 0.d0
C_A(3) = 0.d0
print *, 'XXXX', irp_here
ao_coef_normalized = 0.d0
do i=1,ao_num
@ -67,7 +65,6 @@ BEGIN_PROVIDER [ double precision, ao_coef_normalization_libint_factor, (ao_num)
integer :: l, powA(3), nz
integer :: i,j,k
nz=100
print *, 'XXXX', irp_here
C_A(1) = 0.d0
C_A(2) = 0.d0
C_A(3) = 0.d0
@ -102,7 +99,6 @@ END_PROVIDER
integer :: iorder(ao_prim_num_max)
double precision :: d(ao_prim_num_max,2)
integer :: i,j
print *, 'XXXX', irp_here
do i=1,ao_num
do j=1,ao_prim_num(i)
iorder(j) = j
@ -125,7 +121,6 @@ BEGIN_PROVIDER [ double precision, ao_coef_normalized_ordered_transp, (ao_prim_n
! Transposed :c:data:`ao_coef_normalized_ordered`
END_DOC
integer :: i,j
print *, 'XXXX', irp_here
do j=1, ao_num
do i=1, ao_prim_num_max
ao_coef_normalized_ordered_transp(i,j) = ao_coef_normalized_ordered(j,i)
@ -140,7 +135,6 @@ BEGIN_PROVIDER [ double precision, ao_expo_ordered_transp, (ao_prim_num_max,ao_n
! Transposed :c:data:`ao_expo_ordered`
END_DOC
integer :: i,j
print *, 'XXXX', irp_here
do j=1, ao_num
do i=1, ao_prim_num_max
ao_expo_ordered_transp(i,j) = ao_expo_ordered(j,i)
@ -157,7 +151,6 @@ END_PROVIDER
! :math:`l` value of the |AO|: :math`a+b+c` in :math:`x^a y^b z^c`
END_DOC
integer :: i
print *, 'XXXX', irp_here
do i=1,ao_num
ao_l(i) = ao_power(i,1) + ao_power(i,2) + ao_power(i,3)
ao_l_char(i) = l_to_character(ao_l(i))
@ -174,7 +167,6 @@ integer function ao_power_index(nx,ny,nz)
! :math:`\frac{1}{2} (l-n_x) (l-n_x+1) + n_z + 1`
END_DOC
integer :: l
print *, 'XXXX', irp_here
l = nx + ny + nz
ao_power_index = ((l-nx)*(l-nx+1))/2 + nz + 1
end
@ -185,7 +177,6 @@ BEGIN_PROVIDER [ character*(128), l_to_character, (0:7)]
! Character corresponding to the "l" value of an |AO|
END_DOC
implicit none
print *, 'XXXX', irp_here
l_to_character(0)='s'
l_to_character(1)='p'
l_to_character(2)='d'
@ -204,7 +195,6 @@ END_PROVIDER
! Number of |AOs| per atom
END_DOC
integer :: i
print *, 'XXXX', irp_here
Nucl_N_Aos = 0
do i = 1, ao_num
Nucl_N_Aos(ao_nucl(i)) +=1
@ -219,7 +209,6 @@ END_PROVIDER
END_DOC
integer :: i
integer, allocatable :: nucl_tmp(:)
print *, 'XXXX', irp_here
allocate(nucl_tmp(nucl_num))
nucl_tmp = 0
Nucl_Aos = 0
@ -240,7 +229,6 @@ END_PROVIDER
! By convention, for p,d,f and g |AOs|, we take the index
! of the |AO| with the the corresponding power in the x axis
END_DOC
print *, 'XXXX', irp_here
do i = 1, nucl_num
Nucl_num_shell_Aos(i) = 0
@ -288,7 +276,6 @@ BEGIN_PROVIDER [ character*(4), ao_l_char_space, (ao_num) ]
END_DOC
integer :: i
character*(4) :: give_ao_character_space
print *, 'XXXX', irp_here
do i=1,ao_num
if(ao_l(i)==0)then

1
src/ao_one_e_ints/ao_overlap.irp.f
View File

@ -109,7 +109,6 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num,ao_num) ]
double precision :: A_center(3), B_center(3)
integer :: power_A(3), power_B(3)
double precision :: lower_exp_val, dx
print *, "XXX---", irp_here
if (is_periodic) then
do j=1,ao_num
do i= 1,ao_num

10
src/ao_one_e_ints/pot_ao_ints.irp.f
View File

@ -3,6 +3,8 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)]
! Nucleus-electron interaction, in the |AO| basis set.
!
! :math:`\langle \chi_i | -\sum_A \frac{1}{|r-R_A|} | \chi_j \rangle`
!
! These integrals also contain the pseudopotential integrals.
END_DOC
implicit none
double precision :: alpha, beta, gama, delta
@ -75,11 +77,11 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)]
!$OMP END DO
!$OMP END PARALLEL
endif
IF (DO_PSEUDO) THEN
ao_integrals_n_e += ao_pseudo_integrals
ENDIF
IF (DO_PSEUDO) THEN
ao_integrals_n_e += ao_pseudo_integrals
ENDIF
endif
if (write_ao_integrals_n_e) then

5
src/ao_one_e_ints/screening.irp.f
View File

@ -3,14 +3,11 @@ logical function ao_one_e_integral_zero(i,k)
integer, intent(in) :: i,k
ao_one_e_integral_zero = .False.
if (.not.(read_ao_one_e_integrals.or.is_periodic)) then
if (.not.((io_ao_integrals_overlap/='None').or.is_periodic)) then
if (ao_overlap_abs(i,k) < ao_integrals_threshold) then
ao_one_e_integral_zero = .True.
return
endif
endif
if (ao_two_e_integral_schwartz(i,k) < ao_integrals_threshold) then
ao_one_e_integral_zero = .True.
endif
end

6
src/ao_two_e_ints/screening.irp.f
View File

@ -8,8 +8,8 @@ logical function ao_two_e_integral_zero(i,j,k,l)
ao_two_e_integral_zero = .True.
return
endif
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.
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

253
src/ao_two_e_ints/two_e_integrals.irp.f
View File

@ -18,89 +18,89 @@ double precision function ao_two_e_integral(i,j,k,l)
if (ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then
ao_two_e_integral = ao_two_e_integral_schwartz_accel(i,j,k,l)
return
endif
else
dim1 = n_pt_max_integrals
dim1 = n_pt_max_integrals
num_i = ao_nucl(i)
num_j = ao_nucl(j)
num_k = ao_nucl(k)
num_l = ao_nucl(l)
ao_two_e_integral = 0.d0
num_i = ao_nucl(i)
num_j = ao_nucl(j)
num_k = ao_nucl(k)
num_l = ao_nucl(l)
ao_two_e_integral = 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_power(i,p)
J_power(p) = ao_power(j,p)
K_power(p) = ao_power(k,p)
L_power(p) = ao_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
if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k)then
do p = 1, 3
I_power(p) = ao_power(i,p)
J_power(p) = ao_power(j,p)
K_power(p) = ao_power(k,p)
L_power(p) = ao_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
double precision :: coef1, coef2, coef3, coef4
double precision :: p_inv,q_inv
double precision :: general_primitive_integral
do p = 1, ao_prim_num(i)
coef1 = ao_coef_normalized_ordered_transp(p,i)
do q = 1, ao_prim_num(j)
coef2 = coef1*ao_coef_normalized_ordered_transp(q,j)
call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,&
ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j), &
I_power,J_power,I_center,J_center,dim1)
p_inv = 1.d0/pp
do r = 1, ao_prim_num(k)
coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
do s = 1, ao_prim_num(l)
coef4 = coef3*ao_coef_normalized_ordered_transp(s,l)
call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,&
ao_expo_ordered_transp(r,k),ao_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_two_e_integral = ao_two_e_integral + coef4 * integral
enddo ! s
enddo ! r
enddo ! q
enddo ! p
do p = 1, ao_prim_num(i)
coef1 = ao_coef_normalized_ordered_transp(p,i)
do q = 1, ao_prim_num(j)
coef2 = coef1*ao_coef_normalized_ordered_transp(q,j)
call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,&
ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j), &
I_power,J_power,I_center,J_center,dim1)
p_inv = 1.d0/pp
do r = 1, ao_prim_num(k)
coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
do s = 1, ao_prim_num(l)
coef4 = coef3*ao_coef_normalized_ordered_transp(s,l)
call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,&
ao_expo_ordered_transp(r,k),ao_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_two_e_integral = ao_two_e_integral + coef4 * integral
enddo ! s
enddo ! r
enddo ! q
enddo ! p
else
else
do p = 1, 3
I_power(p) = ao_power(i,p)
J_power(p) = ao_power(j,p)
K_power(p) = ao_power(k,p)
L_power(p) = ao_power(l,p)
enddo
double precision :: ERI
do p = 1, 3
I_power(p) = ao_power(i,p)
J_power(p) = ao_power(j,p)
K_power(p) = ao_power(k,p)
L_power(p) = ao_power(l,p)
enddo
double precision :: ERI
do p = 1, ao_prim_num(i)
coef1 = ao_coef_normalized_ordered_transp(p,i)
do q = 1, ao_prim_num(j)
coef2 = coef1*ao_coef_normalized_ordered_transp(q,j)
do r = 1, ao_prim_num(k)
coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
do s = 1, ao_prim_num(l)
coef4 = coef3*ao_coef_normalized_ordered_transp(s,l)
integral = ERI( &
ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_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_two_e_integral = ao_two_e_integral + coef4 * integral
enddo ! s
enddo ! r
enddo ! q
enddo ! p
do p = 1, ao_prim_num(i)
coef1 = ao_coef_normalized_ordered_transp(p,i)
do q = 1, ao_prim_num(j)
coef2 = coef1*ao_coef_normalized_ordered_transp(q,j)
do r = 1, ao_prim_num(k)
coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
do s = 1, ao_prim_num(l)
coef4 = coef3*ao_coef_normalized_ordered_transp(s,l)
integral = ERI( &
ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_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_two_e_integral = ao_two_e_integral + coef4 * integral
enddo ! s
enddo ! r
enddo ! q
enddo ! p
endif
endif
end
double precision function ao_two_e_integral_schwartz_accel(i,j,k,l)
@ -350,71 +350,72 @@ BEGIN_PROVIDER [ logical, ao_two_e_integrals_in_map ]
call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map)
print*, 'AO integrals provided'
ao_two_e_integrals_in_map = .True.
return
endif
else
print*, 'Providing the AO integrals'
call wall_time(wall_0)
call wall_time(wall_1)
call cpu_time(cpu_1)
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_two_e_integral(1,1,1,1)
endif
integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_socket_pull
call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'ao_integrals')
character(len=:), allocatable :: task
allocate(character(len=ao_num*12) :: task)
write(fmt,*) '(', ao_num, '(I5,X,I5,''|''))'
do l=1,ao_num
write(task,fmt) (i,l, i=1,l)
integer, external :: add_task_to_taskserver
if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task)) == -1) then
stop 'Unable to add task to server'
if (.True.) then
! Avoid openMP
integral = ao_two_e_integral(1,1,1,1)
endif
enddo
deallocate(task)
integer, external :: zmq_set_running
if (zmq_set_running(zmq_to_qp_run_socket) == -1) then
print *, irp_here, ': Failed in zmq_set_running'
endif
integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_socket_pull
call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'ao_integrals')
PROVIDE nproc
!$OMP PARALLEL DEFAULT(shared) private(i) num_threads(nproc+1)
i = omp_get_thread_num()
if (i==0) then
call ao_two_e_integrals_in_map_collector(zmq_socket_pull)
else
call ao_two_e_integrals_in_map_slave_inproc(i)
character(len=:), allocatable :: task
allocate(character(len=ao_num*12) :: task)
write(fmt,*) '(', ao_num, '(I5,X,I5,''|''))'
do l=1,ao_num
write(task,fmt) (i,l, i=1,l)
integer, external :: add_task_to_taskserver
if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task)) == -1) then
stop 'Unable to add task to server'
endif
!$OMP END PARALLEL
enddo
deallocate(task)
call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'ao_integrals')
integer, external :: zmq_set_running
if (zmq_set_running(zmq_to_qp_run_socket) == -1) then
print *, irp_here, ': Failed in zmq_set_running'
endif
PROVIDE nproc
!$OMP PARALLEL DEFAULT(shared) private(i) num_threads(nproc+1)
i = omp_get_thread_num()
if (i==0) then
call ao_two_e_integrals_in_map_collector(zmq_socket_pull)
else
call ao_two_e_integrals_in_map_slave_inproc(i)
endif
!$OMP END PARALLEL
call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'ao_integrals')
print*, 'Sorting the map'
call map_sort(ao_integrals_map)
call cpu_time(cpu_2)
call wall_time(wall_2)
integer(map_size_kind) :: get_ao_map_size, ao_map_size
ao_map_size = get_ao_map_size()
print*, 'Sorting the map'
call map_sort(ao_integrals_map)
call cpu_time(cpu_2)
call wall_time(wall_2)
integer(map_size_kind) :: get_ao_map_size, ao_map_size
ao_map_size = get_ao_map_size()
print*, 'AO integrals provided:'
print*, ' Size of AO map : ', map_mb(ao_integrals_map) ,'MB'
print*, ' Number of AO integrals :', ao_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)), ' )'
print*, 'AO integrals provided:'
print*, ' Size of AO map : ', map_mb(ao_integrals_map) ,'MB'
print*, ' Number of AO integrals :', ao_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_two_e_integrals_in_map = .True.
ao_two_e_integrals_in_map = .True.
if (write_ao_two_e_integrals.and.mpi_master) then
call ezfio_set_work_empty(.False.)
call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map)
call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read')
endif
if (write_ao_two_e_integrals.and.mpi_master) then
call ezfio_set_work_empty(.False.)
call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map)
call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read')
endif
END_PROVIDER

2
src/mo_two_e_ints/mo_bi_integrals.irp.f
View File

@ -189,7 +189,6 @@ subroutine add_integrals_to_map(mask_ijkl)
two_e_tmp_2 = 0.d0
do j1 = 1,ao_num
call get_ao_two_e_integrals(j1,k1,l1,ao_num,two_e_tmp_0(1,j1))
! call compute_ao_two_e_integrals(j1,k1,l1,ao_num,two_e_tmp_0(1,j1))
enddo
do j1 = 1,ao_num
kmax = 0
@ -747,7 +746,6 @@ subroutine add_integrals_to_map_no_exit_34(mask_ijkl)
two_e_tmp_2 = 0.d0
do j1 = 1,ao_num
call get_ao_two_e_integrals(j1,k1,l1,ao_num,two_e_tmp_0(1,j1))
! call compute_ao_two_e_integrals(j1,k1,l1,ao_num,two_e_tmp_0(1,j1))
enddo
do j1 = 1,ao_num
kmax = 0