diff --git a/src/bi_ort_ints/bi_ort_ints.irp.f b/src/bi_ort_ints/bi_ort_ints.irp.f index eae0affe..5653a2e2 100644 --- a/src/bi_ort_ints/bi_ort_ints.irp.f +++ b/src/bi_ort_ints/bi_ort_ints.irp.f @@ -60,60 +60,60 @@ subroutine test_5idx do l = 1, mo_num do m = 1, mo_num -! new = three_e_5_idx_direct_bi_ort(m,l,j,k,i) -! ref = three_e_5_idx_direct_bi_ort_old(m,l,j,k,i) -! contrib = dabs(new - ref) -! accu += contrib -! if(contrib .gt. 1.d-10)then -! print*,'direct' -! print*,i,k,j,l,m -! print*,ref,new,contrib -! stop -! endif + new = three_e_5_idx_direct_bi_ort(m,l,j,k,i) + ref = three_e_5_idx_direct_bi_ort_old(m,l,j,k,i) + contrib = dabs(new - ref) + accu += contrib + if(contrib .gt. 1.d-10)then + print*,'direct' + print*,i,k,j,l,m + print*,ref,new,contrib + stop + endif + + new = three_e_5_idx_exch12_bi_ort(m,l,j,k,i) + ref = three_e_5_idx_exch12_bi_ort_old(m,l,j,k,i) + contrib = dabs(new - ref) + accu += contrib + if(contrib .gt. 1.d-10)then + print*,'exch12' + print*,i,k,j,l,m + print*,ref,new,contrib + stop + endif ! -! new = three_e_5_idx_exch12_bi_ort(m,l,j,k,i) -! ref = three_e_5_idx_exch12_bi_ort_old(m,l,j,k,i) -! contrib = dabs(new - ref) -! accu += contrib -! if(contrib .gt. 1.d-10)then -! print*,'exch12' -! print*,i,k,j,l,m -! print*,ref,new,contrib -! stop -! endif + new = three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) + ref = three_e_5_idx_cycle_1_bi_ort_old(m,l,j,k,i) + contrib = dabs(new - ref) + accu += contrib + if(contrib .gt. 1.d-10)then + print*,'cycle1' + print*,i,k,j,l,m + print*,ref,new,contrib + stop + endif -! new = three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) -! ref = three_e_5_idx_cycle_1_bi_ort_old(m,l,j,k,i) -! contrib = dabs(new - ref) -! accu += contrib -! if(contrib .gt. 1.d-10)then -! print*,'cycle1' -! print*,i,k,j,l,m -! print*,ref,new,contrib -! stop -! endif + new = three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) + ref = three_e_5_idx_cycle_2_bi_ort_old(m,l,j,k,i) + contrib = dabs(new - ref) + accu += contrib + if(contrib .gt. 1.d-10)then + print*,'cycle2' + print*,i,k,j,l,m + print*,ref,new,contrib + stop + endif -! new = three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) -! ref = three_e_5_idx_cycle_2_bi_ort_old(m,l,j,k,i) -! contrib = dabs(new - ref) -! accu += contrib -! if(contrib .gt. 1.d-10)then -! print*,'cycle2' -! print*,i,k,j,l,m -! print*,ref,new,contrib -! stop -! endif - -! new = three_e_5_idx_exch23_bi_ort(m,l,j,k,i) -! ref = three_e_5_idx_exch23_bi_ort_old(m,l,j,k,i) -! contrib = dabs(new - ref) -! accu += contrib -! if(contrib .gt. 1.d-10)then -! print*,'exch23' -! print*,i,k,j,l,m -! print*,ref,new,contrib -! stop -! endif + new = three_e_5_idx_exch23_bi_ort(m,l,j,k,i) + ref = three_e_5_idx_exch23_bi_ort_old(m,l,j,k,i) + contrib = dabs(new - ref) + accu += contrib + if(contrib .gt. 1.d-10)then + print*,'exch23' + print*,i,k,j,l,m + print*,ref,new,contrib + stop + endif new = three_e_5_idx_exch13_bi_ort(m,l,j,k,i) ref = three_e_5_idx_exch13_bi_ort_old(m,l,j,k,i) diff --git a/src/bi_ort_ints/three_body_ijmkl.irp.f b/src/bi_ort_ints/three_body_ijmkl.irp.f index 9f316771..c9e88ab9 100644 --- a/src/bi_ort_ints/three_body_ijmkl.irp.f +++ b/src/bi_ort_ints/three_body_ijmkl.irp.f @@ -1,7 +1,11 @@ ! --- - BEGIN_PROVIDER [ double precision, three_e_5_idx_direct_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)] -&BEGIN_PROVIDER [ double precision, three_e_5_idx_exch12_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)] + BEGIN_PROVIDER [ double precision, three_e_5_idx_direct_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)] +&BEGIN_PROVIDER [ double precision, three_e_5_idx_exch12_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)] +&BEGIN_PROVIDER [ double precision, three_e_5_idx_exch23_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)] +&BEGIN_PROVIDER [ double precision, three_e_5_idx_exch13_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)] +&BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_1_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)] +&BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_2_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)] BEGIN_DOC ! @@ -16,18 +20,22 @@ integer :: i, j, k, m, l double precision :: wall1, wall0 integer :: ipoint - double precision, allocatable :: grad_mli(:,:,:), m2grad_r(:,:,:,:), m2grad_l(:,:,:,:) - double precision, allocatable :: tmp_mat(:,:,:,:), orb_mat(:,:,:) - allocate(m2grad_r(n_points_final_grid,3,mo_num,mo_num)) - allocate(m2grad_l(n_points_final_grid,3,mo_num,mo_num)) - allocate(tmp_mat(mo_num,mo_num,mo_num,mo_num)) + double precision, allocatable :: grad_mli(:,:,:), orb_mat(:,:,:) + double precision, allocatable :: lk_grad_mi(:,:,:,:), rk_grad_im(:,:,:,:) + double precision, allocatable :: lm_grad_ik(:,:,:,:), rm_grad_ik(:,:,:,:) + double precision, allocatable :: tmp_mat(:,:,:,:) + allocate(rm_grad_ik(n_points_final_grid,3,mo_num,mo_num)) + allocate(lm_grad_ik(n_points_final_grid,3,mo_num,mo_num)) + allocate(lk_grad_mi(n_points_final_grid,3,mo_num,mo_num)) + allocate(rk_grad_im(n_points_final_grid,3,mo_num,mo_num)) allocate(grad_mli(n_points_final_grid,mo_num,mo_num)) allocate(orb_mat(n_points_final_grid,mo_num,mo_num)) + allocate(tmp_mat(mo_num,mo_num,mo_num,mo_num)) provide mos_r_in_r_array_transp mos_l_in_r_array_transp PROVIDE mo_l_coef mo_r_coef int2_grad1_u12_bimo_t - print *, ' Providing the three_e_5_idx_direct_bi_ort ...' + print *, ' Providing the three_e_5_idx_bi_ort ...' call wall_time(wall0) do m = 1, mo_num @@ -38,7 +46,8 @@ !$OMP SHARED (m,mo_num,n_points_final_grid, & !$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, & !$OMP int2_grad1_u12_bimo_t, final_weight_at_r_vector, & - !$OMP m2grad_r, m2grad_l, grad_mli, tmp_mat, orb_mat) + !$OMP rm_grad_ik, lm_grad_ik, rk_grad_im, lk_grad_mi, & + !$OMP grad_mli, tmp_mat, orb_mat) !$OMP DO COLLAPSE(2) do i=1,mo_num do l=1,mo_num @@ -51,13 +60,21 @@ orb_mat(ipoint,l,i) = mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,i) - m2grad_l(ipoint,1,l,i) = mos_l_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,1,l,i) * final_weight_at_r_vector(ipoint) - m2grad_l(ipoint,2,l,i) = mos_l_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i) * final_weight_at_r_vector(ipoint) - m2grad_l(ipoint,3,l,i) = mos_l_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i) * final_weight_at_r_vector(ipoint) + lm_grad_ik(ipoint,1,l,i) = mos_l_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,1,l,i) * final_weight_at_r_vector(ipoint) + lm_grad_ik(ipoint,2,l,i) = mos_l_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i) * final_weight_at_r_vector(ipoint) + lm_grad_ik(ipoint,3,l,i) = mos_l_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i) * final_weight_at_r_vector(ipoint) - m2grad_r(ipoint,1,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,1,l,i) - m2grad_r(ipoint,2,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i) - m2grad_r(ipoint,3,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i) + rm_grad_ik(ipoint,1,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,1,l,i) + rm_grad_ik(ipoint,2,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i) + rm_grad_ik(ipoint,3,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i) + + lk_grad_mi(ipoint,1,l,i) = mos_l_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,1,m,i) * final_weight_at_r_vector(ipoint) + lk_grad_mi(ipoint,2,l,i) = mos_l_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,2,m,i) * final_weight_at_r_vector(ipoint) + lk_grad_mi(ipoint,3,l,i) = mos_l_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,3,m,i) * final_weight_at_r_vector(ipoint) + + rk_grad_im(ipoint,1,l,i) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,1,i,m) + rk_grad_im(ipoint,2,l,i) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,2,i,m) + rk_grad_im(ipoint,3,l,i) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,3,i,m) enddo enddo @@ -84,8 +101,8 @@ !$OMP END PARALLEL DO call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0, & - m2grad_l, 3*n_points_final_grid, & - m2grad_r, 3*n_points_final_grid, 0.d0, & + lm_grad_ik, 3*n_points_final_grid, & + rm_grad_ik, 3*n_points_final_grid, 0.d0, & tmp_mat, mo_num*mo_num) !$OMP PARALLEL DO PRIVATE(i,j,k,l) @@ -101,83 +118,6 @@ enddo !$OMP END PARALLEL DO - enddo - - call wall_time(wall1) - print *, ' wall time for three_e_5_idx_direct_bi_ort', wall1 - wall0 - -END_PROVIDER - -! --- - - BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_1_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)] -&BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_2_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)] -&BEGIN_PROVIDER [ double precision, three_e_5_idx_exch23_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)] -&BEGIN_PROVIDER [ double precision, three_e_5_idx_exch13_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)] - - BEGIN_DOC - ! - ! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs - ! - ! three_e_5_idx_direct_bi_ort(m,l,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO - ! - ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign - END_DOC - - implicit none - integer :: i, j, k, m, l - double precision :: wall1, wall0 - integer :: ipoint - double precision, allocatable :: lk_grad_mi(:,:,:,:), rk_grad_im(:,:,:,:) - double precision, allocatable :: lm_grad_ik(:,:,:,:), rm_grad_ik(:,:,:,:) - double precision, allocatable :: tmp_mat(:,:,:,:) - allocate(lk_grad_mi(n_points_final_grid,3,mo_num,mo_num)) - allocate(lm_grad_ik(n_points_final_grid,3,mo_num,mo_num)) - allocate(rk_grad_im(n_points_final_grid,3,mo_num,mo_num)) - allocate(rm_grad_ik(n_points_final_grid,3,mo_num,mo_num)) - allocate(tmp_mat(mo_num,mo_num,mo_num,mo_num)) - - provide mos_r_in_r_array_transp mos_l_in_r_array_transp - PROVIDE mo_l_coef mo_r_coef int2_grad1_u12_bimo_t - - print *, ' Providing the three_e_5_idx_cycle_bi_ort ...' - call wall_time(wall0) - - do m = 1, mo_num - - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,l,ipoint) & - !$OMP SHARED (m,mo_num,n_points_final_grid, & - !$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, & - !$OMP int2_grad1_u12_bimo_t, final_weight_at_r_vector, & - !$OMP rk_grad_im, rm_grad_ik, lk_grad_mi, lm_grad_ik, tmp_mat) - !$OMP DO COLLAPSE(2) - do i=1,mo_num - do l=1,mo_num - do ipoint=1, n_points_final_grid - lk_grad_mi(ipoint,1,l,i) = mos_l_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,1,m,i) * final_weight_at_r_vector(ipoint) - lk_grad_mi(ipoint,2,l,i) = mos_l_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,2,m,i) * final_weight_at_r_vector(ipoint) - lk_grad_mi(ipoint,3,l,i) = mos_l_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,3,m,i) * final_weight_at_r_vector(ipoint) - - lm_grad_ik(ipoint,1,l,i) = mos_l_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,1,l,i) * final_weight_at_r_vector(ipoint) - lm_grad_ik(ipoint,2,l,i) = mos_l_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i) * final_weight_at_r_vector(ipoint) - lm_grad_ik(ipoint,3,l,i) = mos_l_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i) * final_weight_at_r_vector(ipoint) - - rm_grad_ik(ipoint,1,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,1,l,i) - rm_grad_ik(ipoint,2,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i) - rm_grad_ik(ipoint,3,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i) - - rk_grad_im(ipoint,1,l,i) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,1,i,m) - rk_grad_im(ipoint,2,l,i) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,2,i,m) - rk_grad_im(ipoint,3,l,i) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,3,i,m) - - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0, & lk_grad_mi, 3*n_points_final_grid, & rm_grad_ik, 3*n_points_final_grid, 0.d0, & @@ -237,14 +177,11 @@ END_PROVIDER enddo enddo !$OMP END PARALLEL DO + enddo call wall_time(wall1) - print *, ' wall time for three_e_5_idx_cycle_bi_ort', wall1 - wall0 + print *, ' wall time for three_e_5_idx_bi_ort', wall1 - wall0 END_PROVIDER -! --- - - -