diff --git a/src/ao_basis/aos.irp.f b/src/ao_basis/aos.irp.f index f747b1de..e95ac711 100644 --- a/src/ao_basis/aos.irp.f +++ b/src/ao_basis/aos.irp.f @@ -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 diff --git a/src/ao_one_e_ints/ao_overlap.irp.f b/src/ao_one_e_ints/ao_overlap.irp.f index d7a4d640..d7300936 100644 --- a/src/ao_one_e_ints/ao_overlap.irp.f +++ b/src/ao_one_e_ints/ao_overlap.irp.f @@ -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 diff --git a/src/ao_one_e_ints/pot_ao_ints.irp.f b/src/ao_one_e_ints/pot_ao_ints.irp.f index 486ff534..4108ce71 100644 --- a/src/ao_one_e_ints/pot_ao_ints.irp.f +++ b/src/ao_one_e_ints/pot_ao_ints.irp.f @@ -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 diff --git a/src/ao_one_e_ints/screening.irp.f b/src/ao_one_e_ints/screening.irp.f index aa23a9b8..1bbe3c73 100644 --- a/src/ao_one_e_ints/screening.irp.f +++ b/src/ao_one_e_ints/screening.irp.f @@ -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 diff --git a/src/ao_two_e_ints/screening.irp.f b/src/ao_two_e_ints/screening.irp.f index 5095deec..d3230370 100644 --- a/src/ao_two_e_ints/screening.irp.f +++ b/src/ao_two_e_ints/screening.irp.f @@ -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 diff --git a/src/ao_two_e_ints/two_e_integrals.irp.f b/src/ao_two_e_ints/two_e_integrals.irp.f index 973a9011..b6e959d7 100644 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ b/src/ao_two_e_ints/two_e_integrals.irp.f @@ -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 diff --git a/src/mo_two_e_ints/mo_bi_integrals.irp.f b/src/mo_two_e_ints/mo_bi_integrals.irp.f index a9983e51..b926d688 100644 --- a/src/mo_two_e_ints/mo_bi_integrals.irp.f +++ b/src/mo_two_e_ints/mo_bi_integrals.irp.f @@ -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