From 84445aa591668092ac3f0a8c80508fb72f099acc Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Thu, 8 Aug 2024 11:39:42 +0200 Subject: [PATCH] Combine calculation of Left & RIGHT MOs in r --- plugins/local/bi_ort_ints/bi_ort_ints.irp.f | 56 +++++- plugins/local/bi_ort_ints/no_dressing.irp.f | 2 + .../local/bi_ortho_mos/bi_ort_mos_in_r.irp.f | 165 ++++++------------ .../bi_ortho_mos/bi_ort_mos_in_r_old.irp.f | 137 +++++++++++++++ plugins/local/tc_bi_ortho/tc_utils.irp.f | 2 + 5 files changed, 246 insertions(+), 116 deletions(-) create mode 100644 plugins/local/bi_ortho_mos/bi_ort_mos_in_r_old.irp.f diff --git a/plugins/local/bi_ort_ints/bi_ort_ints.irp.f b/plugins/local/bi_ort_ints/bi_ort_ints.irp.f index 0349c731..b691e47e 100644 --- a/plugins/local/bi_ort_ints/bi_ort_ints.irp.f +++ b/plugins/local/bi_ort_ints/bi_ort_ints.irp.f @@ -17,12 +17,14 @@ program bi_ort_ints ! call test_3e ! call test_5idx ! call test_5idx2 - call test_4idx() +! call test_4idx() !call test_4idx_n4() !call test_4idx2() !call test_5idx2 !call test_5idx + call test_mos_in_r() + end subroutine test_5idx2 @@ -472,4 +474,56 @@ subroutine test_4idx() return end +! --- + +subroutine test_mos_in_r() + + implicit none + + integer :: i, j + double precision :: err_tot, nrm_tot, err_loc, acc_thr + + PROVIDE mos_l_in_r_array_transp_old mos_r_in_r_array_transp_old + PROVIDE mos_l_in_r_array_transp mos_r_in_r_array_transp + + acc_thr = 1d-12 + + err_tot = 0.d0 + nrm_tot = 0.d0 + do i = 1, mo_num + do j = 1, n_points_final_grid + err_loc = dabs(mos_l_in_r_array_transp_old(j,i) - mos_l_in_r_array_transp(j,i)) + if(err_loc > acc_thr) then + print*, " error on", j, i + print*, " old res", mos_l_in_r_array_transp_old(j,i) + print*, " new res", mos_l_in_r_array_transp (j,i) + stop + endif + err_tot = err_tot + err_loc + nrm_tot = nrm_tot + dabs(mos_l_in_r_array_transp_old(j,i)) + enddo + enddo + print *, ' absolute accuracy on mos_l_in_r_array_transp (%) =', 100.d0 * err_tot / nrm_tot + + err_tot = 0.d0 + nrm_tot = 0.d0 + do i = 1, mo_num + do j = 1, n_points_final_grid + err_loc = dabs(mos_r_in_r_array_transp_old(j,i) - mos_r_in_r_array_transp(j,i)) + if(err_loc > acc_thr) then + print*, " error on", j, i + print*, " old res", mos_r_in_r_array_transp_old(j,i) + print*, " new res", mos_r_in_r_array_transp (j,i) + stop + endif + err_tot = err_tot + err_loc + nrm_tot = nrm_tot + dabs(mos_r_in_r_array_transp_old(j,i)) + enddo + enddo + print *, ' absolute accuracy on mos_r_in_r_array_transp (%) =', 100.d0 * err_tot / nrm_tot + + return +end + +! --- diff --git a/plugins/local/bi_ort_ints/no_dressing.irp.f b/plugins/local/bi_ort_ints/no_dressing.irp.f index fd2c6285..3bbf31f9 100644 --- a/plugins/local/bi_ort_ints/no_dressing.irp.f +++ b/plugins/local/bi_ort_ints/no_dressing.irp.f @@ -710,6 +710,8 @@ BEGIN_PROVIDER [double precision, noL_0e] endif + print*, " noL_0e =", noL_0e + END_PROVIDER ! --- diff --git a/plugins/local/bi_ortho_mos/bi_ort_mos_in_r.irp.f b/plugins/local/bi_ortho_mos/bi_ort_mos_in_r.irp.f index 25572854..15ed2ce4 100644 --- a/plugins/local/bi_ortho_mos/bi_ort_mos_in_r.irp.f +++ b/plugins/local/bi_ortho_mos/bi_ort_mos_in_r.irp.f @@ -1,135 +1,70 @@ -! TODO: left & right MO without duplicate AO calculation - -! --- - -BEGIN_PROVIDER[double precision, mos_r_in_r_array, (mo_num, n_points_final_grid)] + BEGIN_PROVIDER[double precision, mos_l_in_r_array_transp, (n_points_final_grid, mo_num)] +&BEGIN_PROVIDER[double precision, mos_r_in_r_array_transp, (n_points_final_grid, mo_num)] BEGIN_DOC - ! mos_in_r_array(i,j) = value of the ith RIGHT mo on the jth grid point + ! + ! mos_l_in_r_array_transp(i,j) = value of the jth left-mo on the ith grid point + ! mos_r_in_r_array_transp(i,j) = value of the jth right-mo on the ith grid point + ! END_DOC implicit none - integer :: i, j - double precision :: mos_array(mo_num), r(3) - !$OMP PARALLEL DO & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i, j, r, mos_array) & - !$OMP SHARED (mos_r_in_r_array, n_points_final_grid, mo_num, final_grid_points) + integer :: i + double precision :: tt0, tt1, tt2, tt3 + double precision :: r(3) + double precision, allocatable :: aos_r(:,:) + + call cpu_time(tt0) + + allocate(aos_r(ao_num,n_points_final_grid)) + + ! provide everything required before OpenMP + r(1) = final_grid_points(1,1) + r(2) = final_grid_points(2,1) + r(3) = final_grid_points(3,1) + call give_all_aos_at_r(r, aos_r(1,1)) + + + call cpu_time(tt2) + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i, r) & + !$OMP SHARED(n_points_final_grid, final_grid_points, aos_r) + !$OMP DO do i = 1, n_points_final_grid r(1) = final_grid_points(1,i) r(2) = final_grid_points(2,i) r(3) = final_grid_points(3,i) - call give_all_mos_r_at_r(r, mos_array) - do j = 1, mo_num - mos_r_in_r_array(j,i) = mos_array(j) - enddo + call give_all_aos_at_r(r, aos_r(1,i)) enddo - !$OMP END PARALLEL DO - -END_PROVIDER + !$OMP END DO + !$OMP END PARALLEL -! --- + call cpu_time(tt3) + write(*,"(A,2X,F15.7)") ' wall time for AOs on r (sec) = ', (tt3 - tt2) -BEGIN_PROVIDER[double precision, mos_r_in_r_array_transp, (n_points_final_grid, mo_num)] - BEGIN_DOC - ! mos_r_in_r_array_transp(i,j) = value of the jth mo on the ith grid point - END_DOC + call dgemm("T", "N", n_points_final_grid, mo_num, ao_num, & + 1.d0, & + aos_r(1,1), ao_num, & + mo_l_coef(1,1), ao_num, & + 0.d0, & + mos_l_in_r_array_transp(1,1), n_points_final_grid) - implicit none - integer :: i,j + call dgemm("T", "N", n_points_final_grid, mo_num, ao_num, & + 1.d0, & + aos_r(1,1), ao_num, & + mo_r_coef(1,1), ao_num, & + 0.d0, & + mos_r_in_r_array_transp(1,1), n_points_final_grid) - do i = 1, n_points_final_grid - do j = 1, mo_num - mos_r_in_r_array_transp(i,j) = mos_r_in_r_array(j,i) - enddo - enddo - -END_PROVIDER - -! --- - -subroutine give_all_mos_r_at_r(r, mos_r_array) - - BEGIN_DOC - ! mos_r_array(i) = ith RIGHT MO function evaluated at "r" - END_DOC - - implicit none - double precision, intent(in) :: r(3) - double precision, intent(out) :: mos_r_array(mo_num) - double precision :: aos_array(ao_num) - - call give_all_aos_at_r(r, aos_array) - call dgemv('N', mo_num, ao_num, 1.d0, mo_r_coef_transp, mo_num, aos_array, 1, 0.d0, mos_r_array, 1) - -end subroutine give_all_mos_r_at_r - -! --- - -BEGIN_PROVIDER[double precision, mos_l_in_r_array, (mo_num, n_points_final_grid)] - - BEGIN_DOC - ! mos_in_r_array(i,j) = value of the ith LEFT mo on the jth grid point - END_DOC - - implicit none - integer :: i, j - double precision :: mos_array(mo_num), r(3) - - !$OMP PARALLEL DO & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,r,mos_array,j) & - !$OMP SHARED(mos_l_in_r_array,n_points_final_grid,mo_num,final_grid_points) - do i = 1, n_points_final_grid - r(1) = final_grid_points(1,i) - r(2) = final_grid_points(2,i) - r(3) = final_grid_points(3,i) - call give_all_mos_l_at_r(r, mos_array) - do j = 1, mo_num - mos_l_in_r_array(j,i) = mos_array(j) - enddo - enddo - !$OMP END PARALLEL DO - -END_PROVIDER - -! --- - -subroutine give_all_mos_l_at_r(r, mos_l_array) - - BEGIN_DOC - ! mos_l_array(i) = ith LEFT MO function evaluated at "r" - END_DOC - - implicit none - double precision, intent(in) :: r(3) - double precision, intent(out) :: mos_l_array(mo_num) - double precision :: aos_array(ao_num) - - call give_all_aos_at_r(r, aos_array) - call dgemv('N', mo_num, ao_num, 1.d0, mo_l_coef_transp, mo_num, aos_array, 1, 0.d0, mos_l_array, 1) - -end subroutine give_all_mos_l_at_r - -! --- - -BEGIN_PROVIDER[double precision, mos_l_in_r_array_transp, (n_points_final_grid,mo_num)] - - BEGIN_DOC - ! mos_l_in_r_array_transp(i,j) = value of the jth mo on the ith grid point - END_DOC - - implicit none - integer :: i, j - - do i = 1, n_points_final_grid - do j = 1, mo_num - mos_l_in_r_array_transp(i,j) = mos_l_in_r_array(j,i) - enddo - enddo + deallocate(aos_r) + + call cpu_time(tt1) + write(*,"(A,2X,F15.7)") ' wall time for mos_l_in_r_array_transp & mos_r_in_r_array_transp (sec) = ', (tt1 - tt0) END_PROVIDER diff --git a/plugins/local/bi_ortho_mos/bi_ort_mos_in_r_old.irp.f b/plugins/local/bi_ortho_mos/bi_ort_mos_in_r_old.irp.f new file mode 100644 index 00000000..9fd671f8 --- /dev/null +++ b/plugins/local/bi_ortho_mos/bi_ort_mos_in_r_old.irp.f @@ -0,0 +1,137 @@ + +! TODO: left & right MO without duplicate AO calculation + +! --- + +BEGIN_PROVIDER[double precision, mos_r_in_r_array, (mo_num, n_points_final_grid)] + + BEGIN_DOC + ! mos_in_r_array(i,j) = value of the ith RIGHT mo on the jth grid point + END_DOC + + implicit none + integer :: i, j + double precision :: mos_array(mo_num), r(3) + + !$OMP PARALLEL DO & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i, j, r, mos_array) & + !$OMP SHARED (mos_r_in_r_array, n_points_final_grid, mo_num, final_grid_points) + do i = 1, n_points_final_grid + r(1) = final_grid_points(1,i) + r(2) = final_grid_points(2,i) + r(3) = final_grid_points(3,i) + call give_all_mos_r_at_r(r, mos_array) + do j = 1, mo_num + mos_r_in_r_array(j,i) = mos_array(j) + enddo + enddo + !$OMP END PARALLEL DO + +END_PROVIDER + +! --- + +BEGIN_PROVIDER[double precision, mos_r_in_r_array_transp_old, (n_points_final_grid, mo_num)] + + BEGIN_DOC + ! mos_r_in_r_array_transp_old(i,j) = value of the jth mo on the ith grid point + END_DOC + + implicit none + integer :: i,j + + do i = 1, n_points_final_grid + do j = 1, mo_num + mos_r_in_r_array_transp_old(i,j) = mos_r_in_r_array(j,i) + enddo + enddo + +END_PROVIDER + +! --- + +subroutine give_all_mos_r_at_r(r, mos_r_array) + + BEGIN_DOC + ! mos_r_array(i) = ith RIGHT MO function evaluated at "r" + END_DOC + + implicit none + double precision, intent(in) :: r(3) + double precision, intent(out) :: mos_r_array(mo_num) + double precision :: aos_array(ao_num) + + call give_all_aos_at_r(r, aos_array) + call dgemv('N', mo_num, ao_num, 1.d0, mo_r_coef_transp, mo_num, aos_array, 1, 0.d0, mos_r_array, 1) + +end subroutine give_all_mos_r_at_r + +! --- + +BEGIN_PROVIDER[double precision, mos_l_in_r_array, (mo_num, n_points_final_grid)] + + BEGIN_DOC + ! mos_in_r_array(i,j) = value of the ith LEFT mo on the jth grid point + END_DOC + + implicit none + integer :: i, j + double precision :: mos_array(mo_num), r(3) + + !$OMP PARALLEL DO & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,r,mos_array,j) & + !$OMP SHARED(mos_l_in_r_array,n_points_final_grid,mo_num,final_grid_points) + do i = 1, n_points_final_grid + r(1) = final_grid_points(1,i) + r(2) = final_grid_points(2,i) + r(3) = final_grid_points(3,i) + call give_all_mos_l_at_r(r, mos_array) + do j = 1, mo_num + mos_l_in_r_array(j,i) = mos_array(j) + enddo + enddo + !$OMP END PARALLEL DO + +END_PROVIDER + +! --- + +subroutine give_all_mos_l_at_r(r, mos_l_array) + + BEGIN_DOC + ! mos_l_array(i) = ith LEFT MO function evaluated at "r" + END_DOC + + implicit none + double precision, intent(in) :: r(3) + double precision, intent(out) :: mos_l_array(mo_num) + double precision :: aos_array(ao_num) + + call give_all_aos_at_r(r, aos_array) + call dgemv('N', mo_num, ao_num, 1.d0, mo_l_coef_transp, mo_num, aos_array, 1, 0.d0, mos_l_array, 1) + +end subroutine give_all_mos_l_at_r + +! --- + +BEGIN_PROVIDER[double precision, mos_l_in_r_array_transp_old, (n_points_final_grid,mo_num)] + + BEGIN_DOC + ! mos_l_in_r_array_transp_old(i,j) = value of the jth mo on the ith grid point + END_DOC + + implicit none + integer :: i, j + + do i = 1, n_points_final_grid + do j = 1, mo_num + mos_l_in_r_array_transp_old(i,j) = mos_l_in_r_array(j,i) + enddo + enddo + +END_PROVIDER + +! --- + diff --git a/plugins/local/tc_bi_ortho/tc_utils.irp.f b/plugins/local/tc_bi_ortho/tc_utils.irp.f index 2aa148a3..4dfd4316 100644 --- a/plugins/local/tc_bi_ortho/tc_utils.irp.f +++ b/plugins/local/tc_bi_ortho/tc_utils.irp.f @@ -8,6 +8,8 @@ subroutine write_tc_energy() double precision :: E_1e, E_2e, E_3e double precision, allocatable :: E_TC_tmp(:), E_1e_tmp(:), E_2e_tmp(:), E_3e_tmp(:) + call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,1), psi_det(1,1,1), N_int, hmono, htwoe, hthree, htot) + ! GS ! ---