diff --git a/src/tc_bi_ortho/normal_ordered.irp.f b/src/tc_bi_ortho/normal_ordered.irp.f index a092762b..59e78b92 100644 --- a/src/tc_bi_ortho/normal_ordered.irp.f +++ b/src/tc_bi_ortho/normal_ordered.irp.f @@ -22,8 +22,6 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_ print*,' Providing normal_two_body_bi_orth ...' call wall_time(wall0) - PROVIDE N_int - if(read_tc_norm_ord) then open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/normal_two_body_bi_orth', action="read") @@ -48,12 +46,13 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_ endif ! opposite spin double excitations : s1 /= s2 - normal_two_body_bi_orth(:,:,:,:) = no_aba_contraction(:,:,:,:) + PROVIDE no_aba_contraction - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (hh1, h1, hh2, h2, pp1, p1, pp2, p2, hthree_aab, hthree_aaa) & - !$OMP SHARED (N_int, n_act_orb, list_act, Ne, occ, normal_two_body_bi_orth) + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (hh1, h1, hh2, h2, pp1, p1, pp2, p2, hthree_aab, hthree_aaa) & + !$OMP SHARED (N_int, n_act_orb, list_act, Ne, occ, normal_two_body_bi_orth, & + !$OMP no_aba_contraction) !$OMP DO SCHEDULE (static) do hh1 = 1, n_act_orb h1 = list_act(hh1) @@ -97,7 +96,7 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_ endif - normal_two_body_bi_orth(p2,h2,p1,h1) = 0.5d0*(hthree_aab + hthree_aaa) + normal_two_body_bi_orth(p2,h2,p1,h1) = no_aba_contraction(p2,h2,p1,h1) + 0.5d0*(hthree_aab + hthree_aaa) enddo enddo enddo @@ -124,103 +123,103 @@ END_PROVIDER ! --- -subroutine give_aaa_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree) - - BEGIN_DOC - ! pure same spin contribution to same spin double excitation s1=h1,p1, s2=h2,p2, with s1==s2 - END_DOC - - use bitmasks ! you need to include the bitmasks_module.f90 features - - implicit none - integer, intent(in) :: Nint, h1, h2, p1, p2 - integer, intent(in) :: Ne(2), occ(Nint*bit_kind_size,2) - double precision, intent(out) :: hthree - integer :: ii,i - double precision :: int_direct,int_exc_12,int_exc_13,int_exc_23 - double precision :: integral,int_exc_l,int_exc_ll - - hthree = 0.d0 - do ii = 1, Ne(2) ! purely closed shell part - i = occ(ii,2) - - call give_integrals_3_body_bi_ort(i, p2, p1, i, h2, h1, integral) - int_direct = -1.d0 * integral - - call give_integrals_3_body_bi_ort(p2, p1, i, i, h2, h1, integral) - int_exc_l = -1.d0 * integral - - call give_integrals_3_body_bi_ort(p1, i, p2, i, h2, h1, integral) - int_exc_ll= -1.d0 * integral - - call give_integrals_3_body_bi_ort(p2, i, p1, i, h2, h1, integral) - int_exc_12= -1.d0 * integral - - call give_integrals_3_body_bi_ort(p1, p2, i, i, h2, h1, integral) - int_exc_13= -1.d0 * integral - - call give_integrals_3_body_bi_ort(i, p1, p2, i, h2, h1, integral) - int_exc_23= -1.d0 * integral - - hthree += 1.d0 * int_direct + int_exc_l + int_exc_ll - (int_exc_12 + int_exc_13 + int_exc_23) - enddo - - do ii = Ne(2)+1,Ne(1) ! purely open-shell part - i = occ(ii,1) - - call give_integrals_3_body_bi_ort(i, p2, p1, i, h2, h1, integral) - int_direct = -1.d0 * integral - - call give_integrals_3_body_bi_ort(p2, p1, i , i, h2, h1, integral) - int_exc_l = -1.d0 * integral - - call give_integrals_3_body_bi_ort(p1, i, p2, i, h2, h1, integral) - int_exc_ll = -1.d0 * integral - - call give_integrals_3_body_bi_ort(p2, i, p1, i, h2, h1, integral) - int_exc_12 = -1.d0 * integral - - call give_integrals_3_body_bi_ort(p1, p2, i, i, h2, h1, integral) - int_exc_13 = -1.d0 * integral - - call give_integrals_3_body_bi_ort(i, p1, p2, i, h2, h1, integral) - int_exc_23 = -1.d0 * integral - - hthree += 1.d0 * int_direct + 0.5d0 * (int_exc_l + int_exc_ll - (int_exc_12 + int_exc_13 + int_exc_23)) - enddo - - return -end +!subroutine give_aaa_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree) +! +! BEGIN_DOC +! ! pure same spin contribution to same spin double excitation s1=h1,p1, s2=h2,p2, with s1==s2 +! END_DOC +! +! use bitmasks ! you need to include the bitmasks_module.f90 features +! +! implicit none +! integer, intent(in) :: Nint, h1, h2, p1, p2 +! integer, intent(in) :: Ne(2), occ(Nint*bit_kind_size,2) +! double precision, intent(out) :: hthree +! integer :: ii,i +! double precision :: int_direct,int_exc_12,int_exc_13,int_exc_23 +! double precision :: integral,int_exc_l,int_exc_ll +! +! hthree = 0.d0 +! do ii = 1, Ne(2) ! purely closed shell part +! i = occ(ii,2) +! +! call give_integrals_3_body_bi_ort(i, p2, p1, i, h2, h1, integral) +! int_direct = -1.d0 * integral +! +! call give_integrals_3_body_bi_ort(p2, p1, i, i, h2, h1, integral) +! int_exc_l = -1.d0 * integral +! +! call give_integrals_3_body_bi_ort(p1, i, p2, i, h2, h1, integral) +! int_exc_ll= -1.d0 * integral +! +! call give_integrals_3_body_bi_ort(p2, i, p1, i, h2, h1, integral) +! int_exc_12= -1.d0 * integral +! +! call give_integrals_3_body_bi_ort(p1, p2, i, i, h2, h1, integral) +! int_exc_13= -1.d0 * integral +! +! call give_integrals_3_body_bi_ort(i, p1, p2, i, h2, h1, integral) +! int_exc_23= -1.d0 * integral +! +! hthree += 1.d0 * int_direct + int_exc_l + int_exc_ll - (int_exc_12 + int_exc_13 + int_exc_23) +! enddo +! +! do ii = Ne(2)+1,Ne(1) ! purely open-shell part +! i = occ(ii,1) +! +! call give_integrals_3_body_bi_ort(i, p2, p1, i, h2, h1, integral) +! int_direct = -1.d0 * integral +! +! call give_integrals_3_body_bi_ort(p2, p1, i , i, h2, h1, integral) +! int_exc_l = -1.d0 * integral +! +! call give_integrals_3_body_bi_ort(p1, i, p2, i, h2, h1, integral) +! int_exc_ll = -1.d0 * integral +! +! call give_integrals_3_body_bi_ort(p2, i, p1, i, h2, h1, integral) +! int_exc_12 = -1.d0 * integral +! +! call give_integrals_3_body_bi_ort(p1, p2, i, i, h2, h1, integral) +! int_exc_13 = -1.d0 * integral +! +! call give_integrals_3_body_bi_ort(i, p1, p2, i, h2, h1, integral) +! int_exc_23 = -1.d0 * integral +! +! hthree += 1.d0 * int_direct + 0.5d0 * (int_exc_l + int_exc_ll - (int_exc_12 + int_exc_13 + int_exc_23)) +! enddo +! +! return +!end ! --- -subroutine give_aab_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree) - - use bitmasks ! you need to include the bitmasks_module.f90 features - - implicit none - integer, intent(in) :: Nint, h1, h2, p1, p2 - integer, intent(in) :: Ne(2), occ(Nint*bit_kind_size,2) - double precision, intent(out) :: hthree - integer :: ii, i - double precision :: int_direct, int_exc_12, int_exc_13, int_exc_23 - double precision :: integral, int_exc_l, int_exc_ll - - hthree = 0.d0 - do ii = 1, Ne(2) ! purely closed shell part - i = occ(ii,2) - - call give_integrals_3_body_bi_ort(p2, p1, i, h2, h1, i, integral) - int_direct = -1.d0 * integral - - call give_integrals_3_body_bi_ort(p1, p2, i, h2, h1, i, integral) - int_exc_23= -1.d0 * integral - - hthree += 1.d0 * int_direct - int_exc_23 - enddo - - return -end +!subroutine give_aab_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree) +! +! use bitmasks ! you need to include the bitmasks_module.f90 features +! +! implicit none +! integer, intent(in) :: Nint, h1, h2, p1, p2 +! integer, intent(in) :: Ne(2), occ(Nint*bit_kind_size,2) +! double precision, intent(out) :: hthree +! integer :: ii, i +! double precision :: int_direct, int_exc_12, int_exc_13, int_exc_23 +! double precision :: integral, int_exc_l, int_exc_ll +! +! hthree = 0.d0 +! do ii = 1, Ne(2) ! purely closed shell part +! i = occ(ii,2) +! +! call give_integrals_3_body_bi_ort(p2, p1, i, h2, h1, i, integral) +! int_direct = -1.d0 * integral +! +! call give_integrals_3_body_bi_ort(p1, p2, i, h2, h1, i, integral) +! int_exc_23= -1.d0 * integral +! +! hthree += 1.d0 * int_direct - int_exc_23 +! enddo +! +! return +!end ! --- @@ -264,6 +263,10 @@ BEGIN_PROVIDER [ double precision, no_aba_contraction, (mo_num,mo_num,mo_num,mo_ allocate(tmpvec_1(n_points_final_grid,3)) allocate(tmpvec_2(n_points_final_grid,3)) + double precision, allocatable :: tmp_2d(:,:) + allocate(tmp_2d(mo_num,mo_num)) + + ! purely closed shell part do ii = 1, Ne(2) i = occ(ii,2) @@ -313,9 +316,10 @@ BEGIN_PROVIDER [ double precision, no_aba_contraction, (mo_num,mo_num,mo_num,mo_ !$OMP END DO !$OMP END PARALLEL - call dgemm( 'T', 'N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 1.d0 & - , int2_grad1_u12_bimo_t, 3*n_points_final_grid, tmp1, 3*n_points_final_grid & - , 0.d0, tmp_3d, mo_num) + call dgemm( 'T', 'N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 1.d0 & + , int2_grad1_u12_bimo_t(1,1,1,1), 3*n_points_final_grid & + , tmp1(1,1,1), 3*n_points_final_grid & + , 0.d0, tmp_3d(1,1,1), mo_num*mo_num) !$OMP PARALLEL DO PRIVATE(p1,h2,p2) do p1 = 1, mo_num @@ -364,38 +368,163 @@ BEGIN_PROVIDER [ double precision, no_aba_contraction, (mo_num,mo_num,mo_num,mo_ !$OMP END DO !$OMP END PARALLEL - call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 1.d0 & - , mos_l_in_r_array_transp, n_points_final_grid, tmp2, n_points_final_grid & - , 1.d0, no_aba_contraction(p2,h2,1,1), mo_num*mo_num) + call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 1.d0 & + , mos_l_in_r_array_transp(1,1), n_points_final_grid & + , tmp2(1,1), n_points_final_grid & + , 0.d0, tmp_2d(1,1), mo_num) + + !$OMP PARALLEL DO PRIVATE(h2,p2) + do h2 = 1, mo_num + do p2 = 1, mo_num + no_aba_contraction(p2,h2,p1,h1) = no_aba_contraction(p2,h2,p1,h1) + tmp_2d(p2,h2) + enddo + enddo + !$OMP END PARALLEL DO enddo ! p1 enddo ! h1 enddo ! i - double precision :: integral, int_direct, int_exc_13, int_exc_12 - ! TODO + + + + + + ! purely open-shell part if(Ne(2) < Ne(1)) then - do ii = Ne(2) + 1, Ne(1) i = occ(ii,1) - call give_integrals_3_body_bi_ort(i, p2, p1, i, h2, h1, integral) - int_direct = -1.d0 * integral + do h1 = 1, mo_num - call give_integrals_3_body_bi_ort(p1, p2, i, i, h2, h1, integral) - int_exc_13 = -1.d0 * integral + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint) & + !$OMP SHARED (n_points_final_grid, i, h1, & + !$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 tmpval_1, tmpval_2, tmpvec_1, tmpvec_2) + !$OMP DO + do ipoint = 1, n_points_final_grid + tmpval_1(ipoint) = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint, i) + tmpval_2(ipoint) = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,h1) + tmpvec_1(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,i, i) * mos_r_in_r_array_transp(ipoint,h1) + tmpvec_1(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,i, i) * mos_r_in_r_array_transp(ipoint,h1) + tmpvec_1(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,i, i) * mos_r_in_r_array_transp(ipoint,h1) + tmpvec_2(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,i,h1) * mos_r_in_r_array_transp(ipoint, i) + tmpvec_2(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,i,h1) * mos_r_in_r_array_transp(ipoint, i) + tmpvec_2(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,i,h1) * mos_r_in_r_array_transp(ipoint, i) + enddo + !$OMP END DO + !$OMP END PARALLEL - call give_integrals_3_body_bi_ort(p2, i, p1, i, h2, h1, integral) - int_exc_12 = -1.d0 * integral + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (p1, ipoint) & + !$OMP SHARED (mo_num, n_points_final_grid, h1, i, & + !$OMP mos_l_in_r_array_transp, int2_grad1_u12_bimo_t, & + !$OMP tmpval_1, tmpval_2, tmpvec_1, tmpvec_2, tmp1) + !$OMP DO + do p1 = 1, mo_num + do ipoint = 1, n_points_final_grid + tmp1(ipoint,1,p1) = mos_l_in_r_array_transp(ipoint,p1) * (tmpvec_1(ipoint,1) - tmpvec_2(ipoint,1)) & + + tmpval_1(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p1,h1) - tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p1,i) + tmp1(ipoint,2,p1) = mos_l_in_r_array_transp(ipoint,p1) * (tmpvec_1(ipoint,2) - tmpvec_2(ipoint,2)) & + + tmpval_1(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p1,h1) - tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p1,i) + tmp1(ipoint,3,p1) = mos_l_in_r_array_transp(ipoint,p1) * (tmpvec_1(ipoint,3) - tmpvec_2(ipoint,3)) & + + tmpval_1(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p1,h1) - tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p1,i) + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL - no_aba_contraction(p2,h2,p1,h1) += 1.d0 * int_direct - 0.5d0 * (int_exc_13 + int_exc_12) - enddo + call dgemm( 'T', 'N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 0.5d0 & + , int2_grad1_u12_bimo_t(1,1,1,1), 3*n_points_final_grid & + , tmp1(1,1,1), 3*n_points_final_grid & + , 0.d0, tmp_3d(1,1,1), mo_num*mo_num) + + !$OMP PARALLEL DO PRIVATE(p1,h2,p2) + do p1 = 1, mo_num + do h2 = 1, mo_num + do p2 = 1, mo_num + no_aba_contraction(p2,h2,p1,h1) = no_aba_contraction(p2,h2,p1,h1) + tmp_3d(p2,h2,p1) + enddo + enddo + enddo + !$OMP END PARALLEL DO + + do p1 = 1, mo_num + + ! to minimize the number of operations + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint) & + !$OMP SHARED (n_points_final_grid, i, h1, p1, & + !$OMP int2_grad1_u12_bimo_t, final_weight_at_r_vector, & + !$OMP tmpval_1) + !$OMP DO + do ipoint = 1, n_points_final_grid + tmpval_1(ipoint) = final_weight_at_r_vector(ipoint) * ( int2_grad1_u12_bimo_t(ipoint,1, i,i) * int2_grad1_u12_bimo_t(ipoint,1,p1,h1) & + + int2_grad1_u12_bimo_t(ipoint,2, i,i) * int2_grad1_u12_bimo_t(ipoint,2,p1,h1) & + + int2_grad1_u12_bimo_t(ipoint,3, i,i) * int2_grad1_u12_bimo_t(ipoint,3,p1,h1) & + - int2_grad1_u12_bimo_t(ipoint,1,p1,i) * int2_grad1_u12_bimo_t(ipoint,1, i,h1) & + - int2_grad1_u12_bimo_t(ipoint,2,p1,i) * int2_grad1_u12_bimo_t(ipoint,2, i,h1) & + - int2_grad1_u12_bimo_t(ipoint,3,p1,i) * int2_grad1_u12_bimo_t(ipoint,3, i,h1) ) + enddo + !$OMP END DO + !$OMP END PARALLEL + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (h2, ipoint) & + !$OMP SHARED (mo_num, n_points_final_grid, & + !$OMP mos_r_in_r_array_transp, & + !$OMP tmpval_1, tmp2) + !$OMP DO + do h2 = 1, mo_num + do ipoint = 1, n_points_final_grid + tmp2(ipoint,h2) = mos_r_in_r_array_transp(ipoint,h2) * tmpval_1(ipoint) + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 0.5d0 & + , mos_l_in_r_array_transp(1,1), n_points_final_grid & + , tmp2(1,1), n_points_final_grid & + , 0.d0, tmp_2d(1,1), mo_num) + + !$OMP PARALLEL DO PRIVATE(h2,p2) + do h2 = 1, mo_num + do p2 = 1, mo_num + no_aba_contraction(p2,h2,p1,h1) = no_aba_contraction(p2,h2,p1,h1) + tmp_2d(p2,h2) + enddo + enddo + !$OMP END PARALLEL DO + + enddo ! p1 + enddo ! h1 + enddo !i endif - ! --- + + + + + + + + + + + + + + + deallocate(tmp_3d) deallocate(tmp1, tmp2) @@ -403,17 +532,121 @@ BEGIN_PROVIDER [ double precision, no_aba_contraction, (mo_num,mo_num,mo_num,mo_ deallocate(tmpvec_1, tmpvec_2) - !$OMP PARALLEL DO PRIVATE(h1,h2,p1,p2) - do h1 = 1, mo_num - do p1 = 1, mo_num - do h2 = 1, mo_num - do p2 = 1, mo_num - no_aba_contraction(p2,h2,p1,h1) = -0.5d0 * (no_aba_contraction(p2,h2,p1,h1) + no_aba_contraction(p1,h1,p2,h2)) - enddo - enddo - enddo - enddo - !$OMP END PARALLEL DO + + + + + + + + no_aba_contraction = -0.5d0 * no_aba_contraction + call sum_A_At(no_aba_contraction(1,1,1,1), mo_num*mo_num) + +! do h1 = 1, mo_num +! do p1 = 1, mo_num +! do h2 = 1, mo_num +! do p2 = 1, mo_num +! no_aba_contraction(p2,h2,p1,h1) = -0.5d0 * (tmp_4d(p2,h2,p1,h1) + tmp_4d(p1,h1,p2,h2)) +! enddo +! enddo +! enddo +! enddo + + + ! --- + + double precision :: integral, int_direct, int_exc_13, int_exc_12 + +! no_aba_contraction = 0.d0 +! +! ! purely closed shell part +! do ii = 1, Ne(2) +! i = occ(ii,1) +! +! !$OMP PARALLEL & +! !$OMP DEFAULT (NONE) & +! !$OMP PRIVATE (h1, h2, p1, p2, int_direct, int_exc_13, int_exc_12, integral) & +! !$OMP SHARED (mo_num, i, no_aba_contraction) +! !$OMP DO SCHEDULE (static) +! do h1 = 1, mo_num +! do p1 = 1, mo_num +! do h2 = 1, mo_num +! do p2 = 1, mo_num +! +! call give_integrals_3_body_bi_ort(i, p2, p1, i, h2, h1, integral) +! int_direct = -1.d0 * integral +! +! call give_integrals_3_body_bi_ort(p1, p2, i, i, h2, h1, integral) +! int_exc_13 = -1.d0 * integral +! +! call give_integrals_3_body_bi_ort(p2, i, p1, i, h2, h1, integral) +! int_exc_12 = -1.d0 * integral +! +! !no_aba_contraction(p2,h2,p1,h1) += 1.d0 * int_direct - 0.5d0 * (int_exc_13 + int_exc_12) +! enddo +! enddo +! enddo +! enddo +! !$OMP END DO +! !$OMP END PARALLEL +! enddo + +! ! purely open-shell part +! if(Ne(2) < Ne(1)) then +! +! do ii = Ne(2) + 1, Ne(1) +! i = occ(ii,1) +! +! !$OMP PARALLEL & +! !$OMP DEFAULT (NONE) & +! !$OMP PRIVATE (h1, h2, p1, p2, int_direct, int_exc_13, int_exc_12, integral) & +! !$OMP SHARED (mo_num, i, no_aba_contraction) +! !$OMP DO SCHEDULE (static) +! do h1 = 1, mo_num +! do p1 = 1, mo_num +! do h2 = 1, mo_num +! do p2 = 1, mo_num +! +! call give_integrals_3_body_bi_ort(i, p2, p1, i, h2, h1, integral) +! int_direct = -1.d0 * integral +! +! call give_integrals_3_body_bi_ort(p1, p2, i, i, h2, h1, integral) +! int_exc_13 = -1.d0 * integral +! +! call give_integrals_3_body_bi_ort(p2, i, p1, i, h2, h1, integral) +! int_exc_12 = -1.d0 * integral +! +! no_aba_contraction(p2,h2,p1,h1) += 0.5d0 * int_direct - 0.25d0 * (int_exc_13 + int_exc_12) +! enddo +! enddo +! enddo +! enddo +! !$OMP END DO +! !$OMP END PARALLEL +! enddo +! endif + + ! --- + +! !$OMP PARALLEL & +! !$OMP DEFAULT (NONE) & +! !$OMP PRIVATE (h1, h2, p1, p2, integral) & +! !$OMP SHARED (mo_num, N_int,Ne, occ, no_aba_contraction) +! !$OMP DO SCHEDULE (static) +! do h1 = 1, mo_num +! do p1 = 1, mo_num +! do h2 = 1, mo_num +! do p2 = 1, mo_num +! call give_aba_contraction(N_int, h1, h2, p1, p2, Ne, occ, integral) +! no_aba_contraction(p2,h2,p1,h1) = 0.5d0 * integral +! enddo +! enddo +! enddo +! enddo +! !$OMP END DO +! !$OMP END PARALLEL + + END_PROVIDER diff --git a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f index df86ea65..33b5c5aa 100644 --- a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f +++ b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f @@ -11,12 +11,14 @@ program tc_bi_ortho touch read_wf touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid - call test_h_u0 +! call test_h_u0 ! call test_slater_tc_opt ! call timing_tot ! call timing_diag ! call timing_single ! call timing_double + + call test_no() end subroutine test_h_u0 @@ -252,3 +254,47 @@ subroutine timing_double end +! --- + +subroutine test_no() + + implicit none + integer :: i, j, k, l + double precision :: accu, contrib, new, ref, thr + + print*, ' testing normal_two_body_bi_orth ...' + + thr = 1d-8 + + PROVIDE normal_two_body_bi_orth_old + PROVIDE normal_two_body_bi_orth + + accu = 0.d0 + do i = 1, mo_num + do j = 1, mo_num + do k = 1, mo_num + do l = 1, mo_num + + new = normal_two_body_bi_orth (l,k,j,i) + ref = normal_two_body_bi_orth_old(l,k,j,i) + contrib = dabs(new - ref) + accu += contrib + if(contrib .gt. thr) then + print*, ' problem on normal_two_body_bi_orth' + print*, l, k, j, i + print*, ref, new, contrib + stop + endif + + enddo + enddo + enddo + enddo + print*, ' accu on normal_two_body_bi_orth = ', accu / dble(mo_num)**4 + + return +end + +! --- + +