diff --git a/src/becke_numerical_grid/grid_becke_vector.irp.f b/src/becke_numerical_grid/grid_becke_vector.irp.f index e1147ec7..0386f3c6 100644 --- a/src/becke_numerical_grid/grid_becke_vector.irp.f +++ b/src/becke_numerical_grid/grid_becke_vector.irp.f @@ -20,8 +20,8 @@ BEGIN_PROVIDER [integer, n_points_final_grid] enddo enddo - !print*,' n_points_final_grid = ', n_points_final_grid - !print*,' n max point = ', n_points_integration_angular*(n_points_radial_grid*nucl_num - 1) + print*,' n_points_final_grid = ', n_points_final_grid + print*,' n max point = ', n_points_integration_angular*(n_points_radial_grid*nucl_num - 1) call ezfio_set_becke_numerical_grid_n_points_final_grid(n_points_final_grid) END_PROVIDER diff --git a/src/tc_bi_ortho/normal_ordered.irp.f b/src/tc_bi_ortho/normal_ordered.irp.f index 1a9cddda..e65df450 100644 --- a/src/tc_bi_ortho/normal_ordered.irp.f +++ b/src/tc_bi_ortho/normal_ordered.irp.f @@ -13,13 +13,18 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_ implicit none - integer :: i, ii, h1, p1, h2, p2, ipoint + integer :: i, ii, ipoint + integer :: h1, p1, h2, p2 integer :: hh1, hh2, pp1, pp2 integer :: Ne(2) double precision :: wall0, wall1, walli, wallf integer, allocatable :: occ(:,:) integer(bit_kind), allocatable :: key_i_core(:,:) + PROVIDE mo_class + PROVIDE list_act n_act_orb + PROVIDE N_int + print*,' Providing normal_two_body_bi_orth ...' call wall_time(walli) @@ -31,8 +36,8 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_ else - double precision, allocatable :: tmp_2d(:,:), tmp_3d(:,:,:) - double precision, allocatable :: tmp1(:,:,:), tmp2(:,:), tmp3(:,:,:) + double precision, allocatable :: tmp_3d(:,:,:) + double precision, allocatable :: tmp1(:,:,:), tmp2(:,:) double precision, allocatable :: tmpval_1(:), tmpval_2(:), tmpvec_1(:,:), tmpvec_2(:,:), tmpvec_3(:,:) double precision, allocatable :: tmp(:,:,:,:) double precision, allocatable :: int2_grad1_u12_bimo_t_tmp(:,:,:,:), mos_l_in_r_array_transp_tmp(:,:), mos_r_in_r_array_transp_tmp(:,:) @@ -72,10 +77,6 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_ !$OMP END DO !$OMP END PARALLEL - - PROVIDE mo_class - PROVIDE N_int - allocate( occ(N_int*bit_kind_size,2) ) allocate( key_i_core(N_int,2) ) @@ -97,120 +98,98 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_ print*,' Providing aba_contraction ...' call wall_time(wall0) + tmp = 0.d0 + call set_multiple_levels_omp(.false.) !$OMP PARALLEL & !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, h1, p1, h2, p2, i, ii, & - !$OMP tmp_3d, tmp_2d, tmp1, tmp2, & + !$OMP PRIVATE (ipoint, hh1, pp1, hh2, pp2, h1, p1, h2, p2, i, ii, tmp1, tmp2, & !$OMP tmpval_1, tmpval_2, tmpvec_1, tmpvec_2) & - !$OMP SHARED (n_points_final_grid, Ne, occ, mo_num, mo_class, & + !$OMP SHARED (n_points_final_grid, Ne, occ, mo_num, n_act_orb, list_act, & !$OMP mos_l_in_r_array_transp_tmp, mos_r_in_r_array_transp_tmp, & - !$OMP int2_grad1_u12_bimo_t_tmp, final_weight_at_r_vector, & + !$OMP final_weight_at_r_vector, int2_grad1_u12_bimo_t_tmp, & !$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, int2_grad1_u12_bimo_t, & !$OMP tmp) - allocate(tmp_3d(mo_num,mo_num,mo_num), tmp_2d(mo_num,mo_num)) allocate(tmp1(n_points_final_grid,3,mo_num), tmp2(n_points_final_grid,mo_num)) allocate(tmpval_1(n_points_final_grid), tmpval_2(n_points_final_grid)) allocate(tmpvec_1(n_points_final_grid,3), tmpvec_2(n_points_final_grid,3)) - tmp_3d = 0.d0 - tmp_2d = 0.d0 - tmp1 = 0.d0 - tmp2 = 0.d0 - tmpval_1 = 0.d0 - tmpval_2 = 0.d0 - tmpvec_1 = 0.d0 - tmpvec_2 = 0.d0 - !$OMP DO + do hh1 = 1, n_act_orb + h1 = list_act(hh1) - do h1 = 1, mo_num - tmp(:,:,:,h1) = 0.d0 - if(mo_class(h1) .ne. "Active") cycle - + tmp1 = 0.d0 do ii = 1, Ne(2) i = occ(ii,2) 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_tmp(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_tmp(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_tmp(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_tmp(ipoint,h1) + 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 - do p1 = 1, mo_num - tmp1(:,:,p1) = 0.d0 - if(mo_class(p1) .ne. "Active") cycle - + do pp1 = 1, n_act_orb + p1 = list_act(pp1) do ipoint = 1, n_points_final_grid - tmp1(ipoint,1,p1) = mos_l_in_r_array_transp_tmp(ipoint,p1) * (tmpvec_1(ipoint,1) - tmpvec_2(ipoint,1)) & - + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(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_tmp(ipoint,p1) * (tmpvec_1(ipoint,2) - tmpvec_2(ipoint,2)) & - + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(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_tmp(ipoint,p1) * (tmpvec_1(ipoint,3) - tmpvec_2(ipoint,3)) & - + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(ipoint,3,p1,h1) - tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p1,i) + tmp1(ipoint,1,p1) = 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) = 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) = 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 - call dgemm( 'T', 'N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 1.d0 & - , int2_grad1_u12_bimo_t_tmp(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) + enddo ! ii - do p1 = 1, mo_num - do h2 = 1, mo_num - do p2 = 1, mo_num - tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_3d(p2,h2,p1) - enddo - enddo - enddo + call dgemm( 'T', 'N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 1.d0 & + , int2_grad1_u12_bimo_t_tmp(1,1,1,1), 3*n_points_final_grid & + , tmp1(1,1,1), 3*n_points_final_grid & + , 1.d0, tmp(1,1,1,h1), mo_num*mo_num) - do p1 = 1, mo_num - if(mo_class(p1) .ne. "Active") cycle + do pp1 = 1, n_act_orb + p1 = list_act(pp1) + + tmp2 = 0.d0 + do ii = 1, Ne(2) + i = occ(ii,2) 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_tmp(ipoint,1,p1,h1) & - + int2_grad1_u12_bimo_t(ipoint,2,i,i) * int2_grad1_u12_bimo_t_tmp(ipoint,2,p1,h1) & - + int2_grad1_u12_bimo_t(ipoint,3,i,i) * int2_grad1_u12_bimo_t_tmp(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) & + 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 - do h2 = 1, mo_num - tmp2(:,h2) = 0.d0 - if(mo_class(h2) .ne. "Active") cycle - + do hh2 = 1, n_act_orb + h2 = list_act(hh2) do ipoint = 1, n_points_final_grid - tmp2(ipoint,h2) = mos_r_in_r_array_transp_tmp(ipoint,h2) * tmpval_1(ipoint) + tmp2(ipoint,h2) = tmp2(ipoint,h2) + mos_r_in_r_array_transp(ipoint,h2) * tmpval_1(ipoint) enddo enddo - call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 1.d0 & - , mos_l_in_r_array_transp_tmp(1,1), n_points_final_grid & - , tmp2(1,1), n_points_final_grid & - , 0.d0, tmp_2d(1,1), mo_num) + enddo ! ii - do h2 = 1, mo_num - do p2 = 1, mo_num - tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_2d(p2,h2) - enddo - enddo + call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 1.d0 & + , mos_l_in_r_array_transp_tmp(1,1), n_points_final_grid & + , tmp2(1,1), n_points_final_grid & + , 1.d0, tmp(1,1,p1,h1), mo_num) - enddo ! p1 - enddo ! i + enddo ! p1 enddo ! h1 !$OMP END DO - deallocate(tmp_3d, tmp_2d) deallocate(tmp1, tmp2) deallocate(tmpval_1, tmpval_2) deallocate(tmpvec_1, tmpvec_2) @@ -225,114 +204,93 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_ !$OMP PARALLEL & !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, h1, p1, h2, p2, i, ii, & - !$OMP tmp_3d, tmp_2d, tmp1, tmp2, & + !$OMP PRIVATE (ipoint, hh1, pp1, hh2, pp2, h1, p1, h2, p2, i, ii, tmp1, tmp2, & !$OMP tmpval_1, tmpval_2, tmpvec_1, tmpvec_2) & - !$OMP SHARED (n_points_final_grid, Ne, occ, mo_num, mo_class, & + !$OMP SHARED (n_points_final_grid, Ne, occ, mo_num, n_act_orb, list_act, & !$OMP mos_l_in_r_array_transp_tmp, mos_r_in_r_array_transp_tmp, & - !$OMP int2_grad1_u12_bimo_t_tmp, final_weight_at_r_vector, & + !$OMP final_weight_at_r_vector, int2_grad1_u12_bimo_t_tmp, & !$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, int2_grad1_u12_bimo_t, & !$OMP tmp) - Allocate(tmp_3d(mo_num,mo_num,mo_num), tmp_2d(mo_num,mo_num)) Allocate(tmp1(n_points_final_grid,3,mo_num), tmp2(n_points_final_grid,mo_num)) Allocate(tmpval_1(n_points_final_grid), tmpval_2(n_points_final_grid)) Allocate(tmpvec_1(n_points_final_grid,3), tmpvec_2(n_points_final_grid,3)) - Tmp_3d = 0.d0 - Tmp_2d = 0.d0 - Tmp1 = 0.d0 - Tmp2 = 0.d0 - Tmpval_1 = 0.d0 - Tmpval_2 = 0.d0 - Tmpvec_1 = 0.d0 - Tmpvec_2 = 0.d0 - !$OMP DO - do h1 = 1, mo_num - if(mo_class(h1) .ne. "Active") cycle + do hh1 = 1, n_act_orb + h1 = list_act(hh1) + tmp1 = 0.d0 do ii = Ne(2) + 1, Ne(1) i = occ(ii,1) 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_tmp(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_tmp(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_tmp(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_tmp(ipoint,h1) + 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 - do p1 = 1, mo_num - tmp1(:,:,p1) = 0.d0 - if(mo_class(p1) .ne. "Active") cycle + do pp1 = 1, n_act_orb + p1 = list_act(pp1) do ipoint = 1, n_points_final_grid - tmp1(ipoint,1,p1) = mos_l_in_r_array_transp_tmp(ipoint,p1) * (tmpvec_1(ipoint,1) - tmpvec_2(ipoint,1)) & - + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(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_tmp(ipoint,p1) * (tmpvec_1(ipoint,2) - tmpvec_2(ipoint,2)) & - + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(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_tmp(ipoint,p1) * (tmpvec_1(ipoint,3) - tmpvec_2(ipoint,3)) & - + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(ipoint,3,p1,h1) - tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p1,i) + tmp1(ipoint,1,p1) = 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) = 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) = 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 - call dgemm( 'T', 'N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 0.5d0 & - , int2_grad1_u12_bimo_t_tmp(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) + enddo ! ii - do p1 = 1, mo_num - do h2 = 1, mo_num - do p2 = 1, mo_num - tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_3d(p2,h2,p1) - enddo - enddo - enddo + call dgemm( 'T', 'N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 0.5d0 & + , int2_grad1_u12_bimo_t_tmp(1,1,1,1), 3*n_points_final_grid & + , tmp1(1,1,1), 3*n_points_final_grid & + , 1.d0, tmp(1,1,1,h1), mo_num*mo_num) - do p1 = 1, mo_num - if(mo_class(p1) .ne. "Active") cycle + do pp1 = 1, n_act_orb + p1 = list_act(pp1) + + tmp2 = 0.d0 + do ii = Ne(2) + 1, Ne(1) + i = occ(ii,1) 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_tmp(ipoint,1,p1,h1) & - + int2_grad1_u12_bimo_t(ipoint,2,i,i) * int2_grad1_u12_bimo_t_tmp(ipoint,2,p1,h1) & - + int2_grad1_u12_bimo_t(ipoint,3,i,i) * int2_grad1_u12_bimo_t_tmp(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) & + 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 - do h2 = 1, mo_num - tmp2(:,h2) = 0.d0 - if(mo_class(h2) .ne. "Active") cycle - + do hh2 = 1, n_act_orb + h2 = list_act(hh2) do ipoint = 1, n_points_final_grid - tmp2(ipoint,h2) = mos_r_in_r_array_transp_tmp(ipoint,h2) * tmpval_1(ipoint) + tmp2(ipoint,h2) = tmp2(ipoint,h2) + mos_r_in_r_array_transp(ipoint,h2) * tmpval_1(ipoint) enddo enddo - call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 0.5d0 & - , mos_l_in_r_array_transp_tmp(1,1), n_points_final_grid & - , tmp2(1,1), n_points_final_grid & - , 0.d0, tmp_2d(1,1), mo_num) + enddo ! ii - do h2 = 1, mo_num - do p2 = 1, mo_num - tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_2d(p2,h2) - enddo - enddo + call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 0.5d0 & + , mos_l_in_r_array_transp_tmp(1,1), n_points_final_grid & + , tmp2(1,1), n_points_final_grid & + , 1.d0, tmp(1,1,p1,h1), mo_num) - enddo ! p1 - enddo ! i + enddo ! p1 enddo ! h1 !$OMP END DO - deallocate(tmp_3d, tmp_2d) deallocate(tmp1, tmp2) deallocate(tmpval_1, tmpval_2) deallocate(tmpvec_1, tmpvec_2) @@ -354,109 +312,89 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_ print*,' Providing aab_contraction ...' call wall_time(wall0) + tmp = 0.d0 + call set_multiple_levels_omp(.false.) !$OMP PARALLEL & !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, ii, i, h1, p1, h2, p2, & - !$OMP tmp_2d, tmp_3d, tmp1, tmp2, & + !$OMP PRIVATE (ipoint, hh1, pp1, hh2, pp2, ii, i, h1, p1, h2, p2, tmp1, tmp2, & !$OMP tmpval_1, tmpvec_1) & - !$OMP SHARED (n_points_final_grid, mo_num, Ne, occ, mo_class, & + !$OMP SHARED (n_points_final_grid, mo_num, Ne, occ, n_act_orb, list_act, & !$OMP mos_l_in_r_array_transp_tmp, mos_r_in_r_array_transp_tmp, & - !$OMP int2_grad1_u12_bimo_t_tmp, final_weight_at_r_vector, & + !$OMP final_weight_at_r_vector, int2_grad1_u12_bimo_t_tmp, & !$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, int2_grad1_u12_bimo_t, & !$OMP tmp) - allocate(tmp_2d(mo_num,mo_num)) - allocate(tmp_3d(mo_num,mo_num,mo_num)) allocate(tmp1(n_points_final_grid,3,mo_num)) allocate(tmp2(n_points_final_grid,mo_num)) allocate(tmpval_1(n_points_final_grid)) allocate(tmpvec_1(n_points_final_grid,3)) - tmp_2d = 0.d0 - tmp_3d = 0.d0 - tmp1 = 0.d0 - tmp2 = 0.d0 - tmpval_1 = 0.d0 - tmpvec_1 = 0.d0 - !$OMP DO - do h1 = 1, mo_num - tmp(:,:,:,h1) = 0.d0 - if(mo_class(h1) .ne. "Active") cycle + do hh1 = 1, n_act_orb + h1 = list_act(hh1) + tmp1 = 0.d0 do ii = 1, Ne(2) i = occ(ii,2) 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) - 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_tmp(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_tmp(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_tmp(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) enddo - do p1 = 1, mo_num - tmp1(:,:,p1) = 0.d0 - if(mo_class(p1) .ne. "Active") cycle + do pp1 = 1, n_act_orb + p1 = list_act(pp1) + do ipoint = 1, n_points_final_grid + tmp1(ipoint,1,p1) = tmp1(ipoint,1,p1) + mos_l_in_r_array_transp(ipoint,p1) * tmpvec_1(ipoint,1) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p1,h1) + tmp1(ipoint,2,p1) = tmp1(ipoint,2,p1) + mos_l_in_r_array_transp(ipoint,p1) * tmpvec_1(ipoint,2) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p1,h1) + tmp1(ipoint,3,p1) = tmp1(ipoint,3,p1) + mos_l_in_r_array_transp(ipoint,p1) * tmpvec_1(ipoint,3) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p1,h1) + enddo + enddo + + enddo ! ii + + call dgemm( 'T', 'N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 1.d0 & + , int2_grad1_u12_bimo_t_tmp(1,1,1,1), 3*n_points_final_grid & + , tmp1(1,1,1), 3*n_points_final_grid & + , 1.d0, tmp(1,1,1,h1), mo_num*mo_num) + + do pp1 = 1, n_act_orb + p1 = list_act(pp1) + + tmp2 = 0.d0 + do ii = 1, Ne(2) + i = occ(ii,2) do ipoint = 1, n_points_final_grid - tmp1(ipoint,1,p1) = mos_l_in_r_array_transp_tmp(ipoint,p1) * tmpvec_1(ipoint,1) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(ipoint,1,p1,h1) - tmp1(ipoint,2,p1) = mos_l_in_r_array_transp_tmp(ipoint,p1) * tmpvec_1(ipoint,2) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(ipoint,2,p1,h1) - tmp1(ipoint,3,p1) = mos_l_in_r_array_transp_tmp(ipoint,p1) * tmpvec_1(ipoint,3) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(ipoint,3,p1,h1) - enddo - enddo - - call dgemm( 'T', 'N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 1.d0 & - , int2_grad1_u12_bimo_t_tmp(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) - - do p1 = 1, mo_num - do h2 = 1, mo_num - do p2 = 1, mo_num - tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_3d(p2,h2,p1) - enddo - enddo - enddo - - do p1 = 1, mo_num - if(mo_class(p1) .ne. "Active") cycle - - 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_tmp(ipoint,1,p1,h1) & - + int2_grad1_u12_bimo_t(ipoint,2,i,i) * int2_grad1_u12_bimo_t_tmp(ipoint,2,p1,h1) & - + int2_grad1_u12_bimo_t(ipoint,3,i,i) * int2_grad1_u12_bimo_t_tmp(ipoint,3,p1,h1) ) + 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) ) enddo - do h2 = 1, mo_num - if(mo_class(h2) .ne. "Active") cycle - tmp2(:,h2) = 0.d0 - + do hh2 = 1, n_act_orb + h2 = list_act(hh2) do ipoint = 1, n_points_final_grid - tmp2(ipoint,h2) = mos_r_in_r_array_transp_tmp(ipoint,h2) * tmpval_1(ipoint) + tmp2(ipoint,h2) = tmp2(ipoint,h2) + mos_r_in_r_array_transp(ipoint,h2) * tmpval_1(ipoint) enddo enddo - call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 1.d0 & - , mos_l_in_r_array_transp_tmp(1,1), n_points_final_grid & - , tmp2(1,1), n_points_final_grid & - , 0.d0, tmp_2d(1,1), mo_num) + enddo ! ii - do h2 = 1, mo_num - do p2 = 1, mo_num - tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_2d(p2,h2) - enddo - enddo + call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 1.d0 & + , mos_l_in_r_array_transp_tmp(1,1), n_points_final_grid & + , tmp2(1,1), n_points_final_grid & + , 1.d0, tmp(1,1,p1,h1), mo_num) - enddo ! p1 - enddo ! i + enddo ! p1 enddo ! h1 !$OMP END DO - deallocate(tmp_3d) deallocate(tmp1, tmp2) deallocate(tmpval_1) deallocate(tmpvec_1) @@ -520,212 +458,189 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_ print*,' Providing aaa_contraction ...' call wall_time(wall0) + tmp = 0.d0 + call set_multiple_levels_omp(.false.) !$OMP PARALLEL & !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, i, ii, h1, h2, p1, p2, & - !$OMP tmp_2d, tmp_3d, tmp1, tmp2, tmp3, & + !$OMP PRIVATE (ipoint, i, ii, hh1, hh2, pp1, pp2, h1, h2, p1, p2, tmp1, tmp2, tmp_3d, & !$OMP tmpval_1, tmpval_2, & !$OMP tmpvec_1, tmpvec_2, tmpvec_3) & - !$OMP SHARED (n_points_final_grid, Ne, occ, mo_num, mo_class, & + !$OMP SHARED (n_points_final_grid, Ne, occ, mo_num, n_act_orb, list_act, & !$OMP mos_l_in_r_array_transp_tmp, mos_r_in_r_array_transp_tmp, & - !$OMP int2_grad1_u12_bimo_t_tmp, final_weight_at_r_vector, & + !$OMP final_weight_at_r_vector, int2_grad1_u12_bimo_t_tmp, & !$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, int2_grad1_u12_bimo_t, & !$OMP tmp) - allocate(tmp_2d(mo_num,mo_num)) allocate(tmp_3d(mo_num,mo_num,mo_num)) allocate(tmp1(n_points_final_grid,3,mo_num)) allocate(tmp2(n_points_final_grid,mo_num)) - allocate(tmp3(n_points_final_grid,3,mo_num)) allocate(tmpval_1(n_points_final_grid)) allocate(tmpval_2(n_points_final_grid)) allocate(tmpvec_1(n_points_final_grid,3)) allocate(tmpvec_2(n_points_final_grid,3)) allocate(tmpvec_3(n_points_final_grid,3)) - tmp_2d = 0.d0 - tmp_3d = 0.d0 - tmp1 = 0.d0 - tmp2 = 0.d0 - tmp3 = 0.d0 - tmpval_1 = 0.d0 - tmpval_2 = 0.d0 - tmpvec_1 = 0.d0 - tmpvec_2 = 0.d0 - tmpvec_3 = 0.d0 - !$OMP DO - do h1 = 1, mo_num - tmp(:,:,:,h1) = 0.d0 - if(mo_class(h1) .ne. "Active") cycle + do hh1 = 1, n_act_orb + h1 = list_act(hh1) + tmp1 = 0.d0 + do ii = 1, Ne(2) + i = occ(ii,2) + + 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) + + 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) + enddo + + do pp1 = 1, n_act_orb + p1 = list_act(pp1) + do ipoint = 1, n_points_final_grid + tmp1(ipoint,1,p1) = tmp1(ipoint,1,p1) + mos_l_in_r_array_transp(ipoint,p1) * tmpvec_1(ipoint,1) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p1,h1) + tmp1(ipoint,2,p1) = tmp1(ipoint,2,p1) + mos_l_in_r_array_transp(ipoint,p1) * tmpvec_1(ipoint,2) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p1,h1) + tmp1(ipoint,3,p1) = tmp1(ipoint,3,p1) + mos_l_in_r_array_transp(ipoint,p1) * tmpvec_1(ipoint,3) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p1,h1) + enddo + enddo + + enddo ! ii + + call dgemm( 'T', 'N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 1.d0 & + , int2_grad1_u12_bimo_t_tmp(1,1,1,1), 3*n_points_final_grid & + , tmp1(1,1,1), 3*n_points_final_grid & + , 1.d0, tmp(1,1,1,h1), mo_num*mo_num) + + tmp1 = 0.d0 do ii = 1, Ne(2) i = occ(ii,2) 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_tmp(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_tmp(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_tmp(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_tmp(ipoint,h1) + 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_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 - do p1 = 1, mo_num - tmp1(:,:,p1) = 0.d0 - if(mo_class(p1) .ne. "Active") cycle - + do pp2 = 1, n_act_orb + p2 = list_act(pp2) do ipoint = 1, n_points_final_grid - tmp1(ipoint,1,p1) = mos_l_in_r_array_transp_tmp(ipoint,p1) * tmpvec_1(ipoint,1) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(ipoint,1,p1,h1) - tmp1(ipoint,2,p1) = mos_l_in_r_array_transp_tmp(ipoint,p1) * tmpvec_1(ipoint,2) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(ipoint,2,p1,h1) - tmp1(ipoint,3,p1) = mos_l_in_r_array_transp_tmp(ipoint,p1) * tmpvec_1(ipoint,3) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(ipoint,3,p1,h1) + tmp1(ipoint,1,p2) = tmp1(ipoint,1,p2) + tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p2,i) + mos_l_in_r_array_transp(ipoint,p2) * tmpvec_2(ipoint,1) + tmp1(ipoint,2,p2) = tmp1(ipoint,2,p2) + tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p2,i) + mos_l_in_r_array_transp(ipoint,p2) * tmpvec_2(ipoint,2) + tmp1(ipoint,3,p2) = tmp1(ipoint,3,p2) + tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p2,i) + mos_l_in_r_array_transp(ipoint,p2) * tmpvec_2(ipoint,3) enddo enddo - call dgemm( 'T', 'N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 1.d0 & - , int2_grad1_u12_bimo_t_tmp(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) + enddo ! ii - do p1 = 1, mo_num - do h2 = 1, mo_num - do p2 = 1, mo_num - tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_3d(p2,h2,p1) + call dgemm( 'T', 'N', mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0 & + , tmp1(1,1,1), 3*n_points_final_grid & + , int2_grad1_u12_bimo_t_tmp(1,1,1,1), 3*n_points_final_grid & + , 0.d0, tmp_3d(1,1,1), mo_num) + + do p1 = 1, mo_num + do h2 = 1, mo_num + do p2 = 1, mo_num + tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_3d(p2,p1,h2) + enddo + enddo + enddo + + do pp1 = 1, n_act_orb + p1 = list_act(pp1) + + tmp2 = 0.d0 + do ii = 1, Ne(2) + i = occ(ii,2) + + 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) ) + + tmpvec_1(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p1,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,p1,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,p1,i) * mos_r_in_r_array_transp(ipoint,h1) + enddo + + do hh2 = 1, n_act_orb + h2 = list_act(hh2) + + do ipoint = 1, n_points_final_grid + + tmp2(ipoint,h2) = tmp2(ipoint,h2) + mos_r_in_r_array_transp(ipoint,h2) * tmpval_1(ipoint) & + + int2_grad1_u12_bimo_t(ipoint,1,i,h2) * tmpvec_1(ipoint,1) & + + int2_grad1_u12_bimo_t(ipoint,2,i,h2) * tmpvec_1(ipoint,2) & + + int2_grad1_u12_bimo_t(ipoint,3,i,h2) * tmpvec_1(ipoint,3) enddo enddo - enddo - do p2 = 1, mo_num - tmp1(:,:,p2) = 0.d0 - if(mo_class(p2) .ne. "Active") cycle + enddo ! ii - do ipoint = 1, n_points_final_grid - tmp1(ipoint,1,p2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p2,i) + mos_l_in_r_array_transp_tmp(ipoint,p2) * tmpvec_2(ipoint,1) - tmp1(ipoint,2,p2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p2,i) + mos_l_in_r_array_transp_tmp(ipoint,p2) * tmpvec_2(ipoint,2) - tmp1(ipoint,3,p2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p2,i) + mos_l_in_r_array_transp_tmp(ipoint,p2) * tmpvec_2(ipoint,3) - enddo - enddo + call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 1.d0 & + , mos_l_in_r_array_transp_tmp(1,1), n_points_final_grid & + , tmp2(1,1), n_points_final_grid & + , 1.d0, tmp(1,1,p1,h1), mo_num) - call dgemm( 'T', 'N', mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0 & - , tmp1(1,1,1), 3*n_points_final_grid & - , int2_grad1_u12_bimo_t_tmp(1,1,1,1), 3*n_points_final_grid & - , 0.d0, tmp_3d(1,1,1), mo_num) - - do p1 = 1, mo_num - do h2 = 1, mo_num - do p2 = 1, mo_num - tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_3d(p2,p1,h2) - enddo - enddo - enddo - - do p1 = 1, mo_num - if(mo_class(p1) .ne. "Active") cycle + tmp1 = 0.d0 + tmp2 = 0.d0 + do ii = 1, Ne(2) + i = occ(ii,2) 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_tmp(ipoint,1,p1,h1) & - + int2_grad1_u12_bimo_t(ipoint,2,i,i) * int2_grad1_u12_bimo_t_tmp(ipoint,2,p1,h1) & - + int2_grad1_u12_bimo_t(ipoint,3,i,i) * int2_grad1_u12_bimo_t_tmp(ipoint,3,p1,h1) ) + tmpval_2(ipoint) = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,p1) * mos_r_in_r_array_transp(ipoint,i) - tmpval_2(ipoint) = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp_tmp(ipoint,p1) * mos_r_in_r_array_transp(ipoint,i) - - tmpvec_1(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p1,i) * mos_r_in_r_array_transp_tmp(ipoint,h1) - tmpvec_1(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p1,i) * mos_r_in_r_array_transp_tmp(ipoint,h1) - tmpvec_1(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p1,i) * mos_r_in_r_array_transp_tmp(ipoint,h1) - - tmpvec_2(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,i,h1) * mos_l_in_r_array_transp_tmp(ipoint,p1) - tmpvec_2(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,i,h1) * mos_l_in_r_array_transp_tmp(ipoint,p1) - tmpvec_2(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,i,h1) * mos_l_in_r_array_transp_tmp(ipoint,p1) + tmpvec_2(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,i,h1) * mos_l_in_r_array_transp(ipoint,p1) + tmpvec_2(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,i,h1) * mos_l_in_r_array_transp(ipoint,p1) + tmpvec_2(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,i,h1) * mos_l_in_r_array_transp(ipoint,p1) tmpvec_3(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p1,i) * mos_l_in_r_array_transp(ipoint,i) tmpvec_3(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p1,i) * mos_l_in_r_array_transp(ipoint,i) tmpvec_3(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p1,i) * mos_l_in_r_array_transp(ipoint,i) enddo - do h2 = 1, mo_num - tmp2( :,h2) = 0.d0 - tmp1(:,:,h2) = 0.d0 - if(mo_class(h2) .ne. "Active") cycle - + do pp2 = 1, n_act_orb + p2 = list_act(pp2) do ipoint = 1, n_points_final_grid - tmp2(ipoint,h2) = mos_r_in_r_array_transp_tmp(ipoint,h2) * tmpval_1(ipoint) & - + int2_grad1_u12_bimo_t(ipoint,1,i,h2) * tmpvec_1(ipoint,1) & - + int2_grad1_u12_bimo_t(ipoint,2,i,h2) * tmpvec_1(ipoint,2) & - + int2_grad1_u12_bimo_t(ipoint,3,i,h2) * tmpvec_1(ipoint,3) - - tmp1(ipoint,1,h2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,i,h2) - tmp1(ipoint,2,h2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,i,h2) - tmp1(ipoint,3,h2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,i,h2) + tmp2(ipoint,p2) = tmp2(ipoint,p2) + int2_grad1_u12_bimo_t(ipoint,1,p2,i) * tmpvec_2(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,p2,h1) * tmpvec_3(ipoint,1) & + + int2_grad1_u12_bimo_t(ipoint,2,p2,i) * tmpvec_2(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,p2,h1) * tmpvec_3(ipoint,2) & + + int2_grad1_u12_bimo_t(ipoint,3,p2,i) * tmpvec_2(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,p2,h1) * tmpvec_3(ipoint,3) + tmp1(ipoint,1,p2) = tmp1(ipoint,1,p2) + tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,i,p2) + tmp1(ipoint,2,p2) = tmp1(ipoint,2,p2) + tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,i,p2) + tmp1(ipoint,3,p2) = tmp1(ipoint,3,p2) + tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,i,p2) enddo enddo - call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 1.d0 & - , mos_l_in_r_array_transp_tmp(1,1), n_points_final_grid & - , tmp2(1,1), n_points_final_grid & - , 0.d0, tmp_2d(1,1), mo_num) + enddo ! ii - do h2 = 1, mo_num - do p2 = 1, mo_num - tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_2d(p2,h2) - enddo - enddo + call dgemm( 'T', 'N', mo_num, mo_num, 3*n_points_final_grid, 1.d0 & + , int2_grad1_u12_bimo_t_tmp(1,1,1,h1), 3*n_points_final_grid & + , tmp1(1,1,1), 3*n_points_final_grid & + , 1.d0, tmp(1,1,p1,h1), mo_num) - do p2 = 1, mo_num - tmp2( :,p2) = 0.d0 - tmp3(:,:,p2) = 0.d0 - if(mo_class(p2) .ne. "Active") cycle + call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 1.d0 & + , tmp2(1,1), n_points_final_grid & + , mos_r_in_r_array_transp_tmp(1,1), n_points_final_grid & + , 1.d0, tmp(1,1,p1,h1), mo_num) - do ipoint = 1, n_points_final_grid - - tmp2(ipoint,p2) = int2_grad1_u12_bimo_t(ipoint,1,p2,i) * tmpvec_2(ipoint,1) + int2_grad1_u12_bimo_t_tmp(ipoint,1,p2,h1) * tmpvec_3(ipoint,1) & - + int2_grad1_u12_bimo_t(ipoint,2,p2,i) * tmpvec_2(ipoint,2) + int2_grad1_u12_bimo_t_tmp(ipoint,2,p2,h1) * tmpvec_3(ipoint,2) & - + int2_grad1_u12_bimo_t(ipoint,3,p2,i) * tmpvec_2(ipoint,3) + int2_grad1_u12_bimo_t_tmp(ipoint,3,p2,h1) * tmpvec_3(ipoint,3) - - tmp3(ipoint,1,p2) = int2_grad1_u12_bimo_t_tmp(ipoint,1,p2,h1) - tmp3(ipoint,2,p2) = int2_grad1_u12_bimo_t_tmp(ipoint,2,p2,h1) - tmp3(ipoint,3,p2) = int2_grad1_u12_bimo_t_tmp(ipoint,3,p2,h1) - enddo - enddo - - call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 1.d0 & - , tmp2(1,1), n_points_final_grid & - , mos_r_in_r_array_transp_tmp(1,1), n_points_final_grid & - , 0.d0, tmp_2d(1,1), mo_num) - - call dgemm( 'T', 'N', mo_num, mo_num, 3*n_points_final_grid, 1.d0 & - , tmp3(1,1,1), 3*n_points_final_grid & - , tmp1(1,1,1), 3*n_points_final_grid & - , 1.d0, tmp_2d(1,1), mo_num) - - do h2 = 1, mo_num - do p2 = 1, mo_num - tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_2d(p2,h2) - enddo - enddo - - enddo ! p1 - enddo ! i + enddo ! p1 enddo ! h1 !$OMP END DO - deallocate(tmp_2d) deallocate(tmp_3d) deallocate(tmp1) deallocate(tmp2) - deallocate(tmp3) deallocate(tmpval_1) deallocate(tmpval_2) deallocate(tmpvec_1) @@ -741,41 +656,29 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_ !$OMP PARALLEL & !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, i, ii, h1, h2, p1, p2, tmp_2d, tmp_3d, tmp1, tmp2, tmp3, & + !$OMP PRIVATE (ipoint, i, ii, hh1, hh2, pp1, pp2, h1, h2, p1, p2, tmp_3d, tmp1, tmp2, & !$OMP tmpval_1, tmpval_2, tmpvec_1, tmpvec_2, tmpvec_3) & - !$OMP SHARED (n_points_final_grid, Ne, occ, mo_num, mo_class, & + !$OMP SHARED (n_points_final_grid, Ne, occ, mo_num, n_act_orb, list_act, & !$OMP mos_l_in_r_array_transp_tmp, mos_r_in_r_array_transp_tmp, & - !$OMP int2_grad1_u12_bimo_t_tmp, final_weight_at_r_vector, & + !$OMP final_weight_at_r_vector, int2_grad1_u12_bimo_t_tmp, & !$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, int2_grad1_u12_bimo_t, & !$OMP tmp) - allocate(tmp_2d(mo_num,mo_num)) allocate(tmp_3d(mo_num,mo_num,mo_num)) allocate(tmp1(n_points_final_grid,3,mo_num)) allocate(tmp2(n_points_final_grid,mo_num)) - allocate(tmp3(n_points_final_grid,3,mo_num)) allocate(tmpval_1(n_points_final_grid)) allocate(tmpval_2(n_points_final_grid)) allocate(tmpvec_1(n_points_final_grid,3)) allocate(tmpvec_2(n_points_final_grid,3)) allocate(tmpvec_3(n_points_final_grid,3)) - tmp_2d = 0.d0 - tmp_3d = 0.d0 - tmp1 = 0.d0 - tmp2 = 0.d0 - tmp3 = 0.d0 - tmpval_1 = 0.d0 - tmpval_2 = 0.d0 - tmpvec_1 = 0.d0 - tmpvec_2 = 0.d0 - tmpvec_3 = 0.d0 - !$OMP DO - do h1 = 1, mo_num - if(mo_class(h1) .ne. "Active") cycle + do hh1 = 1, n_act_orb + h1 = list_act(hh1) + tmp1 = 0.d0 do ii = Ne(2) + 1, Ne(1) i = occ(ii,1) @@ -783,163 +686,150 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_ 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_tmp(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) + enddo - 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_tmp(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_tmp(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_tmp(ipoint,h1) + do pp1 = 1, n_act_orb + p1 = list_act(pp1) + do ipoint = 1, n_points_final_grid + tmp1(ipoint,1,p1) = tmp1(ipoint,1,p1) + mos_l_in_r_array_transp(ipoint,p1) * tmpvec_1(ipoint,1) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p1,h1) + tmp1(ipoint,2,p1) = tmp1(ipoint,2,p1) + mos_l_in_r_array_transp(ipoint,p1) * tmpvec_1(ipoint,2) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p1,h1) + tmp1(ipoint,3,p1) = tmp1(ipoint,3,p1) + mos_l_in_r_array_transp(ipoint,p1) * tmpvec_1(ipoint,3) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p1,h1) + enddo + enddo + + enddo ! ii + + call dgemm( 'T', 'N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 0.5d0 & + , int2_grad1_u12_bimo_t_tmp(1,1,1,1), 3*n_points_final_grid & + , tmp1(1,1,1), 3*n_points_final_grid & + , 1.d0, tmp(1,1,1,h1), mo_num*mo_num) + + tmp1 = 0.d0 + do ii = Ne(2) + 1, Ne(1) + i = occ(ii,1) + + do ipoint = 1, n_points_final_grid + + 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_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 - do p1 = 1, mo_num - tmp1(:,:,p1) = 0.d0 - if(mo_class(p1) .ne. "Active") cycle - + do pp2 = 1, n_act_orb + p2 = list_act(pp2) do ipoint = 1, n_points_final_grid - tmp1(ipoint,1,p1) = mos_l_in_r_array_transp_tmp(ipoint,p1) * tmpvec_1(ipoint,1) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(ipoint,1,p1,h1) - tmp1(ipoint,2,p1) = mos_l_in_r_array_transp_tmp(ipoint,p1) * tmpvec_1(ipoint,2) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(ipoint,2,p1,h1) - tmp1(ipoint,3,p1) = mos_l_in_r_array_transp_tmp(ipoint,p1) * tmpvec_1(ipoint,3) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(ipoint,3,p1,h1) + tmp1(ipoint,1,p2) = tmp1(ipoint,1,p2) + tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p2,i) + mos_l_in_r_array_transp(ipoint,p2) * tmpvec_2(ipoint,1) + tmp1(ipoint,2,p2) = tmp1(ipoint,2,p2) + tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p2,i) + mos_l_in_r_array_transp(ipoint,p2) * tmpvec_2(ipoint,2) + tmp1(ipoint,3,p2) = tmp1(ipoint,3,p2) + tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p2,i) + mos_l_in_r_array_transp(ipoint,p2) * tmpvec_2(ipoint,3) enddo enddo - call dgemm( 'T', 'N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 0.5d0 & - , int2_grad1_u12_bimo_t_tmp(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) + enddo ! ii - do p1 = 1, mo_num - do h2 = 1, mo_num - do p2 = 1, mo_num - tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_3d(p2,h2,p1) + call dgemm( 'T', 'N', mo_num, mo_num*mo_num, 3*n_points_final_grid, 0.5d0 & + , tmp1(1,1,1), 3*n_points_final_grid & + , int2_grad1_u12_bimo_t_tmp(1,1,1,1), 3*n_points_final_grid & + , 0.d0, tmp_3d(1,1,1), mo_num) + + do p1 = 1, mo_num + do h2 = 1, mo_num + do p2 = 1, mo_num + tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_3d(p2,p1,h2) + enddo + enddo + enddo + + do pp1 = 1, n_act_orb + p1 = list_act(pp1) + + tmp2 = 0.d0 + do ii = Ne(2) + 1, Ne(1) + i = occ(ii,1) + + 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) ) + + tmpvec_1(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p1,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,p1,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,p1,i) * mos_r_in_r_array_transp(ipoint,h1) + enddo + + do hh2 = 1, n_act_orb + h2 = list_act(hh2) + do ipoint = 1, n_points_final_grid + tmp2(ipoint,h2) = tmp2(ipoint,h2) + mos_r_in_r_array_transp(ipoint,h2) * tmpval_1(ipoint) & + + int2_grad1_u12_bimo_t(ipoint,1,i,h2) * tmpvec_1(ipoint,1) & + + int2_grad1_u12_bimo_t(ipoint,2,i,h2) * tmpvec_1(ipoint,2) & + + int2_grad1_u12_bimo_t(ipoint,3,i,h2) * tmpvec_1(ipoint,3) enddo enddo - enddo - do p2 = 1, mo_num - tmp1(:,:,p2) = 0.d0 - if(mo_class(p2) .ne. "Active") cycle + enddo ! ii - do ipoint = 1, n_points_final_grid - tmp1(ipoint,1,p2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p2,i) + mos_l_in_r_array_transp_tmp(ipoint,p2) * tmpvec_2(ipoint,1) - tmp1(ipoint,2,p2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p2,i) + mos_l_in_r_array_transp_tmp(ipoint,p2) * tmpvec_2(ipoint,2) - tmp1(ipoint,3,p2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p2,i) + mos_l_in_r_array_transp_tmp(ipoint,p2) * tmpvec_2(ipoint,3) - enddo - enddo + call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 0.5d0 & + , mos_l_in_r_array_transp_tmp(1,1), n_points_final_grid & + , tmp2(1,1), n_points_final_grid & + , 1.d0, tmp(1,1,p1,h1), mo_num) - call dgemm( 'T', 'N', mo_num, mo_num*mo_num, 3*n_points_final_grid, 0.5d0 & - , tmp1(1,1,1), 3*n_points_final_grid & - , int2_grad1_u12_bimo_t_tmp(1,1,1,1), 3*n_points_final_grid & - , 0.d0, tmp_3d(1,1,1), mo_num) - do p1 = 1, mo_num - do h2 = 1, mo_num - do p2 = 1, mo_num - tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_3d(p2,p1,h2) - enddo - enddo - enddo - - do p1 = 1, mo_num - if(mo_class(p1) .ne. "Active") cycle + tmp1 = 0.d0 + tmp2 = 0.d0 + do ii = Ne(2) + 1, Ne(1) + i = occ(ii,1) 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_tmp(ipoint,1,p1,h1) & - + int2_grad1_u12_bimo_t(ipoint,2,i,i) * int2_grad1_u12_bimo_t_tmp(ipoint,2,p1,h1) & - + int2_grad1_u12_bimo_t(ipoint,3,i,i) * int2_grad1_u12_bimo_t_tmp(ipoint,3,p1,h1) ) + tmpval_2(ipoint) = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,p1) * mos_r_in_r_array_transp(ipoint,i) - tmpval_2(ipoint) = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp_tmp(ipoint,p1) * mos_r_in_r_array_transp(ipoint,i) - - tmpvec_1(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p1,i) * mos_r_in_r_array_transp_tmp(ipoint,h1) - tmpvec_1(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p1,i) * mos_r_in_r_array_transp_tmp(ipoint,h1) - tmpvec_1(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p1,i) * mos_r_in_r_array_transp_tmp(ipoint,h1) - - tmpvec_2(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,i,h1) * mos_l_in_r_array_transp_tmp(ipoint,p1) - tmpvec_2(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,i,h1) * mos_l_in_r_array_transp_tmp(ipoint,p1) - tmpvec_2(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,i,h1) * mos_l_in_r_array_transp_tmp(ipoint,p1) + tmpvec_2(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,i,h1) * mos_l_in_r_array_transp(ipoint,p1) + tmpvec_2(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,i,h1) * mos_l_in_r_array_transp(ipoint,p1) + tmpvec_2(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,i,h1) * mos_l_in_r_array_transp(ipoint,p1) tmpvec_3(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p1,i) * mos_l_in_r_array_transp(ipoint,i) tmpvec_3(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p1,i) * mos_l_in_r_array_transp(ipoint,i) tmpvec_3(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p1,i) * mos_l_in_r_array_transp(ipoint,i) enddo - do h2 = 1, mo_num - tmp2( :,h2) = 0.d0 - tmp1(:,:,h2) = 0.d0 - if(mo_class(h2) .ne. "Active") cycle - + do pp2 = 1, n_act_orb + p2 = list_act(pp2) do ipoint = 1, n_points_final_grid - - tmp2(ipoint,h2) = mos_r_in_r_array_transp_tmp(ipoint,h2) * tmpval_1(ipoint) & - + int2_grad1_u12_bimo_t(ipoint,1,i,h2) * tmpvec_1(ipoint,1) & - + int2_grad1_u12_bimo_t(ipoint,2,i,h2) * tmpvec_1(ipoint,2) & - + int2_grad1_u12_bimo_t(ipoint,3,i,h2) * tmpvec_1(ipoint,3) - - tmp1(ipoint,1,h2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,i,h2) - tmp1(ipoint,2,h2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,i,h2) - tmp1(ipoint,3,h2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,i,h2) - + tmp1(ipoint,1,p2) = tmp1(ipoint,1,p2) + tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,i,p2) + tmp1(ipoint,2,p2) = tmp1(ipoint,2,p2) + tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,i,p2) + tmp1(ipoint,3,p2) = tmp1(ipoint,3,p2) + tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,i,p2) + tmp2(ipoint,p2) = tmp2(ipoint,p2) + int2_grad1_u12_bimo_t(ipoint,1,p2,i) * tmpvec_2(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,p2,h1) * tmpvec_3(ipoint,1) & + + int2_grad1_u12_bimo_t(ipoint,2,p2,i) * tmpvec_2(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,p2,h1) * tmpvec_3(ipoint,2) & + + int2_grad1_u12_bimo_t(ipoint,3,p2,i) * tmpvec_2(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,p2,h1) * tmpvec_3(ipoint,3) enddo enddo - call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 0.5d0 & - , mos_l_in_r_array_transp_tmp(1,1), n_points_final_grid & - , tmp2(1,1), n_points_final_grid & - , 0.d0, tmp_2d(1,1), mo_num) + enddo ! ii - do h2 = 1, mo_num - do p2 = 1, mo_num - tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_2d(p2,h2) - enddo - enddo + call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 0.5d0 & + , tmp2(1,1), n_points_final_grid & + , mos_r_in_r_array_transp_tmp(1,1), n_points_final_grid & + , 1.d0, tmp(1,1,p1,h1), mo_num) - do p2 = 1, mo_num - tmp2( :,p2) = 0.d0 - tmp3(:,:,p2) = 0.d0 - if(mo_class(p2) .ne. "Active") cycle + call dgemm( 'T', 'N', mo_num, mo_num, 3*n_points_final_grid, 0.5d0 & + , int2_grad1_u12_bimo_t_tmp(1,1,1,h1), 3*n_points_final_grid & + , tmp1(1,1,1), 3*n_points_final_grid & + , 1.d0, tmp(1,1,p1,h1), mo_num) - do ipoint = 1, n_points_final_grid - - tmp2(ipoint,p2) = int2_grad1_u12_bimo_t(ipoint,1,p2,i) * tmpvec_2(ipoint,1) + int2_grad1_u12_bimo_t_tmp(ipoint,1,p2,h1) * tmpvec_3(ipoint,1) & - + int2_grad1_u12_bimo_t(ipoint,2,p2,i) * tmpvec_2(ipoint,2) + int2_grad1_u12_bimo_t_tmp(ipoint,2,p2,h1) * tmpvec_3(ipoint,2) & - + int2_grad1_u12_bimo_t(ipoint,3,p2,i) * tmpvec_2(ipoint,3) + int2_grad1_u12_bimo_t_tmp(ipoint,3,p2,h1) * tmpvec_3(ipoint,3) - - tmp3(ipoint,1,p2) = int2_grad1_u12_bimo_t_tmp(ipoint,1,p2,h1) - tmp3(ipoint,2,p2) = int2_grad1_u12_bimo_t_tmp(ipoint,2,p2,h1) - tmp3(ipoint,3,p2) = int2_grad1_u12_bimo_t_tmp(ipoint,3,p2,h1) - enddo - enddo - - call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 0.5d0 & - , tmp2(1,1), n_points_final_grid & - , mos_r_in_r_array_transp_tmp(1,1), n_points_final_grid & - , 0.d0, tmp_2d(1,1), mo_num) - - call dgemm( 'T', 'N', mo_num, mo_num, 3*n_points_final_grid, 0.5d0 & - , tmp3(1,1,1), 3*n_points_final_grid & - , tmp1(1,1,1), 3*n_points_final_grid & - , 1.d0, tmp_2d(1,1), mo_num) - - do h2 = 1, mo_num - do p2 = 1, mo_num - tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_2d(p2,h2) - enddo - enddo - - enddo ! p1 - enddo ! i + enddo ! p1 enddo ! h1 !$OMP END DO - deallocate(tmp_2d) deallocate(tmp_3d) deallocate(tmp1) deallocate(tmp2) - deallocate(tmp3) deallocate(tmpval_1) deallocate(tmpval_2) deallocate(tmpvec_1) diff --git a/src/tc_bi_ortho/normal_ordered_v0.irp.f b/src/tc_bi_ortho/normal_ordered_v0.irp.f new file mode 100644 index 00000000..784af9db --- /dev/null +++ b/src/tc_bi_ortho/normal_ordered_v0.irp.f @@ -0,0 +1,1022 @@ + +! --- + +BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth_v0, (mo_num, mo_num, mo_num, mo_num)] + + BEGIN_DOC + ! + ! Normal ordering of the three body interaction on the HF density + ! + END_DOC + + use bitmasks ! you need to include the bitmasks_module.f90 features + + implicit none + + integer :: i, ii, h1, p1, h2, p2, ipoint + integer :: hh1, hh2, pp1, pp2 + integer :: Ne(2) + double precision :: wall0, wall1, walli, wallf + integer, allocatable :: occ(:,:) + integer(bit_kind), allocatable :: key_i_core(:,:) + + PROVIDE mo_class + PROVIDE N_int + + print*,' Providing normal_two_body_bi_orth_v0 ...' + call wall_time(walli) + + if(read_tc_norm_ord) then + + open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/normal_two_body_bi_orth', action="read") + read(11) normal_two_body_bi_orth_v0 + close(11) + + else + + double precision, allocatable :: tmp_2d(:,:), tmp_3d(:,:,:) + double precision, allocatable :: tmp1(:,:,:), tmp2(:,:), tmp3(:,:,:) + double precision, allocatable :: tmpval_1(:), tmpval_2(:), tmpvec_1(:,:), tmpvec_2(:,:), tmpvec_3(:,:) + double precision, allocatable :: tmp(:,:,:,:) + double precision, allocatable :: int2_grad1_u12_bimo_t_tmp(:,:,:,:), mos_l_in_r_array_transp_tmp(:,:), mos_r_in_r_array_transp_tmp(:,:) + + PROVIDE int2_grad1_u12_bimo_t + PROVIDE mos_l_in_r_array_transp mos_r_in_r_array_transp + + allocate(int2_grad1_u12_bimo_t_tmp(n_points_final_grid,3,mo_num,mo_num)) + allocate(mos_l_in_r_array_transp_tmp(n_points_final_grid,mo_num)) + allocate(mos_r_in_r_array_transp_tmp(n_points_final_grid,mo_num)) + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (h1, p1) & + !$OMP SHARED (mo_num, mo_class, & + !$OMP int2_grad1_u12_bimo_t, int2_grad1_u12_bimo_t_tmp, & + !$OMP mos_l_in_r_array_transp, mos_l_in_r_array_transp_tmp, & + !$OMP mos_r_in_r_array_transp, mos_r_in_r_array_transp_tmp) + !$OMP DO + do h1 = 1, mo_num + + mos_l_in_r_array_transp_tmp(:,h1) = 0.d0 + mos_r_in_r_array_transp_tmp(:,h1) = 0.d0 + + if(mo_class(h1) .ne. "Active") cycle + + mos_l_in_r_array_transp_tmp(:,h1) = mos_l_in_r_array_transp(:,h1) + mos_r_in_r_array_transp_tmp(:,h1) = mos_r_in_r_array_transp(:,h1) + + do p1 = 1, mo_num + int2_grad1_u12_bimo_t_tmp(:,:,p1,h1) = 0.d0 + if(mo_class(p1) .ne. "Active") cycle + + int2_grad1_u12_bimo_t_tmp(:,:,p1,h1) = int2_grad1_u12_bimo_t(:,:,p1,h1) + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + allocate( occ(N_int*bit_kind_size,2) ) + allocate( key_i_core(N_int,2) ) + + if(core_tc_op) then + do i = 1, N_int + key_i_core(i,1) = xor(ref_bitmask(i,1), core_bitmask(i,1)) + key_i_core(i,2) = xor(ref_bitmask(i,2), core_bitmask(i,2)) + enddo + call bitstring_to_list_ab(key_i_core, occ, Ne, N_int) + else + call bitstring_to_list_ab(ref_bitmask, occ, Ne, N_int) + endif + + allocate(tmp(mo_num,mo_num,mo_num,mo_num)) + + ! --- + ! aba contraction + + print*,' Providing aba_contraction_v0 ...' + call wall_time(wall0) + + call set_multiple_levels_omp(.false.) + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, h1, p1, h2, p2, i, ii, & + !$OMP tmp_3d, tmp_2d, tmp1, tmp2, & + !$OMP tmpval_1, tmpval_2, tmpvec_1, tmpvec_2) & + !$OMP SHARED (n_points_final_grid, Ne, occ, mo_num, mo_class, & + !$OMP mos_l_in_r_array_transp_tmp, mos_r_in_r_array_transp_tmp, & + !$OMP int2_grad1_u12_bimo_t_tmp, final_weight_at_r_vector, & + !$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, int2_grad1_u12_bimo_t, & + !$OMP tmp) + + allocate(tmp_3d(mo_num,mo_num,mo_num), tmp_2d(mo_num,mo_num)) + allocate(tmp1(n_points_final_grid,3,mo_num), tmp2(n_points_final_grid,mo_num)) + allocate(tmpval_1(n_points_final_grid), tmpval_2(n_points_final_grid)) + allocate(tmpvec_1(n_points_final_grid,3), tmpvec_2(n_points_final_grid,3)) + + tmp_3d = 0.d0 + tmp_2d = 0.d0 + tmp1 = 0.d0 + tmp2 = 0.d0 + tmpval_1 = 0.d0 + tmpval_2 = 0.d0 + tmpvec_1 = 0.d0 + tmpvec_2 = 0.d0 + + ! TODO: active electrons + + !$OMP DO + + do h1 = 1, mo_num + tmp(:,:,:,h1) = 0.d0 + if(mo_class(h1) .ne. "Active") cycle + + do ii = 1, Ne(2) + i = occ(ii,2) + + 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_tmp(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_tmp(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_tmp(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_tmp(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 + + do p1 = 1, mo_num + tmp1(:,:,p1) = 0.d0 + if(mo_class(p1) .ne. "Active") cycle + + do ipoint = 1, n_points_final_grid + tmp1(ipoint,1,p1) = mos_l_in_r_array_transp_tmp(ipoint,p1) * (tmpvec_1(ipoint,1) - tmpvec_2(ipoint,1)) & + + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(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_tmp(ipoint,p1) * (tmpvec_1(ipoint,2) - tmpvec_2(ipoint,2)) & + + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(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_tmp(ipoint,p1) * (tmpvec_1(ipoint,3) - tmpvec_2(ipoint,3)) & + + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(ipoint,3,p1,h1) - tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p1,i) + enddo + enddo + + call dgemm( 'T', 'N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 1.d0 & + , int2_grad1_u12_bimo_t_tmp(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) + + do p1 = 1, mo_num + do h2 = 1, mo_num + do p2 = 1, mo_num + tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_3d(p2,h2,p1) + enddo + enddo + enddo + + do p1 = 1, mo_num + if(mo_class(p1) .ne. "Active") cycle + + 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_tmp(ipoint,1,p1,h1) & + + int2_grad1_u12_bimo_t(ipoint,2,i,i) * int2_grad1_u12_bimo_t_tmp(ipoint,2,p1,h1) & + + int2_grad1_u12_bimo_t(ipoint,3,i,i) * int2_grad1_u12_bimo_t_tmp(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 + + do h2 = 1, mo_num + tmp2(:,h2) = 0.d0 + if(mo_class(h2) .ne. "Active") cycle + + do ipoint = 1, n_points_final_grid + tmp2(ipoint,h2) = mos_r_in_r_array_transp_tmp(ipoint,h2) * tmpval_1(ipoint) + enddo + enddo + + call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 1.d0 & + , mos_l_in_r_array_transp_tmp(1,1), n_points_final_grid & + , tmp2(1,1), n_points_final_grid & + , 0.d0, tmp_2d(1,1), mo_num) + + do h2 = 1, mo_num + do p2 = 1, mo_num + tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_2d(p2,h2) + enddo + enddo + + enddo ! p1 + enddo ! i + enddo ! h1 + + !$OMP END DO + + deallocate(tmp_3d, tmp_2d) + deallocate(tmp1, tmp2) + deallocate(tmpval_1, tmpval_2) + deallocate(tmpvec_1, tmpvec_2) + + !$OMP END PARALLEL + + + ! purely open-shell part + if(Ne(2) < Ne(1)) then + + call set_multiple_levels_omp(.false.) + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, h1, p1, h2, p2, i, ii, & + !$OMP tmp_3d, tmp_2d, tmp1, tmp2, & + !$OMP tmpval_1, tmpval_2, tmpvec_1, tmpvec_2) & + !$OMP SHARED (n_points_final_grid, Ne, occ, mo_num, mo_class, & + !$OMP mos_l_in_r_array_transp_tmp, mos_r_in_r_array_transp_tmp, & + !$OMP int2_grad1_u12_bimo_t_tmp, final_weight_at_r_vector, & + !$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, int2_grad1_u12_bimo_t, & + !$OMP tmp) + + Allocate(tmp_3d(mo_num,mo_num,mo_num), tmp_2d(mo_num,mo_num)) + Allocate(tmp1(n_points_final_grid,3,mo_num), tmp2(n_points_final_grid,mo_num)) + Allocate(tmpval_1(n_points_final_grid), tmpval_2(n_points_final_grid)) + Allocate(tmpvec_1(n_points_final_grid,3), tmpvec_2(n_points_final_grid,3)) + + Tmp_3d = 0.d0 + Tmp_2d = 0.d0 + Tmp1 = 0.d0 + Tmp2 = 0.d0 + Tmpval_1 = 0.d0 + Tmpval_2 = 0.d0 + Tmpvec_1 = 0.d0 + Tmpvec_2 = 0.d0 + + !$OMP DO + + do h1 = 1, mo_num + if(mo_class(h1) .ne. "Active") cycle + + do ii = Ne(2) + 1, Ne(1) + i = occ(ii,1) + + 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_tmp(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_tmp(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_tmp(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_tmp(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 + + do p1 = 1, mo_num + tmp1(:,:,p1) = 0.d0 + if(mo_class(p1) .ne. "Active") cycle + + do ipoint = 1, n_points_final_grid + tmp1(ipoint,1,p1) = mos_l_in_r_array_transp_tmp(ipoint,p1) * (tmpvec_1(ipoint,1) - tmpvec_2(ipoint,1)) & + + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(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_tmp(ipoint,p1) * (tmpvec_1(ipoint,2) - tmpvec_2(ipoint,2)) & + + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(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_tmp(ipoint,p1) * (tmpvec_1(ipoint,3) - tmpvec_2(ipoint,3)) & + + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(ipoint,3,p1,h1) - tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p1,i) + enddo + enddo + + call dgemm( 'T', 'N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 0.5d0 & + , int2_grad1_u12_bimo_t_tmp(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) + + do p1 = 1, mo_num + do h2 = 1, mo_num + do p2 = 1, mo_num + tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_3d(p2,h2,p1) + enddo + enddo + enddo + + do p1 = 1, mo_num + if(mo_class(p1) .ne. "Active") cycle + + 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_tmp(ipoint,1,p1,h1) & + + int2_grad1_u12_bimo_t(ipoint,2,i,i) * int2_grad1_u12_bimo_t_tmp(ipoint,2,p1,h1) & + + int2_grad1_u12_bimo_t(ipoint,3,i,i) * int2_grad1_u12_bimo_t_tmp(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 + + do h2 = 1, mo_num + tmp2(:,h2) = 0.d0 + if(mo_class(h2) .ne. "Active") cycle + + do ipoint = 1, n_points_final_grid + tmp2(ipoint,h2) = mos_r_in_r_array_transp_tmp(ipoint,h2) * tmpval_1(ipoint) + enddo + enddo + + call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 0.5d0 & + , mos_l_in_r_array_transp_tmp(1,1), n_points_final_grid & + , tmp2(1,1), n_points_final_grid & + , 0.d0, tmp_2d(1,1), mo_num) + + do h2 = 1, mo_num + do p2 = 1, mo_num + tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_2d(p2,h2) + enddo + enddo + + enddo ! p1 + enddo ! i + enddo ! h1 + !$OMP END DO + + deallocate(tmp_3d, tmp_2d) + deallocate(tmp1, tmp2) + deallocate(tmpval_1, tmpval_2) + deallocate(tmpvec_1, tmpvec_2) + + !$OMP END PARALLEL + endif + + tmp = -0.5d0 * tmp + call sum_A_At(tmp(1,1,1,1), mo_num*mo_num) + + call wall_time(wall1) + print*,' Wall time for aba_contraction_v0', wall1-wall0 + + normal_two_body_bi_orth_v0 = tmp + + ! --- + ! aab contraction + + print*,' Providing aab_contraction_v0 ...' + call wall_time(wall0) + + call set_multiple_levels_omp(.false.) + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, ii, i, h1, p1, h2, p2, & + !$OMP tmp_2d, tmp_3d, tmp1, tmp2, & + !$OMP tmpval_1, tmpvec_1) & + !$OMP SHARED (n_points_final_grid, mo_num, Ne, occ, mo_class, & + !$OMP mos_l_in_r_array_transp_tmp, mos_r_in_r_array_transp_tmp, & + !$OMP int2_grad1_u12_bimo_t_tmp, final_weight_at_r_vector, & + !$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, int2_grad1_u12_bimo_t, & + !$OMP tmp) + + allocate(tmp_2d(mo_num,mo_num)) + allocate(tmp_3d(mo_num,mo_num,mo_num)) + allocate(tmp1(n_points_final_grid,3,mo_num)) + allocate(tmp2(n_points_final_grid,mo_num)) + allocate(tmpval_1(n_points_final_grid)) + allocate(tmpvec_1(n_points_final_grid,3)) + + tmp_2d = 0.d0 + tmp_3d = 0.d0 + tmp1 = 0.d0 + tmp2 = 0.d0 + tmpval_1 = 0.d0 + tmpvec_1 = 0.d0 + + !$OMP DO + + do h1 = 1, mo_num + tmp(:,:,:,h1) = 0.d0 + if(mo_class(h1) .ne. "Active") cycle + + do ii = 1, Ne(2) + i = occ(ii,2) + + 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) + 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_tmp(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_tmp(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_tmp(ipoint,h1) + enddo + + do p1 = 1, mo_num + tmp1(:,:,p1) = 0.d0 + if(mo_class(p1) .ne. "Active") cycle + + do ipoint = 1, n_points_final_grid + tmp1(ipoint,1,p1) = mos_l_in_r_array_transp_tmp(ipoint,p1) * tmpvec_1(ipoint,1) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(ipoint,1,p1,h1) + tmp1(ipoint,2,p1) = mos_l_in_r_array_transp_tmp(ipoint,p1) * tmpvec_1(ipoint,2) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(ipoint,2,p1,h1) + tmp1(ipoint,3,p1) = mos_l_in_r_array_transp_tmp(ipoint,p1) * tmpvec_1(ipoint,3) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(ipoint,3,p1,h1) + enddo + enddo + + call dgemm( 'T', 'N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 1.d0 & + , int2_grad1_u12_bimo_t_tmp(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) + + do p1 = 1, mo_num + do h2 = 1, mo_num + do p2 = 1, mo_num + tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_3d(p2,h2,p1) + enddo + enddo + enddo + + do p1 = 1, mo_num + if(mo_class(p1) .ne. "Active") cycle + + 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_tmp(ipoint,1,p1,h1) & + + int2_grad1_u12_bimo_t(ipoint,2,i,i) * int2_grad1_u12_bimo_t_tmp(ipoint,2,p1,h1) & + + int2_grad1_u12_bimo_t(ipoint,3,i,i) * int2_grad1_u12_bimo_t_tmp(ipoint,3,p1,h1) ) + enddo + + do h2 = 1, mo_num + if(mo_class(h2) .ne. "Active") cycle + tmp2(:,h2) = 0.d0 + + do ipoint = 1, n_points_final_grid + tmp2(ipoint,h2) = mos_r_in_r_array_transp_tmp(ipoint,h2) * tmpval_1(ipoint) + enddo + enddo + + call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 1.d0 & + , mos_l_in_r_array_transp_tmp(1,1), n_points_final_grid & + , tmp2(1,1), n_points_final_grid & + , 0.d0, tmp_2d(1,1), mo_num) + + do h2 = 1, mo_num + do p2 = 1, mo_num + tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_2d(p2,h2) + enddo + enddo + + enddo ! p1 + enddo ! i + enddo ! h1 + + !$OMP END DO + + deallocate(tmp_3d) + deallocate(tmp1, tmp2) + deallocate(tmpval_1) + deallocate(tmpvec_1) + + !$OMP END PARALLEL + + tmp = -0.5d0 * tmp + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (h1, h2, p1, p2) & + !$OMP SHARED (tmp, mo_num) + + !$OMP DO COLLAPSE(2) + do h1 = 1, mo_num + do h2 = 1, mo_num + do p1 = 1, mo_num + do p2 = p1, mo_num + tmp(p2,h2,p1,h1) -= tmp(p1,h2,p2,h1) + enddo + enddo + enddo + enddo + !$OMP END DO + + !$OMP DO COLLAPSE(2) + do h1 = 1, mo_num + do h2 = 1, mo_num + do p1 = 2, mo_num + do p2 = 1, p1-1 + tmp(p2,h2,p1,h1) = -tmp(p1,h2,p2,h1) + enddo + enddo + enddo + enddo + !$OMP END DO + + !$OMP DO + do h1 = 1, mo_num-1 + do h2 = h1+1, mo_num + do p1 = 2, mo_num + do p2 = 1, p1-1 + tmp(p2,h2,p1,h1) *= -1.d0 + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(wall1) + print*,' Wall time for aab_contraction_v0', wall1-wall0 + + normal_two_body_bi_orth_v0 += tmp + + ! --- + ! aaa contraction + + if(Ne(2) .ge. 3) then + + print*,' Providing aaa_contraction_v0 ...' + call wall_time(wall0) + + call set_multiple_levels_omp(.false.) + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, ii, h1, h2, p1, p2, & + !$OMP tmp_2d, tmp_3d, tmp1, tmp2, tmp3, & + !$OMP tmpval_1, tmpval_2, & + !$OMP tmpvec_1, tmpvec_2, tmpvec_3) & + !$OMP SHARED (n_points_final_grid, Ne, occ, mo_num, mo_class, & + !$OMP mos_l_in_r_array_transp_tmp, mos_r_in_r_array_transp_tmp, & + !$OMP int2_grad1_u12_bimo_t_tmp, final_weight_at_r_vector, & + !$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, int2_grad1_u12_bimo_t, & + !$OMP tmp) + + allocate(tmp_2d(mo_num,mo_num)) + allocate(tmp_3d(mo_num,mo_num,mo_num)) + allocate(tmp1(n_points_final_grid,3,mo_num)) + allocate(tmp2(n_points_final_grid,mo_num)) + allocate(tmp3(n_points_final_grid,3,mo_num)) + allocate(tmpval_1(n_points_final_grid)) + allocate(tmpval_2(n_points_final_grid)) + allocate(tmpvec_1(n_points_final_grid,3)) + allocate(tmpvec_2(n_points_final_grid,3)) + allocate(tmpvec_3(n_points_final_grid,3)) + + tmp_2d = 0.d0 + tmp_3d = 0.d0 + tmp1 = 0.d0 + tmp2 = 0.d0 + tmp3 = 0.d0 + tmpval_1 = 0.d0 + tmpval_2 = 0.d0 + tmpvec_1 = 0.d0 + tmpvec_2 = 0.d0 + tmpvec_3 = 0.d0 + + !$OMP DO + + do h1 = 1, mo_num + tmp(:,:,:,h1) = 0.d0 + if(mo_class(h1) .ne. "Active") cycle + + do ii = 1, Ne(2) + i = occ(ii,2) + + 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_tmp(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_tmp(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_tmp(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_tmp(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 + + do p1 = 1, mo_num + tmp1(:,:,p1) = 0.d0 + if(mo_class(p1) .ne. "Active") cycle + + do ipoint = 1, n_points_final_grid + tmp1(ipoint,1,p1) = mos_l_in_r_array_transp_tmp(ipoint,p1) * tmpvec_1(ipoint,1) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(ipoint,1,p1,h1) + tmp1(ipoint,2,p1) = mos_l_in_r_array_transp_tmp(ipoint,p1) * tmpvec_1(ipoint,2) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(ipoint,2,p1,h1) + tmp1(ipoint,3,p1) = mos_l_in_r_array_transp_tmp(ipoint,p1) * tmpvec_1(ipoint,3) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(ipoint,3,p1,h1) + enddo + enddo + + call dgemm( 'T', 'N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 1.d0 & + , int2_grad1_u12_bimo_t_tmp(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) + + do p1 = 1, mo_num + do h2 = 1, mo_num + do p2 = 1, mo_num + tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_3d(p2,h2,p1) + enddo + enddo + enddo + + do p2 = 1, mo_num + tmp1(:,:,p2) = 0.d0 + if(mo_class(p2) .ne. "Active") cycle + + do ipoint = 1, n_points_final_grid + tmp1(ipoint,1,p2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p2,i) + mos_l_in_r_array_transp_tmp(ipoint,p2) * tmpvec_2(ipoint,1) + tmp1(ipoint,2,p2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p2,i) + mos_l_in_r_array_transp_tmp(ipoint,p2) * tmpvec_2(ipoint,2) + tmp1(ipoint,3,p2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p2,i) + mos_l_in_r_array_transp_tmp(ipoint,p2) * tmpvec_2(ipoint,3) + enddo + enddo + + call dgemm( 'T', 'N', mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0 & + , tmp1(1,1,1), 3*n_points_final_grid & + , int2_grad1_u12_bimo_t_tmp(1,1,1,1), 3*n_points_final_grid & + , 0.d0, tmp_3d(1,1,1), mo_num) + + do p1 = 1, mo_num + do h2 = 1, mo_num + do p2 = 1, mo_num + tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_3d(p2,p1,h2) + enddo + enddo + enddo + + do p1 = 1, mo_num + if(mo_class(p1) .ne. "Active") cycle + + 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_tmp(ipoint,1,p1,h1) & + + int2_grad1_u12_bimo_t(ipoint,2,i,i) * int2_grad1_u12_bimo_t_tmp(ipoint,2,p1,h1) & + + int2_grad1_u12_bimo_t(ipoint,3,i,i) * int2_grad1_u12_bimo_t_tmp(ipoint,3,p1,h1) ) + + tmpval_2(ipoint) = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp_tmp(ipoint,p1) * mos_r_in_r_array_transp(ipoint,i) + + tmpvec_1(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p1,i) * mos_r_in_r_array_transp_tmp(ipoint,h1) + tmpvec_1(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p1,i) * mos_r_in_r_array_transp_tmp(ipoint,h1) + tmpvec_1(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p1,i) * mos_r_in_r_array_transp_tmp(ipoint,h1) + + tmpvec_2(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,i,h1) * mos_l_in_r_array_transp_tmp(ipoint,p1) + tmpvec_2(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,i,h1) * mos_l_in_r_array_transp_tmp(ipoint,p1) + tmpvec_2(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,i,h1) * mos_l_in_r_array_transp_tmp(ipoint,p1) + + tmpvec_3(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p1,i) * mos_l_in_r_array_transp(ipoint,i) + tmpvec_3(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p1,i) * mos_l_in_r_array_transp(ipoint,i) + tmpvec_3(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p1,i) * mos_l_in_r_array_transp(ipoint,i) + enddo + + do h2 = 1, mo_num + tmp2( :,h2) = 0.d0 + tmp1(:,:,h2) = 0.d0 + if(mo_class(h2) .ne. "Active") cycle + + do ipoint = 1, n_points_final_grid + + tmp2(ipoint,h2) = mos_r_in_r_array_transp_tmp(ipoint,h2) * tmpval_1(ipoint) & + + int2_grad1_u12_bimo_t(ipoint,1,i,h2) * tmpvec_1(ipoint,1) & + + int2_grad1_u12_bimo_t(ipoint,2,i,h2) * tmpvec_1(ipoint,2) & + + int2_grad1_u12_bimo_t(ipoint,3,i,h2) * tmpvec_1(ipoint,3) + + tmp1(ipoint,1,h2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,i,h2) + tmp1(ipoint,2,h2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,i,h2) + tmp1(ipoint,3,h2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,i,h2) + + enddo + enddo + + call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 1.d0 & + , mos_l_in_r_array_transp_tmp(1,1), n_points_final_grid & + , tmp2(1,1), n_points_final_grid & + , 0.d0, tmp_2d(1,1), mo_num) + + do h2 = 1, mo_num + do p2 = 1, mo_num + tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_2d(p2,h2) + enddo + enddo + + do p2 = 1, mo_num + tmp2( :,p2) = 0.d0 + tmp3(:,:,p2) = 0.d0 + if(mo_class(p2) .ne. "Active") cycle + + do ipoint = 1, n_points_final_grid + + tmp2(ipoint,p2) = int2_grad1_u12_bimo_t(ipoint,1,p2,i) * tmpvec_2(ipoint,1) + int2_grad1_u12_bimo_t_tmp(ipoint,1,p2,h1) * tmpvec_3(ipoint,1) & + + int2_grad1_u12_bimo_t(ipoint,2,p2,i) * tmpvec_2(ipoint,2) + int2_grad1_u12_bimo_t_tmp(ipoint,2,p2,h1) * tmpvec_3(ipoint,2) & + + int2_grad1_u12_bimo_t(ipoint,3,p2,i) * tmpvec_2(ipoint,3) + int2_grad1_u12_bimo_t_tmp(ipoint,3,p2,h1) * tmpvec_3(ipoint,3) + + tmp3(ipoint,1,p2) = int2_grad1_u12_bimo_t_tmp(ipoint,1,p2,h1) + tmp3(ipoint,2,p2) = int2_grad1_u12_bimo_t_tmp(ipoint,2,p2,h1) + tmp3(ipoint,3,p2) = int2_grad1_u12_bimo_t_tmp(ipoint,3,p2,h1) + enddo + enddo + + call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 1.d0 & + , tmp2(1,1), n_points_final_grid & + , mos_r_in_r_array_transp_tmp(1,1), n_points_final_grid & + , 0.d0, tmp_2d(1,1), mo_num) + + call dgemm( 'T', 'N', mo_num, mo_num, 3*n_points_final_grid, 1.d0 & + , tmp3(1,1,1), 3*n_points_final_grid & + , tmp1(1,1,1), 3*n_points_final_grid & + , 1.d0, tmp_2d(1,1), mo_num) + + do h2 = 1, mo_num + do p2 = 1, mo_num + tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_2d(p2,h2) + enddo + enddo + + enddo ! p1 + enddo ! i + enddo ! h1 + !$OMP END DO + + deallocate(tmp_2d) + deallocate(tmp_3d) + deallocate(tmp1) + deallocate(tmp2) + deallocate(tmp3) + deallocate(tmpval_1) + deallocate(tmpval_2) + deallocate(tmpvec_1) + deallocate(tmpvec_2) + deallocate(tmpvec_3) + + !$OMP END PARALLEL + + ! purely open-shell part + if(Ne(2) < Ne(1)) then + + call set_multiple_levels_omp(.false.) + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, ii, h1, h2, p1, p2, tmp_2d, tmp_3d, tmp1, tmp2, tmp3, & + !$OMP tmpval_1, tmpval_2, tmpvec_1, tmpvec_2, tmpvec_3) & + !$OMP SHARED (n_points_final_grid, Ne, occ, mo_num, mo_class, & + !$OMP mos_l_in_r_array_transp_tmp, mos_r_in_r_array_transp_tmp, & + !$OMP int2_grad1_u12_bimo_t_tmp, final_weight_at_r_vector, & + !$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, int2_grad1_u12_bimo_t, & + !$OMP tmp) + + allocate(tmp_2d(mo_num,mo_num)) + allocate(tmp_3d(mo_num,mo_num,mo_num)) + allocate(tmp1(n_points_final_grid,3,mo_num)) + allocate(tmp2(n_points_final_grid,mo_num)) + allocate(tmp3(n_points_final_grid,3,mo_num)) + allocate(tmpval_1(n_points_final_grid)) + allocate(tmpval_2(n_points_final_grid)) + allocate(tmpvec_1(n_points_final_grid,3)) + allocate(tmpvec_2(n_points_final_grid,3)) + allocate(tmpvec_3(n_points_final_grid,3)) + + tmp_2d = 0.d0 + tmp_3d = 0.d0 + tmp1 = 0.d0 + tmp2 = 0.d0 + tmp3 = 0.d0 + tmpval_1 = 0.d0 + tmpval_2 = 0.d0 + tmpvec_1 = 0.d0 + tmpvec_2 = 0.d0 + tmpvec_3 = 0.d0 + + !$OMP DO + + do h1 = 1, mo_num + if(mo_class(h1) .ne. "Active") cycle + + do ii = Ne(2) + 1, Ne(1) + i = occ(ii,1) + + 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_tmp(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_tmp(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_tmp(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_tmp(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 + + do p1 = 1, mo_num + tmp1(:,:,p1) = 0.d0 + if(mo_class(p1) .ne. "Active") cycle + + do ipoint = 1, n_points_final_grid + tmp1(ipoint,1,p1) = mos_l_in_r_array_transp_tmp(ipoint,p1) * tmpvec_1(ipoint,1) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(ipoint,1,p1,h1) + tmp1(ipoint,2,p1) = mos_l_in_r_array_transp_tmp(ipoint,p1) * tmpvec_1(ipoint,2) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(ipoint,2,p1,h1) + tmp1(ipoint,3,p1) = mos_l_in_r_array_transp_tmp(ipoint,p1) * tmpvec_1(ipoint,3) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t_tmp(ipoint,3,p1,h1) + enddo + enddo + + call dgemm( 'T', 'N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 0.5d0 & + , int2_grad1_u12_bimo_t_tmp(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) + + do p1 = 1, mo_num + do h2 = 1, mo_num + do p2 = 1, mo_num + tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_3d(p2,h2,p1) + enddo + enddo + enddo + + do p2 = 1, mo_num + tmp1(:,:,p2) = 0.d0 + if(mo_class(p2) .ne. "Active") cycle + + do ipoint = 1, n_points_final_grid + tmp1(ipoint,1,p2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p2,i) + mos_l_in_r_array_transp_tmp(ipoint,p2) * tmpvec_2(ipoint,1) + tmp1(ipoint,2,p2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p2,i) + mos_l_in_r_array_transp_tmp(ipoint,p2) * tmpvec_2(ipoint,2) + tmp1(ipoint,3,p2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p2,i) + mos_l_in_r_array_transp_tmp(ipoint,p2) * tmpvec_2(ipoint,3) + enddo + enddo + + call dgemm( 'T', 'N', mo_num, mo_num*mo_num, 3*n_points_final_grid, 0.5d0 & + , tmp1(1,1,1), 3*n_points_final_grid & + , int2_grad1_u12_bimo_t_tmp(1,1,1,1), 3*n_points_final_grid & + , 0.d0, tmp_3d(1,1,1), mo_num) + + do p1 = 1, mo_num + do h2 = 1, mo_num + do p2 = 1, mo_num + tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_3d(p2,p1,h2) + enddo + enddo + enddo + + do p1 = 1, mo_num + if(mo_class(p1) .ne. "Active") cycle + + 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_tmp(ipoint,1,p1,h1) & + + int2_grad1_u12_bimo_t(ipoint,2,i,i) * int2_grad1_u12_bimo_t_tmp(ipoint,2,p1,h1) & + + int2_grad1_u12_bimo_t(ipoint,3,i,i) * int2_grad1_u12_bimo_t_tmp(ipoint,3,p1,h1) ) + + tmpval_2(ipoint) = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp_tmp(ipoint,p1) * mos_r_in_r_array_transp(ipoint,i) + + tmpvec_1(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p1,i) * mos_r_in_r_array_transp_tmp(ipoint,h1) + tmpvec_1(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p1,i) * mos_r_in_r_array_transp_tmp(ipoint,h1) + tmpvec_1(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p1,i) * mos_r_in_r_array_transp_tmp(ipoint,h1) + + tmpvec_2(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,i,h1) * mos_l_in_r_array_transp_tmp(ipoint,p1) + tmpvec_2(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,i,h1) * mos_l_in_r_array_transp_tmp(ipoint,p1) + tmpvec_2(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,i,h1) * mos_l_in_r_array_transp_tmp(ipoint,p1) + + tmpvec_3(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p1,i) * mos_l_in_r_array_transp(ipoint,i) + tmpvec_3(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p1,i) * mos_l_in_r_array_transp(ipoint,i) + tmpvec_3(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p1,i) * mos_l_in_r_array_transp(ipoint,i) + enddo + + do h2 = 1, mo_num + tmp2( :,h2) = 0.d0 + tmp1(:,:,h2) = 0.d0 + if(mo_class(h2) .ne. "Active") cycle + + do ipoint = 1, n_points_final_grid + + tmp2(ipoint,h2) = mos_r_in_r_array_transp_tmp(ipoint,h2) * tmpval_1(ipoint) & + + int2_grad1_u12_bimo_t(ipoint,1,i,h2) * tmpvec_1(ipoint,1) & + + int2_grad1_u12_bimo_t(ipoint,2,i,h2) * tmpvec_1(ipoint,2) & + + int2_grad1_u12_bimo_t(ipoint,3,i,h2) * tmpvec_1(ipoint,3) + + tmp1(ipoint,1,h2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,i,h2) + tmp1(ipoint,2,h2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,i,h2) + tmp1(ipoint,3,h2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,i,h2) + + enddo + enddo + + call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 0.5d0 & + , mos_l_in_r_array_transp_tmp(1,1), n_points_final_grid & + , tmp2(1,1), n_points_final_grid & + , 0.d0, tmp_2d(1,1), mo_num) + + do h2 = 1, mo_num + do p2 = 1, mo_num + tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_2d(p2,h2) + enddo + enddo + + do p2 = 1, mo_num + tmp2( :,p2) = 0.d0 + tmp3(:,:,p2) = 0.d0 + if(mo_class(p2) .ne. "Active") cycle + + do ipoint = 1, n_points_final_grid + + tmp2(ipoint,p2) = int2_grad1_u12_bimo_t(ipoint,1,p2,i) * tmpvec_2(ipoint,1) + int2_grad1_u12_bimo_t_tmp(ipoint,1,p2,h1) * tmpvec_3(ipoint,1) & + + int2_grad1_u12_bimo_t(ipoint,2,p2,i) * tmpvec_2(ipoint,2) + int2_grad1_u12_bimo_t_tmp(ipoint,2,p2,h1) * tmpvec_3(ipoint,2) & + + int2_grad1_u12_bimo_t(ipoint,3,p2,i) * tmpvec_2(ipoint,3) + int2_grad1_u12_bimo_t_tmp(ipoint,3,p2,h1) * tmpvec_3(ipoint,3) + + tmp3(ipoint,1,p2) = int2_grad1_u12_bimo_t_tmp(ipoint,1,p2,h1) + tmp3(ipoint,2,p2) = int2_grad1_u12_bimo_t_tmp(ipoint,2,p2,h1) + tmp3(ipoint,3,p2) = int2_grad1_u12_bimo_t_tmp(ipoint,3,p2,h1) + enddo + enddo + + call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 0.5d0 & + , tmp2(1,1), n_points_final_grid & + , mos_r_in_r_array_transp_tmp(1,1), n_points_final_grid & + , 0.d0, tmp_2d(1,1), mo_num) + + call dgemm( 'T', 'N', mo_num, mo_num, 3*n_points_final_grid, 0.5d0 & + , tmp3(1,1,1), 3*n_points_final_grid & + , tmp1(1,1,1), 3*n_points_final_grid & + , 1.d0, tmp_2d(1,1), mo_num) + + do h2 = 1, mo_num + do p2 = 1, mo_num + tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_2d(p2,h2) + enddo + enddo + + enddo ! p1 + enddo ! i + enddo ! h1 + !$OMP END DO + + deallocate(tmp_2d) + deallocate(tmp_3d) + deallocate(tmp1) + deallocate(tmp2) + deallocate(tmp3) + deallocate(tmpval_1) + deallocate(tmpval_2) + deallocate(tmpvec_1) + deallocate(tmpvec_2) + deallocate(tmpvec_3) + + !$OMP END PARALLEL + endif + + tmp = -0.5d0 * tmp + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (h1, h2, p1, p2) & + !$OMP SHARED (tmp, mo_num) + + !$OMP DO COLLAPSE(2) + do h1 = 1, mo_num + do h2 = 1, mo_num + do p1 = 1, mo_num + do p2 = p1, mo_num + tmp(p2,h2,p1,h1) -= tmp(p1,h2,p2,h1) + enddo + enddo + enddo + enddo + !$OMP END DO + + !$OMP DO COLLAPSE(2) + do h1 = 1, mo_num + do h2 = 1, mo_num + do p1 = 2, mo_num + do p2 = 1, p1-1 + tmp(p2,h2,p1,h1) = -tmp(p1,h2,p2,h1) + enddo + enddo + enddo + enddo + !$OMP END DO + + !$OMP DO + do h1 = 1, mo_num-1 + do h2 = h1+1, mo_num + do p1 = 2, mo_num + do p2 = 1, p1-1 + tmp(p2,h2,p1,h1) *= -1.d0 + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(wall1) + print*,' Wall time for aaa_contraction_v0', wall1-wall0 + + normal_two_body_bi_orth_v0 += tmp + endif ! Ne(2) .ge. 3 + + deallocate(tmp) + deallocate(int2_grad1_u12_bimo_t_tmp, mos_l_in_r_array_transp_tmp, mos_r_in_r_array_transp_tmp) + + endif ! read_tc_norm_ord + + if(write_tc_norm_ord.and.mpi_master) then + open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/normal_two_body_bi_orth', action="write") + call ezfio_set_work_empty(.False.) + write(11) normal_two_body_bi_orth_v0 + close(11) + call ezfio_set_tc_keywords_io_tc_integ('Read') + endif + + call wall_time(wallf) + print*,' Wall time for normal_two_body_bi_orth_v0 ', wallf-walli + +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 3301fcc2..d25a1f70 100644 --- a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f +++ b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f @@ -34,7 +34,8 @@ program tc_bi_ortho !call test_no_aab() !call test_no_aaa() - call test_no() + !call test_no() + call test_no_v0() end @@ -273,6 +274,52 @@ end ! --- +subroutine test_no_v0() + + implicit none + integer :: i, j, k, l + double precision :: accu, contrib, new, ref, thr, norm + + print*, ' test_no_v0 ...' + + thr = 1d-8 + + PROVIDE normal_two_body_bi_orth_v0 + PROVIDE normal_two_body_bi_orth + + accu = 0.d0 + norm = 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_v0(l,k,j,i) + + contrib = dabs(new - ref) + if(contrib .gt. thr) then + print*, ' problem on normal_two_body_bi_orth' + print*, l, k, j, i + print*, ref, new, contrib + stop + endif + + accu += contrib + norm += dabs(ref) + enddo + enddo + enddo + enddo + + print*, ' accu (%) = ', 100.d0*accu/norm + + return +end subroutine test_no + +! --- + + subroutine test_no() implicit none