diff --git a/src/ao_many_one_e_ints/ao_erf_gauss.irp.f b/src/ao_many_one_e_ints/ao_erf_gauss.irp.f index d9c35a8c..fe25e9c0 100644 --- a/src/ao_many_one_e_ints/ao_erf_gauss.irp.f +++ b/src/ao_many_one_e_ints/ao_erf_gauss.irp.f @@ -298,10 +298,16 @@ subroutine NAI_pol_x_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_cen double precision, intent(out) :: ints(3) integer :: i, j, power_Ai(3), power_Aj(3), n_pt_in, power_xA(3), m - double precision :: Ai_center(3), Aj_center(3), integral, alphai, alphaj, coef + double precision :: Ai_center(3), Aj_center(3), integral, alphai, alphaj, coef, coefi double precision, external :: NAI_pol_mult_erf_with1s + ASSERT(beta .ge. 0.d0) + if(beta .lt. 1d-10) then + call NAI_pol_x_mult_erf_ao(i_ao, j_ao, mu_in, C_center, ints) + return + endif + ints = 0.d0 if(ao_overlap_abs(j_ao,i_ao) .lt. 1.d-12) then return @@ -316,26 +322,27 @@ subroutine NAI_pol_x_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_cen n_pt_in = n_pt_max_integrals do i = 1, ao_prim_num(i_ao) - alphai = ao_expo_ordered_transp(i,i_ao) + alphai = ao_expo_ordered_transp (i,i_ao) + coefi = ao_coef_normalized_ordered_transp(i,i_ao) do m = 1, 3 - power_xA = power_Ai ! x * phi_i(r) = x * (x-Ax)**ax = (x-Ax)**(ax+1) + Ax * (x-Ax)**ax + power_xA = power_Ai power_xA(m) += 1 do j = 1, ao_prim_num(j_ao) - alphaj = ao_expo_ordered_transp(j,j_ao) - coef = ao_coef_normalized_ordered_transp(j,j_ao) * ao_coef_normalized_ordered_transp(i,i_ao) + alphaj = ao_expo_ordered_transp (j,j_ao) + coef = coefi * ao_coef_normalized_ordered_transp(j,j_ao) ! First term = (x-Ax)**(ax+1) - integral = NAI_pol_mult_erf_with1s( Ai_center, Aj_center, power_xA, power_Aj, alphai, alphaj & - , beta, b_center, c_center, n_pt_in, mu_in ) + integral = NAI_pol_mult_erf_with1s( Ai_center, Aj_center, power_xA, power_Aj, alphai, alphaj & + , beta, B_center, C_center, n_pt_in, mu_in ) ints(m) += integral * coef ! Second term = Ax * (x-Ax)**(ax) - integral = NAI_pol_mult_erf_with1s( Ai_center, Aj_center, power_Ai, power_Aj, alphai, alphaj & - , beta, b_center, c_center, n_pt_in, mu_in ) + integral = NAI_pol_mult_erf_with1s( Ai_center, Aj_center, power_Ai, power_Aj, alphai, alphaj & + , beta, B_center, C_center, n_pt_in, mu_in ) ints(m) += Ai_center(m) * integral * coef enddo diff --git a/src/ao_many_one_e_ints/ao_gaus_gauss.irp.f b/src/ao_many_one_e_ints/ao_gaus_gauss.irp.f index fadec343..c058d0d8 100644 --- a/src/ao_many_one_e_ints/ao_gaus_gauss.irp.f +++ b/src/ao_many_one_e_ints/ao_gaus_gauss.irp.f @@ -116,7 +116,7 @@ double precision function overlap_gauss_r12_ao(D_center, delta, i, j) double precision, intent(in) :: D_center(3), delta integer :: power_A(3), power_B(3), l, k - double precision :: A_center(3), B_center(3), alpha, beta, coef, analytical_j + double precision :: A_center(3), B_center(3), alpha, beta, coef, coef1, analytical_j double precision, external :: overlap_gauss_r12 @@ -133,10 +133,12 @@ double precision function overlap_gauss_r12_ao(D_center, delta, i, j) B_center(1:3) = nucl_coord(ao_nucl(j),1:3) do l = 1, ao_prim_num(i) - alpha = ao_expo_ordered_transp(l,i) + alpha = ao_expo_ordered_transp (l,i) + coef1 = ao_coef_normalized_ordered_transp(l,i) + do k = 1, ao_prim_num(j) beta = ao_expo_ordered_transp(k,j) - coef = ao_coef_normalized_ordered_transp(l,i) * ao_coef_normalized_ordered_transp(k,j) + coef = coef1 * ao_coef_normalized_ordered_transp(k,j) if(dabs(coef) .lt. 1d-12) cycle @@ -153,7 +155,9 @@ end function overlap_gauss_r12_ao double precision function overlap_gauss_r12_ao_with1s(B_center, beta, D_center, delta, i, j) BEGIN_DOC + ! ! \int dr AO_i(r) AO_j(r) e^{-beta |r-B_center^2|} e^{-delta |r-D_center|^2} + ! END_DOC implicit none @@ -161,7 +165,7 @@ double precision function overlap_gauss_r12_ao_with1s(B_center, beta, D_center, double precision, intent(in) :: B_center(3), beta, D_center(3), delta integer :: power_A1(3), power_A2(3), l, k - double precision :: A1_center(3), A2_center(3), alpha1, alpha2, coef12, analytical_j + double precision :: A1_center(3), A2_center(3), alpha1, alpha2, coef1, coef12, analytical_j double precision :: G_center(3), gama, fact_g, gama_inv double precision, external :: overlap_gauss_r12, overlap_gauss_r12_ao @@ -188,8 +192,8 @@ double precision function overlap_gauss_r12_ao_with1s(B_center, beta, D_center, fact_g = beta * delta * gama_inv * ( (B_center(1) - D_center(1)) * (B_center(1) - D_center(1)) & + (B_center(2) - D_center(2)) * (B_center(2) - D_center(2)) & + (B_center(3) - D_center(3)) * (B_center(3) - D_center(3)) ) - fact_g = dexp(-fact_g) - if(fact_g .lt. 1.d-12) return + if(fact_g .gt. 80d0) return + fact_g = dexp(-fact_g) ! --- @@ -200,11 +204,13 @@ double precision function overlap_gauss_r12_ao_with1s(B_center, beta, D_center, A2_center(1:3) = nucl_coord(ao_nucl(j),1:3) do l = 1, ao_prim_num(i) - alpha1 = ao_expo_ordered_transp(l,i) - do k = 1, ao_prim_num(j) - alpha2 = ao_expo_ordered_transp(k,j) - coef12 = fact_g * ao_coef_normalized_ordered_transp(l,i) * ao_coef_normalized_ordered_transp(k,j) + alpha1 = ao_expo_ordered_transp (l,i) + coef1 = fact_g * ao_coef_normalized_ordered_transp(l,i) + !if(dabs(coef1) .lt. 1d-12) cycle + do k = 1, ao_prim_num(j) + alpha2 = ao_expo_ordered_transp (k,j) + coef12 = coef1 * ao_coef_normalized_ordered_transp(k,j) if(dabs(coef12) .lt. 1d-12) cycle analytical_j = overlap_gauss_r12(G_center, gama, A1_center, A2_center, power_A1, power_A2, alpha1, alpha2) diff --git a/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f b/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f index 7e08bd97..50b4bf96 100644 --- a/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f +++ b/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f @@ -13,8 +13,8 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n integer :: i, j, ipoint, i_1s, i_fit double precision :: r(3), int_fit, expo_fit, coef_fit double precision :: coef, beta, B_center(3) + double precision :: tmp double precision :: wall0, wall1 - double precision, allocatable :: tmp(:,:,:) double precision, external :: overlap_gauss_r12_ao_with1s @@ -31,19 +31,17 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n !$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, & !$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, & !$OMP List_all_comb_b3_cent, int2_grad1u2_grad2u2_j1b2) - - allocate( tmp(ao_num,ao_num,n_points_final_grid) ) - tmp = 0.d0 - !$OMP DO !do ipoint = 1, 10 do ipoint = 1, n_points_final_grid + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + do i = 1, ao_num do j = i, ao_num - r(1) = final_grid_points(1,ipoint) - r(2) = final_grid_points(2,ipoint) - r(3) = final_grid_points(3,ipoint) + tmp = 0.d0 do i_1s = 1, List_all_comb_b3_size coef = List_all_comb_b3_coef (i_1s) @@ -58,29 +56,19 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n coef_fit = coef_gauss_1_erf_x_2(i_fit) int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) - tmp(j,i,ipoint) += -0.25d0 * coef * coef_fit * int_fit + tmp += -0.25d0 * coef * coef_fit * int_fit enddo enddo + + int2_grad1u2_grad2u2_j1b2(j,i,ipoint) = tmp enddo enddo enddo !$OMP END DO - - !$OMP CRITICAL - do ipoint = 1, n_points_final_grid - do i = 1, ao_num - do j = i, ao_num - int2_grad1u2_grad2u2_j1b2(j,i,ipoint) += tmp(j,i,ipoint) - enddo - enddo - enddo - !$OMP END CRITICAL - - deallocate( tmp ) !$OMP END PARALLEL do ipoint = 1, n_points_final_grid - do i = 1, ao_num + do i = 2, ao_num do j = 1, i-1 int2_grad1u2_grad2u2_j1b2(j,i,ipoint) = int2_grad1u2_grad2u2_j1b2(i,j,ipoint) enddo @@ -105,9 +93,8 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final implicit none integer :: i, j, ipoint, i_1s, i_fit double precision :: r(3), int_fit, expo_fit, coef_fit - double precision :: coef, beta, B_center(3) + double precision :: coef, beta, B_center(3), tmp double precision :: wall0, wall1 - double precision, allocatable :: tmp(:,:,:) double precision, external :: overlap_gauss_r12_ao_with1s @@ -124,19 +111,17 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final !$OMP expo_gauss_j_mu_x_2, coef_gauss_j_mu_x_2, & !$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, & !$OMP List_all_comb_b3_cent, int2_u2_j1b2) - - allocate( tmp(ao_num,ao_num,n_points_final_grid) ) - tmp = 0.d0 - !$OMP DO !do ipoint = 1, 10 do ipoint = 1, n_points_final_grid + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + do i = 1, ao_num do j = i, ao_num - r(1) = final_grid_points(1,ipoint) - r(2) = final_grid_points(2,ipoint) - r(3) = final_grid_points(3,ipoint) + tmp = 0.d0 do i_1s = 1, List_all_comb_b3_size coef = List_all_comb_b3_coef (i_1s) @@ -151,29 +136,19 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final coef_fit = coef_gauss_j_mu_x_2(i_fit) int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) - tmp(j,i,ipoint) += coef * coef_fit * int_fit + tmp += coef * coef_fit * int_fit enddo enddo + + int2_u2_j1b2(j,i,ipoint) = tmp enddo enddo enddo !$OMP END DO - - !$OMP CRITICAL - do ipoint = 1, n_points_final_grid - do i = 1, ao_num - do j = i, ao_num - int2_u2_j1b2(j,i,ipoint) += tmp(j,i,ipoint) - enddo - enddo - enddo - !$OMP END CRITICAL - - deallocate( tmp ) !$OMP END PARALLEL do ipoint = 1, n_points_final_grid - do i = 1, ao_num + do i = 2, ao_num do j = 1, i-1 int2_u2_j1b2(j,i,ipoint) = int2_u2_j1b2(i,j,ipoint) enddo @@ -187,7 +162,7 @@ END_PROVIDER ! --- -BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b, (3, ao_num, ao_num, n_points_final_grid)] +BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2, (3, ao_num, ao_num, n_points_final_grid)] BEGIN_DOC ! @@ -196,39 +171,40 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b, (3, ao_num, ao_num, n_po END_DOC implicit none - integer :: i, j, ipoint, i_1s, i_fit - double precision :: r(3), int_fit(3), expo_fit, coef_fit - double precision :: coef, beta, B_center(3) - double precision :: alpha_1s, alpha_1s_inv, centr_1s(3), expo_coef_1s, coeff_1s - double precision :: wall0, wall1 - double precision, allocatable :: tmp(:,:,:,:) + integer :: i, j, ipoint, i_1s, i_fit + double precision :: r(3), int_fit(3), expo_fit, coef_fit + double precision :: coef, beta, B_center(3), dist + double precision :: alpha_1s, alpha_1s_inv, centr_1s(3), expo_coef_1s, coef_tmp + double precision :: tmp_x, tmp_y, tmp_z + double precision :: wall0, wall1 provide mu_erf final_grid_points j1b_pen call wall_time(wall0) - int2_u_grad1u_x_j1b = 0.d0 + int2_u_grad1u_x_j1b2 = 0.d0 !$OMP PARALLEL DEFAULT (NONE) & !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & - !$OMP coef_fit, expo_fit, int_fit, tmp, alpha_1s, & - !$OMP alpha_1s_inv, centr_1s, expo_coef_1s, coeff_1s) & + !$OMP coef_fit, expo_fit, int_fit, alpha_1s, dist, & + !$OMP alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp, & + !$OMP tmp_x, tmp_y, tmp_z) & !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, & !$OMP final_grid_points, n_max_fit_slat, & !$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, & !$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, & - !$OMP List_all_comb_b3_cent, int2_u_grad1u_x_j1b) - - allocate( tmp(3,ao_num,ao_num,n_points_final_grid) ) - tmp = 0.d0 - + !$OMP List_all_comb_b3_cent, int2_u_grad1u_x_j1b2) !$OMP DO do ipoint = 1, n_points_final_grid + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + do i = 1, ao_num do j = i, ao_num - r(1) = final_grid_points(1,ipoint) - r(2) = final_grid_points(2,ipoint) - r(3) = final_grid_points(3,ipoint) + tmp_x = 0.d0 + tmp_y = 0.d0 + tmp_z = 0.d0 do i_1s = 1, List_all_comb_b3_size coef = List_all_comb_b3_coef (i_1s) @@ -236,6 +212,9 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b, (3, ao_num, ao_num, n_po B_center(1) = List_all_comb_b3_cent(1,i_1s) B_center(2) = List_all_comb_b3_cent(2,i_1s) B_center(3) = List_all_comb_b3_cent(3,i_1s) + dist = (B_center(1) - r(1)) * (B_center(1) - r(1)) & + + (B_center(2) - r(2)) * (B_center(2) - r(2)) & + + (B_center(3) - r(3)) * (B_center(3) - r(3)) do i_fit = 1, n_max_fit_slat @@ -244,56 +223,45 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b, (3, ao_num, ao_num, n_po alpha_1s = beta + expo_fit alpha_1s_inv = 1.d0 / alpha_1s + centr_1s(1) = alpha_1s_inv * (beta * B_center(1) + expo_fit * r(1)) centr_1s(2) = alpha_1s_inv * (beta * B_center(2) + expo_fit * r(2)) centr_1s(3) = alpha_1s_inv * (beta * B_center(3) + expo_fit * r(3)) - expo_coef_1s = -beta * expo_fit * alpha_1s_inv & - * ( (B_center(1) - r(1)) * (B_center(1) - r(1)) & - + (B_center(2) - r(2)) * (B_center(2) - r(2)) & - + (B_center(3) - r(3)) * (B_center(3) - r(3)) ) - if(expo_coef_1s .gt. 80.d0) cycle - coeff_1s = dexp(-expo_coef_1s) + + expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist + !if(expo_coef_1s .gt. 80.d0) cycle + coef_tmp = coef * coef_fit * dexp(-expo_coef_1s) + !if(dabs(coef_tmp) .lt. 1d-10) cycle - call NAI_pol_x_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r, int_fit) + call NAI_pol_x_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r, int_fit) - - tmp(1,j,i,ipoint) += coef * coef_fit * coeff_1s * int_fit(1) - tmp(2,j,i,ipoint) += coef * coef_fit * coeff_1s * int_fit(2) - tmp(3,j,i,ipoint) += coef * coef_fit * coeff_1s * int_fit(3) + tmp_x += coef_tmp * int_fit(1) + tmp_y += coef_tmp * int_fit(2) + tmp_z += coef_tmp * int_fit(3) enddo enddo + + int2_u_grad1u_x_j1b2(1,j,i,ipoint) = tmp_x + int2_u_grad1u_x_j1b2(2,j,i,ipoint) = tmp_y + int2_u_grad1u_x_j1b2(3,j,i,ipoint) = tmp_z enddo enddo enddo !$OMP END DO - - !$OMP CRITICAL - do ipoint = 1, n_points_final_grid - do i = 1, ao_num - do j = i, ao_num - int2_u_grad1u_x_j1b(1,j,i,ipoint) += tmp(1,j,i,ipoint) - int2_u_grad1u_x_j1b(2,j,i,ipoint) += tmp(2,j,i,ipoint) - int2_u_grad1u_x_j1b(3,j,i,ipoint) += tmp(3,j,i,ipoint) - enddo - enddo - enddo - !$OMP END CRITICAL - - deallocate( tmp ) !$OMP END PARALLEL do ipoint = 1, n_points_final_grid - do i = 1, ao_num + do i = 2, ao_num do j = 1, i-1 - int2_u_grad1u_x_j1b(1,j,i,ipoint) = int2_u_grad1u_x_j1b(1,i,j,ipoint) - int2_u_grad1u_x_j1b(2,j,i,ipoint) = int2_u_grad1u_x_j1b(2,i,j,ipoint) - int2_u_grad1u_x_j1b(3,j,i,ipoint) = int2_u_grad1u_x_j1b(3,i,j,ipoint) + int2_u_grad1u_x_j1b2(1,j,i,ipoint) = int2_u_grad1u_x_j1b2(1,i,j,ipoint) + int2_u_grad1u_x_j1b2(2,j,i,ipoint) = int2_u_grad1u_x_j1b2(2,i,j,ipoint) + int2_u_grad1u_x_j1b2(3,j,i,ipoint) = int2_u_grad1u_x_j1b2(3,i,j,ipoint) enddo enddo enddo call wall_time(wall1) - print*, ' wall time for int2_u_grad1u_x_j1b', wall1 - wall0 + print*, ' wall time for int2_u_grad1u_x_j1b2', wall1 - wall0 END_PROVIDER @@ -309,11 +277,10 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points implicit none integer :: i, j, ipoint, i_1s, i_fit - double precision :: r(3), int_fit, expo_fit, coef_fit - double precision :: coef, beta, B_center(3) - double precision :: alpha_1s, alpha_1s_inv, centr_1s(3), expo_coef_1s, coeff_1s + double precision :: r(3), int_fit, expo_fit, coef_fit, coef_tmp + double precision :: coef, beta, B_center(3), dist + double precision :: alpha_1s, alpha_1s_inv, centr_1s(3), expo_coef_1s, tmp double precision :: wall0, wall1 - double precision, allocatable :: tmp(:,:,:) double precision, external :: NAI_pol_mult_erf_ao_with1s provide mu_erf final_grid_points j1b_pen @@ -323,17 +290,13 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points !$OMP PARALLEL DEFAULT (NONE) & !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & - !$OMP coef_fit, expo_fit, int_fit, tmp, alpha_1s, & - !$OMP alpha_1s_inv, centr_1s, expo_coef_1s, coeff_1s) & + !$OMP coef_fit, expo_fit, int_fit, tmp, alpha_1s, dist, & + !$OMP alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp) & !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, & !$OMP final_grid_points, n_max_fit_slat, & !$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, & !$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, & !$OMP List_all_comb_b3_cent, int2_u_grad1u_j1b2) - - allocate( tmp(ao_num,ao_num,n_points_final_grid) ) - tmp = 0.d0 - !$OMP DO do ipoint = 1, n_points_final_grid do i = 1, ao_num @@ -342,6 +305,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points r(2) = final_grid_points(2,ipoint) r(3) = final_grid_points(3,ipoint) + tmp = 0.d0 do i_1s = 1, List_all_comb_b3_size coef = List_all_comb_b3_coef (i_1s) @@ -349,6 +313,9 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points B_center(1) = List_all_comb_b3_cent(1,i_1s) B_center(2) = List_all_comb_b3_cent(2,i_1s) B_center(3) = List_all_comb_b3_cent(3,i_1s) + dist = (B_center(1) - r(1)) * (B_center(1) - r(1)) & + + (B_center(2) - r(2)) * (B_center(2) - r(2)) & + + (B_center(3) - r(3)) * (B_center(3) - r(3)) do i_fit = 1, n_max_fit_slat @@ -360,39 +327,27 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points centr_1s(1) = alpha_1s_inv * (beta * B_center(1) + expo_fit * r(1)) centr_1s(2) = alpha_1s_inv * (beta * B_center(2) + expo_fit * r(2)) centr_1s(3) = alpha_1s_inv * (beta * B_center(3) + expo_fit * r(3)) - expo_coef_1s = -beta * expo_fit * alpha_1s_inv & - * ( (B_center(1) - r(1)) * (B_center(1) - r(1)) & - + (B_center(2) - r(2)) * (B_center(2) - r(2)) & - + (B_center(3) - r(3)) * (B_center(3) - r(3)) ) - if(expo_coef_1s .gt. 80.d0) cycle - coeff_1s = dexp(-expo_coef_1s) + + expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist + !if(expo_coef_1s .gt. 80.d0) cycle + coef_tmp = coef * coef_fit * dexp(-expo_coef_1s) + !if(dabs(coef_tmp) .lt. 1d-10) cycle int_fit = NAI_pol_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r) - - tmp(j,i,ipoint) += coef * coef_fit * coeff_1s * int_fit + tmp += coef_tmp * int_fit enddo enddo + + int2_u_grad1u_j1b2(j,i,ipoint) = tmp enddo enddo enddo !$OMP END DO - - !$OMP CRITICAL - do ipoint = 1, n_points_final_grid - do i = 1, ao_num - do j = i, ao_num - int2_u_grad1u_j1b2(j,i,ipoint) += tmp(j,i,ipoint) - enddo - enddo - enddo - !$OMP END CRITICAL - - deallocate( tmp ) !$OMP END PARALLEL do ipoint = 1, n_points_final_grid - do i = 1, ao_num + do i = 2, ao_num do j = 1, i-1 int2_u_grad1u_j1b2(j,i,ipoint) = int2_u_grad1u_j1b2(i,j,ipoint) enddo @@ -405,3 +360,4 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points END_PROVIDER ! --- + diff --git a/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f b/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f index fab50805..552e7069 100644 --- a/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f +++ b/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f @@ -10,13 +10,12 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po END_DOC implicit none - integer :: i, j, ipoint, i_1s - double precision :: r(3), int_mu, int_coulomb - double precision :: coef, beta, B_center(3) - double precision :: wall0, wall1 - double precision, allocatable :: tmp(:,:,:) - - double precision, external :: NAI_pol_mult_erf_ao_with1s + integer :: i, j, ipoint, i_1s + double precision :: r(3), int_mu, int_coulomb + double precision :: coef, beta, B_center(3) + double precision :: tmp + double precision :: wall0, wall1 + double precision, external :: NAI_pol_mult_erf_ao_with1s provide mu_erf final_grid_points j1b_pen call wall_time(wall0) @@ -28,19 +27,17 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, final_grid_points, & !$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, List_all_comb_b2_cent, & !$OMP v_ij_erf_rk_cst_mu_j1b, mu_erf) - - allocate( tmp(ao_num,ao_num,n_points_final_grid) ) - tmp = 0.d0 - !$OMP DO !do ipoint = 1, 10 do ipoint = 1, n_points_final_grid + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + do i = 1, ao_num do j = i, ao_num - r(1) = final_grid_points(1,ipoint) - r(2) = final_grid_points(2,ipoint) - r(3) = final_grid_points(3,ipoint) + tmp = 0.d0 do i_1s = 1, List_all_comb_b2_size coef = List_all_comb_b2_coef (i_1s) @@ -52,28 +49,18 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po int_mu = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r) int_coulomb = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r) - tmp(j,i,ipoint) += coef * (int_mu - int_coulomb) + tmp += coef * (int_mu - int_coulomb) enddo + + v_ij_erf_rk_cst_mu_j1b(j,i,ipoint) = tmp enddo enddo enddo !$OMP END DO - - !$OMP CRITICAL - do ipoint = 1, n_points_final_grid - do i = 1, ao_num - do j = i, ao_num - v_ij_erf_rk_cst_mu_j1b(j,i,ipoint) += tmp(j,i,ipoint) - enddo - enddo - enddo - !$OMP END CRITICAL - - deallocate( tmp ) !$OMP END PARALLEL do ipoint = 1, n_points_final_grid - do i = 1, ao_num + do i = 2, ao_num do j = 1, i-1 v_ij_erf_rk_cst_mu_j1b(j,i,ipoint) = v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) enddo @@ -123,33 +110,34 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp_j1b, (3, ao_num, ao_ END_DOC implicit none - integer :: i, j, ipoint, i_1s - double precision :: coef, beta, B_center(3), r(3), ints(3), ints_coulomb(3) - double precision :: wall0, wall1 - double precision, allocatable :: tmp(:,:,:,:) + integer :: i, j, ipoint, i_1s + double precision :: coef, beta, B_center(3), r(3), ints(3), ints_coulomb(3) + double precision :: tmp_x, tmp_y, tmp_z + double precision :: wall0, wall1 call wall_time(wall0) x_v_ij_erf_rk_cst_mu_tmp_j1b = 0.d0 !$OMP PARALLEL DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, ints, ints_coulomb, tmp) & + !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, ints, ints_coulomb, & + !$OMP tmp_x, tmp_y, tmp_z) & !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, final_grid_points,& !$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, List_all_comb_b2_cent, & !$OMP x_v_ij_erf_rk_cst_mu_tmp_j1b, mu_erf) - - allocate( tmp(3,ao_num,ao_num,n_points_final_grid) ) - tmp = 0.d0 - !$OMP DO !do ipoint = 1, 10 do ipoint = 1, n_points_final_grid + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + do i = 1, ao_num do j = i, ao_num - r(1) = final_grid_points(1,ipoint) - r(2) = final_grid_points(2,ipoint) - r(3) = final_grid_points(3,ipoint) + tmp_x = 0.d0 + tmp_y = 0.d0 + tmp_z = 0.d0 do i_1s = 1, List_all_comb_b2_size coef = List_all_comb_b2_coef (i_1s) @@ -161,32 +149,22 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp_j1b, (3, ao_num, ao_ call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, ints ) call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, ints_coulomb) - tmp(1,j,i,ipoint) += coef * (ints(1) - ints_coulomb(1)) - tmp(2,j,i,ipoint) += coef * (ints(2) - ints_coulomb(2)) - tmp(3,j,i,ipoint) += coef * (ints(3) - ints_coulomb(3)) + tmp_x += coef * (ints(1) - ints_coulomb(1)) + tmp_y += coef * (ints(2) - ints_coulomb(2)) + tmp_z += coef * (ints(3) - ints_coulomb(3)) enddo + + x_v_ij_erf_rk_cst_mu_tmp_j1b(1,j,i,ipoint) = tmp_x + x_v_ij_erf_rk_cst_mu_tmp_j1b(2,j,i,ipoint) = tmp_y + x_v_ij_erf_rk_cst_mu_tmp_j1b(3,j,i,ipoint) = tmp_z enddo enddo enddo !$OMP END DO - - !$OMP CRITICAL - do ipoint = 1, n_points_final_grid - do i = 1, ao_num - do j = i, ao_num - x_v_ij_erf_rk_cst_mu_tmp_j1b(1,j,i,ipoint) += tmp(1,j,i,ipoint) - x_v_ij_erf_rk_cst_mu_tmp_j1b(2,j,i,ipoint) += tmp(2,j,i,ipoint) - x_v_ij_erf_rk_cst_mu_tmp_j1b(3,j,i,ipoint) += tmp(3,j,i,ipoint) - enddo - enddo - enddo - !$OMP END CRITICAL - - deallocate( tmp ) !$OMP END PARALLEL do ipoint = 1, n_points_final_grid - do i = 1, ao_num + do i = 2, ao_num do j = 1, i-1 x_v_ij_erf_rk_cst_mu_tmp_j1b(1,j,i,ipoint) = x_v_ij_erf_rk_cst_mu_tmp_j1b(1,i,j,ipoint) x_v_ij_erf_rk_cst_mu_tmp_j1b(2,j,i,ipoint) = x_v_ij_erf_rk_cst_mu_tmp_j1b(2,i,j,ipoint) @@ -202,6 +180,7 @@ END_PROVIDER ! --- +! TODO analytically BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_final_grid)] BEGIN_DOC @@ -211,13 +190,13 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_ END_DOC implicit none - integer :: i, j, ipoint, i_1s, i_fit - double precision :: r(3), int_fit, expo_fit, coef_fit - double precision :: coef, beta, B_center(3) - double precision :: wall0, wall1 - double precision, allocatable :: tmp(:,:,:) + integer :: i, j, ipoint, i_1s, i_fit + double precision :: r(3), int_fit, expo_fit, coef_fit + double precision :: coef, beta, B_center(3) + double precision :: tmp + double precision :: wall0, wall1 - double precision, external :: overlap_gauss_r12_ao_with1s + double precision, external :: overlap_gauss_r12_ao_with1s provide mu_erf final_grid_points j1b_pen call wall_time(wall0) @@ -232,19 +211,17 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_ !$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, & !$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, & !$OMP List_all_comb_b2_cent, v_ij_u_cst_mu_j1b) - - allocate( tmp(ao_num,ao_num,n_points_final_grid) ) - tmp = 0.d0 - !$OMP DO !do ipoint = 1, 10 do ipoint = 1, n_points_final_grid + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + do i = 1, ao_num do j = i, ao_num - r(1) = final_grid_points(1,ipoint) - r(2) = final_grid_points(2,ipoint) - r(3) = final_grid_points(3,ipoint) + tmp = 0.d0 do i_1s = 1, List_all_comb_b2_size coef = List_all_comb_b2_coef (i_1s) @@ -259,29 +236,19 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_ coef_fit = coef_gauss_j_mu_x(i_fit) int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) - tmp(j,i,ipoint) += coef * coef_fit * int_fit + tmp += coef * coef_fit * int_fit enddo enddo + + v_ij_u_cst_mu_j1b(j,i,ipoint) = tmp enddo enddo enddo !$OMP END DO - - !$OMP CRITICAL - do ipoint = 1, n_points_final_grid - do i = 1, ao_num - do j = i, ao_num - v_ij_u_cst_mu_j1b(j,i,ipoint) += tmp(j,i,ipoint) - enddo - enddo - enddo - !$OMP END CRITICAL - - deallocate( tmp ) !$OMP END PARALLEL do ipoint = 1, n_points_final_grid - do i = 1, ao_num + do i = 2, ao_num do j = 1, i-1 v_ij_u_cst_mu_j1b(j,i,ipoint) = v_ij_u_cst_mu_j1b(i,j,ipoint) enddo @@ -294,3 +261,4 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_ END_PROVIDER ! --- + diff --git a/src/ao_many_one_e_ints/grad_related_ints.irp.f b/src/ao_many_one_e_ints/grad_related_ints.irp.f index 7b183d83..67fb0fe7 100644 --- a/src/ao_many_one_e_ints/grad_related_ints.irp.f +++ b/src/ao_many_one_e_ints/grad_related_ints.irp.f @@ -28,11 +28,12 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu, (ao_num, ao_num, n_points !$OMP SHARED (ao_num, n_points_final_grid, v_ij_erf_rk_cst_mu, final_grid_points, mu_erf) !$OMP DO SCHEDULE (dynamic) do ipoint = 1, n_points_final_grid + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + do i = 1, ao_num do j = i, ao_num - r(1) = final_grid_points(1,ipoint) - r(2) = final_grid_points(2,ipoint) - r(3) = final_grid_points(3,ipoint) int_mu = NAI_pol_mult_erf_ao(i, j, mu_erf, r) int_coulomb = NAI_pol_mult_erf_ao(i, j, 1.d+9, r) @@ -45,7 +46,7 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu, (ao_num, ao_num, n_points !$OMP END PARALLEL do ipoint = 1, n_points_final_grid - do i = 1, ao_num + do i = 2, ao_num do j = 1, i-1 v_ij_erf_rk_cst_mu(j,i,ipoint) = v_ij_erf_rk_cst_mu(i,j,ipoint) enddo @@ -53,54 +54,61 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu, (ao_num, ao_num, n_points enddo call wall_time(wall1) - print*, 'wall time for v_ij_erf_rk_cst_mu ', wall1 - wall0 + print*, ' wall time for v_ij_erf_rk_cst_mu ', wall1 - wall0 END_PROVIDER ! --- BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_transp, (n_points_final_grid, ao_num, ao_num)] - implicit none - BEGIN_DOC -! int dr phi_i(r) phi_j(r) (erf(mu(R) |r - R| - 1)/|r - R| - END_DOC - integer :: i,j,ipoint - double precision :: mu,r(3),NAI_pol_mult_erf_ao - double precision :: int_mu, int_coulomb - provide mu_erf final_grid_points - double precision :: wall0, wall1 - call wall_time(wall0) - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,ipoint,mu,r,int_mu,int_coulomb) & + + BEGIN_DOC + ! int dr phi_i(r) phi_j(r) (erf(mu(R) |r - R| - 1)/|r - R| + END_DOC + + implicit none + integer :: i, j, ipoint + double precision :: r(3) + double precision :: int_mu, int_coulomb + double precision :: wall0, wall1 + double precision :: NAI_pol_mult_erf_ao + + provide mu_erf final_grid_points + call wall_time(wall0) + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,ipoint,r,int_mu,int_coulomb) & !$OMP SHARED (ao_num,n_points_final_grid,v_ij_erf_rk_cst_mu_transp,final_grid_points,mu_erf) !$OMP DO SCHEDULE (dynamic) - do i = 1, ao_num - do j = i, ao_num - do ipoint = 1, n_points_final_grid - mu = mu_erf - r(1) = final_grid_points(1,ipoint) - r(2) = final_grid_points(2,ipoint) - r(3) = final_grid_points(3,ipoint) - int_mu = NAI_pol_mult_erf_ao(i,j,mu,r) - int_coulomb = NAI_pol_mult_erf_ao(i,j,1.d+9,r) - v_ij_erf_rk_cst_mu_transp(ipoint,j,i)= (int_mu - int_coulomb ) - enddo + do ipoint = 1, n_points_final_grid + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + + do i = 1, ao_num + do j = i, ao_num + int_mu = NAI_pol_mult_erf_ao(i, j, mu_erf, r) + int_coulomb = NAI_pol_mult_erf_ao(i, j, 1.d+9, r) + + v_ij_erf_rk_cst_mu_transp(ipoint,j,i) = int_mu - int_coulomb + enddo + enddo enddo - enddo !$OMP END DO !$OMP END PARALLEL - do i = 1, ao_num - do j = 1, i-1 - do ipoint = 1, n_points_final_grid - v_ij_erf_rk_cst_mu_transp(ipoint,j,i)= v_ij_erf_rk_cst_mu_transp(ipoint,i,j) + do i = 2, ao_num + do j = 1, i-1 + do ipoint = 1, n_points_final_grid + v_ij_erf_rk_cst_mu_transp(ipoint,j,i) = v_ij_erf_rk_cst_mu_transp(ipoint,i,j) + enddo enddo - enddo enddo - call wall_time(wall1) - print*,'wall time for v_ij_erf_rk_cst_mu_transp ',wall1 - wall0 + call wall_time(wall1) + print *, ' wall time for v_ij_erf_rk_cst_mu_transp ', wall1 - wall0 + END_PROVIDER ! --- @@ -112,30 +120,31 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp, (3, ao_num, ao_num, END_DOC implicit none - integer :: i, j, ipoint, m + integer :: i, j, ipoint double precision :: r(3), ints(3), ints_coulomb(3) double precision :: wall0, wall1 call wall_time(wall0) - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,ipoint,r,ints,m,ints_coulomb) & + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,ipoint,r,ints,ints_coulomb) & !$OMP SHARED (ao_num,n_points_final_grid,x_v_ij_erf_rk_cst_mu_tmp,final_grid_points,mu_erf) !$OMP DO SCHEDULE (dynamic) do ipoint = 1, n_points_final_grid + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + do i = 1, ao_num do j = i, ao_num - r(1) = final_grid_points(1,ipoint) - r(2) = final_grid_points(2,ipoint) - r(3) = final_grid_points(3,ipoint) call NAI_pol_x_mult_erf_ao(i, j, mu_erf, r, ints ) call NAI_pol_x_mult_erf_ao(i, j, 1.d+9 , r, ints_coulomb) - do m = 1, 3 - x_v_ij_erf_rk_cst_mu_tmp(m,j,i,ipoint) = (ints(m) - ints_coulomb(m)) - enddo + x_v_ij_erf_rk_cst_mu_tmp(1,j,i,ipoint) = ints(1) - ints_coulomb(1) + x_v_ij_erf_rk_cst_mu_tmp(2,j,i,ipoint) = ints(2) - ints_coulomb(2) + x_v_ij_erf_rk_cst_mu_tmp(3,j,i,ipoint) = ints(3) - ints_coulomb(3) enddo enddo enddo @@ -143,11 +152,11 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp, (3, ao_num, ao_num, !$OMP END PARALLEL do ipoint = 1, n_points_final_grid - do i = 1, ao_num + do i = 2, ao_num do j = 1, i-1 - do m = 1, 3 - x_v_ij_erf_rk_cst_mu_tmp(m,j,i,ipoint) = x_v_ij_erf_rk_cst_mu_tmp(m,i,j,ipoint) - enddo + x_v_ij_erf_rk_cst_mu_tmp(1,j,i,ipoint) = x_v_ij_erf_rk_cst_mu_tmp(1,i,j,ipoint) + x_v_ij_erf_rk_cst_mu_tmp(2,j,i,ipoint) = x_v_ij_erf_rk_cst_mu_tmp(2,i,j,ipoint) + x_v_ij_erf_rk_cst_mu_tmp(3,j,i,ipoint) = x_v_ij_erf_rk_cst_mu_tmp(3,i,j,ipoint) enddo enddo enddo @@ -160,208 +169,249 @@ END_PROVIDER ! --- BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu, (ao_num, ao_num,n_points_final_grid,3)] - implicit none - BEGIN_DOC -! int dr x * phi_i(r) phi_j(r) (erf(mu(R) |r - R|) - 1)/|r - R| - END_DOC - integer :: i,j,ipoint,m - double precision :: mu,r(3),ints,ints_coulomb - double precision :: wall0, wall1 - call wall_time(wall0) - do ipoint = 1, n_points_final_grid - do i = 1, ao_num - do j = 1, ao_num - do m = 1, 3 - x_v_ij_erf_rk_cst_mu(j,i,ipoint,m)= x_v_ij_erf_rk_cst_mu_tmp(m,j,i,ipoint) - enddo - enddo - enddo - enddo - call wall_time(wall1) - print*,'wall time for x_v_ij_erf_rk_cst_mu',wall1 - wall0 + BEGIN_DOC + ! int dr x * phi_i(r) phi_j(r) (erf(mu(R) |r - R|) - 1)/|r - R| + END_DOC + + implicit none + integer :: i, j, ipoint + double precision :: wall0, wall1 + + call wall_time(wall0) + + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + do j = 1, ao_num + x_v_ij_erf_rk_cst_mu(j,i,ipoint,1) = x_v_ij_erf_rk_cst_mu_tmp(1,j,i,ipoint) + x_v_ij_erf_rk_cst_mu(j,i,ipoint,2) = x_v_ij_erf_rk_cst_mu_tmp(2,j,i,ipoint) + x_v_ij_erf_rk_cst_mu(j,i,ipoint,3) = x_v_ij_erf_rk_cst_mu_tmp(3,j,i,ipoint) + enddo + enddo + enddo + + call wall_time(wall1) + print *, ' wall time for x_v_ij_erf_rk_cst_mu', wall1 - wall0 END_PROVIDER +! --- BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_transp, (ao_num, ao_num,3,n_points_final_grid)] - implicit none - BEGIN_DOC -! int dr x * phi_i(r) phi_j(r) (erf(mu(R) |r - R|) - 1)/|r - R| - END_DOC - integer :: i,j,ipoint,m - double precision :: mu,r(3),ints,ints_coulomb - double precision :: wall0, wall1 - call wall_time(wall0) - do ipoint = 1, n_points_final_grid - do m = 1, 3 - do i = 1, ao_num - do j = 1, ao_num - x_v_ij_erf_rk_cst_mu_transp(j,i,m,ipoint)= x_v_ij_erf_rk_cst_mu_tmp(m,j,i,ipoint) - enddo - enddo - enddo - enddo - call wall_time(wall1) - print*,'wall time for x_v_ij_erf_rk_cst_mu_transp',wall1 - wall0 + BEGIN_DOC + ! int dr x * phi_i(r) phi_j(r) (erf(mu(R) |r - R|) - 1)/|r - R| + END_DOC + + implicit none + integer :: i, j, ipoint + double precision :: wall0, wall1 + + call wall_time(wall0) + + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + do j = 1, ao_num + x_v_ij_erf_rk_cst_mu_transp(j,i,1,ipoint) = x_v_ij_erf_rk_cst_mu_tmp(1,j,i,ipoint) + x_v_ij_erf_rk_cst_mu_transp(j,i,2,ipoint) = x_v_ij_erf_rk_cst_mu_tmp(2,j,i,ipoint) + x_v_ij_erf_rk_cst_mu_transp(j,i,3,ipoint) = x_v_ij_erf_rk_cst_mu_tmp(3,j,i,ipoint) + enddo + enddo + enddo + + call wall_time(wall1) + print *, ' wall time for x_v_ij_erf_rk_cst_mu_transp', wall1 - wall0 END_PROVIDER +! --- + BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_transp_bis, (n_points_final_grid,ao_num, ao_num,3)] - implicit none - BEGIN_DOC -! int dr x * phi_i(r) phi_j(r) (erf(mu(R) |r - R|) - 1)/|r - R| - END_DOC - integer :: i,j,ipoint,m - double precision :: mu,r(3),ints,ints_coulomb - double precision :: wall0, wall1 - call wall_time(wall0) - do m = 1, 3 - do i = 1, ao_num - do j = 1, ao_num - do ipoint = 1, n_points_final_grid - x_v_ij_erf_rk_cst_mu_transp_bis(ipoint,j,i,m)= x_v_ij_erf_rk_cst_mu_tmp(m,j,i,ipoint) - enddo - enddo - enddo - enddo - call wall_time(wall1) - print*,'wall time for x_v_ij_erf_rk_cst_mu_transp',wall1 - wall0 + BEGIN_DOC + ! int dr x * phi_i(r) phi_j(r) (erf(mu(R) |r - R|) - 1)/|r - R| + END_DOC + + implicit none + integer :: i, j, ipoint + double precision :: wall0, wall1 + + call wall_time(wall0) + + do i = 1, ao_num + do j = 1, ao_num + do ipoint = 1, n_points_final_grid + x_v_ij_erf_rk_cst_mu_transp_bis(ipoint,j,i,1) = x_v_ij_erf_rk_cst_mu_tmp(1,j,i,ipoint) + x_v_ij_erf_rk_cst_mu_transp_bis(ipoint,j,i,2) = x_v_ij_erf_rk_cst_mu_tmp(2,j,i,ipoint) + x_v_ij_erf_rk_cst_mu_transp_bis(ipoint,j,i,3) = x_v_ij_erf_rk_cst_mu_tmp(3,j,i,ipoint) + enddo + enddo + enddo + + call wall_time(wall1) + print *, ' wall time for x_v_ij_erf_rk_cst_mu_transp_bis', wall1 - wall0 END_PROVIDER +! --- + +BEGIN_PROVIDER [ double precision, d_dx_v_ij_erf_rk_cst_mu_tmp, (3, n_points_final_grid, ao_num, ao_num)] + + BEGIN_DOC + ! d_dx_v_ij_erf_rk_cst_mu_tmp(m,R,j,i) = int dr phi_j(r)) (erf(mu(R) |r - R|) - 1)/|r - R| d/dx (phi_i(r) + ! + ! with m == 1 -> d/dx , m == 2 -> d/dy , m == 3 -> d/dz + END_DOC -BEGIN_PROVIDER [ double precision, d_dx_v_ij_erf_rk_cst_mu_tmp, (3,n_points_final_grid,ao_num, ao_num)] implicit none - BEGIN_DOC -! d_dx_v_ij_erf_rk_cst_mu_tmp(m,R,j,i) = int dr phi_j(r)) (erf(mu(R) |r - R|) - 1)/|r - R| d/dx (phi_i(r) -! -! with m == 1 -> d/dx , m == 2 -> d/dy , m == 3 -> d/dz - END_DOC - integer :: i,j,ipoint,m - double precision :: mu,r(3),ints(3),ints_coulomb(3) + integer :: i, j, ipoint + double precision :: r(3), ints(3), ints_coulomb(3) double precision :: wall0, wall1 + call wall_time(wall0) - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,ipoint,mu,r,ints,m,ints_coulomb) & + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,ipoint,r,ints,ints_coulomb) & !$OMP SHARED (ao_num,n_points_final_grid,d_dx_v_ij_erf_rk_cst_mu_tmp,final_grid_points,mu_erf) !$OMP DO SCHEDULE (dynamic) - do i = 1, ao_num - do j = 1, ao_num - do ipoint = 1, n_points_final_grid - mu = mu_erf + do ipoint = 1, n_points_final_grid r(1) = final_grid_points(1,ipoint) r(2) = final_grid_points(2,ipoint) r(3) = final_grid_points(3,ipoint) - call phi_j_erf_mu_r_dxyz_phi(j,i,mu, r, ints) - call phi_j_erf_mu_r_dxyz_phi(j,i,1.d+9, r, ints_coulomb) - do m = 1, 3 - d_dx_v_ij_erf_rk_cst_mu_tmp(m,ipoint,j,i) = ( ints(m) - ints_coulomb(m)) + + do i = 1, ao_num + do j = 1, ao_num + call phi_j_erf_mu_r_dxyz_phi(j, i, mu_erf, r, ints) + call phi_j_erf_mu_r_dxyz_phi(j, i, 1.d+9, r, ints_coulomb) + + d_dx_v_ij_erf_rk_cst_mu_tmp(1,ipoint,j,i) = ints(1) - ints_coulomb(1) + d_dx_v_ij_erf_rk_cst_mu_tmp(2,ipoint,j,i) = ints(2) - ints_coulomb(2) + d_dx_v_ij_erf_rk_cst_mu_tmp(3,ipoint,j,i) = ints(3) - ints_coulomb(3) + enddo enddo - enddo enddo - enddo !$OMP END DO !$OMP END PARALLEL - call wall_time(wall1) - print*,'wall time for d_dx_v_ij_erf_rk_cst_mu_tmp',wall1 - wall0 - + call wall_time(wall1) + print *, ' wall time for d_dx_v_ij_erf_rk_cst_mu_tmp', wall1 - wall0 END_PROVIDER -BEGIN_PROVIDER [ double precision, d_dx_v_ij_erf_rk_cst_mu, (n_points_final_grid,ao_num, ao_num,3)] - implicit none - BEGIN_DOC -! d_dx_v_ij_erf_rk_cst_mu_tmp(j,i,R,m) = int dr phi_j(r)) (erf(mu(R) |r - R|) - 1)/|r - R| d/dx (phi_i(r) -! -! with m == 1 -> d/dx , m == 2 -> d/dy , m == 3 -> d/dz - END_DOC - integer :: i,j,ipoint,m - double precision :: mu,r(3),ints,ints_coulomb - double precision :: wall0, wall1 - call wall_time(wall0) - do i = 1, ao_num - do j = 1, ao_num - do m = 1, 3 - do ipoint = 1, n_points_final_grid - d_dx_v_ij_erf_rk_cst_mu(ipoint,j,i,m)= d_dx_v_ij_erf_rk_cst_mu_tmp(m,ipoint,j,i) +! --- + +BEGIN_PROVIDER [ double precision, d_dx_v_ij_erf_rk_cst_mu, (n_points_final_grid, ao_num, ao_num, 3)] + + BEGIN_DOC + ! + ! d_dx_v_ij_erf_rk_cst_mu_tmp(j,i,R,m) = int dr phi_j(r)) (erf(mu(R) |r - R|) - 1)/|r - R| d/dx (phi_i(r) + ! + ! with m == 1 -> d/dx , m == 2 -> d/dy , m == 3 -> d/dz + ! + END_DOC + + implicit none + integer :: i, j, ipoint + double precision :: wall0, wall1 + + call wall_time(wall0) + do i = 1, ao_num + do j = 1, ao_num + do ipoint = 1, n_points_final_grid + d_dx_v_ij_erf_rk_cst_mu(ipoint,j,i,1) = d_dx_v_ij_erf_rk_cst_mu_tmp(1,ipoint,j,i) + d_dx_v_ij_erf_rk_cst_mu(ipoint,j,i,2) = d_dx_v_ij_erf_rk_cst_mu_tmp(2,ipoint,j,i) + d_dx_v_ij_erf_rk_cst_mu(ipoint,j,i,3) = d_dx_v_ij_erf_rk_cst_mu_tmp(3,ipoint,j,i) + enddo enddo - enddo enddo - enddo - call wall_time(wall1) - print*,'wall time for d_dx_v_ij_erf_rk_cst_mu',wall1 - wall0 + call wall_time(wall1) + print *, ' wall time for d_dx_v_ij_erf_rk_cst_mu', wall1 - wall0 END_PROVIDER -BEGIN_PROVIDER [ double precision, x_d_dx_v_ij_erf_rk_cst_mu_tmp, (3,n_points_final_grid,ao_num, ao_num)] - implicit none - BEGIN_DOC -! x_d_dx_v_ij_erf_rk_cst_mu_tmp(m,j,i,R) = int dr x phi_j(r)) (erf(mu(R) |r - R|) - 1)/|r - R| d/dx (phi_i(r) -! -! with m == 1 -> d/dx , m == 2 -> d/dy , m == 3 -> d/dz - END_DOC - integer :: i,j,ipoint,m - double precision :: mu,r(3),ints(3),ints_coulomb(3) - double precision :: wall0, wall1 - call wall_time(wall0) - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,ipoint,mu,r,ints,m,ints_coulomb) & +! --- + +BEGIN_PROVIDER [ double precision, x_d_dx_v_ij_erf_rk_cst_mu_tmp, (3, n_points_final_grid, ao_num, ao_num)] + + BEGIN_DOC + ! + ! x_d_dx_v_ij_erf_rk_cst_mu_tmp(m,j,i,R) = int dr x phi_j(r)) (erf(mu(R) |r - R|) - 1)/|r - R| d/dx (phi_i(r) + ! + ! with m == 1 -> d/dx , m == 2 -> d/dy , m == 3 -> d/dz + ! + END_DOC + + implicit none + integer :: i, j, ipoint + double precision :: r(3), ints(3), ints_coulomb(3) + double precision :: wall0, wall1 + + call wall_time(wall0) + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,ipoint,r,ints,ints_coulomb) & !$OMP SHARED (ao_num,n_points_final_grid,x_d_dx_v_ij_erf_rk_cst_mu_tmp,final_grid_points,mu_erf) !$OMP DO SCHEDULE (dynamic) - do i = 1, ao_num - do j = 1, ao_num - do ipoint = 1, n_points_final_grid - mu = mu_erf + do ipoint = 1, n_points_final_grid r(1) = final_grid_points(1,ipoint) r(2) = final_grid_points(2,ipoint) r(3) = final_grid_points(3,ipoint) - call phi_j_erf_mu_r_xyz_dxyz_phi(j,i,mu, r, ints) - call phi_j_erf_mu_r_xyz_dxyz_phi(j,i,1.d+9, r, ints_coulomb) - do m = 1, 3 - x_d_dx_v_ij_erf_rk_cst_mu_tmp(m,ipoint,j,i) = ( ints(m) - ints_coulomb(m)) + + do i = 1, ao_num + do j = 1, ao_num + call phi_j_erf_mu_r_xyz_dxyz_phi(j, i, mu_erf, r, ints) + call phi_j_erf_mu_r_xyz_dxyz_phi(j, i, 1.d+9, r, ints_coulomb) + + x_d_dx_v_ij_erf_rk_cst_mu_tmp(1,ipoint,j,i) = ints(1) - ints_coulomb(1) + x_d_dx_v_ij_erf_rk_cst_mu_tmp(2,ipoint,j,i) = ints(2) - ints_coulomb(2) + x_d_dx_v_ij_erf_rk_cst_mu_tmp(3,ipoint,j,i) = ints(3) - ints_coulomb(3) + enddo enddo - enddo enddo - enddo !$OMP END DO !$OMP END PARALLEL - call wall_time(wall1) - print*,'wall time for x_d_dx_v_ij_erf_rk_cst_mu_tmp',wall1 - wall0 - + call wall_time(wall1) + print *, ' wall time for x_d_dx_v_ij_erf_rk_cst_mu_tmp', wall1 - wall0 END_PROVIDER +! --- + BEGIN_PROVIDER [ double precision, x_d_dx_v_ij_erf_rk_cst_mu, (n_points_final_grid,ao_num, ao_num,3)] - implicit none - BEGIN_DOC -! x_d_dx_v_ij_erf_rk_cst_mu_tmp(j,i,R,m) = int dr x phi_j(r)) (erf(mu(R) |r - R|) - 1)/|r - R| d/dx (phi_i(r) -! -! with m == 1 -> d/dx , m == 2 -> d/dy , m == 3 -> d/dz - END_DOC - integer :: i,j,ipoint,m - double precision :: mu,r(3),ints,ints_coulomb - double precision :: wall0, wall1 - call wall_time(wall0) - do i = 1, ao_num - do j = 1, ao_num - do ipoint = 1, n_points_final_grid - do m = 1, 3 - x_d_dx_v_ij_erf_rk_cst_mu(ipoint,j,i,m)= x_d_dx_v_ij_erf_rk_cst_mu_tmp(m,ipoint,j,i) - enddo - enddo - enddo - enddo - call wall_time(wall1) - print*,'wall time for x_d_dx_v_ij_erf_rk_cst_mu',wall1 - wall0 + BEGIN_DOC + ! + ! x_d_dx_v_ij_erf_rk_cst_mu_tmp(j,i,R,m) = int dr x phi_j(r)) (erf(mu(R) |r - R|) - 1)/|r - R| d/dx (phi_i(r) + ! + ! with m == 1 -> d/dx , m == 2 -> d/dy , m == 3 -> d/dz + ! + END_DOC + + implicit none + integer :: i, j, ipoint + double precision :: wall0, wall1 + + call wall_time(wall0) + + do i = 1, ao_num + do j = 1, ao_num + do ipoint = 1, n_points_final_grid + x_d_dx_v_ij_erf_rk_cst_mu(ipoint,j,i,1) = x_d_dx_v_ij_erf_rk_cst_mu_tmp(1,ipoint,j,i) + x_d_dx_v_ij_erf_rk_cst_mu(ipoint,j,i,2) = x_d_dx_v_ij_erf_rk_cst_mu_tmp(2,ipoint,j,i) + x_d_dx_v_ij_erf_rk_cst_mu(ipoint,j,i,3) = x_d_dx_v_ij_erf_rk_cst_mu_tmp(3,ipoint,j,i) + enddo + enddo + enddo + + call wall_time(wall1) + print *, ' wall time for x_d_dx_v_ij_erf_rk_cst_mu', wall1 - wall0 END_PROVIDER +! --- + + diff --git a/src/ao_one_e_ints/pot_ao_erf_ints.irp.f b/src/ao_one_e_ints/pot_ao_erf_ints.irp.f index 1d2d8faf..263e9845 100644 --- a/src/ao_one_e_ints/pot_ao_erf_ints.irp.f +++ b/src/ao_one_e_ints/pot_ao_erf_ints.irp.f @@ -78,7 +78,7 @@ double precision function NAI_pol_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, double precision, intent(in) :: mu_in, C_center(3) integer :: i, j, power_A1(3), power_A2(3), n_pt_in - double precision :: A1_center(3), A2_center(3), alpha1, alpha2, coef12, integral + double precision :: A1_center(3), A2_center(3), alpha1, alpha2, coef12, coef1, integral double precision, external :: NAI_pol_mult_erf_with1s, NAI_pol_mult_erf_ao @@ -98,11 +98,12 @@ double precision function NAI_pol_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, NAI_pol_mult_erf_ao_with1s = 0.d0 do i = 1, ao_prim_num(i_ao) - alpha1 = ao_expo_ordered_transp(i,i_ao) + alpha1 = ao_expo_ordered_transp (i,i_ao) + coef1 = ao_coef_normalized_ordered_transp(i,i_ao) + do j = 1, ao_prim_num(j_ao) alpha2 = ao_expo_ordered_transp(j,j_ao) - - coef12 = ao_coef_normalized_ordered_transp(j,j_ao) * ao_coef_normalized_ordered_transp(i,i_ao) + coef12 = coef1 * ao_coef_normalized_ordered_transp(j,j_ao) if(dabs(coef12) .lt. 1d-14) cycle integral = NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A1, power_A2, alpha1, alpha2 & @@ -242,9 +243,9 @@ double precision function NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A A12_center(1) = (alpha1 * A1_center(1) + alpha2 * A2_center(1)) * alpha12_inv A12_center(2) = (alpha1 * A1_center(2) + alpha2 * A2_center(2)) * alpha12_inv A12_center(3) = (alpha1 * A1_center(3) + alpha2 * A2_center(3)) * alpha12_inv - dist12 = ( (A1_center(1) - A2_center(1)) * (A1_center(1) - A2_center(1)) & - + (A1_center(2) - A2_center(2)) * (A1_center(2) - A2_center(2)) & - + (A1_center(3) - A2_center(3)) * (A1_center(3) - A2_center(3)) ) + dist12 = (A1_center(1) - A2_center(1)) * (A1_center(1) - A2_center(1)) & + + (A1_center(2) - A2_center(2)) * (A1_center(2) - A2_center(2)) & + + (A1_center(3) - A2_center(3)) * (A1_center(3) - A2_center(3)) const_factor12 = dist12 * rho12 if(const_factor12 > 80.d0) then @@ -262,9 +263,9 @@ double precision function NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A P_center(1) = (alpha12 * A12_center(1) + beta * B_center(1)) * p_inv P_center(2) = (alpha12 * A12_center(2) + beta * B_center(2)) * p_inv P_center(3) = (alpha12 * A12_center(3) + beta * B_center(3)) * p_inv - dist = ( (A12_center(1) - B_center(1)) * (A12_center(1) - B_center(1)) & - + (A12_center(2) - B_center(2)) * (A12_center(2) - B_center(2)) & - + (A12_center(3) - B_center(3)) * (A12_center(3) - B_center(3)) ) + dist = (A12_center(1) - B_center(1)) * (A12_center(1) - B_center(1)) & + + (A12_center(2) - B_center(2)) * (A12_center(2) - B_center(2)) & + + (A12_center(3) - B_center(3)) * (A12_center(3) - B_center(3)) const_factor = const_factor12 + dist * rho if(const_factor > 80.d0) then @@ -272,11 +273,9 @@ double precision function NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A return endif - dist_integral = 0.d0 - do i = 1, 3 - dist_integral += (P_center(i) - C_center(i)) * (P_center(i) - C_center(i)) - enddo - + dist_integral = (P_center(1) - C_center(1)) * (P_center(1) - C_center(1)) & + + (P_center(2) - C_center(2)) * (P_center(2) - C_center(2)) & + + (P_center(3) - C_center(3)) * (P_center(3) - C_center(3)) ! --- diff --git a/src/ao_tc_eff_map/compute_ints_eff_pot.irp.f b/src/ao_tc_eff_map/compute_ints_eff_pot.irp.f index 2e7e21c0..7a567979 100644 --- a/src/ao_tc_eff_map/compute_ints_eff_pot.irp.f +++ b/src/ao_tc_eff_map/compute_ints_eff_pot.irp.f @@ -61,7 +61,6 @@ subroutine compute_ao_tc_sym_two_e_pot_jl(j, l, n_integrals, buffer_i, buffer_va integral = integral + j1b_gauss_2e_j2(i, k, j, l) endif - if(abs(integral) < thr) then cycle endif diff --git a/src/non_h_ints_mu/fit_j.irp.f b/src/ao_tc_eff_map/fit_j.irp.f similarity index 100% rename from src/non_h_ints_mu/fit_j.irp.f rename to src/ao_tc_eff_map/fit_j.irp.f diff --git a/src/ao_tc_eff_map/two_e_ints_gauss.irp.f b/src/ao_tc_eff_map/two_e_ints_gauss.irp.f index 988b0b58..51ef73a0 100644 --- a/src/ao_tc_eff_map/two_e_ints_gauss.irp.f +++ b/src/ao_tc_eff_map/two_e_ints_gauss.irp.f @@ -254,6 +254,7 @@ double precision function general_primitive_integral_gauss(dim, & rho_old = (p*q)/(p+q) prefactor = pi_3 * inv_pq_3_2 * fact_p * fact_q do i = 1, n_gauss_eff_pot ! browse the gaussians with different expo/coef + !do i = 1, n_gauss_eff_pot-1 aa = expo_gauss_eff_pot(i) c_a = coef_gauss_eff_pot(i) t_a = dsqrt( aa /(rho_old + aa) ) diff --git a/src/ao_two_e_ints/map_integrals.irp.f b/src/ao_two_e_ints/map_integrals.irp.f index de4195ba..0d34d95e 100644 --- a/src/ao_two_e_ints/map_integrals.irp.f +++ b/src/ao_two_e_ints/map_integrals.irp.f @@ -321,8 +321,9 @@ BEGIN_PROVIDER [ double precision, ao_integrals_cache, (0:64*64*64*64) ] !$OMP END PARALLEL DO END_PROVIDER +! --- -double precision function get_ao_two_e_integral(i,j,k,l,map) result(result) +double precision function get_ao_two_e_integral(i, j, k, l, map) result(result) use map_module implicit none BEGIN_DOC diff --git a/src/bi_ort_ints/biorthog_mo_for_h.irp.f b/src/bi_ort_ints/biorthog_mo_for_h.irp.f index a8e7630b..452c13f1 100644 --- a/src/bi_ort_ints/biorthog_mo_for_h.irp.f +++ b/src/bi_ort_ints/biorthog_mo_for_h.irp.f @@ -1,37 +1,6 @@ ! --- -BEGIN_PROVIDER [double precision, ao_two_e_coul, (ao_num, ao_num, ao_num, ao_num) ] - - BEGIN_DOC - ! - ! ao_two_e_coul(k,i,l,j) = ( k i | 1/r12 | l j ) = < l k | 1/r12 | j i > - ! - END_DOC - - integer :: i, j, k, l - double precision :: integral - double precision, external :: get_ao_two_e_integral - - PROVIDE ao_integrals_map - - do j = 1, ao_num - do l = 1, ao_num - do i = 1, ao_num - do k = 1, ao_num - - integral = get_ao_two_e_integral(i, j, k, l, ao_integrals_map) - - ao_two_e_coul(k,i,l,j) = integral - enddo - enddo - enddo - enddo - -END_PROVIDER - -! --- - double precision function bi_ortho_mo_coul_ints(l, k, j, i) BEGIN_DOC @@ -155,7 +124,7 @@ BEGIN_PROVIDER [double precision, mo_bi_ortho_coul_e, (mo_num, mo_num, mo_num, m do i = 1, mo_num do l = 1, mo_num do k = 1, mo_num - ! < k l | V12 | i j > (k i|l j) + ! < k l | V12 | i j > (k i|l j) mo_bi_ortho_coul_e(k,l,i,j) = mo_bi_ortho_coul_e_chemist(k,i,l,j) enddo enddo @@ -169,13 +138,14 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, mo_bi_ortho_one_e, (mo_num, mo_num)] BEGIN_DOC - ! mo_bi_ortho_one_e(k,i) = + ! + ! mo_bi_ortho_one_e(k,i) = < MO^L_k | h_c | MO^R_i > + ! END_DOC implicit none - call ao_to_mo_bi_ortho( ao_one_e_integrals, ao_num & - , mo_bi_ortho_one_e , mo_num ) + call ao_to_mo_bi_ortho(ao_one_e_integrals, ao_num, mo_bi_ortho_one_e , mo_num) END_PROVIDER diff --git a/src/bi_ort_ints/one_e_bi_ort.irp.f b/src/bi_ort_ints/one_e_bi_ort.irp.f index 5efcb637..8997991d 100644 --- a/src/bi_ort_ints/one_e_bi_ort.irp.f +++ b/src/bi_ort_ints/one_e_bi_ort.irp.f @@ -10,7 +10,7 @@ BEGIN_PROVIDER [double precision, ao_one_e_integrals_tc_tot, (ao_num,ao_num)] provide j1b_type - if(j1b_type .ne. 0) then + if( (j1b_type .eq. 1) .or. (j1b_type .eq. 2) ) then do i = 1, ao_num do j = 1, ao_num diff --git a/src/bi_ort_ints/semi_num_ints_mo.irp.f b/src/bi_ort_ints/semi_num_ints_mo.irp.f index 27fcb7de..e7c1fdd1 100644 --- a/src/bi_ort_ints/semi_num_ints_mo.irp.f +++ b/src/bi_ort_ints/semi_num_ints_mo.irp.f @@ -86,10 +86,8 @@ BEGIN_PROVIDER [ double precision, mo_x_v_ki_bi_ortho_erf_rk_cst_mu, (mo_num, mo call ao_to_mo_bi_ortho( x_v_ij_erf_rk_cst_mu_transp (1,1,1,ipoint), size(x_v_ij_erf_rk_cst_mu_transp, 1) & , mo_x_v_ki_bi_ortho_erf_rk_cst_mu(1,1,1,ipoint), size(mo_x_v_ki_bi_ortho_erf_rk_cst_mu, 1) ) - call ao_to_mo_bi_ortho( x_v_ij_erf_rk_cst_mu_transp (1,1,2,ipoint), size(x_v_ij_erf_rk_cst_mu_transp, 1) & , mo_x_v_ki_bi_ortho_erf_rk_cst_mu(1,1,2,ipoint), size(mo_x_v_ki_bi_ortho_erf_rk_cst_mu, 1) ) - call ao_to_mo_bi_ortho( x_v_ij_erf_rk_cst_mu_transp (1,1,3,ipoint), size(x_v_ij_erf_rk_cst_mu_transp, 1) & , mo_x_v_ki_bi_ortho_erf_rk_cst_mu(1,1,3,ipoint), size(mo_x_v_ki_bi_ortho_erf_rk_cst_mu, 1) ) @@ -103,7 +101,55 @@ END_PROVIDER ! --- -BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo, (3, ao_num, ao_num, n_points_final_grid)] +BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_transp, (ao_num, ao_num, 3, n_points_final_grid)] + + implicit none + integer :: i, j, ipoint + double precision :: wall0, wall1 + + call wall_time(wall0) + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + do j = 1, ao_num + int2_grad1_u12_ao_transp(j,i,1,ipoint) = int2_grad1_u12_ao(1,j,i,ipoint) + int2_grad1_u12_ao_transp(j,i,2,ipoint) = int2_grad1_u12_ao(2,j,i,ipoint) + int2_grad1_u12_ao_transp(j,i,3,ipoint) = int2_grad1_u12_ao(3,j,i,ipoint) + enddo + enddo + enddo + call wall_time(wall1) + print *, ' wall time for int2_grad1_u12_ao_transp ', wall1 - wall0 + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_transp, (mo_num, mo_num, 3, n_points_final_grid)] + + implicit none + integer :: ipoint + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint) & + !$OMP SHARED (n_points_final_grid,int2_grad1_u12_ao_transp,int2_grad1_u12_bimo_transp) + !$OMP DO SCHEDULE (dynamic) + do ipoint = 1, n_points_final_grid + call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp (1,1,1,ipoint), size(int2_grad1_u12_ao_transp , 1) & + , int2_grad1_u12_bimo_transp(1,1,1,ipoint), size(int2_grad1_u12_bimo_transp, 1) ) + call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp (1,1,2,ipoint), size(int2_grad1_u12_ao_transp , 1) & + , int2_grad1_u12_bimo_transp(1,1,2,ipoint), size(int2_grad1_u12_bimo_transp, 1) ) + call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp (1,1,3,ipoint), size(int2_grad1_u12_ao_transp , 1) & + , int2_grad1_u12_bimo_transp(1,1,3,ipoint), size(int2_grad1_u12_bimo_transp, 1) ) + enddo + !$OMP END DO + !$OMP END PARALLEL + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo, (3, mo_num, mo_num, n_points_final_grid)] BEGIN_DOC ! @@ -121,14 +167,12 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo, (3, ao_num, ao_num, n_po !$OMP DO SCHEDULE (dynamic) do ipoint = 1, n_points_final_grid - call ao_to_mo_bi_ortho( int2_grad1_u12_ao (1,1,1,ipoint), size(int2_grad1_u12_ao , 1) & - , int2_grad1_u12_bimo(1,1,1,ipoint), size(int2_grad1_u12_bimo, 1) ) - - call ao_to_mo_bi_ortho( int2_grad1_u12_ao (2,1,1,ipoint), size(int2_grad1_u12_ao , 1) & - , int2_grad1_u12_bimo(2,1,1,ipoint), size(int2_grad1_u12_bimo, 1) ) - - call ao_to_mo_bi_ortho( int2_grad1_u12_ao (3,1,1,ipoint), size(int2_grad1_u12_ao , 1) & - , int2_grad1_u12_bimo(3,1,1,ipoint), size(int2_grad1_u12_bimo, 1) ) + call ao_to_mo_bi_ortho( int2_grad1_u12_ao (1,1,1,ipoint), size(int2_grad1_u12_ao , 2) & + , int2_grad1_u12_bimo(1,1,1,ipoint), size(int2_grad1_u12_bimo, 2) ) + call ao_to_mo_bi_ortho( int2_grad1_u12_ao (2,1,1,ipoint), size(int2_grad1_u12_ao , 2) & + , int2_grad1_u12_bimo(2,1,1,ipoint), size(int2_grad1_u12_bimo, 2) ) + call ao_to_mo_bi_ortho( int2_grad1_u12_ao (3,1,1,ipoint), size(int2_grad1_u12_ao , 2) & + , int2_grad1_u12_bimo(3,1,1,ipoint), size(int2_grad1_u12_bimo, 2) ) enddo !$OMP END DO diff --git a/src/bi_ort_ints/three_body_ijm.irp.f b/src/bi_ort_ints/three_body_ijm.irp.f index 4fd85756..0e42264b 100644 --- a/src/bi_ort_ints/three_body_ijm.irp.f +++ b/src/bi_ort_ints/three_body_ijm.irp.f @@ -1,304 +1,366 @@ + +! --- + BEGIN_PROVIDER [ double precision, three_e_3_idx_direct_bi_ort, (mo_num, mo_num, mo_num)] - implicit none - BEGIN_DOC -! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the direct terms -! -! three_e_3_idx_direct_bi_ort(m,j,i) = -! -! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign - END_DOC - integer :: i,j,m - double precision :: integral, wall1, wall0 - character*(128) :: name_file - three_e_3_idx_direct_bi_ort = 0.d0 - print*,'Providing the three_e_3_idx_direct_bi_ort ...' - call wall_time(wall0) - name_file = 'six_index_tensor' - provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & + + BEGIN_DOC + ! + ! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the direct terms + ! + ! three_e_3_idx_direct_bi_ort(m,j,i) = + ! + ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign + ! + END_DOC + + implicit none + integer :: i, j, m + double precision :: integral, wall1, wall0 + + three_e_3_idx_direct_bi_ort = 0.d0 + print *, ' Providing the three_e_3_idx_direct_bi_ort ...' + call wall_time(wall0) + + provide mos_r_in_r_array_transp mos_l_in_r_array_transp + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & !$OMP PRIVATE (i,j,m,integral) & !$OMP SHARED (mo_num,three_e_3_idx_direct_bi_ort) !$OMP DO SCHEDULE (dynamic) do i = 1, mo_num - do j = 1, mo_num - do m = j, mo_num - call give_integrals_3_body_bi_ort(m,j,i,m,j,i,integral) - three_e_3_idx_direct_bi_ort(m,j,i) = -1.d0 * integral + do j = 1, mo_num + do m = j, mo_num + call give_integrals_3_body_bi_ort(m, j, i, m, j, i, integral) + three_e_3_idx_direct_bi_ort(m,j,i) = -1.d0 * integral + enddo enddo - enddo enddo !$OMP END DO !$OMP END PARALLEL - call wall_time(wall1) - print*,'wall time for three_e_3_idx_direct_bi_ort',wall1 - wall0 - + do i = 1, mo_num - do j = 1, mo_num - do m = 1, j - three_e_3_idx_direct_bi_ort(m,j,i) = three_e_3_idx_direct_bi_ort(j,m,i) + do j = 1, mo_num + do m = 1, j + three_e_3_idx_direct_bi_ort(m,j,i) = three_e_3_idx_direct_bi_ort(j,m,i) + enddo enddo - enddo enddo + call wall_time(wall1) + print *, ' wall time for three_e_3_idx_direct_bi_ort', wall1 - wall0 + END_PROVIDER +! --- + BEGIN_PROVIDER [ double precision, three_e_3_idx_cycle_1_bi_ort, (mo_num, mo_num, mo_num)] - implicit none - BEGIN_DOC -! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the first cyclic permutation -! -! three_e_3_idx_direct_bi_ort(m,j,i) = -! -! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign - END_DOC - integer :: i,j,m - double precision :: integral, wall1, wall0 - character*(128) :: name_file - three_e_3_idx_cycle_1_bi_ort = 0.d0 - print*,'Providing the three_e_3_idx_cycle_1_bi_ort ...' - call wall_time(wall0) - name_file = 'six_index_tensor' - provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & + + BEGIN_DOC + ! + ! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the first cyclic permutation + ! + ! three_e_3_idx_direct_bi_ort(m,j,i) = + ! + ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign + ! + END_DOC + + implicit none + integer :: i, j, m + double precision :: integral, wall1, wall0 + + three_e_3_idx_cycle_1_bi_ort = 0.d0 + print *, ' Providing the three_e_3_idx_cycle_1_bi_ort ...' + call wall_time(wall0) + + provide mos_r_in_r_array_transp mos_l_in_r_array_transp + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & !$OMP PRIVATE (i,j,m,integral) & !$OMP SHARED (mo_num,three_e_3_idx_cycle_1_bi_ort) !$OMP DO SCHEDULE (dynamic) do i = 1, mo_num - do j = 1, mo_num - do m = j, mo_num - call give_integrals_3_body_bi_ort(m,j,i,j,i,m,integral) - three_e_3_idx_cycle_1_bi_ort(m,j,i) = -1.d0 * integral + do j = 1, mo_num + do m = j, mo_num + call give_integrals_3_body_bi_ort(m, j, i, j, i, m, integral) + three_e_3_idx_cycle_1_bi_ort(m,j,i) = -1.d0 * integral + enddo enddo - enddo enddo !$OMP END DO !$OMP END PARALLEL - call wall_time(wall1) do i = 1, mo_num - do j = 1, mo_num - do m = 1, j - three_e_3_idx_cycle_1_bi_ort(m,j,i) = three_e_3_idx_cycle_1_bi_ort(j,m,i) + do j = 1, mo_num + do m = 1, j + three_e_3_idx_cycle_1_bi_ort(m,j,i) = three_e_3_idx_cycle_1_bi_ort(j,m,i) + enddo enddo - enddo enddo - print*,'wall time for three_e_3_idx_cycle_1_bi_ort',wall1 - wall0 + + call wall_time(wall1) + print *, ' wall time for three_e_3_idx_cycle_1_bi_ort', wall1 - wall0 END_PROVIDER +! --- + BEGIN_PROVIDER [ double precision, three_e_3_idx_cycle_2_bi_ort, (mo_num, mo_num, mo_num)] - implicit none - BEGIN_DOC -! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the second cyclic permutation -! -! three_e_3_idx_direct_bi_ort(m,j,i) = -! -! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign - END_DOC - integer :: i,j,m - double precision :: integral, wall1, wall0 - character*(128) :: name_file - three_e_3_idx_cycle_2_bi_ort = 0.d0 - print*,'Providing the three_e_3_idx_cycle_2_bi_ort ...' - call wall_time(wall0) - name_file = 'six_index_tensor' - provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & + + BEGIN_DOC + ! + ! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the second cyclic permutation + ! + ! three_e_3_idx_direct_bi_ort(m,j,i) = + ! + ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign + ! + END_DOC + + implicit none + integer :: i, j, m + double precision :: integral, wall1, wall0 + + three_e_3_idx_cycle_2_bi_ort = 0.d0 + print *, ' Providing the three_e_3_idx_cycle_2_bi_ort ...' + call wall_time(wall0) + + provide mos_r_in_r_array_transp mos_l_in_r_array_transp + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & !$OMP PRIVATE (i,j,m,integral) & !$OMP SHARED (mo_num,three_e_3_idx_cycle_2_bi_ort) !$OMP DO SCHEDULE (dynamic) do i = 1, mo_num - do j = 1, mo_num - do m = j, mo_num - call give_integrals_3_body_bi_ort(m,j,i,i,m,j,integral) - three_e_3_idx_cycle_2_bi_ort(m,j,i) = -1.d0 * integral + do j = 1, mo_num + do m = j, mo_num + call give_integrals_3_body_bi_ort(m, j, i, i, m, j, integral) + three_e_3_idx_cycle_2_bi_ort(m,j,i) = -1.d0 * integral + enddo enddo - enddo enddo !$OMP END DO !$OMP END PARALLEL - call wall_time(wall1) + do i = 1, mo_num - do j = 1, mo_num - do m = 1, j - three_e_3_idx_cycle_2_bi_ort(m,j,i) = three_e_3_idx_cycle_2_bi_ort(j,m,i) + do j = 1, mo_num + do m = 1, j + three_e_3_idx_cycle_2_bi_ort(m,j,i) = three_e_3_idx_cycle_2_bi_ort(j,m,i) + enddo enddo - enddo enddo - print*,'wall time for three_e_3_idx_cycle_2_bi_ort',wall1 - wall0 + + call wall_time(wall1) + print *, ' wall time for three_e_3_idx_cycle_2_bi_ort', wall1 - wall0 END_PROVIDER +! --- + BEGIN_PROVIDER [ double precision, three_e_3_idx_exch23_bi_ort, (mo_num, mo_num, mo_num)] - implicit none - BEGIN_DOC -! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the permutations of particle 2 and 3 -! -! three_e_3_idx_exch23_bi_ort(m,j,i) = -! -! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign - END_DOC - integer :: i,j,m - double precision :: integral, wall1, wall0 - character*(128) :: name_file - three_e_3_idx_exch23_bi_ort = 0.d0 - print*,'Providing the three_e_3_idx_exch23_bi_ort ...' - call wall_time(wall0) - name_file = 'six_index_tensor' - provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & + + BEGIN_DOC + ! + ! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the permutations of particle 2 and 3 + ! + ! three_e_3_idx_exch23_bi_ort(m,j,i) = + ! + ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign + ! + END_DOC + + implicit none + integer :: i, j, m + double precision :: integral, wall1, wall0 + + three_e_3_idx_exch23_bi_ort = 0.d0 + print*,'Providing the three_e_3_idx_exch23_bi_ort ...' + call wall_time(wall0) + + provide mos_r_in_r_array_transp mos_l_in_r_array_transp + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & !$OMP PRIVATE (i,j,m,integral) & !$OMP SHARED (mo_num,three_e_3_idx_exch23_bi_ort) !$OMP DO SCHEDULE (dynamic) do i = 1, mo_num - do j = 1, mo_num - do m = j, mo_num - call give_integrals_3_body_bi_ort(m,j,i,j,m,i,integral) - three_e_3_idx_exch23_bi_ort(m,j,i) = -1.d0 * integral + do j = 1, mo_num + do m = j, mo_num + call give_integrals_3_body_bi_ort(m, j, i, j, m, i, integral) + three_e_3_idx_exch23_bi_ort(m,j,i) = -1.d0 * integral + enddo enddo - enddo enddo !$OMP END DO !$OMP END PARALLEL + do i = 1, mo_num - do j = 1, mo_num - do m = 1, j - three_e_3_idx_exch23_bi_ort(m,j,i) = three_e_3_idx_exch23_bi_ort(j,m,i) + do j = 1, mo_num + do m = 1, j + three_e_3_idx_exch23_bi_ort(m,j,i) = three_e_3_idx_exch23_bi_ort(j,m,i) + enddo enddo - enddo enddo - call wall_time(wall1) - print*,'wall time for three_e_3_idx_exch23_bi_ort',wall1 - wall0 + + call wall_time(wall1) + print *, ' wall time for three_e_3_idx_exch23_bi_ort', wall1 - wall0 END_PROVIDER +! --- + BEGIN_PROVIDER [ double precision, three_e_3_idx_exch13_bi_ort, (mo_num, mo_num, mo_num)] - implicit none - BEGIN_DOC -! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the permutations of particle 1 and 3 -! -! three_e_3_idx_exch13_bi_ort(m,j,i) = -! -! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign - END_DOC - integer :: i,j,m - double precision :: integral, wall1, wall0 - character*(128) :: name_file - three_e_3_idx_exch13_bi_ort = 0.d0 - print*,'Providing the three_e_3_idx_exch13_bi_ort ...' - call wall_time(wall0) - name_file = 'six_index_tensor' - provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & + + BEGIN_DOC + ! + ! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the permutations of particle 1 and 3 + ! + ! three_e_3_idx_exch13_bi_ort(m,j,i) = + ! + ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign + ! + END_DOC + + implicit none + integer :: i,j,m + double precision :: integral, wall1, wall0 + + three_e_3_idx_exch13_bi_ort = 0.d0 + print *, ' Providing the three_e_3_idx_exch13_bi_ort ...' + call wall_time(wall0) + + provide mos_r_in_r_array_transp mos_l_in_r_array_transp + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & !$OMP PRIVATE (i,j,m,integral) & !$OMP SHARED (mo_num,three_e_3_idx_exch13_bi_ort) !$OMP DO SCHEDULE (dynamic) do i = 1, mo_num - do j = 1, mo_num - do m = j, mo_num - call give_integrals_3_body_bi_ort(m,j,i,i,j,m,integral) - three_e_3_idx_exch13_bi_ort(m,j,i) = -1.d0 * integral + do j = 1, mo_num + do m = j, mo_num + call give_integrals_3_body_bi_ort(m, j, i, i, j, m,integral) + three_e_3_idx_exch13_bi_ort(m,j,i) = -1.d0 * integral + enddo enddo - enddo enddo !$OMP END DO !$OMP END PARALLEL + do i = 1, mo_num - do j = 1, mo_num - do m = 1, j - three_e_3_idx_exch13_bi_ort(m,j,i) = three_e_3_idx_exch13_bi_ort(j,m,i) + do j = 1, mo_num + do m = 1, j + three_e_3_idx_exch13_bi_ort(m,j,i) = three_e_3_idx_exch13_bi_ort(j,m,i) + enddo enddo - enddo enddo - call wall_time(wall1) - print*,'wall time for three_e_3_idx_exch13_bi_ort',wall1 - wall0 + + call wall_time(wall1) + print *, ' wall time for three_e_3_idx_exch13_bi_ort', wall1 - wall0 END_PROVIDER +! --- + BEGIN_PROVIDER [ double precision, three_e_3_idx_exch12_bi_ort, (mo_num, mo_num, mo_num)] - implicit none - BEGIN_DOC -! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the permutations of particle 1 and 2 -! -! three_e_3_idx_exch12_bi_ort(m,j,i) = -! -! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign - END_DOC - integer :: i,j,m - double precision :: integral, wall1, wall0 - character*(128) :: name_file - three_e_3_idx_exch12_bi_ort = 0.d0 - print*,'Providing the three_e_3_idx_exch12_bi_ort ...' - call wall_time(wall0) - name_file = 'six_index_tensor' - provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & + + BEGIN_DOC + ! + ! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the permutations of particle 1 and 2 + ! + ! three_e_3_idx_exch12_bi_ort(m,j,i) = + ! + ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign + ! + END_DOC + + implicit none + integer :: i, j, m + double precision :: integral, wall1, wall0 + + three_e_3_idx_exch12_bi_ort = 0.d0 + print *, ' Providing the three_e_3_idx_exch12_bi_ort ...' + call wall_time(wall0) + + provide mos_r_in_r_array_transp mos_l_in_r_array_transp + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & !$OMP PRIVATE (i,j,m,integral) & !$OMP SHARED (mo_num,three_e_3_idx_exch12_bi_ort) !$OMP DO SCHEDULE (dynamic) do i = 1, mo_num - do j = 1, mo_num - do m = 1, mo_num - call give_integrals_3_body_bi_ort(m,j,i,m,i,j,integral) - three_e_3_idx_exch12_bi_ort(m,j,i) = -1.d0 * integral + do j = 1, mo_num + do m = 1, mo_num + call give_integrals_3_body_bi_ort(m, j, i, m, i, j, integral) + three_e_3_idx_exch12_bi_ort(m,j,i) = -1.d0 * integral + enddo enddo - enddo enddo !$OMP END DO !$OMP END PARALLEL - call wall_time(wall1) - print*,'wall time for three_e_3_idx_exch12_bi_ort',wall1 - wall0 + + call wall_time(wall1) + print *, ' wall time for three_e_3_idx_exch12_bi_ort', wall1 - wall0 END_PROVIDER +! --- BEGIN_PROVIDER [ double precision, three_e_3_idx_exch12_bi_ort_new, (mo_num, mo_num, mo_num)] - implicit none - BEGIN_DOC -! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the permutations of particle 1 and 2 -! -! three_e_3_idx_exch12_bi_ort_new(m,j,i) = -! -! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign - END_DOC - integer :: i,j,m - double precision :: integral, wall1, wall0 - character*(128) :: name_file - three_e_3_idx_exch12_bi_ort_new = 0.d0 - print*,'Providing the three_e_3_idx_exch12_bi_ort_new ...' - call wall_time(wall0) - name_file = 'six_index_tensor' - provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & + + BEGIN_DOC + ! + ! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the permutations of particle 1 and 2 + ! + ! three_e_3_idx_exch12_bi_ort_new(m,j,i) = + ! + ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign + ! + END_DOC + + implicit none + integer :: i, j, m + double precision :: integral, wall1, wall0 + + three_e_3_idx_exch12_bi_ort_new = 0.d0 + print *, ' Providing the three_e_3_idx_exch12_bi_ort_new ...' + call wall_time(wall0) + + provide mos_r_in_r_array_transp mos_l_in_r_array_transp + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & !$OMP PRIVATE (i,j,m,integral) & !$OMP SHARED (mo_num,three_e_3_idx_exch12_bi_ort_new) !$OMP DO SCHEDULE (dynamic) do i = 1, mo_num - do j = 1, mo_num - do m = j, mo_num - call give_integrals_3_body_bi_ort(m,j,i,m,i,j,integral) - three_e_3_idx_exch12_bi_ort_new(m,j,i) = -1.d0 * integral + do j = 1, mo_num + do m = j, mo_num + call give_integrals_3_body_bi_ort(m, j, i, m, i, j, integral) + three_e_3_idx_exch12_bi_ort_new(m,j,i) = -1.d0 * integral enddo enddo enddo !$OMP END DO !$OMP END PARALLEL + do i = 1, mo_num - do j = 1, mo_num - do m = 1, j - three_e_3_idx_exch12_bi_ort_new(m,j,i) = three_e_3_idx_exch12_bi_ort_new(j,m,i) + do j = 1, mo_num + do m = 1, j + three_e_3_idx_exch12_bi_ort_new(m,j,i) = three_e_3_idx_exch12_bi_ort_new(j,m,i) + enddo enddo - enddo enddo - call wall_time(wall1) - print*,'wall time for three_e_3_idx_exch12_bi_ort_new',wall1 - wall0 + + call wall_time(wall1) + print *, ' wall time for three_e_3_idx_exch12_bi_ort_new', wall1 - wall0 END_PROVIDER +! --- + diff --git a/src/bi_ort_ints/three_body_ijmk.irp.f b/src/bi_ort_ints/three_body_ijmk.irp.f index 40c34ddf..0d5016ce 100644 --- a/src/bi_ort_ints/three_body_ijmk.irp.f +++ b/src/bi_ort_ints/three_body_ijmk.irp.f @@ -1,228 +1,284 @@ + +! --- + BEGIN_PROVIDER [ double precision, three_e_4_idx_direct_bi_ort, (mo_num, mo_num, mo_num, mo_num)] + + BEGIN_DOC + ! + ! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs + ! + ! three_e_4_idx_direct_bi_ort(m,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO + ! + ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign + ! + END_DOC + implicit none - BEGIN_DOC -! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs -! -!three_e_4_idx_direct_bi_ort(m,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO -! -! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign - END_DOC - integer :: i,j,k,m + integer :: i, j, k, m double precision :: integral, wall1, wall0 - character*(128) :: name_file - three_e_4_idx_direct_bi_ort = 0.d0 - print*,'Providing the three_e_4_idx_direct_bi_ort ...' - call wall_time(wall0) - provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,k,m,integral) & - !$OMP SHARED (mo_num,three_e_4_idx_direct_bi_ort) - !$OMP DO SCHEDULE (dynamic) - do i = 1, mo_num + + three_e_4_idx_direct_bi_ort = 0.d0 + print *, ' Providing the three_e_4_idx_direct_bi_ort ...' + call wall_time(wall0) + + provide mos_r_in_r_array_transp mos_l_in_r_array_transp + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,m,integral) & + !$OMP SHARED (mo_num,three_e_4_idx_direct_bi_ort) + !$OMP DO SCHEDULE (dynamic) + do i = 1, mo_num do k = 1, mo_num - do j = 1, mo_num - do m = 1, mo_num - call give_integrals_3_body_bi_ort(m,j,k,m,j,i,integral) - three_e_4_idx_direct_bi_ort(m,j,k,i) = -1.d0 * integral + do j = 1, mo_num + do m = 1, mo_num + call give_integrals_3_body_bi_ort(m, j, k, m, j, i, integral) + three_e_4_idx_direct_bi_ort(m,j,k,i) = -1.d0 * integral + enddo enddo - enddo enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - call wall_time(wall1) - print*,'wall time for three_e_4_idx_direct_bi_ort',wall1 - wall0 + enddo + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(wall1) + print *, ' wall time for three_e_4_idx_direct_bi_ort', wall1 - wall0 END_PROVIDER +! --- + BEGIN_PROVIDER [ double precision, three_e_4_idx_cycle_1_bi_ort, (mo_num, mo_num, mo_num, mo_num)] - implicit none - BEGIN_DOC -! matrix element of the -L three-body operator FOR THE FIRST CYCLIC PERMUTATION TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs -! -!three_e_4_idx_cycle_1_bi_ort(m,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO -! -! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign - END_DOC - integer :: i,j,k,m - double precision :: integral, wall1, wall0 - character*(128) :: name_file - three_e_4_idx_cycle_1_bi_ort = 0.d0 - print*,'Providing the three_e_4_idx_cycle_1_bi_ort ...' - call wall_time(wall0) - provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,k,m,integral) & - !$OMP SHARED (mo_num,three_e_4_idx_cycle_1_bi_ort) - !$OMP DO SCHEDULE (dynamic) - do i = 1, mo_num + + BEGIN_DOC + ! + ! matrix element of the -L three-body operator FOR THE FIRST CYCLIC PERMUTATION TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs + ! + ! three_e_4_idx_cycle_1_bi_ort(m,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO + ! + ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign + ! + END_DOC + + implicit none + integer :: i, j, k, m + double precision :: integral, wall1, wall0 + + three_e_4_idx_cycle_1_bi_ort = 0.d0 + print *, ' Providing the three_e_4_idx_cycle_1_bi_ort ...' + call wall_time(wall0) + + provide mos_r_in_r_array_transp mos_l_in_r_array_transp + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,m,integral) & + !$OMP SHARED (mo_num,three_e_4_idx_cycle_1_bi_ort) + !$OMP DO SCHEDULE (dynamic) + do i = 1, mo_num do k = 1, mo_num - do j = 1, mo_num - do m = 1, mo_num - call give_integrals_3_body_bi_ort(m,j,k,j,i,m,integral) - three_e_4_idx_cycle_1_bi_ort(m,j,k,i) = -1.d0 * integral + do j = 1, mo_num + do m = 1, mo_num + call give_integrals_3_body_bi_ort(m, j, k, j, i, m, integral) + three_e_4_idx_cycle_1_bi_ort(m,j,k,i) = -1.d0 * integral + enddo enddo - enddo enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - call wall_time(wall1) - print*,'wall time for three_e_4_idx_cycle_1_bi_ort',wall1 - wall0 + enddo + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(wall1) + print *, ' wall time for three_e_4_idx_cycle_1_bi_ort', wall1 - wall0 END_PROVIDER +! -- BEGIN_PROVIDER [ double precision, three_e_4_idx_cycle_2_bi_ort, (mo_num, mo_num, mo_num, mo_num)] - implicit none - BEGIN_DOC -! matrix element of the -L three-body operator FOR THE FIRST CYCLIC PERMUTATION TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs -! -!three_e_4_idx_cycle_2_bi_ort(m,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO -! -! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign - END_DOC - integer :: i,j,k,m - double precision :: integral, wall1, wall0 - character*(128) :: name_file - three_e_4_idx_cycle_2_bi_ort = 0.d0 - print*,'Providing the three_e_4_idx_cycle_2_bi_ort ...' - call wall_time(wall0) - provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,k,m,integral) & - !$OMP SHARED (mo_num,three_e_4_idx_cycle_2_bi_ort) - !$OMP DO SCHEDULE (dynamic) - do i = 1, mo_num + + BEGIN_DOC + ! + ! matrix element of the -L three-body operator FOR THE FIRST CYCLIC PERMUTATION TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs + ! + ! three_e_4_idx_cycle_2_bi_ort(m,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO + ! + ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign + ! + END_DOC + + implicit none + integer :: i, j, k, m + double precision :: integral, wall1, wall0 + + three_e_4_idx_cycle_2_bi_ort = 0.d0 + print *, ' Providing the three_e_4_idx_cycle_2_bi_ort ...' + call wall_time(wall0) + + provide mos_r_in_r_array_transp mos_l_in_r_array_transp + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,m,integral) & + !$OMP SHARED (mo_num,three_e_4_idx_cycle_2_bi_ort) + !$OMP DO SCHEDULE (dynamic) + do i = 1, mo_num do k = 1, mo_num - do j = 1, mo_num - do m = 1, mo_num - call give_integrals_3_body_bi_ort(m,j,k,i,m,j,integral) - three_e_4_idx_cycle_2_bi_ort(m,j,k,i) = -1.d0 * integral + do j = 1, mo_num + do m = 1, mo_num + call give_integrals_3_body_bi_ort(m, j, k, i, m, j, integral) + three_e_4_idx_cycle_2_bi_ort(m,j,k,i) = -1.d0 * integral + enddo enddo - enddo enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - call wall_time(wall1) - print*,'wall time for three_e_4_idx_cycle_2_bi_ort',wall1 - wall0 + enddo + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(wall1) + print *, ' wall time for three_e_4_idx_cycle_2_bi_ort', wall1 - wall0 END_PROVIDER +! --- + BEGIN_PROVIDER [ double precision, three_e_4_idx_exch23_bi_ort, (mo_num, mo_num, mo_num, mo_num)] - implicit none - BEGIN_DOC -! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs -! -!three_e_4_idx_exch23_bi_ort(m,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO -! -! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign - END_DOC - integer :: i,j,k,m - double precision :: integral, wall1, wall0 - character*(128) :: name_file - three_e_4_idx_exch23_bi_ort = 0.d0 - print*,'Providing the three_e_4_idx_exch23_bi_ort ...' - call wall_time(wall0) - provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,k,m,integral) & - !$OMP SHARED (mo_num,three_e_4_idx_exch23_bi_ort) - !$OMP DO SCHEDULE (dynamic) - do i = 1, mo_num - do k = 1, mo_num - do j = 1, mo_num - do m = 1, mo_num - call give_integrals_3_body_bi_ort(m,j,k,j,m,i,integral) - three_e_4_idx_exch23_bi_ort(m,j,k,i) = -1.d0 * integral - enddo - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - call wall_time(wall1) - print*,'wall time for three_e_4_idx_exch23_bi_ort',wall1 - wall0 -END_PROVIDER + BEGIN_DOC + ! + ! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs + ! + ! three_e_4_idx_exch23_bi_ort(m,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO + ! + ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign + ! + END_DOC + + implicit none + integer :: i, j, k, m + double precision :: integral, wall1, wall0 + + three_e_4_idx_exch23_bi_ort = 0.d0 + print *, ' Providing the three_e_4_idx_exch23_bi_ort ...' + call wall_time(wall0) + + provide mos_r_in_r_array_transp mos_l_in_r_array_transp + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,m,integral) & + !$OMP SHARED (mo_num,three_e_4_idx_exch23_bi_ort) + !$OMP DO SCHEDULE (dynamic) + do i = 1, mo_num + do k = 1, mo_num + do j = 1, mo_num + do m = 1, mo_num + call give_integrals_3_body_bi_ort(m, j, k, j, m, i, integral) + three_e_4_idx_exch23_bi_ort(m,j,k,i) = -1.d0 * integral + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(wall1) + print *, ' wall time for three_e_4_idx_exch23_bi_ort', wall1 - wall0 + +END_PROVIDER + +! --- BEGIN_PROVIDER [ double precision, three_e_4_idx_exch13_bi_ort, (mo_num, mo_num, mo_num, mo_num)] - implicit none - BEGIN_DOC -! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs -! -!three_e_4_idx_exch13_bi_ort(m,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO -! -! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign - END_DOC - integer :: i,j,k,m - double precision :: integral, wall1, wall0 - character*(128) :: name_file - three_e_4_idx_exch13_bi_ort = 0.d0 - print*,'Providing the three_e_4_idx_exch13_bi_ort ...' - call wall_time(wall0) - provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,k,m,integral) & - !$OMP SHARED (mo_num,three_e_4_idx_exch13_bi_ort) - !$OMP DO SCHEDULE (dynamic) - do i = 1, mo_num + + BEGIN_DOC + ! + ! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs + ! + ! three_e_4_idx_exch13_bi_ort(m,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO + ! + ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign + END_DOC + + implicit none + integer :: i, j, k, m + double precision :: integral, wall1, wall0 + + three_e_4_idx_exch13_bi_ort = 0.d0 + print *, ' Providing the three_e_4_idx_exch13_bi_ort ...' + call wall_time(wall0) + + provide mos_r_in_r_array_transp mos_l_in_r_array_transp + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,m,integral) & + !$OMP SHARED (mo_num,three_e_4_idx_exch13_bi_ort) + !$OMP DO SCHEDULE (dynamic) + do i = 1, mo_num do k = 1, mo_num - do j = 1, mo_num - do m = 1, mo_num - call give_integrals_3_body_bi_ort(m,j,k,i,j,m,integral) - three_e_4_idx_exch13_bi_ort(m,j,k,i) = -1.d0 * integral + do j = 1, mo_num + do m = 1, mo_num + call give_integrals_3_body_bi_ort(m, j, k, i, j, m, integral) + three_e_4_idx_exch13_bi_ort(m,j,k,i) = -1.d0 * integral + enddo enddo - enddo enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - call wall_time(wall1) - print*,'wall time for three_e_4_idx_exch13_bi_ort',wall1 - wall0 + enddo + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(wall1) + print *, ' wall time for three_e_4_idx_exch13_bi_ort', wall1 - wall0 END_PROVIDER +! --- + BEGIN_PROVIDER [ double precision, three_e_4_idx_exch12_bi_ort, (mo_num, mo_num, mo_num, mo_num)] - implicit none - BEGIN_DOC -! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs -! -!three_e_4_idx_exch12_bi_ort(m,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO -! -! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign - END_DOC - integer :: i,j,k,m - double precision :: integral, wall1, wall0 - character*(128) :: name_file - three_e_4_idx_exch12_bi_ort = 0.d0 - print*,'Providing the three_e_4_idx_exch12_bi_ort ...' - call wall_time(wall0) - provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,k,m,integral) & - !$OMP SHARED (mo_num,three_e_4_idx_exch12_bi_ort) - !$OMP DO SCHEDULE (dynamic) - do i = 1, mo_num + + BEGIN_DOC + ! + ! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs + ! + ! three_e_4_idx_exch12_bi_ort(m,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO + ! + ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign + ! + END_DOC + + implicit none + integer :: i, j, k, m + double precision :: integral, wall1, wall0 + + three_e_4_idx_exch12_bi_ort = 0.d0 + print *, ' Providing the three_e_4_idx_exch12_bi_ort ...' + call wall_time(wall0) + + provide mos_r_in_r_array_transp mos_l_in_r_array_transp + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,m,integral) & + !$OMP SHARED (mo_num,three_e_4_idx_exch12_bi_ort) + !$OMP DO SCHEDULE (dynamic) + do i = 1, mo_num do k = 1, mo_num - do j = 1, mo_num - do m = 1, mo_num - call give_integrals_3_body_bi_ort(m,j,k,m,i,j,integral) - three_e_4_idx_exch12_bi_ort(m,j,k,i) = -1.d0 * integral + do j = 1, mo_num + do m = 1, mo_num + call give_integrals_3_body_bi_ort(m, j, k, m, i, j, integral) + three_e_4_idx_exch12_bi_ort(m,j,k,i) = -1.d0 * integral + enddo enddo - enddo enddo - enddo + enddo !$OMP END DO !$OMP END PARALLEL - call wall_time(wall1) - print*,'wall time for three_e_4_idx_exch12_bi_ort',wall1 - wall0 + + call wall_time(wall1) + print *, ' wall time for three_e_4_idx_exch12_bi_ort', wall1 - wall0 END_PROVIDER + +! --- + diff --git a/src/bi_ort_ints/three_body_ijmkl.irp.f b/src/bi_ort_ints/three_body_ijmkl.irp.f index 72e93955..6287c5a3 100644 --- a/src/bi_ort_ints/three_body_ijmkl.irp.f +++ b/src/bi_ort_ints/three_body_ijmkl.irp.f @@ -1,240 +1,296 @@ + +! --- + BEGIN_PROVIDER [ double precision, three_e_5_idx_direct_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)] - implicit none - BEGIN_DOC -! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs -! -!three_e_5_idx_direct_bi_ort(m,l,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO -! -! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign - END_DOC - integer :: i,j,k,m,l - double precision :: integral, wall1, wall0 - character*(128) :: name_file - three_e_5_idx_direct_bi_ort = 0.d0 - print*,'Providing the three_e_5_idx_direct_bi_ort ...' - call wall_time(wall0) - provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,k,m,l,integral) & - !$OMP SHARED (mo_num,three_e_5_idx_direct_bi_ort) - !$OMP DO SCHEDULE (dynamic) - do i = 1, mo_num + + BEGIN_DOC + ! + ! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs + ! + ! three_e_5_idx_direct_bi_ort(m,l,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO + ! + ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign + END_DOC + + implicit none + integer :: i, j, k, m, l + double precision :: integral, wall1, wall0 + + three_e_5_idx_direct_bi_ort = 0.d0 + print *, ' Providing the three_e_5_idx_direct_bi_ort ...' + call wall_time(wall0) + + provide mos_r_in_r_array_transp mos_l_in_r_array_transp + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,m,l,integral) & + !$OMP SHARED (mo_num,three_e_5_idx_direct_bi_ort) + !$OMP DO SCHEDULE (dynamic) + do i = 1, mo_num do k = 1, mo_num - do j = 1, mo_num - do l = 1, mo_num - do m = 1, mo_num - call give_integrals_3_body_bi_ort(m,l,k,m,j,i,integral) - three_e_5_idx_direct_bi_ort(m,l,j,k,i) = -1.d0 * integral - enddo + do j = 1, mo_num + do l = 1, mo_num + do m = 1, mo_num + call give_integrals_3_body_bi_ort(m, l, k, m, j, i, integral) + three_e_5_idx_direct_bi_ort(m,l,j,k,i) = -1.d0 * integral + enddo + enddo enddo - enddo enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - call wall_time(wall1) - print*,'wall time for three_e_5_idx_direct_bi_ort',wall1 - wall0 + enddo + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(wall1) + print *, ' wall time for three_e_5_idx_direct_bi_ort', wall1 - wall0 END_PROVIDER +! --- + BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_1_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)] - implicit none - BEGIN_DOC -! matrix element of the -L three-body operator FOR THE FIRST CYCLIC PERMUTATION TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs -! -!three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO -! -! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign - END_DOC - integer :: i,j,k,m,l - double precision :: integral, wall1, wall0 - character*(128) :: name_file - three_e_5_idx_cycle_1_bi_ort = 0.d0 - print*,'Providing the three_e_5_idx_cycle_1_bi_ort ...' - call wall_time(wall0) - provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,k,m,l,integral) & - !$OMP SHARED (mo_num,three_e_5_idx_cycle_1_bi_ort) - !$OMP DO SCHEDULE (dynamic) - do i = 1, mo_num + + BEGIN_DOC + ! + ! matrix element of the -L three-body operator FOR THE FIRST CYCLIC PERMUTATION TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs + ! + ! three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO + ! + ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign + ! + END_DOC + + implicit none + integer :: i, j, k, m, l + double precision :: integral, wall1, wall0 + + three_e_5_idx_cycle_1_bi_ort = 0.d0 + print *, ' Providing the three_e_5_idx_cycle_1_bi_ort ...' + call wall_time(wall0) + + provide mos_r_in_r_array_transp mos_l_in_r_array_transp + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,m,l,integral) & + !$OMP SHARED (mo_num,three_e_5_idx_cycle_1_bi_ort) + !$OMP DO SCHEDULE (dynamic) + do i = 1, mo_num do k = 1, mo_num - do j = 1, mo_num - do l = 1, mo_num - do m = 1, mo_num - call give_integrals_3_body_bi_ort(m,l,k,j,i,m,integral) - three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = -1.d0 * integral - enddo + do j = 1, mo_num + do l = 1, mo_num + do m = 1, mo_num + call give_integrals_3_body_bi_ort(m, l, k, j, i, m, integral) + three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = -1.d0 * integral + enddo + enddo enddo - enddo enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - call wall_time(wall1) - print*,'wall time for three_e_5_idx_cycle_1_bi_ort',wall1 - wall0 + enddo + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(wall1) + print *, ' wall time for three_e_5_idx_cycle_1_bi_ort', wall1 - wall0 END_PROVIDER +! --- BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_2_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)] - implicit none - BEGIN_DOC -! matrix element of the -L three-body operator FOR THE FIRST CYCLIC PERMUTATION TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs -! -!three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO -! -! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign - END_DOC - integer :: i,j,k,m,l - double precision :: integral, wall1, wall0 - character*(128) :: name_file - three_e_5_idx_cycle_2_bi_ort = 0.d0 - print*,'Providing the three_e_5_idx_cycle_2_bi_ort ...' - call wall_time(wall0) - provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,k,m,l,integral) & - !$OMP SHARED (mo_num,three_e_5_idx_cycle_2_bi_ort) - !$OMP DO SCHEDULE (dynamic) - do i = 1, mo_num + + BEGIN_DOC + ! + ! matrix element of the -L three-body operator FOR THE FIRST CYCLIC PERMUTATION TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs + ! + ! three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO + ! + ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign + ! + END_DOC + + implicit none + integer :: i, j, k, m, l + double precision :: integral, wall1, wall0 + + three_e_5_idx_cycle_2_bi_ort = 0.d0 + print *, ' Providing the three_e_5_idx_cycle_2_bi_ort ...' + call wall_time(wall0) + + provide mos_r_in_r_array_transp mos_l_in_r_array_transp + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,m,l,integral) & + !$OMP SHARED (mo_num,three_e_5_idx_cycle_2_bi_ort) + !$OMP DO SCHEDULE (dynamic) + do i = 1, mo_num do k = 1, mo_num - do j = 1, mo_num - do m = 1, mo_num - do l = 1, mo_num - call give_integrals_3_body_bi_ort(m,l,k,i,m,j,integral) - three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = -1.d0 * integral - enddo + do j = 1, mo_num + do m = 1, mo_num + do l = 1, mo_num + call give_integrals_3_body_bi_ort(m, l, k, i, m, j, integral) + three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = -1.d0 * integral + enddo + enddo enddo - enddo enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - call wall_time(wall1) - print*,'wall time for three_e_5_idx_cycle_2_bi_ort',wall1 - wall0 + enddo + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(wall1) + print *, ' wall time for three_e_5_idx_cycle_2_bi_ort', wall1 - wall0 END_PROVIDER +! --- + BEGIN_PROVIDER [ double precision, three_e_5_idx_exch23_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)] - implicit none - BEGIN_DOC -! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs -! -!three_e_5_idx_exch23_bi_ort(m,l,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO -! -! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign - END_DOC - integer :: i,j,k,m,l - double precision :: integral, wall1, wall0 - character*(128) :: name_file - three_e_5_idx_exch23_bi_ort = 0.d0 - print*,'Providing the three_e_5_idx_exch23_bi_ort ...' - call wall_time(wall0) - provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,k,m,l,integral) & - !$OMP SHARED (mo_num,three_e_5_idx_exch23_bi_ort) - !$OMP DO SCHEDULE (dynamic) - do i = 1, mo_num + + BEGIN_DOC + ! + ! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs + ! + ! three_e_5_idx_exch23_bi_ort(m,l,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO + ! + ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign + ! + END_DOC + + implicit none + integer :: i, j, k, m, l + double precision :: integral, wall1, wall0 + + three_e_5_idx_exch23_bi_ort = 0.d0 + print *, ' Providing the three_e_5_idx_exch23_bi_ort ...' + call wall_time(wall0) + + provide mos_r_in_r_array_transp mos_l_in_r_array_transp + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,m,l,integral) & + !$OMP SHARED (mo_num,three_e_5_idx_exch23_bi_ort) + !$OMP DO SCHEDULE (dynamic) + do i = 1, mo_num do k = 1, mo_num - do j = 1, mo_num - do l = 1, mo_num - do m = 1, mo_num - call give_integrals_3_body_bi_ort(m,l,k,j,m,i,integral) - three_e_5_idx_exch23_bi_ort(m,l,j,k,i) = -1.d0 * integral - enddo + do j = 1, mo_num + do l = 1, mo_num + do m = 1, mo_num + call give_integrals_3_body_bi_ort(m, l, k, j, m, i, integral) + three_e_5_idx_exch23_bi_ort(m,l,j,k,i) = -1.d0 * integral + enddo + enddo enddo - enddo enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - call wall_time(wall1) - print*,'wall time for three_e_5_idx_exch23_bi_ort',wall1 - wall0 + enddo + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(wall1) + print *, ' wall time for three_e_5_idx_exch23_bi_ort', wall1 - wall0 END_PROVIDER +! --- + BEGIN_PROVIDER [ double precision, three_e_5_idx_exch13_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)] - implicit none - BEGIN_DOC -! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs -! -!three_e_5_idx_exch13_bi_ort(m,l,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO -! -! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign - END_DOC - integer :: i,j,k,m,l - double precision :: integral, wall1, wall0 - character*(128) :: name_file - three_e_5_idx_exch13_bi_ort = 0.d0 - print*,'Providing the three_e_5_idx_exch13_bi_ort ...' - call wall_time(wall0) - provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,k,m,l,integral) & - !$OMP SHARED (mo_num,three_e_5_idx_exch13_bi_ort) - !$OMP DO SCHEDULE (dynamic) - do i = 1, mo_num + + BEGIN_DOC + ! + ! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs + ! + ! three_e_5_idx_exch13_bi_ort(m,l,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO + ! + ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign + ! + END_DOC + + implicit none + integer :: i, j, k, m, l + double precision :: integral, wall1, wall0 + + three_e_5_idx_exch13_bi_ort = 0.d0 + print *, ' Providing the three_e_5_idx_exch13_bi_ort ...' + call wall_time(wall0) + + provide mos_r_in_r_array_transp mos_l_in_r_array_transp + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,m,l,integral) & + !$OMP SHARED (mo_num,three_e_5_idx_exch13_bi_ort) + !$OMP DO SCHEDULE (dynamic) + do i = 1, mo_num do k = 1, mo_num - do j = 1, mo_num - do l = 1, mo_num - do m = 1, mo_num - call give_integrals_3_body_bi_ort(m,l,k,i,j,m,integral) - three_e_5_idx_exch13_bi_ort(m,l,j,k,i) = -1.d0 * integral - enddo + do j = 1, mo_num + do l = 1, mo_num + do m = 1, mo_num + call give_integrals_3_body_bi_ort(m, l, k, i, j, m, integral) + three_e_5_idx_exch13_bi_ort(m,l,j,k,i) = -1.d0 * integral + enddo + enddo enddo - enddo enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - call wall_time(wall1) - print*,'wall time for three_e_5_idx_exch13_bi_ort',wall1 - wall0 + enddo + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(wall1) + print *, ' wall time for three_e_5_idx_exch13_bi_ort', wall1 - wall0 END_PROVIDER +! --- + BEGIN_PROVIDER [ double precision, three_e_5_idx_exch12_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)] - implicit none - BEGIN_DOC -! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs -! -!three_e_5_idx_exch12_bi_ort(m,l,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO -! -! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign - END_DOC - integer :: i,j,k,m,l - double precision :: integral, wall1, wall0 - character*(128) :: name_file - three_e_5_idx_exch12_bi_ort = 0.d0 - print*,'Providing the three_e_5_idx_exch12_bi_ort ...' - call wall_time(wall0) - provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,k,m,l,integral) & - !$OMP SHARED (mo_num,three_e_5_idx_exch12_bi_ort) - !$OMP DO SCHEDULE (dynamic) - do i = 1, mo_num + + BEGIN_DOC + ! + ! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs + ! + ! three_e_5_idx_exch12_bi_ort(m,l,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO + ! + ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign + ! + END_DOC + + implicit none + integer :: i, j, k, m, l + double precision :: integral, wall1, wall0 + + three_e_5_idx_exch12_bi_ort = 0.d0 + print *, ' Providing the three_e_5_idx_exch12_bi_ort ...' + call wall_time(wall0) + + provide mos_r_in_r_array_transp mos_l_in_r_array_transp + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,m,l,integral) & + !$OMP SHARED (mo_num,three_e_5_idx_exch12_bi_ort) + !$OMP DO SCHEDULE (dynamic) + do i = 1, mo_num do k = 1, mo_num - do j = 1, mo_num - do l = 1, mo_num - do m = 1, mo_num - call give_integrals_3_body_bi_ort(m,l,k,m,i,j,integral) - three_e_5_idx_exch12_bi_ort(m,l,j,k,i) = -1.d0 * integral - enddo + do j = 1, mo_num + do l = 1, mo_num + do m = 1, mo_num + call give_integrals_3_body_bi_ort(m, l, k, m, i, j, integral) + three_e_5_idx_exch12_bi_ort(m,l,j,k,i) = -1.d0 * integral + enddo + enddo enddo - enddo enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - call wall_time(wall1) - print*,'wall time for three_e_5_idx_exch12_bi_ort',wall1 - wall0 + enddo + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(wall1) + print *, ' wall time for three_e_5_idx_exch12_bi_ort', wall1 - wall0 END_PROVIDER + +! --- + diff --git a/src/bi_ort_ints/three_body_ints_bi_ort.irp.f b/src/bi_ort_ints/three_body_ints_bi_ort.irp.f index 12361ace..2cca84a2 100644 --- a/src/bi_ort_ints/three_body_ints_bi_ort.irp.f +++ b/src/bi_ort_ints/three_body_ints_bi_ort.irp.f @@ -1,17 +1,24 @@ + +! --- + BEGIN_PROVIDER [ double precision, three_body_ints_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num, mo_num)] - implicit none + BEGIN_DOC ! matrix element of the -L three-body operator ! ! notice the -1 sign: in this way three_body_ints_bi_ort can be directly used to compute Slater rules :) END_DOC - integer :: i,j,k,l,m,n + + implicit none + integer :: i, j, k, l, m, n double precision :: integral, wall1, wall0 - character*(128) :: name_file - three_body_ints_bi_ort = 0.d0 - print*,'Providing the three_body_ints_bi_ort ...' - call wall_time(wall0) - name_file = 'six_index_tensor' + character*(128) :: name_file + + three_body_ints_bi_ort = 0.d0 + print*,'Providing the three_body_ints_bi_ort ...' + call wall_time(wall0) + name_file = 'six_index_tensor' + ! if(read_three_body_ints_bi_ort)then ! call read_fcidump_3_tc(three_body_ints_bi_ort) ! else @@ -19,32 +26,37 @@ BEGIN_PROVIDER [ double precision, three_body_ints_bi_ort, (mo_num, mo_num, mo_n ! print*,'Reading three_body_ints_bi_ort from disk ...' ! call read_array_6_index_tensor(mo_num,three_body_ints_bi_ort,name_file) ! else - provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,k,l,m,n,integral) & - !$OMP SHARED (mo_num,three_body_ints_bi_ort) - !$OMP DO SCHEDULE (dynamic) - do i = 1, mo_num + + !provide x_W_ki_bi_ortho_erf_rk + provide mos_r_in_r_array_transp mos_l_in_r_array_transp + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,l,m,n,integral) & + !$OMP SHARED (mo_num,three_body_ints_bi_ort) + !$OMP DO SCHEDULE (dynamic) + do i = 1, mo_num do j = 1, mo_num - do m = 1, mo_num - do k = 1, mo_num - do l = 1, mo_num - do n = 1, mo_num - call give_integrals_3_body_bi_ort(n,l,k,m,j,i,integral) - three_body_ints_bi_ort(n,l,k,m,j,i) = -1.d0 * integral + do m = 1, mo_num + do k = 1, mo_num + do l = 1, mo_num + do n = 1, mo_num + call give_integrals_3_body_bi_ort(n, l, k, m, j, i, integral) + + three_body_ints_bi_ort(n,l,k,m,j,i) = -1.d0 * integral + enddo + enddo enddo - enddo enddo - enddo enddo - enddo - !$OMP END DO - !$OMP END PARALLEL + enddo + !$OMP END DO + !$OMP END PARALLEL ! endif ! endif - call wall_time(wall1) - print*,'wall time for three_body_ints_bi_ort',wall1 - wall0 + + call wall_time(wall1) + print *, ' wall time for three_body_ints_bi_ort', wall1 - wall0 ! if(write_three_body_ints_bi_ort)then ! print*,'Writing three_body_ints_bi_ort on disk ...' ! call write_array_6_index_tensor(mo_num,three_body_ints_bi_ort,name_file) @@ -64,7 +76,7 @@ subroutine give_integrals_3_body_bi_ort(n, l, k, m, j, i, integral) END_DOC implicit none - integer, intent(in) :: n,l,k,m,j,i + integer, intent(in) :: n, l, k, m, j, i double precision, intent(out) :: integral integer :: ipoint double precision :: weight @@ -86,18 +98,31 @@ subroutine give_integrals_3_body_bi_ort(n, l, k, m, j, i, integral) ! + x_W_ki_bi_ortho_erf_rk(ipoint,2,l,j) * x_W_ki_bi_ortho_erf_rk(ipoint,2,k,i) & ! + x_W_ki_bi_ortho_erf_rk(ipoint,3,l,j) * x_W_ki_bi_ortho_erf_rk(ipoint,3,k,i) ) - integral += weight * mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) & - * ( int2_grad1_u12_bimo(1,ipoint,n,m) * int2_grad1_u12_bimo(1,ipoint,l,j) & - + int2_grad1_u12_bimo(2,ipoint,n,m) * int2_grad1_u12_bimo(2,ipoint,l,j) & - + int2_grad1_u12_bimo(3,ipoint,n,m) * int2_grad1_u12_bimo(3,ipoint,l,j) ) - integral += weight * mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) & - * ( int2_grad1_u12_bimo(1,ipoint,n,m) * int2_grad1_u12_bimo(1,ipoint,k,i) & - + int2_grad1_u12_bimo(2,ipoint,n,m) * int2_grad1_u12_bimo(2,ipoint,k,i) & - + int2_grad1_u12_bimo(3,ipoint,n,m) * int2_grad1_u12_bimo(3,ipoint,k,i) ) - integral += weight * mos_l_in_r_array_transp(ipoint,n) * mos_r_in_r_array_transp(ipoint,m) & - * ( int2_grad1_u12_bimo(1,ipoint,l,j) * int2_grad1_u12_bimo(1,ipoint,k,i) & - + int2_grad1_u12_bimo(2,ipoint,l,j) * int2_grad1_u12_bimo(2,ipoint,k,i) & - + int2_grad1_u12_bimo(3,ipoint,l,j) * int2_grad1_u12_bimo(3,ipoint,k,i) ) +! integral += weight * mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) & +! * ( int2_grad1_u12_bimo(1,n,m,ipoint) * int2_grad1_u12_bimo(1,l,j,ipoint) & +! + int2_grad1_u12_bimo(2,n,m,ipoint) * int2_grad1_u12_bimo(2,l,j,ipoint) & +! + int2_grad1_u12_bimo(3,n,m,ipoint) * int2_grad1_u12_bimo(3,l,j,ipoint) ) +! integral += weight * mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) & +! * ( int2_grad1_u12_bimo(1,n,m,ipoint) * int2_grad1_u12_bimo(1,k,i,ipoint) & +! + int2_grad1_u12_bimo(2,n,m,ipoint) * int2_grad1_u12_bimo(2,k,i,ipoint) & +! + int2_grad1_u12_bimo(3,n,m,ipoint) * int2_grad1_u12_bimo(3,k,i,ipoint) ) +! integral += weight * mos_l_in_r_array_transp(ipoint,n) * mos_r_in_r_array_transp(ipoint,m) & +! * ( int2_grad1_u12_bimo(1,l,j,ipoint) * int2_grad1_u12_bimo(1,k,i,ipoint) & +! + int2_grad1_u12_bimo(2,l,j,ipoint) * int2_grad1_u12_bimo(2,k,i,ipoint) & +! + int2_grad1_u12_bimo(3,l,j,ipoint) * int2_grad1_u12_bimo(3,k,i,ipoint) ) + + integral += weight * mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) & + * ( int2_grad1_u12_bimo_transp(n,m,1,ipoint) * int2_grad1_u12_bimo_transp(l,j,1,ipoint) & + + int2_grad1_u12_bimo_transp(n,m,2,ipoint) * int2_grad1_u12_bimo_transp(l,j,2,ipoint) & + + int2_grad1_u12_bimo_transp(n,m,3,ipoint) * int2_grad1_u12_bimo_transp(l,j,3,ipoint) ) + integral += weight * mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) & + * ( int2_grad1_u12_bimo_transp(n,m,1,ipoint) * int2_grad1_u12_bimo_transp(k,i,1,ipoint) & + + int2_grad1_u12_bimo_transp(n,m,2,ipoint) * int2_grad1_u12_bimo_transp(k,i,2,ipoint) & + + int2_grad1_u12_bimo_transp(n,m,3,ipoint) * int2_grad1_u12_bimo_transp(k,i,3,ipoint) ) + integral += weight * mos_l_in_r_array_transp(ipoint,n) * mos_r_in_r_array_transp(ipoint,m) & + * ( int2_grad1_u12_bimo_transp(l,j,1,ipoint) * int2_grad1_u12_bimo_transp(k,i,1,ipoint) & + + int2_grad1_u12_bimo_transp(l,j,2,ipoint) * int2_grad1_u12_bimo_transp(k,i,2,ipoint) & + + int2_grad1_u12_bimo_transp(l,j,3,ipoint) * int2_grad1_u12_bimo_transp(k,i,3,ipoint) ) enddo diff --git a/src/bi_ort_ints/total_twoe_pot.irp.f b/src/bi_ort_ints/total_twoe_pot.irp.f index 72ded7cf..89f46a05 100644 --- a/src/bi_ort_ints/total_twoe_pot.irp.f +++ b/src/bi_ort_ints/total_twoe_pot.irp.f @@ -17,23 +17,42 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n double precision :: integral_sym, integral_nsym double precision, external :: get_ao_tc_sym_two_e_pot - PROVIDE ao_tc_sym_two_e_pot_in_map + provide j1b_type - do j = 1, ao_num - do l = 1, ao_num - do i = 1, ao_num - do k = 1, ao_num + if(j1b_type .eq. 3) then - integral_sym = get_ao_tc_sym_two_e_pot(i,j,k,l,ao_tc_sym_two_e_pot_map) - - ! ao_non_hermit_term_chemist(k,i,l,j) = < k l | [erf( mu r12) - 1] d/d_r12 | i j > on the AO basis - integral_nsym = ao_non_hermit_term_chemist(k,i,l,j) - - ao_two_e_tc_tot(k,i,l,j) = integral_sym + integral_nsym + do j = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do k = 1, ao_num + ao_two_e_tc_tot(k,i,l,j) = ao_tc_int_chemist(k,i,l,j) + !write(222,*) ao_two_e_tc_tot(k,i,l,j) + enddo enddo enddo enddo - enddo + + else + + PROVIDE ao_tc_sym_two_e_pot_in_map + + do j = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do k = 1, ao_num + + integral_sym = get_ao_tc_sym_two_e_pot(i, j, k, l, ao_tc_sym_two_e_pot_map) + ! ao_non_hermit_term_chemist(k,i,l,j) = < k l | [erf( mu r12) - 1] d/d_r12 | i j > on the AO basis + integral_nsym = ao_non_hermit_term_chemist(k,i,l,j) + + ao_two_e_tc_tot(k,i,l,j) = integral_sym + integral_nsym + !write(111,*) ao_two_e_tc_tot(k,i,l,j) + enddo + enddo + enddo + enddo + + endif END_PROVIDER @@ -42,9 +61,11 @@ END_PROVIDER double precision function bi_ortho_mo_ints(l, k, j, i) BEGIN_DOC + ! ! ! ! WARNING :: very naive, super slow, only used to DEBUG. + ! END_DOC implicit none diff --git a/src/bi_ortho_mos/mos_rl.irp.f b/src/bi_ortho_mos/mos_rl.irp.f index b6e93c17..034a436e 100644 --- a/src/bi_ortho_mos/mos_rl.irp.f +++ b/src/bi_ortho_mos/mos_rl.irp.f @@ -1,33 +1,37 @@ + +! --- + subroutine ao_to_mo_bi_ortho(A_ao, LDA_ao, A_mo, LDA_mo) BEGIN_DOC + ! ! Transform A from the |AO| basis to the BI ORTHONORMAL MOS ! ! $C_L^\dagger.A_{ao}.C_R$ where C_L and C_R are the LEFT and RIGHT MO coefs + ! END_DOC implicit none - integer, intent(in) :: LDA_ao,LDA_mo + integer, intent(in) :: LDA_ao, LDA_mo double precision, intent(in) :: A_ao(LDA_ao,ao_num) double precision, intent(out) :: A_mo(LDA_mo,mo_num) double precision, allocatable :: T(:,:) allocate ( T(ao_num,mo_num) ) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T - integer :: i,j,p,q - call dgemm('N', 'N', ao_num, mo_num, ao_num, & - 1.d0, A_ao, LDA_ao, & - mo_r_coef, size(mo_r_coef, 1), & - 0.d0, T, size(T, 1)) + ! T = A_ao x mo_r_coef + call dgemm( 'N', 'N', ao_num, mo_num, ao_num, 1.d0 & + , A_ao, LDA_ao, mo_r_coef, size(mo_r_coef, 1) & + , 0.d0, T, size(T, 1) ) - call dgemm('T', 'N', mo_num, mo_num, ao_num, & - 1.d0, mo_l_coef, size(mo_l_coef, 1), & - T, ao_num, & - 0.d0, A_mo, size(A_mo, 1)) + ! A_mo = mo_l_coef.T x T + call dgemm( 'T', 'N', mo_num, mo_num, ao_num, 1.d0 & + , mo_l_coef, size(mo_l_coef, 1), T, size(T, 1) & + , 0.d0, A_mo, LDA_mo ) ! call restore_symmetry(mo_num,mo_num,A_mo,size(A_mo,1),1.d-12) - deallocate(T) + deallocate(T) end subroutine ao_to_mo_bi_ortho @@ -131,7 +135,7 @@ BEGIN_PROVIDER [ double precision, mo_l_coef, (ao_num, mo_num) ] IRP_ENDIF else - print*, 'mo_r_coef are mo_coef' + print*, 'mo_l_coef are mo_coef' do i = 1, mo_num do j = 1, ao_num mo_l_coef(j,i) = mo_coef(j,i) diff --git a/src/non_h_ints_mu/debug_fit.irp.f b/src/non_h_ints_mu/debug_fit.irp.f new file mode 100644 index 00000000..af441335 --- /dev/null +++ b/src/non_h_ints_mu/debug_fit.irp.f @@ -0,0 +1,512 @@ + +! -- + +program debug_fit + + implicit none + + my_grid_becke = .True. + + my_n_pt_r_grid = 30 + my_n_pt_a_grid = 50 + !my_n_pt_r_grid = 100 + !my_n_pt_a_grid = 170 + !my_n_pt_r_grid = 150 + !my_n_pt_a_grid = 194 + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + + PROVIDE mu_erf j1b_pen + + !call test_j1b_nucl() + call test_grad_j1b_nucl() + !call test_lapl_j1b_nucl() + + !call test_list_b2() + !call test_list_b3() + + call test_fit_u() + !call test_fit_u2() + !call test_fit_ugradu() + +end + +! --- + +subroutine test_j1b_nucl() + + implicit none + integer :: ipoint + double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz + double precision :: r(3) + double precision, external :: j1b_nucl + + print*, ' test_j1b_nucl ...' + + PROVIDE v_1b + + eps_ij = 1d-7 + acc_tot = 0.d0 + normalz = 0.d0 + + do ipoint = 1, n_points_final_grid + + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + + i_exc = v_1b(ipoint) + i_num = j1b_nucl(r) + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in v_1b on', ipoint + print *, ' analyt = ', i_exc + print *, ' numeri = ', i_num + print *, ' diff = ', acc_ij + endif + + acc_tot += acc_ij + normalz += dabs(i_num) + enddo + + print*, ' acc_tot = ', acc_tot + print*, ' normalz = ', normalz + + return +end subroutine test_j1b_nucl + +! --- + +subroutine test_grad_j1b_nucl() + + implicit none + integer :: ipoint + double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz + double precision :: r(3) + double precision, external :: grad_x_j1b_nucl + double precision, external :: grad_y_j1b_nucl + double precision, external :: grad_z_j1b_nucl + + print*, ' test_grad_j1b_nucl ...' + + PROVIDE v_1b_grad + + eps_ij = 1d-7 + acc_tot = 0.d0 + normalz = 0.d0 + + do ipoint = 1, n_points_final_grid + + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + + i_exc = v_1b_grad(1,ipoint) + i_num = grad_x_j1b_nucl(r) + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in x of v_1b_grad on', ipoint + print *, ' analyt = ', i_exc + print *, ' numeri = ', i_num + print *, ' diff = ', acc_ij + endif + + i_exc = v_1b_grad(2,ipoint) + i_num = grad_y_j1b_nucl(r) + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in y of v_1b_grad on', ipoint + print *, ' analyt = ', i_exc + print *, ' numeri = ', i_num + print *, ' diff = ', acc_ij + endif + + i_exc = v_1b_grad(3,ipoint) + i_num = grad_z_j1b_nucl(r) + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in z of v_1b_grad on', ipoint + print *, ' analyt = ', i_exc + print *, ' numeri = ', i_num + print *, ' diff = ', acc_ij + endif + + acc_tot += acc_ij + normalz += dabs(i_num) + enddo + + print*, ' acc_tot = ', acc_tot + print*, ' normalz = ', normalz + + return +end subroutine test_grad_j1b_nucl + +! --- + +subroutine test_lapl_j1b_nucl() + + implicit none + integer :: ipoint + double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz + double precision :: r(3) + double precision, external :: lapl_j1b_nucl + + print*, ' test_lapl_j1b_nucl ...' + + PROVIDE v_1b_lapl + + eps_ij = 1d-5 + acc_tot = 0.d0 + normalz = 0.d0 + + do ipoint = 1, n_points_final_grid + + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + + i_exc = v_1b_lapl(ipoint) + i_num = lapl_j1b_nucl(r) + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in v_1b_lapl on', ipoint + print *, ' analyt = ', i_exc + print *, ' numeri = ', i_num + print *, ' diff = ', acc_ij + endif + + acc_tot += acc_ij + normalz += dabs(i_num) + enddo + + print*, ' acc_tot = ', acc_tot + print*, ' normalz = ', normalz + + return +end subroutine test_lapl_j1b_nucl + +! --- + +subroutine test_list_b2() + + implicit none + integer :: ipoint + double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz + double precision :: r(3) + double precision, external :: j1b_nucl + + print*, ' test_list_b2 ...' + + PROVIDE v_1b_list_b2 + + eps_ij = 1d-7 + acc_tot = 0.d0 + normalz = 0.d0 + + do ipoint = 1, n_points_final_grid + + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + + i_exc = v_1b_list_b2(ipoint) + i_num = j1b_nucl(r) + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in list_b2 on', ipoint + print *, ' analyt = ', i_exc + print *, ' numeri = ', i_num + print *, ' diff = ', acc_ij + endif + + acc_tot += acc_ij + normalz += dabs(i_num) + enddo + + print*, ' acc_tot = ', acc_tot + print*, ' normalz = ', normalz + + return +end subroutine test_list_b2 + +! --- + +subroutine test_list_b3() + + implicit none + integer :: ipoint + double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_tmp, i_num, normalz + double precision :: r(3) + double precision, external :: j1b_nucl + + print*, ' test_list_b3 ...' + + PROVIDE v_1b_list_b3 + + eps_ij = 1d-7 + acc_tot = 0.d0 + normalz = 0.d0 + + do ipoint = 1, n_points_final_grid + + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + + i_exc = v_1b_list_b3(ipoint) + i_tmp = j1b_nucl(r) + i_num = i_tmp * i_tmp + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in list_b3 on', ipoint + print *, ' analyt = ', i_exc + print *, ' numeri = ', i_num + print *, ' diff = ', acc_ij + endif + + acc_tot += acc_ij + normalz += dabs(i_num) + enddo + + print*, ' acc_tot = ', acc_tot + print*, ' normalz = ', normalz + + return +end subroutine test_list_b3 + +! --- + +subroutine test_fit_ugradu() + + implicit none + + integer :: jpoint, ipoint, i + double precision :: i_exc, i_fit, i_num, x2, tmp, dx, dy, dz + double precision :: r1(3), r2(3), grad(3) + double precision :: eps_ij, acc_tot, acc_ij, normalz, coef, expo + + double precision, external :: j12_mu + + print*, ' test_fit_ugradu ...' + + eps_ij = 1d-3 + + do jpoint = 1, n_points_final_grid + r2(1) = final_grid_points(1,jpoint) + r2(2) = final_grid_points(2,jpoint) + r2(3) = final_grid_points(3,jpoint) + + acc_tot = 0.d0 + normalz = 0.d0 + do ipoint = 1, n_points_final_grid + r1(1) = final_grid_points(1,ipoint) + r1(2) = final_grid_points(2,ipoint) + r1(3) = final_grid_points(3,ipoint) + + dx = r1(1) - r2(1) + dy = r1(2) - r2(2) + dz = r1(3) - r2(3) + x2 = dx * dx + dy * dy + dz * dz + if(x2 .lt. 1d-10) cycle + + i_fit = 0.d0 + do i = 1, n_max_fit_slat + expo = expo_gauss_j_mu_1_erf(i) + coef = coef_gauss_j_mu_1_erf(i) + i_fit += coef * dexp(-expo*x2) + enddo + i_fit = i_fit / dsqrt(x2) + + tmp = j12_mu(r1, r2) + call grad1_j12_mu_exc(r1, r2, grad) + + ! --- + + i_exc = tmp * grad(1) + i_num = i_fit * dx + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem on x in test_fit_ugradu on', ipoint + print *, ' analyt = ', i_exc + print *, ' numeri = ', i_num + print *, ' diff = ', acc_ij + endif + acc_tot += acc_ij + normalz += dabs(i_exc) + + ! --- + + i_exc = tmp * grad(2) + i_num = i_fit * dy + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem on y in test_fit_ugradu on', ipoint + print *, ' analyt = ', i_exc + print *, ' numeri = ', i_num + print *, ' diff = ', acc_ij + endif + acc_tot += acc_ij + normalz += dabs(i_exc) + + ! --- + + i_exc = tmp * grad(3) + i_num = i_fit * dz + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem on z in test_fit_ugradu on', ipoint + print *, ' analyt = ', i_exc + print *, ' numeri = ', i_num + print *, ' diff = ', acc_ij + endif + acc_tot += acc_ij + normalz += dabs(i_exc) + + ! --- + + enddo + + if( (acc_tot/normalz) .gt. 1d-3 ) then + print*, ' acc_tot = ', acc_tot + print*, ' normalz = ', normalz + endif + enddo + + return +end subroutine test_fit_ugradu + +! --- + +subroutine test_fit_u() + + implicit none + + integer :: jpoint, ipoint, i + double precision :: i_exc, i_fit, i_num, x2 + double precision :: r1(3), r2(3), dx, dy, dz + double precision :: eps_ij, acc_tot, acc_ij, normalz, coef, expo + + double precision, external :: j12_mu + + print*, ' test_fit_u ...' + + eps_ij = 1d-3 + + do jpoint = 1, n_points_final_grid + r2(1) = final_grid_points(1,jpoint) + r2(2) = final_grid_points(2,jpoint) + r2(3) = final_grid_points(3,jpoint) + + acc_tot = 0.d0 + normalz = 0.d0 + do ipoint = 1, n_points_final_grid + + r1(1) = final_grid_points(1,ipoint) + r1(2) = final_grid_points(2,ipoint) + r1(3) = final_grid_points(3,ipoint) + + dx = r1(1) - r2(1) + dy = r1(2) - r2(2) + dz = r1(3) - r2(3) + x2 = dx * dx + dy * dy + dz * dz + if(x2 .lt. 1d-10) cycle + + i_fit = 0.d0 + do i = 1, n_max_fit_slat + expo = expo_gauss_j_mu_x(i) + coef = coef_gauss_j_mu_x(i) + i_fit += coef * dexp(-expo*x2) + enddo + + i_exc = j12_mu(r1, r2) + i_num = i_fit + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in test_fit_u on', ipoint + print *, ' analyt = ', i_exc + print *, ' numeri = ', i_num + print *, ' diff = ', acc_ij + endif + + acc_tot += acc_ij + normalz += dabs(i_exc) + enddo + + if( (acc_tot/normalz) .gt. 1d-3 ) then + print*, ' acc_tot = ', acc_tot + print*, ' normalz = ', normalz + endif + enddo + + return +end subroutine test_fit_u + +! --- + +subroutine test_fit_u2() + + implicit none + + integer :: jpoint, ipoint, i + double precision :: i_exc, i_fit, i_num, x2 + double precision :: r1(3), r2(3), dx, dy, dz, tmp + double precision :: eps_ij, acc_tot, acc_ij, normalz, coef, expo + + double precision, external :: j12_mu + + print*, ' test_fit_u2 ...' + + eps_ij = 1d-3 + + do jpoint = 1, n_points_final_grid + r2(1) = final_grid_points(1,jpoint) + r2(2) = final_grid_points(2,jpoint) + r2(3) = final_grid_points(3,jpoint) + + acc_tot = 0.d0 + normalz = 0.d0 + do ipoint = 1, n_points_final_grid + + r1(1) = final_grid_points(1,ipoint) + r1(2) = final_grid_points(2,ipoint) + r1(3) = final_grid_points(3,ipoint) + + dx = r1(1) - r2(1) + dy = r1(2) - r2(2) + dz = r1(3) - r2(3) + x2 = dx * dx + dy * dy + dz * dz + if(x2 .lt. 1d-10) cycle + + i_fit = 0.d0 + do i = 1, n_max_fit_slat + expo = expo_gauss_j_mu_x_2(i) + coef = coef_gauss_j_mu_x_2(i) + i_fit += coef * dexp(-expo*x2) + enddo + + tmp = j12_mu(r1, r2) + i_exc = tmp * tmp + i_num = i_fit + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in test_fit_u2 on', ipoint + print *, ' analyt = ', i_exc + print *, ' numeri = ', i_num + print *, ' diff = ', acc_ij + endif + + acc_tot += acc_ij + normalz += dabs(i_exc) + enddo + + if( (acc_tot/normalz) .gt. 1d-3 ) then + print*, ' acc_tot = ', acc_tot + print*, ' normalz = ', normalz + endif + enddo + + return +end subroutine test_fit_u2 + +! --- + + diff --git a/src/non_h_ints_mu/debug_integ_jmu_modif.irp.f b/src/non_h_ints_mu/debug_integ_jmu_modif.irp.f index 7b99cc91..bb585f63 100644 --- a/src/non_h_ints_mu/debug_integ_jmu_modif.irp.f +++ b/src/non_h_ints_mu/debug_integ_jmu_modif.irp.f @@ -17,25 +17,19 @@ program debug_integ_jmu_modif PROVIDE mu_erf j1b_pen - !call test_j1b_nucl() - !call test_grad_j1b_nucl() - !call test_lapl_j1b_nucl() - - !call test_list_b2() - !call test_list_b3() - - !call test_fit_u() - call test_fit_u2() - !call test_fit_ugradu() - - !call test_v_ij_u_cst_mu_j1b() - !call test_v_ij_erf_rk_cst_mu_j1b() - !call test_x_v_ij_erf_rk_cst_mu_j1b() - !call test_int2_u2_j1b2() - !call test_int2_grad1u2_grad2u2_j1b2() - - !call test_int2_grad1_u12_ao() - !call test_gradu_squared_u_ij_mu() + call test_v_ij_u_cst_mu_j1b() +! call test_v_ij_erf_rk_cst_mu_j1b() +! call test_x_v_ij_erf_rk_cst_mu_j1b() +! call test_int2_u2_j1b2() +! call test_int2_grad1u2_grad2u2_j1b2() +! call test_int2_u_grad1u_total_j1b2() +! +! call test_int2_grad1_u12_ao() +! +! call test_grad12_j12() +! call test_u12sq_j1bsq() +! call test_u12_grad1_u12_j1b_grad1_j1b() +! !call test_gradu_squared_u_ij_mu() end @@ -52,8 +46,9 @@ subroutine test_v_ij_u_cst_mu_j1b() PROVIDE v_ij_u_cst_mu_j1b - eps_ij = 1d-8 + eps_ij = 1d-3 acc_tot = 0.d0 + normalz = 0.d0 !do ipoint = 1, 10 do ipoint = 1, n_points_final_grid @@ -76,9 +71,8 @@ subroutine test_v_ij_u_cst_mu_j1b() enddo enddo - acc_tot = acc_tot / normalz - print*, ' normalized acc = ', acc_tot - print*, ' normalz = ', normalz + print*, ' acc_tot = ', acc_tot + print*, ' normalz = ', normalz return end subroutine test_v_ij_u_cst_mu_j1b @@ -96,8 +90,9 @@ subroutine test_v_ij_erf_rk_cst_mu_j1b() PROVIDE v_ij_erf_rk_cst_mu_j1b - eps_ij = 1d-8 + eps_ij = 1d-3 acc_tot = 0.d0 + normalz = 0.d0 !do ipoint = 1, 10 do ipoint = 1, n_points_final_grid @@ -120,9 +115,8 @@ subroutine test_v_ij_erf_rk_cst_mu_j1b() enddo enddo - acc_tot = acc_tot / normalz - print*, ' normalized acc = ', acc_tot - print*, ' normalz = ', normalz + print*, ' acc_tot = ', acc_tot + print*, ' normalz = ', normalz return end subroutine test_v_ij_erf_rk_cst_mu_j1b @@ -140,8 +134,9 @@ subroutine test_x_v_ij_erf_rk_cst_mu_j1b() PROVIDE x_v_ij_erf_rk_cst_mu_j1b - eps_ij = 1d-8 + eps_ij = 1d-3 acc_tot = 0.d0 + normalz = 0.d0 !do ipoint = 1, 10 do ipoint = 1, n_points_final_grid @@ -190,9 +185,8 @@ subroutine test_x_v_ij_erf_rk_cst_mu_j1b() enddo enddo - acc_tot = acc_tot / normalz - print*, ' normalized acc = ', acc_tot - print*, ' normalz = ', normalz + print*, ' acc_tot = ', acc_tot + print*, ' normalz = ', normalz return end subroutine test_x_v_ij_erf_rk_cst_mu_j1b @@ -210,8 +204,9 @@ subroutine test_int2_u2_j1b2() PROVIDE int2_u2_j1b2 - eps_ij = 1d-8 + eps_ij = 1d-3 acc_tot = 0.d0 + normalz = 0.d0 !do ipoint = 1, 10 do ipoint = 1, n_points_final_grid @@ -235,8 +230,8 @@ subroutine test_int2_u2_j1b2() enddo acc_tot = acc_tot / normalz - print*, ' normalized acc = ', acc_tot - print*, ' normalz = ', normalz + print*, ' acc_tot = ', acc_tot + print*, ' normalz = ', normalz return end subroutine test_int2_u2_j1b2 @@ -254,8 +249,9 @@ subroutine test_int2_grad1u2_grad2u2_j1b2() PROVIDE int2_grad1u2_grad2u2_j1b2 - eps_ij = 1d-8 + eps_ij = 1d-3 acc_tot = 0.d0 + normalz = 0.d0 !do ipoint = 1, 10 do ipoint = 1, n_points_final_grid @@ -278,9 +274,8 @@ subroutine test_int2_grad1u2_grad2u2_j1b2() enddo enddo - acc_tot = acc_tot / normalz - print*, ' normalized acc = ', acc_tot - print*, ' normalz = ', normalz + print*, ' acc_tot = ', acc_tot + print*, ' normalz = ', normalz return end subroutine test_int2_grad1u2_grad2u2_j1b2 @@ -298,8 +293,9 @@ subroutine test_int2_grad1_u12_ao() PROVIDE int2_grad1_u12_ao - eps_ij = 1d-6 + eps_ij = 1d-3 acc_tot = 0.d0 + normalz = 0.d0 do ipoint = 1, n_points_final_grid do j = 1, ao_num @@ -347,15 +343,90 @@ subroutine test_int2_grad1_u12_ao() enddo enddo - acc_tot = acc_tot / normalz - print*, ' normalized acc = ', acc_tot - print*, ' normalz = ', normalz + print*, ' acc_tot = ', acc_tot + print*, ' normalz = ', normalz return end subroutine test_int2_grad1_u12_ao ! --- +subroutine test_int2_u_grad1u_total_j1b2() + + implicit none + integer :: i, j, ipoint + double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz + double precision :: x, y, z + double precision :: integ(3) + + print*, ' test_int2_u_grad1u_total_j1b2 ...' + + PROVIDE int2_u_grad1u_j1b2 + PROVIDE int2_u_grad1u_x_j1b2 + + eps_ij = 1d-3 + acc_tot = 0.d0 + normalz = 0.d0 + + !do ipoint = 1, 10 + do ipoint = 1, n_points_final_grid + x = final_grid_points(1,ipoint) + y = final_grid_points(2,ipoint) + z = final_grid_points(3,ipoint) + + do j = 1, ao_num + do i = 1, ao_num + + call num_int2_u_grad1u_total_j1b2(i, j, ipoint, integ) + + i_exc = x * int2_u_grad1u_j1b2(i,j,ipoint) - int2_u_grad1u_x_j1b2(1,i,j,ipoint) + i_num = integ(1) + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in x part of int2_u_grad1u_total_j1b2 on', i, j, ipoint + print *, ' analyt integ = ', i_exc + print *, ' numeri integ = ', i_num + print *, ' diff = ', acc_ij + endif + acc_tot += acc_ij + normalz += dabs(i_num) + + i_exc = y * int2_u_grad1u_j1b2(i,j,ipoint) - int2_u_grad1u_x_j1b2(2,i,j,ipoint) + i_num = integ(2) + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in y part of int2_u_grad1u_total_j1b2 on', i, j, ipoint + print *, ' analyt integ = ', i_exc + print *, ' numeri integ = ', i_num + print *, ' diff = ', acc_ij + endif + acc_tot += acc_ij + normalz += dabs(i_num) + + i_exc = z * int2_u_grad1u_j1b2(i,j,ipoint) - int2_u_grad1u_x_j1b2(3,i,j,ipoint) + i_num = integ(3) + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in z part of int2_u_grad1u_total_j1b2 on', i, j, ipoint + print *, ' analyt integ = ', i_exc + print *, ' numeri integ = ', i_num + print *, ' diff = ', acc_ij + endif + acc_tot += acc_ij + normalz += dabs(i_num) + + enddo + enddo + enddo + + print*, ' acc_tot = ', acc_tot + print*, ' normalz = ', normalz + + return +end subroutine test_int2_u_grad1u_total_j1b2 + +! --- + subroutine test_gradu_squared_u_ij_mu() implicit none @@ -367,8 +438,9 @@ subroutine test_gradu_squared_u_ij_mu() PROVIDE gradu_squared_u_ij_mu - eps_ij = 1d-6 + eps_ij = 1d-3 acc_tot = 0.d0 + normalz = 0.d0 do ipoint = 1, n_points_final_grid do j = 1, ao_num @@ -390,458 +462,140 @@ subroutine test_gradu_squared_u_ij_mu() enddo enddo - acc_tot = acc_tot / normalz - print*, ' normalized acc = ', acc_tot - print*, ' normalz = ', normalz + print*, ' acc_tot = ', acc_tot + print*, ' normalz = ', normalz return end subroutine test_gradu_squared_u_ij_mu ! --- -subroutine test_j1b_nucl() +subroutine test_grad12_j12() implicit none - integer :: ipoint + integer :: i, j, ipoint double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz - double precision :: r(3) - double precision, external :: j1b_nucl + double precision, external :: num_grad12_j12 - print*, ' test_j1b_nucl ...' + print*, ' test_grad12_j12 ...' - PROVIDE v_1b + PROVIDE grad12_j12 - eps_ij = 1d-7 + eps_ij = 1d-3 acc_tot = 0.d0 + normalz = 0.d0 do ipoint = 1, n_points_final_grid - - r(1) = final_grid_points(1,ipoint) - r(2) = final_grid_points(2,ipoint) - r(3) = final_grid_points(3,ipoint) - - i_exc = v_1b(ipoint) - i_num = j1b_nucl(r) - acc_ij = dabs(i_exc - i_num) - if(acc_ij .gt. eps_ij) then - print *, ' problem in v_1b on', ipoint - print *, ' analyt = ', i_exc - print *, ' numeri = ', i_num - print *, ' diff = ', acc_ij - endif - - acc_tot += acc_ij - normalz += dabs(i_num) - enddo - - acc_tot = acc_tot / normalz - print*, ' normalized acc = ', acc_tot - print*, ' normalz = ', normalz - - return -end subroutine test_j1b_nucl - -! --- - -subroutine test_grad_j1b_nucl() - - implicit none - integer :: ipoint - double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz - double precision :: r(3) - double precision, external :: grad_x_j1b_nucl - double precision, external :: grad_y_j1b_nucl - double precision, external :: grad_z_j1b_nucl - - print*, ' test_grad_j1b_nucl ...' - - PROVIDE v_1b_grad - - eps_ij = 1d-6 - acc_tot = 0.d0 - - do ipoint = 1, n_points_final_grid - - r(1) = final_grid_points(1,ipoint) - r(2) = final_grid_points(2,ipoint) - r(3) = final_grid_points(3,ipoint) - - i_exc = v_1b_grad(1,ipoint) - i_num = grad_x_j1b_nucl(r) - acc_ij = dabs(i_exc - i_num) - if(acc_ij .gt. eps_ij) then - print *, ' problem in x of v_1b_grad on', ipoint - print *, ' analyt = ', i_exc - print *, ' numeri = ', i_num - print *, ' diff = ', acc_ij - endif - - i_exc = v_1b_grad(2,ipoint) - i_num = grad_y_j1b_nucl(r) - acc_ij = dabs(i_exc - i_num) - if(acc_ij .gt. eps_ij) then - print *, ' problem in y of v_1b_grad on', ipoint - print *, ' analyt = ', i_exc - print *, ' numeri = ', i_num - print *, ' diff = ', acc_ij - endif - - i_exc = v_1b_grad(3,ipoint) - i_num = grad_z_j1b_nucl(r) - acc_ij = dabs(i_exc - i_num) - if(acc_ij .gt. eps_ij) then - print *, ' problem in z of v_1b_grad on', ipoint - print *, ' analyt = ', i_exc - print *, ' numeri = ', i_num - print *, ' diff = ', acc_ij - endif - - acc_tot += acc_ij - normalz += dabs(i_num) - enddo - - acc_tot = acc_tot / normalz - print*, ' normalized acc = ', acc_tot - print*, ' normalz = ', normalz - - return -end subroutine test_grad_j1b_nucl - -! --- - -subroutine test_lapl_j1b_nucl() - - implicit none - integer :: ipoint - double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz - double precision :: r(3) - double precision, external :: lapl_j1b_nucl - - print*, ' test_lapl_j1b_nucl ...' - - PROVIDE v_1b_lapl - - eps_ij = 1d-5 - acc_tot = 0.d0 - - do ipoint = 1, n_points_final_grid - - r(1) = final_grid_points(1,ipoint) - r(2) = final_grid_points(2,ipoint) - r(3) = final_grid_points(3,ipoint) - - i_exc = v_1b_lapl(ipoint) - i_num = lapl_j1b_nucl(r) - acc_ij = dabs(i_exc - i_num) - if(acc_ij .gt. eps_ij) then - print *, ' problem in v_1b_lapl on', ipoint - print *, ' analyt = ', i_exc - print *, ' numeri = ', i_num - print *, ' diff = ', acc_ij - endif - - acc_tot += acc_ij - normalz += dabs(i_num) - enddo - - acc_tot = acc_tot / normalz - print*, ' normalized acc = ', acc_tot - print*, ' normalz = ', normalz - - return -end subroutine test_lapl_j1b_nucl - -! --- - -subroutine test_list_b2() - - implicit none - integer :: ipoint - double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz - double precision :: r(3) - double precision, external :: j1b_nucl - - print*, ' test_list_b2 ...' - - PROVIDE v_1b_list_b2 - - eps_ij = 1d-7 - acc_tot = 0.d0 - - do ipoint = 1, n_points_final_grid - - r(1) = final_grid_points(1,ipoint) - r(2) = final_grid_points(2,ipoint) - r(3) = final_grid_points(3,ipoint) - - i_exc = v_1b_list_b2(ipoint) - i_num = j1b_nucl(r) - acc_ij = dabs(i_exc - i_num) - if(acc_ij .gt. eps_ij) then - print *, ' problem in list_b2 on', ipoint - print *, ' analyt = ', i_exc - print *, ' numeri = ', i_num - print *, ' diff = ', acc_ij - endif - - acc_tot += acc_ij - normalz += dabs(i_num) - enddo - - acc_tot = acc_tot / normalz - print*, ' normalized acc = ', acc_tot - print*, ' normalz = ', normalz - - return -end subroutine test_list_b2 - -! --- - -subroutine test_list_b3() - - implicit none - integer :: ipoint - double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_tmp, i_num, normalz - double precision :: r(3) - double precision, external :: j1b_nucl - - print*, ' test_list_b3 ...' - - PROVIDE v_1b_list_b3 - - eps_ij = 1d-7 - acc_tot = 0.d0 - - do ipoint = 1, n_points_final_grid - - r(1) = final_grid_points(1,ipoint) - r(2) = final_grid_points(2,ipoint) - r(3) = final_grid_points(3,ipoint) - - i_exc = v_1b_list_b3(ipoint) - i_tmp = j1b_nucl(r) - i_num = i_tmp * i_tmp - acc_ij = dabs(i_exc - i_num) - if(acc_ij .gt. eps_ij) then - print *, ' problem in list_b3 on', ipoint - print *, ' analyt = ', i_exc - print *, ' numeri = ', i_num - print *, ' diff = ', acc_ij - endif - - acc_tot += acc_ij - normalz += dabs(i_num) - enddo - - acc_tot = acc_tot / normalz - print*, ' normalized acc = ', acc_tot - print*, ' normalz = ', normalz - - return -end subroutine test_list_b3 - -! --- - -subroutine test_fit_ugradu() - - implicit none - - integer :: ipoint, i - double precision :: i_exc, i_fit, i_num, x2 - double precision :: r1(3), r2(3), grad(3) - double precision :: eps_ij, acc_tot, acc_ij, normalz, coef, expo - - double precision, external :: j12_mu - - print*, ' test_fit_ugradu ...' - - eps_ij = 1d-7 - acc_tot = 0.d0 - - r2 = 0.d0 - - do ipoint = 1, n_points_final_grid - - r1(1) = final_grid_points(1,ipoint) - r1(2) = final_grid_points(2,ipoint) - r1(3) = final_grid_points(3,ipoint) - x2 = r1(1) * r1(1) + r1(2) * r1(2) + r1(3) * r1(3) - if(x2 .lt. 1d-10) cycle - - i_fit = 0.d0 - do i = 1, n_max_fit_slat - expo = expo_gauss_j_mu_1_erf(i) - coef = coef_gauss_j_mu_1_erf(i) - i_fit += coef * dexp(-expo*x2) + do j = 1, ao_num + do i = 1, ao_num + + i_exc = grad12_j12(i,j,ipoint) + i_num = num_grad12_j12(i, j, ipoint) + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in grad12_j12 on', i, j, ipoint + print *, ' analyt integ = ', i_exc + print *, ' numeri integ = ', i_num + print *, ' diff = ', acc_ij + endif + + acc_tot += acc_ij + normalz += dabs(i_num) + enddo enddo - i_fit = i_fit / dsqrt(x2) - - call grad1_j12_mu_exc(r1, r2, grad) - - ! --- - - i_exc = j12_mu(r1, r2) * grad(1) - i_num = i_fit * r1(1) - acc_ij = dabs(i_exc - i_num) - if(acc_ij .gt. eps_ij) then - print *, ' problem on x in test_fit_ugradu on', ipoint - print *, ' analyt = ', i_exc - print *, ' numeri = ', i_num - print *, ' diff = ', acc_ij - endif - acc_tot += acc_ij - normalz += dabs(i_exc) - - ! --- - - i_exc = j12_mu(r1, r2) * grad(2) - i_num = i_fit * r1(2) - acc_ij = dabs(i_exc - i_num) - if(acc_ij .gt. eps_ij) then - print *, ' problem on y in test_fit_ugradu on', ipoint - print *, ' analyt = ', i_exc - print *, ' numeri = ', i_num - print *, ' diff = ', acc_ij - endif - acc_tot += acc_ij - normalz += dabs(i_exc) - - ! --- - - i_exc = j12_mu(r1, r2) * grad(3) - i_num = i_fit * r1(3) - acc_ij = dabs(i_exc - i_num) - if(acc_ij .gt. eps_ij) then - print *, ' problem on z in test_fit_ugradu on', ipoint - print *, ' analyt = ', i_exc - print *, ' numeri = ', i_num - print *, ' diff = ', acc_ij - endif - acc_tot += acc_ij - normalz += dabs(i_exc) - - ! --- - enddo - acc_tot = acc_tot / normalz - print*, ' normalized acc = ', acc_tot - print*, ' normalz = ', normalz + print*, ' acc_tot = ', acc_tot + print*, ' normalz = ', normalz return -end subroutine test_fit_ugradu +end subroutine test_grad12_j12 ! --- -subroutine test_fit_u() +subroutine test_u12sq_j1bsq() implicit none + integer :: i, j, ipoint + double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz + double precision, external :: num_u12sq_j1bsq - integer :: ipoint, i - double precision :: i_exc, i_fit, i_num, x2 - double precision :: r1(3), r2(3) - double precision :: eps_ij, acc_tot, acc_ij, normalz, coef, expo + print*, ' test_u12sq_j1bsq ...' - double precision, external :: j12_mu + PROVIDE u12sq_j1bsq - print*, ' test_fit_u ...' - - eps_ij = 1d-7 + eps_ij = 1d-3 acc_tot = 0.d0 - - r2 = 0.d0 + normalz = 0.d0 do ipoint = 1, n_points_final_grid + do j = 1, ao_num + do i = 1, ao_num - r1(1) = final_grid_points(1,ipoint) - r1(2) = final_grid_points(2,ipoint) - r1(3) = final_grid_points(3,ipoint) - x2 = r1(1) * r1(1) + r1(2) * r1(2) + r1(3) * r1(3) - if(x2 .lt. 1d-10) cycle + i_exc = u12sq_j1bsq(i,j,ipoint) + i_num = num_u12sq_j1bsq(i, j, ipoint) + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in u12sq_j1bsq on', i, j, ipoint + print *, ' analyt integ = ', i_exc + print *, ' numeri integ = ', i_num + print *, ' diff = ', acc_ij + endif - i_fit = 0.d0 - do i = 1, n_max_fit_slat - expo = expo_gauss_j_mu_x(i) - coef = coef_gauss_j_mu_x(i) - i_fit += coef * dexp(-expo*x2) + acc_tot += acc_ij + normalz += dabs(i_num) + enddo enddo - - i_exc = j12_mu(r1, r2) - i_num = i_fit - acc_ij = dabs(i_exc - i_num) - if(acc_ij .gt. eps_ij) then - print *, ' problem in test_fit_u on', ipoint - print *, ' analyt = ', i_exc - print *, ' numeri = ', i_num - print *, ' diff = ', acc_ij - endif - - acc_tot += acc_ij - normalz += dabs(i_exc) enddo - acc_tot = acc_tot / normalz - print*, ' normalized acc = ', acc_tot - print*, ' normalz = ', normalz + print*, ' acc_tot = ', acc_tot + print*, ' normalz = ', normalz return -end subroutine test_fit_u +end subroutine test_u12sq_j1bsq ! --- -subroutine test_fit_u2() +subroutine test_u12_grad1_u12_j1b_grad1_j1b() implicit none + integer :: i, j, ipoint + double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz + double precision, external :: num_u12_grad1_u12_j1b_grad1_j1b - integer :: ipoint, i - double precision :: i_exc, i_fit, i_num, x2 - double precision :: r1(3), r2(3) - double precision :: eps_ij, acc_tot, acc_ij, normalz, coef, expo + print*, ' test_u12_grad1_u12_j1b_grad1_j1b ...' - double precision, external :: j12_mu + PROVIDE u12_grad1_u12_j1b_grad1_j1b - print*, ' test_fit_u2 ...' - - eps_ij = 1d-7 + eps_ij = 1d-3 acc_tot = 0.d0 - - r2 = 0.d0 + normalz = 0.d0 do ipoint = 1, n_points_final_grid + do j = 1, ao_num + do i = 1, ao_num - r1(1) = final_grid_points(1,ipoint) - r1(2) = final_grid_points(2,ipoint) - r1(3) = final_grid_points(3,ipoint) - x2 = r1(1) * r1(1) + r1(2) * r1(2) + r1(3) * r1(3) - if(x2 .lt. 1d-10) cycle + i_exc = u12_grad1_u12_j1b_grad1_j1b(i,j,ipoint) + i_num = num_u12_grad1_u12_j1b_grad1_j1b(i, j, ipoint) + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in u12_grad1_u12_j1b_grad1_j1b on', i, j, ipoint + print *, ' analyt integ = ', i_exc + print *, ' numeri integ = ', i_num + print *, ' diff = ', acc_ij + endif - i_fit = 0.d0 - do i = 1, n_max_fit_slat - expo = expo_gauss_j_mu_x_2(i) - coef = coef_gauss_j_mu_x_2(i) - i_fit += coef * dexp(-expo*x2) + acc_tot += acc_ij + normalz += dabs(i_num) + enddo enddo - - i_exc = j12_mu(r1, r2) * j12_mu(r1, r2) - i_num = i_fit - acc_ij = dabs(i_exc - i_num) - if(acc_ij .gt. eps_ij) then - print *, ' problem in test_fit_u2 on', ipoint - print *, ' analyt = ', i_exc - print *, ' numeri = ', i_num - print *, ' diff = ', acc_ij - endif - - acc_tot += acc_ij - normalz += dabs(i_exc) enddo - acc_tot = acc_tot / normalz - print*, ' normalized acc = ', acc_tot - print*, ' normalz = ', normalz + print*, ' acc_tot = ', acc_tot + print*, ' normalz = ', normalz return -end subroutine test_fit_u2 +end subroutine test_u12_grad1_u12_j1b_grad1_j1b, ! --- diff --git a/src/non_h_ints_mu/grad_squared.irp.f b/src/non_h_ints_mu/grad_squared.irp.f index bf37c551..4e70bc5c 100644 --- a/src/non_h_ints_mu/grad_squared.irp.f +++ b/src/non_h_ints_mu/grad_squared.irp.f @@ -23,52 +23,63 @@ BEGIN_PROVIDER [ double precision, gradu_squared_u_ij_mu, (ao_num, ao_num, n_poi ! + -1.00 x v1 (grad_1 v1) \int r2 \phi_i(2) \phi_j(2) (grad_1 u12) v2^2 ! = v1^2 x int2_grad1u2_grad2u2_j1b2 ! + -0.5 x (grad_1 v1)^2 x int2_u2_j1b2 - ! + -1.0 X V1 x (grad_1 v1) \cdot int2_u_grad1u_x_j1b + ! + -1.0 X V1 x (grad_1 v1) \cdot [ int2_u_grad1u_j1b2 x r - int2_u_grad1u_x_j1b ] ! ! END_DOC implicit none integer :: ipoint, i, j, m, igauss - double precision :: r(3), delta, coef - double precision :: tmp_v, tmp_x, tmp_y, tmp_z, tmp1, tmp2, tmp3, tmp4, tmp5 + double precision :: x, y, z, r(3), delta, coef + double precision :: tmp_v, tmp_x, tmp_y, tmp_z + double precision :: tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9 double precision :: time0, time1 double precision, external :: overlap_gauss_r12_ao print*, ' providing gradu_squared_u_ij_mu ...' call wall_time(time0) - PROVIDE j1b_type j1b_pen + PROVIDE j1b_type if(j1b_type .eq. 3) then do ipoint = 1, n_points_final_grid + x = final_grid_points(1,ipoint) + y = final_grid_points(2,ipoint) + z = final_grid_points(3,ipoint) tmp_v = v_1b (ipoint) tmp_x = v_1b_grad(1,ipoint) tmp_y = v_1b_grad(2,ipoint) tmp_z = v_1b_grad(3,ipoint) tmp1 = tmp_v * tmp_v - tmp2 = 0.5d0 * (tmp_x * tmp_x + tmp_y * tmp_y + tmp_z * tmp_z) + tmp2 = -0.5d0 * (tmp_x * tmp_x + tmp_y * tmp_y + tmp_z * tmp_z) tmp3 = tmp_v * tmp_x tmp4 = tmp_v * tmp_y tmp5 = tmp_v * tmp_z + tmp6 = -x * tmp3 + tmp7 = -y * tmp4 + tmp8 = -z * tmp5 + do j = 1, ao_num do i = 1, ao_num - gradu_squared_u_ij_mu(j,i,ipoint) += tmp1 * int2_grad1u2_grad2u2_j1b2(i,j,ipoint) & - - tmp2 * int2_u2_j1b2 (i,j,ipoint) & - - tmp3 * int2_u_grad1u_x_j1b (1,i,j,ipoint) & - - tmp4 * int2_u_grad1u_x_j1b (2,i,j,ipoint) & - - tmp5 * int2_u_grad1u_x_j1b (3,i,j,ipoint) + tmp9 = int2_u_grad1u_j1b2(i,j,ipoint) + + gradu_squared_u_ij_mu(i,j,ipoint) = tmp1 * int2_grad1u2_grad2u2_j1b2(i,j,ipoint) & + + tmp2 * int2_u2_j1b2 (i,j,ipoint) & + + tmp6 * tmp9 + tmp3 * int2_u_grad1u_x_j1b2(1,i,j,ipoint) & + + tmp7 * tmp9 + tmp4 * int2_u_grad1u_x_j1b2(2,i,j,ipoint) & + + tmp8 * tmp9 + tmp5 * int2_u_grad1u_x_j1b2(3,i,j,ipoint) enddo enddo enddo else + gradu_squared_u_ij_mu = 0.d0 do ipoint = 1, n_points_final_grid r(1) = final_grid_points(1,ipoint) r(2) = final_grid_points(2,ipoint) @@ -78,7 +89,7 @@ BEGIN_PROVIDER [ double precision, gradu_squared_u_ij_mu, (ao_num, ao_num, n_poi do igauss = 1, n_max_fit_slat delta = expo_gauss_1_erf_x_2(igauss) coef = coef_gauss_1_erf_x_2(igauss) - gradu_squared_u_ij_mu(j,i,ipoint) += -0.25d0 * coef * overlap_gauss_r12_ao(r, delta, i, j) + gradu_squared_u_ij_mu(i,j,ipoint) += -0.25d0 * coef * overlap_gauss_r12_ao(r, delta, i, j) enddo enddo enddo @@ -93,6 +104,57 @@ END_PROVIDER ! --- +!BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao_num)] +! +! BEGIN_DOC +! ! +! ! tc_grad_square_ao(k,i,l,j) = -1/2 +! ! +! END_DOC +! +! implicit none +! integer :: ipoint, i, j, k, l +! double precision :: weight1, ao_ik_r, ao_i_r +! double precision, allocatable :: ac_mat(:,:,:,:) +! +! allocate(ac_mat(ao_num,ao_num,ao_num,ao_num)) +! ac_mat = 0.d0 +! +! do ipoint = 1, n_points_final_grid +! weight1 = final_weight_at_r_vector(ipoint) +! +! do i = 1, ao_num +! ao_i_r = weight1 * aos_in_r_array_transp(ipoint,i) +! +! do k = 1, ao_num +! ao_ik_r = ao_i_r * aos_in_r_array_transp(ipoint,k) +! +! do j = 1, ao_num +! do l = 1, ao_num +! ac_mat(k,i,l,j) += ao_ik_r * gradu_squared_u_ij_mu(l,j,ipoint) +! enddo +! enddo +! enddo +! enddo +! enddo +! +! do j = 1, ao_num +! do l = 1, ao_num +! do i = 1, ao_num +! do k = 1, ao_num +! tc_grad_square_ao(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i) +! !write(11,*) tc_grad_square_ao(k,i,l,j) +! enddo +! enddo +! enddo +! enddo +! +! deallocate(ac_mat) +! +!END_PROVIDER + +! --- + BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao_num)] BEGIN_DOC @@ -103,27 +165,27 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao implicit none integer :: ipoint, i, j, k, l - double precision :: contrib, weight1, ao_k_r, ao_i_r - double precision, allocatable :: ac_mat(:,:,:,:) + double precision :: weight1, ao_ik_r, ao_i_r + double precision, allocatable :: ac_mat(:,:,:,:), bc_mat(:,:,:,:) allocate(ac_mat(ao_num,ao_num,ao_num,ao_num)) ac_mat = 0.d0 + allocate(bc_mat(ao_num,ao_num,ao_num,ao_num)) + bc_mat = 0.d0 do ipoint = 1, n_points_final_grid - weight1 = 0.5d0 * final_weight_at_r_vector(ipoint) + weight1 = final_weight_at_r_vector(ipoint) do i = 1, ao_num - ao_i_r = aos_in_r_array_transp(ipoint,i) + ao_i_r = weight1 * aos_in_r_array_transp(ipoint,i) do k = 1, ao_num - ao_k_r = aos_in_r_array_transp(ipoint,k) + ao_ik_r = ao_i_r * aos_in_r_array_transp(ipoint,k) do j = 1, ao_num do l = 1, ao_num - - contrib = gradu_squared_u_ij_mu(l,j,ipoint) * ao_k_r * ao_i_r - - ac_mat(k,i,l,j) += weight1 * contrib + ac_mat(k,i,l,j) += ao_ik_r * ( u12sq_j1bsq(l,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b(l,j,ipoint) ) + bc_mat(k,i,l,j) += ao_ik_r * grad12_j12(l,j,ipoint) enddo enddo enddo @@ -134,13 +196,147 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao do l = 1, ao_num do i = 1, ao_num do k = 1, ao_num - tc_grad_square_ao(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i) + tc_grad_square_ao(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i) + bc_mat(k,i,l,j) enddo enddo enddo enddo deallocate(ac_mat) + deallocate(bc_mat) + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, grad12_j12, (ao_num, ao_num, n_points_final_grid) ] + + implicit none + integer :: ipoint, i, j, m, igauss + double precision :: r(3), delta, coef + double precision :: tmp1 + double precision :: time0, time1 + double precision, external :: overlap_gauss_r12_ao + + print*, ' providing grad12_j12 ...' + call wall_time(time0) + + PROVIDE j1b_type + + if(j1b_type .eq. 3) then + + do ipoint = 1, n_points_final_grid + tmp1 = v_1b(ipoint) + tmp1 = tmp1 * tmp1 + do j = 1, ao_num + do i = 1, ao_num + grad12_j12(i,j,ipoint) = tmp1 * int2_grad1u2_grad2u2_j1b2(i,j,ipoint) + enddo + enddo + enddo + + else + + grad12_j12 = 0.d0 + do ipoint = 1, n_points_final_grid + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + do j = 1, ao_num + do i = 1, ao_num + do igauss = 1, n_max_fit_slat + delta = expo_gauss_1_erf_x_2(igauss) + coef = coef_gauss_1_erf_x_2(igauss) + grad12_j12(i,j,ipoint) += -0.25d0 * coef * overlap_gauss_r12_ao(r, delta, i, j) + enddo + enddo + enddo + enddo + + endif + + call wall_time(time1) + print*, ' Wall time for grad12_j12 = ', time1 - time0 + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, u12sq_j1bsq, (ao_num, ao_num, n_points_final_grid) ] + + implicit none + integer :: ipoint, i, j + double precision :: tmp_x, tmp_y, tmp_z + double precision :: tmp1 + double precision :: time0, time1 + + print*, ' providing u12sq_j1bsq ...' + call wall_time(time0) + + do ipoint = 1, n_points_final_grid + tmp_x = v_1b_grad(1,ipoint) + tmp_y = v_1b_grad(2,ipoint) + tmp_z = v_1b_grad(3,ipoint) + tmp1 = -0.5d0 * (tmp_x * tmp_x + tmp_y * tmp_y + tmp_z * tmp_z) + do j = 1, ao_num + do i = 1, ao_num + u12sq_j1bsq(i,j,ipoint) = tmp1 * int2_u2_j1b2(i,j,ipoint) + enddo + enddo + enddo + + call wall_time(time1) + print*, ' Wall time for u12sq_j1bsq = ', time1 - time0 + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, u12_grad1_u12_j1b_grad1_j1b, (ao_num, ao_num, n_points_final_grid) ] + + implicit none + integer :: ipoint, i, j, m, igauss + double precision :: x, y, z + double precision :: tmp_v, tmp_x, tmp_y, tmp_z + double precision :: tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9 + double precision :: time0, time1 + double precision, external :: overlap_gauss_r12_ao + + print*, ' providing u12_grad1_u12_j1b_grad1_j1b ...' + call wall_time(time0) + + do ipoint = 1, n_points_final_grid + + x = final_grid_points(1,ipoint) + y = final_grid_points(2,ipoint) + z = final_grid_points(3,ipoint) + tmp_v = v_1b (ipoint) + tmp_x = v_1b_grad(1,ipoint) + tmp_y = v_1b_grad(2,ipoint) + tmp_z = v_1b_grad(3,ipoint) + + tmp3 = tmp_v * tmp_x + tmp4 = tmp_v * tmp_y + tmp5 = tmp_v * tmp_z + + tmp6 = -x * tmp3 + tmp7 = -y * tmp4 + tmp8 = -z * tmp5 + + do j = 1, ao_num + do i = 1, ao_num + + tmp9 = int2_u_grad1u_j1b2(i,j,ipoint) + + u12_grad1_u12_j1b_grad1_j1b(i,j,ipoint) = tmp6 * tmp9 + tmp3 * int2_u_grad1u_x_j1b2(1,i,j,ipoint) & + + tmp7 * tmp9 + tmp4 * int2_u_grad1u_x_j1b2(2,i,j,ipoint) & + + tmp8 * tmp9 + tmp5 * int2_u_grad1u_x_j1b2(3,i,j,ipoint) + enddo + enddo + enddo + + call wall_time(time1) + print*, ' Wall time for u12_grad1_u12_j1b_grad1_j1b = ', time1 - time0 END_PROVIDER diff --git a/src/non_h_ints_mu/grad_tc_int.irp.f b/src/non_h_ints_mu/grad_tc_int.irp.f index 40600335..cb3b71a3 100644 --- a/src/non_h_ints_mu/grad_tc_int.irp.f +++ b/src/non_h_ints_mu/grad_tc_int.irp.f @@ -1,67 +1,75 @@ + +! --- + BEGIN_PROVIDER [double precision, ao_non_hermit_term_chemist, (ao_num, ao_num, ao_num, ao_num)] - implicit none -BEGIN_DOC -! 1 1 2 2 1 2 1 2 -! -! ao_non_hermit_term_chemist(k,i,l,j) = < k l | [erf( mu r12) - 1] d/d_r12 | i j > on the AO basis -END_DOC - integer :: i,j,k,l,ipoint,m - double precision :: weight1,thr,r(3) - thr = 1.d-8 - double precision, allocatable :: b_mat(:,:,:,:),ac_mat(:,:,:,:) -! provide v_ij_erf_rk_cst_mu - provide v_ij_erf_rk_cst_mu x_v_ij_erf_rk_cst_mu -! ao_non_hermit_term_chemist = non_h_ints -! return + + BEGIN_DOC + ! 1 1 2 2 1 2 1 2 + ! + ! ao_non_hermit_term_chemist(k,i,l,j) = < k l | [erf( mu r12) - 1] d/d_r12 | i j > on the AO basis + ! + END_DOC + + implicit none + integer :: i, j, k, l, ipoint, m + double precision :: weight1, r(3) + double precision :: wall1, wall0 + double precision, allocatable :: b_mat(:,:,:,:), ac_mat(:,:,:,:) + + provide v_ij_erf_rk_cst_mu x_v_ij_erf_rk_cst_mu + call wall_time(wall0) - allocate(b_mat(n_points_final_grid,ao_num,ao_num,3),ac_mat(ao_num, ao_num, ao_num, ao_num)) - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & + allocate(b_mat(n_points_final_grid,ao_num,ao_num,3), ac_mat(ao_num,ao_num,ao_num,ao_num)) + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & !$OMP PRIVATE (i,k,m,ipoint,r,weight1) & !$OMP SHARED (aos_in_r_array_transp,aos_grad_in_r_array_transp_bis,b_mat)& !$OMP SHARED (ao_num,n_points_final_grid,final_grid_points,final_weight_at_r_vector) !$OMP DO SCHEDULE (static) - do m = 1, 3 - do i = 1, ao_num - do k = 1, ao_num - do ipoint = 1, n_points_final_grid - r(1) = final_grid_points(1,ipoint) - r(2) = final_grid_points(2,ipoint) - r(3) = final_grid_points(3,ipoint) - weight1 = final_weight_at_r_vector(ipoint) - b_mat(ipoint,k,i,m) = 0.5d0 * aos_in_r_array_transp(ipoint,k) * r(m) * weight1 * aos_grad_in_r_array_transp_bis(ipoint,i,m) + do m = 1, 3 + do i = 1, ao_num + do k = 1, ao_num + do ipoint = 1, n_points_final_grid + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + weight1 = final_weight_at_r_vector(ipoint) + b_mat(ipoint,k,i,m) = 0.5d0 * aos_in_r_array_transp(ipoint,k) * r(m) * weight1 * aos_grad_in_r_array_transp_bis(ipoint,i,m) + enddo + enddo enddo - enddo enddo - enddo !$OMP END DO !$OMP END PARALLEL - - ! (A) b_mat(ipoint,k,i,m) X v_ij_erf_rk_cst_mu(j,l,r1) - ! 1/2 \int dr1 x1 phi_k(1) d/dx1 phi_i(1) \int dr2 (1 - erf(mu_r12))/r12 phi_j(2) phi_l(2) + ! (A) b_mat(ipoint,k,i,m) X v_ij_erf_rk_cst_mu(j,l,r1) + ! 1/2 \int dr1 x1 phi_k(1) d/dx1 phi_i(1) \int dr2 (1 - erf(mu_r12))/r12 phi_j(2) phi_l(2) ac_mat = 0.d0 do m = 1, 3 - ! A B^T dim(A,1) dim(B,2) dim(A,2) alpha * A LDA - call dgemm("N","N",ao_num*ao_num,ao_num*ao_num,n_points_final_grid,1.d0,v_ij_erf_rk_cst_mu(1,1,1),ao_num*ao_num & - ,b_mat(1,1,1,m),n_points_final_grid,1.d0,ac_mat,ao_num*ao_num) + ! A B^T dim(A,1) dim(B,2) dim(A,2) alpha * A LDA + + call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 & + , v_ij_erf_rk_cst_mu(1,1,1), ao_num*ao_num, b_mat(1,1,1,m), n_points_final_grid & + , 1.d0, ac_mat, ao_num*ao_num) + enddo - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & !$OMP PRIVATE (i,k,m,ipoint,weight1) & !$OMP SHARED (aos_in_r_array_transp,aos_grad_in_r_array_transp_bis,b_mat,ao_num,n_points_final_grid,final_weight_at_r_vector) !$OMP DO SCHEDULE (static) - do m = 1, 3 - do i = 1, ao_num - do k = 1, ao_num - do ipoint = 1, n_points_final_grid - weight1 = final_weight_at_r_vector(ipoint) - b_mat(ipoint,k,i,m) = 0.5d0 * aos_in_r_array_transp(ipoint,k) * weight1 * aos_grad_in_r_array_transp_bis(ipoint,i,m) + do m = 1, 3 + do i = 1, ao_num + do k = 1, ao_num + do ipoint = 1, n_points_final_grid + weight1 = final_weight_at_r_vector(ipoint) + b_mat(ipoint,k,i,m) = 0.5d0 * aos_in_r_array_transp(ipoint,k) * weight1 * aos_grad_in_r_array_transp_bis(ipoint,i,m) + enddo + enddo enddo - enddo enddo - enddo !$OMP END DO !$OMP END PARALLEL @@ -69,117 +77,141 @@ END_DOC ! 1/2 \int dr1 phi_k(1) d/dx1 phi_i(1) \int dr2 x2(1 - erf(mu_r12))/r12 phi_j(2) phi_l(2) do m = 1, 3 ! A B^T dim(A,1) dim(B,2) dim(A,2) alpha * A LDA - call dgemm("N","N",ao_num*ao_num,ao_num*ao_num,n_points_final_grid,-1.d0,x_v_ij_erf_rk_cst_mu(1,1,1,m),ao_num*ao_num & - ,b_mat(1,1,1,m),n_points_final_grid,1.d0,ac_mat,ao_num*ao_num) + + call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, -1.d0 & + , x_v_ij_erf_rk_cst_mu(1,1,1,m), ao_num*ao_num, b_mat(1,1,1,m), n_points_final_grid & + , 1.d0, ac_mat, ao_num*ao_num) enddo - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & !$OMP PRIVATE (i,k,j,l) & !$OMP SHARED (ac_mat,ao_non_hermit_term_chemist,ao_num) !$OMP DO SCHEDULE (static) - do j = 1, ao_num - do l = 1, ao_num - do i = 1, ao_num - do k = 1, ao_num - ! (ki|lj) (ki|lj) (lj|ki) - ao_non_hermit_term_chemist(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i) + do j = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do k = 1, ao_num + ! (ki|lj) (ki|lj) (lj|ki) + ao_non_hermit_term_chemist(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i) + enddo + enddo enddo - enddo enddo - enddo !$OMP END DO !$OMP END PARALLEL - double precision :: wall1, wall0 + call wall_time(wall1) - print*,'wall time dgemm ',wall1 - wall0 + print *, ' wall time dgemm ', wall1 - wall0 + END_PROVIDER +! --- + +! TODO :: optimization :: transform into DGEM + BEGIN_PROVIDER [double precision, mo_non_hermit_term_chemist, (mo_num, mo_num, mo_num, mo_num)] - implicit none -BEGIN_DOC -! 1 1 2 2 1 2 1 2 -! -! mo_non_hermit_term_chemist(k,i,l,j) = < k l | [erf( mu r12) - 1] d/d_r12 | i j > on the MO basis -END_DOC - integer :: i,j,k,l,m,n,p,q - double precision, allocatable :: mo_tmp_1(:,:,:,:),mo_tmp_2(:,:,:,:),mo_tmp_3(:,:,:,:) + + BEGIN_DOC + ! 1 1 2 2 1 2 1 2 + ! + ! mo_non_hermit_term_chemist(k,i,l,j) = < k l | [erf( mu r12) - 1] d/d_r12 | i j > on the MO basis + END_DOC + + implicit none + integer :: i, j, k, l, m, n, p, q + double precision, allocatable :: mo_tmp_1(:,:,:,:), mo_tmp_2(:,:,:,:) - allocate(mo_tmp_1(mo_num,ao_num,ao_num,ao_num)) - ! TODO :: optimization :: transform into DGEM - mo_tmp_1 = 0.d0 - do m = 1, ao_num - do p = 1, ao_num - do n = 1, ao_num - do q = 1, ao_num - do k = 1, mo_num - ! (k n|p m) = sum_q c_qk * (q n|p m) - mo_tmp_1(k,n,p,m) += mo_coef_transp(k,q) * ao_non_hermit_term_chemist(q,n,p,m) - enddo + allocate(mo_tmp_1(mo_num,ao_num,ao_num,ao_num)) + mo_tmp_1 = 0.d0 + + do m = 1, ao_num + do p = 1, ao_num + do n = 1, ao_num + do q = 1, ao_num + do k = 1, mo_num + ! (k n|p m) = sum_q c_qk * (q n|p m) + mo_tmp_1(k,n,p,m) += mo_coef_transp(k,q) * ao_non_hermit_term_chemist(q,n,p,m) + enddo + enddo + enddo enddo - enddo enddo - enddo - free ao_non_hermit_term_chemist - allocate(mo_tmp_2(mo_num,mo_num,ao_num,ao_num)) - mo_tmp_2 = 0.d0 - do m = 1, ao_num - do p = 1, ao_num - do n = 1, ao_num - do i = 1, mo_num - do k = 1, mo_num - ! (k i|p m) = sum_n c_ni * (k n|p m) - mo_tmp_2(k,i,p,m) += mo_coef_transp(i,n) * mo_tmp_1(k,n,p,m) - enddo + free ao_non_hermit_term_chemist + + allocate(mo_tmp_2(mo_num,mo_num,ao_num,ao_num)) + mo_tmp_2 = 0.d0 + + do m = 1, ao_num + do p = 1, ao_num + do n = 1, ao_num + do i = 1, mo_num + do k = 1, mo_num + ! (k i|p m) = sum_n c_ni * (k n|p m) + mo_tmp_2(k,i,p,m) += mo_coef_transp(i,n) * mo_tmp_1(k,n,p,m) + enddo + enddo + enddo enddo - enddo enddo - enddo - deallocate(mo_tmp_1) - allocate(mo_tmp_1(mo_num,mo_num,mo_num,ao_num)) - mo_tmp_1 = 0.d0 - do m = 1, ao_num - do p = 1, ao_num - do l = 1, mo_num - do i = 1, mo_num - do k = 1, mo_num - mo_tmp_1(k,i,l,m) += mo_coef_transp(l,p) * mo_tmp_2(k,i,p,m) - enddo + deallocate(mo_tmp_1) + + allocate(mo_tmp_1(mo_num,mo_num,mo_num,ao_num)) + mo_tmp_1 = 0.d0 + + do m = 1, ao_num + do p = 1, ao_num + do l = 1, mo_num + do i = 1, mo_num + do k = 1, mo_num + mo_tmp_1(k,i,l,m) += mo_coef_transp(l,p) * mo_tmp_2(k,i,p,m) + enddo + enddo + enddo enddo - enddo enddo - enddo - deallocate(mo_tmp_2) - mo_non_hermit_term_chemist = 0.d0 - do m = 1, ao_num - do j = 1, mo_num - do l = 1, mo_num - do i = 1, mo_num - do k = 1, mo_num - mo_non_hermit_term_chemist(k,i,l,j) += mo_coef_transp(j,m) * mo_tmp_1(k,i,l,m) - enddo + deallocate(mo_tmp_2) + + mo_non_hermit_term_chemist = 0.d0 + do m = 1, ao_num + do j = 1, mo_num + do l = 1, mo_num + do i = 1, mo_num + do k = 1, mo_num + mo_non_hermit_term_chemist(k,i,l,j) += mo_coef_transp(j,m) * mo_tmp_1(k,i,l,m) + enddo + enddo + enddo enddo - enddo enddo - enddo + deallocate(mo_tmp_1) END_PROVIDER +! --- + BEGIN_PROVIDER [double precision, mo_non_hermit_term, (mo_num, mo_num, mo_num, mo_num)] - implicit none -BEGIN_DOC -! 1 2 1 2 1 2 1 2 -! -! mo_non_hermit_term(k,l,i,j) = < k l | [erf( mu r12) - 1] d/d_r12 | i j > on the MO basis -END_DOC - integer :: i,j,k,l - do j = 1, mo_num - do i = 1, mo_num - do l = 1, mo_num - do k = 1, mo_num - mo_non_hermit_term(k,l,i,j) = mo_non_hermit_term_chemist(k,i,l,j) + + BEGIN_DOC + ! 1 2 1 2 1 2 1 2 + ! + ! mo_non_hermit_term(k,l,i,j) = < k l | [erf( mu r12) - 1] d/d_r12 | i j > on the MO basis + END_DOC + + implicit none + integer :: i, j, k, l + + do j = 1, mo_num + do i = 1, mo_num + do l = 1, mo_num + do k = 1, mo_num + mo_non_hermit_term(k,l,i,j) = mo_non_hermit_term_chemist(k,i,l,j) + enddo + enddo enddo - enddo enddo - enddo + END_PROVIDER + +! --- + diff --git a/src/non_h_ints_mu/j12_nucl_utils.irp.f b/src/non_h_ints_mu/j12_nucl_utils.irp.f index a6dd0939..f3b68f43 100644 --- a/src/non_h_ints_mu/j12_nucl_utils.irp.f +++ b/src/non_h_ints_mu/j12_nucl_utils.irp.f @@ -586,4 +586,38 @@ end subroutine grad1_j12_mu_exc ! --- +subroutine grad1_jmu_modif_num(r1, r2, grad) + + implicit none + + double precision, intent(in) :: r1(3), r2(3) + double precision, intent(out) :: grad(3) + + double precision :: tmp0, tmp1, tmp2, tmp3, tmp4, grad_u12(3) + + double precision, external :: j12_mu + double precision, external :: j1b_nucl + double precision, external :: grad_x_j1b_nucl + double precision, external :: grad_y_j1b_nucl + double precision, external :: grad_z_j1b_nucl + + call grad1_j12_mu_exc(r1, r2, grad_u12) + + tmp0 = j1b_nucl(r1) + tmp1 = j1b_nucl(r2) + tmp2 = j12_mu(r1, r2) + tmp3 = tmp0 * tmp1 + tmp4 = tmp2 * tmp1 + + grad(1) = tmp3 * grad_u12(1) + tmp4 * grad_x_j1b_nucl(r1) + grad(2) = tmp3 * grad_u12(2) + tmp4 * grad_y_j1b_nucl(r1) + grad(3) = tmp3 * grad_u12(3) + tmp4 * grad_z_j1b_nucl(r1) + + return +end subroutine grad1_jmu_modif_num + +! --- + + + diff --git a/src/non_h_ints_mu/new_grad_tc.irp.f b/src/non_h_ints_mu/new_grad_tc.irp.f index db659520..d34e629c 100644 --- a/src/non_h_ints_mu/new_grad_tc.irp.f +++ b/src/non_h_ints_mu/new_grad_tc.irp.f @@ -17,10 +17,10 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao, (3, ao_num, ao_num, n_poin ! if J(r1,r2) = u12 x v1 x v2 ! ! int2_grad1_u12_ao(:,i,j,ipoint) = v1 x [ 0.5 x \int dr2 [(r1 - r2) (erf(mu * r12)-1)r_12] v2 \phi_i(r2) \phi_j(r2) ] - ! + \grad_1 v1 x [ \int dr2 u12 v2 \phi_i(r2) \phi_j(r2) ] + ! - \grad_1 v1 x [ \int dr2 u12 v2 \phi_i(r2) \phi_j(r2) ] ! = 0.5 v_1b(ipoint) * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) * r(:) ! - 0.5 v_1b(ipoint) * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,:) - ! + v_1b_grad[:,ipoint] * v_ij_u_cst_mu_j1b(i,j,ipoint) + ! - v_1b_grad[:,ipoint] * v_ij_u_cst_mu_j1b(i,j,ipoint) ! ! END_DOC @@ -29,7 +29,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao, (3, ao_num, ao_num, n_poin integer :: ipoint, i, j double precision :: x, y, z, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2 - PROVIDE j1b_type j1b_pen + PROVIDE j1b_type if(j1b_type .eq. 3) then @@ -46,12 +46,12 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao, (3, ao_num, ao_num, n_poin do j = 1, ao_num do i = 1, ao_num - tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) - tmp2 = v_ij_u_cst_mu_j1b(i,j,ipoint) + tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) + tmp2 = v_ij_u_cst_mu_j1b(i,j,ipoint) - int2_grad1_u12_ao(1,i,j,ipoint) = tmp1 * x - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1) + tmp_x * tmp2 - int2_grad1_u12_ao(2,i,j,ipoint) = tmp1 * y - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2) + tmp_y * tmp2 - int2_grad1_u12_ao(3,i,j,ipoint) = tmp1 * z - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3) + tmp_z * tmp2 + int2_grad1_u12_ao(1,i,j,ipoint) = tmp1 * x - tmp0 * x_v_ij_erf_rk_cst_mu_tmp_j1b(1,i,j,ipoint) - tmp2 * tmp_x + int2_grad1_u12_ao(2,i,j,ipoint) = tmp1 * y - tmp0 * x_v_ij_erf_rk_cst_mu_tmp_j1b(2,i,j,ipoint) - tmp2 * tmp_y + int2_grad1_u12_ao(3,i,j,ipoint) = tmp1 * z - tmp0 * x_v_ij_erf_rk_cst_mu_tmp_j1b(3,i,j,ipoint) - tmp2 * tmp_z enddo enddo enddo @@ -62,11 +62,14 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao, (3, ao_num, ao_num, n_poin x = final_grid_points(1,ipoint) y = final_grid_points(2,ipoint) z = final_grid_points(3,ipoint) + do j = 1, ao_num do i = 1, ao_num - int2_grad1_u12_ao(1,i,j,ipoint) = v_ij_erf_rk_cst_mu(i,j,ipoint) * x - x_v_ij_erf_rk_cst_mu(i,j,ipoint,1) - int2_grad1_u12_ao(2,i,j,ipoint) = v_ij_erf_rk_cst_mu(i,j,ipoint) * y - x_v_ij_erf_rk_cst_mu(i,j,ipoint,2) - int2_grad1_u12_ao(3,i,j,ipoint) = v_ij_erf_rk_cst_mu(i,j,ipoint) * z - x_v_ij_erf_rk_cst_mu(i,j,ipoint,3) + tmp1 = v_ij_erf_rk_cst_mu(i,j,ipoint) + + int2_grad1_u12_ao(1,i,j,ipoint) = tmp1 * x - x_v_ij_erf_rk_cst_mu_tmp(1,i,j,ipoint) + int2_grad1_u12_ao(2,i,j,ipoint) = tmp1 * y - x_v_ij_erf_rk_cst_mu_tmp(2,i,j,ipoint) + int2_grad1_u12_ao(3,i,j,ipoint) = tmp1 * z - x_v_ij_erf_rk_cst_mu_tmp(3,i,j,ipoint) enddo enddo enddo @@ -93,9 +96,8 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num, implicit none integer :: ipoint, i, j, k, l - double precision :: contrib, weight1, contrib_x, contrib_y, contrib_z - double precision :: ao_k_r, ao_k_dx, ao_k_dy, ao_k_dz - double precision :: ao_i_r, ao_i_dx, ao_i_dy, ao_i_dz + double precision :: weight1, contrib_x, contrib_y, contrib_z, tmp_x, tmp_y, tmp_z + double precision :: ao_k_r, ao_i_r, ao_i_dx, ao_i_dy, ao_i_dz double precision, allocatable :: ac_mat(:,:,:,:) allocate(ac_mat(ao_num,ao_num,ao_num,ao_num)) @@ -105,27 +107,26 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num, weight1 = 0.5d0 * final_weight_at_r_vector(ipoint) do i = 1, ao_num - ao_i_r = aos_in_r_array_transp (ipoint,i) - ao_i_dx = aos_grad_in_r_array_transp_bis(ipoint,i,1) - ao_i_dy = aos_grad_in_r_array_transp_bis(ipoint,i,2) - ao_i_dz = aos_grad_in_r_array_transp_bis(ipoint,i,3) + ao_i_r = weight1 * aos_in_r_array_transp (ipoint,i) + ao_i_dx = weight1 * aos_grad_in_r_array_transp_bis(ipoint,i,1) + ao_i_dy = weight1 * aos_grad_in_r_array_transp_bis(ipoint,i,2) + ao_i_dz = weight1 * aos_grad_in_r_array_transp_bis(ipoint,i,3) do k = 1, ao_num - ao_k_r = aos_in_r_array_transp (ipoint,k) - ao_k_dx = aos_grad_in_r_array_transp_bis(ipoint,k,1) - ao_k_dy = aos_grad_in_r_array_transp_bis(ipoint,k,2) - ao_k_dz = aos_grad_in_r_array_transp_bis(ipoint,k,3) + ao_k_r = aos_in_r_array_transp(ipoint,k) + + tmp_x = ao_k_r * ao_i_dx - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,1) + tmp_y = ao_k_r * ao_i_dy - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,2) + tmp_z = ao_k_r * ao_i_dz - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,3) do j = 1, ao_num do l = 1, ao_num - contrib_x = int2_grad1_u12_ao(1,l,j,ipoint) * ( ao_k_r * ao_i_dx - ao_i_r * ao_k_dx ) - contrib_y = int2_grad1_u12_ao(2,l,j,ipoint) * ( ao_k_r * ao_i_dy - ao_i_r * ao_k_dy ) - contrib_z = int2_grad1_u12_ao(3,l,j,ipoint) * ( ao_k_r * ao_i_dz - ao_i_r * ao_k_dz ) + contrib_x = int2_grad1_u12_ao(1,l,j,ipoint) * tmp_x + contrib_y = int2_grad1_u12_ao(2,l,j,ipoint) * tmp_y + contrib_z = int2_grad1_u12_ao(3,l,j,ipoint) * tmp_z - contrib = weight1 * ( contrib_x + contrib_y + contrib_z ) - - ac_mat(k,i,l,j) += contrib + ac_mat(k,i,l,j) += contrib_x + contrib_y + contrib_z enddo enddo enddo diff --git a/src/non_h_ints_mu/numerical_integ.irp.f b/src/non_h_ints_mu/numerical_integ.irp.f index dae68649..dcd7a52a 100644 --- a/src/non_h_ints_mu/numerical_integ.irp.f +++ b/src/non_h_ints_mu/numerical_integ.irp.f @@ -55,6 +55,7 @@ double precision function num_int2_u2_j1b2(i, j, ipoint) double precision, external :: ao_value double precision, external :: j1b_nucl + double precision, external :: j12_mu r1(1) = final_grid_points(1,ipoint) r1(2) = final_grid_points(2,ipoint) @@ -74,13 +75,14 @@ double precision function num_int2_u2_j1b2(i, j, ipoint) tmp1 = j1b_nucl(r2) tmp2 = tmp1 * tmp1 * ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint) - tmp3 = 0.d0 - do i_fit = 1, n_max_fit_slat - expo = expo_gauss_j_mu_x_2(i_fit) - coef = coef_gauss_j_mu_x_2(i_fit) - - tmp3 += coef * dexp(-expo*x2) - enddo + !tmp3 = 0.d0 + !do i_fit = 1, n_max_fit_slat + ! expo = expo_gauss_j_mu_x_2(i_fit) + ! coef = coef_gauss_j_mu_x_2(i_fit) + ! tmp3 += coef * dexp(-expo*x2) + !enddo + tmp3 = j12_mu(r1, r2) + tmp3 = tmp3 * tmp3 num_int2_u2_j1b2 += tmp2 * tmp3 enddo @@ -127,13 +129,15 @@ double precision function num_int2_grad1u2_grad2u2_j1b2(i, j, ipoint) tmp1 = j1b_nucl(r2) tmp2 = tmp1 * tmp1 * ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint) - tmp3 = 0.d0 - do i_fit = 1, n_max_fit_slat - expo = expo_gauss_1_erf_x_2(i_fit) - coef = coef_gauss_1_erf_x_2(i_fit) + !tmp3 = 0.d0 + !do i_fit = 1, n_max_fit_slat + ! expo = expo_gauss_1_erf_x_2(i_fit) + ! coef = coef_gauss_1_erf_x_2(i_fit) + ! tmp3 += coef * dexp(-expo*x2) + !enddo + tmp3 = derf(mu_erf*r12) - 1.d0 + tmp3 = tmp3 * tmp3 - tmp3 += coef * dexp(-expo*x2) - enddo tmp3 = -0.25d0 * tmp3 num_int2_grad1u2_grad2u2_j1b2 += tmp2 * tmp3 @@ -246,6 +250,12 @@ end subroutine num_x_v_ij_erf_rk_cst_mu_j1b subroutine num_int2_grad1_u12_ao(i, j, ipoint, integ) + BEGIN_DOC + ! + ! \int dr2 [-grad_1 u12] \phi_i(r2) \phi_j(r2) x v12_1b(r1, r2) + ! + END_DOC + implicit none integer, intent(in) :: i, j, ipoint @@ -256,7 +266,6 @@ subroutine num_int2_grad1_u12_ao(i, j, ipoint, integ) double precision :: tmp_x, tmp_y, tmp_z double precision, external :: ao_value - double precision, external :: j12_nucl r1(1) = final_grid_points(1,ipoint) r1(2) = final_grid_points(2,ipoint) @@ -269,9 +278,9 @@ subroutine num_int2_grad1_u12_ao(i, j, ipoint, integ) r2(1) = final_grid_points(1,jpoint) r2(2) = final_grid_points(2,jpoint) r2(3) = final_grid_points(3,jpoint) - tmp = ao_value(i, r2) * ao_value(j, r2) * j12_nucl(r1, r2) * final_weight_at_r_vector(jpoint) + tmp = ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint) - call grad1_j12_mu_exc(r1, r2, grad) + call grad1_jmu_modif_num(r1, r2, grad) tmp_x += tmp * (-1.d0 * grad(1)) tmp_y += tmp * (-1.d0 * grad(2)) @@ -289,16 +298,33 @@ end subroutine num_int2_grad1_u12_ao double precision function num_gradu_squared_u_ij_mu(i, j, ipoint) + BEGIN_DOC + ! + ! -0.50 x \int r2 \phi_i(2) \phi_j(2) x v2^2 + ! [ v1^2 ((grad_1 u12)^2 + (grad_2 u12^2)]) + ! + u12^2 (grad_1 v1)^2 + ! + 2 u12 v1 (grad_1 u12) . (grad_1 v1) + ! + END_DOC + + implicit none integer, intent(in) :: i, j, ipoint integer :: jpoint - double precision :: tmp, r1(3), r2(3), r12 - double precision :: tmp_x, tmp_y, tmp_z, tmp1, tmp2 + double precision :: r1(3), r2(3) + double precision :: tmp_x, tmp_y, tmp_z, r12 + double precision :: dx1_v1, dy1_v1, dz1_v1, grad_u12(3) + double precision :: tmp1, v1_tmp, v2_tmp, u12_tmp + double precision :: fst_term, scd_term, thd_term, tmp double precision, external :: ao_value - double precision, external :: j12_nucl + double precision, external :: j1b_nucl + double precision, external :: j12_mu + double precision, external :: grad_x_j1b_nucl + double precision, external :: grad_y_j1b_nucl + double precision, external :: grad_z_j1b_nucl r1(1) = final_grid_points(1,ipoint) r1(2) = final_grid_points(2,ipoint) @@ -306,16 +332,32 @@ double precision function num_gradu_squared_u_ij_mu(i, j, ipoint) num_gradu_squared_u_ij_mu = 0.d0 do jpoint = 1, n_points_final_grid + r2(1) = final_grid_points(1,jpoint) r2(2) = final_grid_points(2,jpoint) r2(3) = final_grid_points(3,jpoint) + tmp_x = r1(1) - r2(1) tmp_y = r1(2) - r2(2) tmp_z = r1(3) - r2(3) - r12 = dsqrt( tmp_x*tmp_x + tmp_y*tmp_y + tmp_z*tmp_z ) - tmp1 = 1.d0 - derf(mu_erf * r12) - tmp2 = j12_nucl(r1, r2) - tmp = -0.25d0 * tmp1 * tmp1 * tmp2 * tmp2 * ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint) + r12 = dsqrt(tmp_x*tmp_x + tmp_y*tmp_y + tmp_z*tmp_z) + + dx1_v1 = grad_x_j1b_nucl(r1) + dy1_v1 = grad_y_j1b_nucl(r1) + dz1_v1 = grad_z_j1b_nucl(r1) + + call grad1_j12_mu_exc(r1, r2, grad_u12) + + tmp1 = 1.d0 - derf(mu_erf * r12) + v1_tmp = j1b_nucl(r1) + v2_tmp = j1b_nucl(r2) + u12_tmp = j12_mu(r1, r2) + + fst_term = 0.5d0 * tmp1 * tmp1 * v1_tmp * v1_tmp + scd_term = u12_tmp * u12_tmp * (dx1_v1*dx1_v1 + dy1_v1*dy1_v1 + dz1_v1*dz1_v1) + thd_term = 2.d0 * v1_tmp * u12_tmp * (dx1_v1*grad_u12(1) + dy1_v1*grad_u12(2) + dz1_v1*grad_u12(3)) + + tmp = -0.5d0 * ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint) * (fst_term + scd_term + thd_term) * v2_tmp * v2_tmp num_gradu_squared_u_ij_mu += tmp enddo @@ -325,4 +367,257 @@ end function num_gradu_squared_u_ij_mu ! --- +double precision function num_grad12_j12(i, j, ipoint) + BEGIN_DOC + ! + ! -0.50 x \int r2 \phi_i(2) \phi_j(2) x v2^2 [v1^2 ((grad_1 u12)^2 + (grad_2 u12^2)]) ] + ! + END_DOC + + + implicit none + + integer, intent(in) :: i, j, ipoint + + integer :: jpoint + double precision :: r1(3), r2(3) + double precision :: tmp_x, tmp_y, tmp_z, r12 + double precision :: dx1_v1, dy1_v1, dz1_v1, grad_u12(3) + double precision :: tmp1, v1_tmp, v2_tmp, u12_tmp + double precision :: fst_term, scd_term, thd_term, tmp + + double precision, external :: ao_value + double precision, external :: j1b_nucl + double precision, external :: j12_mu + double precision, external :: grad_x_j1b_nucl + double precision, external :: grad_y_j1b_nucl + double precision, external :: grad_z_j1b_nucl + + r1(1) = final_grid_points(1,ipoint) + r1(2) = final_grid_points(2,ipoint) + r1(3) = final_grid_points(3,ipoint) + + num_grad12_j12 = 0.d0 + do jpoint = 1, n_points_final_grid + + r2(1) = final_grid_points(1,jpoint) + r2(2) = final_grid_points(2,jpoint) + r2(3) = final_grid_points(3,jpoint) + + tmp_x = r1(1) - r2(1) + tmp_y = r1(2) - r2(2) + tmp_z = r1(3) - r2(3) + r12 = dsqrt(tmp_x*tmp_x + tmp_y*tmp_y + tmp_z*tmp_z) + + dx1_v1 = grad_x_j1b_nucl(r1) + dy1_v1 = grad_y_j1b_nucl(r1) + dz1_v1 = grad_z_j1b_nucl(r1) + + call grad1_j12_mu_exc(r1, r2, grad_u12) + + tmp1 = 1.d0 - derf(mu_erf * r12) + v1_tmp = j1b_nucl(r1) + v2_tmp = j1b_nucl(r2) + u12_tmp = j12_mu(r1, r2) + + fst_term = 0.5d0 * tmp1 * tmp1 * v1_tmp * v1_tmp + + tmp = -0.5d0 * ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint) * fst_term * v2_tmp * v2_tmp + + num_grad12_j12 += tmp + enddo + + return +end function num_grad12_j12 + +! --- + +double precision function num_u12sq_j1bsq(i, j, ipoint) + + BEGIN_DOC + ! + ! -0.50 x \int r2 \phi_i(2) \phi_j(2) x v2^2 [ u12^2 (grad_1 v1)^2 ] + ! + END_DOC + + + implicit none + + integer, intent(in) :: i, j, ipoint + + integer :: jpoint + double precision :: r1(3), r2(3) + double precision :: tmp_x, tmp_y, tmp_z, r12 + double precision :: dx1_v1, dy1_v1, dz1_v1, grad_u12(3) + double precision :: tmp1, v1_tmp, v2_tmp, u12_tmp + double precision :: fst_term, scd_term, thd_term, tmp + + double precision, external :: ao_value + double precision, external :: j1b_nucl + double precision, external :: j12_mu + double precision, external :: grad_x_j1b_nucl + double precision, external :: grad_y_j1b_nucl + double precision, external :: grad_z_j1b_nucl + + r1(1) = final_grid_points(1,ipoint) + r1(2) = final_grid_points(2,ipoint) + r1(3) = final_grid_points(3,ipoint) + + num_u12sq_j1bsq = 0.d0 + do jpoint = 1, n_points_final_grid + + r2(1) = final_grid_points(1,jpoint) + r2(2) = final_grid_points(2,jpoint) + r2(3) = final_grid_points(3,jpoint) + + tmp_x = r1(1) - r2(1) + tmp_y = r1(2) - r2(2) + tmp_z = r1(3) - r2(3) + r12 = dsqrt(tmp_x*tmp_x + tmp_y*tmp_y + tmp_z*tmp_z) + + dx1_v1 = grad_x_j1b_nucl(r1) + dy1_v1 = grad_y_j1b_nucl(r1) + dz1_v1 = grad_z_j1b_nucl(r1) + + call grad1_j12_mu_exc(r1, r2, grad_u12) + + tmp1 = 1.d0 - derf(mu_erf * r12) + v1_tmp = j1b_nucl(r1) + v2_tmp = j1b_nucl(r2) + u12_tmp = j12_mu(r1, r2) + + scd_term = u12_tmp * u12_tmp * (dx1_v1*dx1_v1 + dy1_v1*dy1_v1 + dz1_v1*dz1_v1) + + tmp = -0.5d0 * ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint) * scd_term * v2_tmp * v2_tmp + + num_u12sq_j1bsq += tmp + enddo + + return +end function num_u12sq_j1bsq + +! --- + +double precision function num_u12_grad1_u12_j1b_grad1_j1b(i, j, ipoint) + + BEGIN_DOC + ! + ! -0.50 x \int r2 \phi_i(2) \phi_j(2) x v2^2 [ 2 u12 v1 (grad_1 u12) . (grad_1 v1) ] + ! + END_DOC + + + implicit none + + integer, intent(in) :: i, j, ipoint + + integer :: jpoint + double precision :: r1(3), r2(3) + double precision :: tmp_x, tmp_y, tmp_z, r12 + double precision :: dx1_v1, dy1_v1, dz1_v1, grad_u12(3) + double precision :: tmp1, v1_tmp, v2_tmp, u12_tmp + double precision :: fst_term, scd_term, thd_term, tmp + + double precision, external :: ao_value + double precision, external :: j1b_nucl + double precision, external :: j12_mu + double precision, external :: grad_x_j1b_nucl + double precision, external :: grad_y_j1b_nucl + double precision, external :: grad_z_j1b_nucl + + r1(1) = final_grid_points(1,ipoint) + r1(2) = final_grid_points(2,ipoint) + r1(3) = final_grid_points(3,ipoint) + + num_u12_grad1_u12_j1b_grad1_j1b = 0.d0 + do jpoint = 1, n_points_final_grid + + r2(1) = final_grid_points(1,jpoint) + r2(2) = final_grid_points(2,jpoint) + r2(3) = final_grid_points(3,jpoint) + + tmp_x = r1(1) - r2(1) + tmp_y = r1(2) - r2(2) + tmp_z = r1(3) - r2(3) + r12 = dsqrt(tmp_x*tmp_x + tmp_y*tmp_y + tmp_z*tmp_z) + + dx1_v1 = grad_x_j1b_nucl(r1) + dy1_v1 = grad_y_j1b_nucl(r1) + dz1_v1 = grad_z_j1b_nucl(r1) + + call grad1_j12_mu_exc(r1, r2, grad_u12) + + tmp1 = 1.d0 - derf(mu_erf * r12) + v1_tmp = j1b_nucl(r1) + v2_tmp = j1b_nucl(r2) + u12_tmp = j12_mu(r1, r2) + + thd_term = 2.d0 * v1_tmp * u12_tmp * (dx1_v1*grad_u12(1) + dy1_v1*grad_u12(2) + dz1_v1*grad_u12(3)) + + tmp = -0.5d0 * ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint) * thd_term * v2_tmp * v2_tmp + + num_u12_grad1_u12_j1b_grad1_j1b += tmp + enddo + + return +end function num_u12_grad1_u12_j1b_grad1_j1b + +! --- + +subroutine num_int2_u_grad1u_total_j1b2(i, j, ipoint, integ) + + BEGIN_DOC + ! + ! \int dr2 u12 (grad_1 u12) \phi_i(r2) \phi_j(r2) x v_1b(r2)^2 + ! + END_DOC + + implicit none + + integer, intent(in) :: i, j, ipoint + double precision, intent(out) :: integ(3) + + integer :: jpoint + double precision :: r1(3), r2(3), grad(3) + double precision :: dx, dy, dz, r12, tmp0, tmp1, tmp2 + double precision :: tmp_x, tmp_y, tmp_z + + double precision, external :: ao_value + double precision, external :: j1b_nucl + double precision, external :: j12_mu + + r1(1) = final_grid_points(1,ipoint) + r1(2) = final_grid_points(2,ipoint) + r1(3) = final_grid_points(3,ipoint) + + tmp_x = 0.d0 + tmp_y = 0.d0 + tmp_z = 0.d0 + do jpoint = 1, n_points_final_grid + r2(1) = final_grid_points(1,jpoint) + r2(2) = final_grid_points(2,jpoint) + r2(3) = final_grid_points(3,jpoint) + dx = r1(1) - r2(1) + dy = r1(2) - r2(2) + dz = r1(3) - r2(3) + r12 = dsqrt( dx * dx + dy * dy + dz * dz ) + if(r12 .lt. 1d-10) cycle + + tmp0 = j1b_nucl(r2) + tmp1 = 0.5d0 * j12_mu(r1, r2) * (1.d0 - derf(mu_erf * r12)) / r12 + tmp2 = tmp0 * tmp0 * tmp1 * ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint) + + tmp_x += tmp2 * dx + tmp_y += tmp2 * dy + tmp_z += tmp2 * dz + enddo + + integ(1) = tmp_x + integ(2) = tmp_y + integ(3) = tmp_z + + return +end subroutine num_int2_u_grad1u_total_j1b2 + +! --- diff --git a/src/non_h_ints_mu/total_tc_int.irp.f b/src/non_h_ints_mu/total_tc_int.irp.f new file mode 100644 index 00000000..979296d1 --- /dev/null +++ b/src/non_h_ints_mu/total_tc_int.irp.f @@ -0,0 +1,60 @@ + +! --- + +BEGIN_PROVIDER [double precision, ao_tc_int_chemist, (ao_num, ao_num, ao_num, ao_num)] + + implicit none + integer :: i, j, k, l + double precision :: wall1, wall0 + + call wall_time(wall0) + + do j = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do k = 1, ao_num + ao_tc_int_chemist(k,i,l,j) = tc_grad_square_ao(k,i,l,j) + tc_grad_and_lapl_ao(k,i,l,j) + ao_two_e_coul(k,i,l,j) + enddo + enddo + enddo + enddo + + call wall_time(wall1) + print *, ' wall time for ao_tc_int_chemist ', wall1 - wall0 + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, ao_two_e_coul, (ao_num, ao_num, ao_num, ao_num) ] + + BEGIN_DOC + ! + ! ao_two_e_coul(k,i,l,j) = ( k i | 1/r12 | l j ) = < l k | 1/r12 | j i > + ! + END_DOC + + integer :: i, j, k, l + double precision :: integral + double precision, external :: get_ao_two_e_integral + + PROVIDE ao_integrals_map + + do j = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do k = 1, ao_num + + ! < 1:k, 2:l | 1:i, 2:j > + integral = get_ao_two_e_integral(i, j, k, l, ao_integrals_map) + + ao_two_e_coul(k,i,l,j) = integral + enddo + enddo + enddo + enddo + +END_PROVIDER + +! --- + diff --git a/src/tc_bi_ortho/compute_deltamu_right.irp.f b/src/tc_bi_ortho/compute_deltamu_right.irp.f index 32566cc8..6464796e 100644 --- a/src/tc_bi_ortho/compute_deltamu_right.irp.f +++ b/src/tc_bi_ortho/compute_deltamu_right.irp.f @@ -34,6 +34,7 @@ subroutine delta_right() !do k = 1, 1 ! get < I_left | H_mu - H | psi_right > + !call get_h_bitc_right(psi_det, psi_r_coef_bi_ortho(:,k), N_det, N_int, delta(:,k)) call get_delta_bitc_right(psi_det, psi_r_coef_bi_ortho(:,k), N_det, N_int, delta(:,k)) ! order as QMCCHEM diff --git a/src/tc_bi_ortho/dressing_vectors_lr.irp.f b/src/tc_bi_ortho/dressing_vectors_lr.irp.f index e69a970b..08913bab 100644 --- a/src/tc_bi_ortho/dressing_vectors_lr.irp.f +++ b/src/tc_bi_ortho/dressing_vectors_lr.irp.f @@ -23,10 +23,12 @@ subroutine get_delta_bitc_right(psidet, psicoef, ndet, Nint, delta) double precision :: htc_mono, htc_twoe, htc_three, htc_tot double precision :: delta_mat + print *, ' get_delta_bitc_right ...' + i = 1 j = 1 call htilde_mu_mat_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot) - call hmat_bi_ortho (psidet(1,1,i), psidet(1,1,j), Nint, h_mono , h_twoe , h_tot) + call hmat_bi_ortho (psidet(1,1,i), psidet(1,1,j), Nint, h_mono, h_twoe, h_tot) delta = 0.d0 !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(dynamic,8) & @@ -39,7 +41,7 @@ subroutine get_delta_bitc_right(psidet, psicoef, ndet, Nint, delta) ! < I | Htilde | J > call htilde_mu_mat_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot) ! < I | H | J > - call hmat_bi_ortho (psidet(1,1,i), psidet(1,1,j), Nint, h_mono , h_twoe , h_tot) + call hmat_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, h_mono, h_twoe, h_tot) delta_mat = htc_tot - h_tot @@ -52,3 +54,102 @@ end subroutine get_delta_bitc_right ! --- +subroutine get_htc_bitc_right(psidet, psicoef, ndet, Nint, delta) + + BEGIN_DOC + ! + ! delta(I) = < I_left | H_TC | Psi_right > + ! + END_DOC + + use bitmasks + + implicit none + + integer, intent(in) :: ndet, Nint + double precision, intent(in) :: psicoef(ndet) + integer(bit_kind), intent(in) :: psidet(Nint,2,ndet) + double precision, intent(out) :: delta(ndet) + + integer :: i, j + double precision :: htc_mono, htc_twoe, htc_three, htc_tot + + print *, ' get_htc_bitc_right ...' + + i = 1 + j = 1 + call htilde_mu_mat_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot) + + delta = 0.d0 + !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(dynamic,8) & + !$OMP SHARED(delta, ndet, psidet, psicoef, Nint) & + !$OMP PRIVATE(i, j, htc_mono, htc_twoe, htc_three, htc_tot) + do i = 1, ndet + do j = 1, ndet + + ! < I | Htilde | J > + call htilde_mu_mat_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot) + + delta(i) = delta(i) + psicoef(j) * htc_tot + enddo + enddo + !$OMP END PARALLEL DO + +end subroutine get_htc_bitc_right + +! --- + +subroutine get_h_bitc_right(psidet, psicoef, ndet, Nint, delta) + + BEGIN_DOC + ! + ! delta(I) = < I_left | H | Psi_right > + ! + END_DOC + + use bitmasks + + implicit none + + integer, intent(in) :: ndet, Nint + double precision, intent(in) :: psicoef(ndet) + integer(bit_kind), intent(in) :: psidet(Nint,2,ndet) + double precision, intent(out) :: delta(ndet) + + integer :: i, j + double precision :: h_mono, h_twoe, h_tot + + print *, ' get_h_bitc_right ...' + + i = 1 + j = 1 + call hmat_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, h_mono, h_twoe, h_tot) + + !double precision :: norm + !norm = 0.d0 + !do i = 1, ndet + ! norm += psicoef(i) * psicoef(i) + !enddo + !print*, ' norm = ', norm + + call hmat_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, h_mono, h_twoe, h_tot) + + delta = 0.d0 +! !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(dynamic,8) & +! !$OMP SHARED(delta, ndet, psidet, psicoef, Nint) & +! !$OMP PRIVATE(i, j, h_mono, h_twoe, h_tot) + do i = 1, ndet + do j = 1, ndet + + ! < I | H | J > + call hmat_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, h_mono, h_twoe, h_tot) + + delta(i) = delta(i) + psicoef(j) * h_tot + enddo + enddo +! !$OMP END PARALLEL DO + +end subroutine get_h_bitc_right + +! --- + diff --git a/src/tc_bi_ortho/h_biortho.irp.f b/src/tc_bi_ortho/h_biortho.irp.f index 0494399f..492e1282 100644 --- a/src/tc_bi_ortho/h_biortho.irp.f +++ b/src/tc_bi_ortho/h_biortho.irp.f @@ -5,7 +5,7 @@ subroutine hmat_bi_ortho(key_j, key_i, Nint, hmono, htwoe, htot) BEGIN_DOC ! - ! where | key_j > is developed on the LEFT basis and | key_i > is developed on the RIGHT basis + ! < key_j | H | key_i > where | key_j > is developed on the LEFT basis and | key_i > is developed on the RIGHT basis ! END_DOC @@ -13,11 +13,11 @@ subroutine hmat_bi_ortho(key_j, key_i, Nint, hmono, htwoe, htot) implicit none - integer, intent(in) :: Nint - integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) - double precision, intent(out) :: hmono, htwoe, htot + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + double precision, intent(out) :: hmono, htwoe, htot - integer :: degree + integer :: degree hmono = 0.d0 htwoe = 0.d0 @@ -31,11 +31,11 @@ subroutine hmat_bi_ortho(key_j, key_i, Nint, hmono, htwoe, htot) call diag_hmat_bi_ortho(Nint, key_i, hmono, htwoe) htot = htot + nuclear_repulsion - else if (degree == 1)then + else if (degree == 1) then call single_hmat_bi_ortho(Nint, key_j, key_i, hmono, htwoe) - else if(degree == 2)then + else if(degree == 2) then call double_hmat_bi_ortho(Nint, key_j, key_i, hmono, htwoe) @@ -59,8 +59,7 @@ subroutine diag_hmat_bi_ortho(Nint, key_i, hmono, htwoe) double precision, intent(out) :: hmono, htwoe integer :: occ(Nint*bit_kind_size,2) - integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk - integer(bit_kind) :: key_i_core(Nint,2) + integer :: Ne(2), i, j, ii, jj, ispin, jspin hmono = 0.d0 htwoe = 0.d0 @@ -125,13 +124,11 @@ subroutine single_hmat_bi_ortho(Nint, key_j, key_i, hmono, htwoe) double precision, intent(out) :: hmono, htwoe integer :: occ(Nint*bit_kind_size,2) - integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk + integer :: Ne(2), i, j, ii, ispin, jspin integer :: degree,exc(0:2,2,2) integer :: h1, p1, h2, p2, s1, s2 integer :: other_spin(2) - integer(bit_kind) :: key_j_core(Nint,2), key_i_core(Nint,2) double precision :: phase - double precision :: direct_int, exchange_int_12, exchange_int_23, exchange_int_13 other_spin(1) = 2 other_spin(2) = 1 @@ -201,11 +198,10 @@ subroutine double_hmat_bi_ortho(Nint, key_j, key_i, hmono, htwoe) double precision, intent(out) :: hmono, htwoe integer :: occ(Nint*bit_kind_size,2) - integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk + integer :: Ne(2), i, j, ii, ispin, jspin integer :: degree,exc(0:2,2,2) integer :: h1, p1, h2, p2, s1, s2 integer :: other_spin(2) - integer(bit_kind) :: key_i_core(Nint,2) double precision :: phase other_spin(1) = 2 @@ -225,7 +221,7 @@ subroutine double_hmat_bi_ortho(Nint, key_j, key_i, hmono, htwoe) call get_double_excitation(key_i, key_j, exc, phase, Nint) call decode_exc(exc, 2, h1, p1, h2, p2, s1, s2) - if(s1.ne.s2) then + if(s1 .ne. s2) then htwoe = mo_bi_ortho_coul_e(p2,p1,h2,h1) @@ -233,10 +229,8 @@ subroutine double_hmat_bi_ortho(Nint, key_j, key_i, hmono, htwoe) ! same spin two-body - ! direct terms - htwoe = mo_bi_ortho_coul_e(p2,p1,h2,h1) - ! exchange terms - htwoe -= mo_bi_ortho_coul_e(p1,p2,h2,h1) + ! direct terms exchange terms + htwoe = mo_bi_ortho_coul_e(p2,p1,h2,h1) - mo_bi_ortho_coul_e(p1,p2,h2,h1) endif diff --git a/src/tc_bi_ortho/psi_r_l_prov.irp.f b/src/tc_bi_ortho/psi_r_l_prov.irp.f index 551da858..2a5887d5 100644 --- a/src/tc_bi_ortho/psi_r_l_prov.irp.f +++ b/src/tc_bi_ortho/psi_r_l_prov.irp.f @@ -41,6 +41,15 @@ BEGIN_PROVIDER [ double precision, psi_l_coef_bi_ortho, (psi_det_size,N_states) enddo deallocate(psi_l_coef_bi_ortho_read) + else + + print*, 'psi_l_coef_bi_ortho are psi_coef' + do k=1,N_states + do i=1,N_det + psi_l_coef_bi_ortho(i,k) = psi_coef(i,k) + enddo + enddo + endif endif endif @@ -100,6 +109,15 @@ BEGIN_PROVIDER [ double precision, psi_r_coef_bi_ortho, (psi_det_size,N_states) enddo deallocate(psi_r_coef_bi_ortho_read) + else + + print*, 'psi_r_coef_bi_ortho are psi_coef' + do k=1,N_states + do i=1,N_det + psi_r_coef_bi_ortho(i,k) = psi_coef(i,k) + enddo + enddo + endif endif endif diff --git a/src/tc_bi_ortho/slater_tc.irp.f b/src/tc_bi_ortho/slater_tc.irp.f index e0a52741..33b738ba 100644 --- a/src/tc_bi_ortho/slater_tc.irp.f +++ b/src/tc_bi_ortho/slater_tc.irp.f @@ -1,4 +1,6 @@ -!!!!!! + +! --- + subroutine htilde_mu_mat_bi_ortho_tot(key_j, key_i, Nint, htot) BEGIN_DOC @@ -15,13 +17,14 @@ subroutine htilde_mu_mat_bi_ortho_tot(key_j, key_i, Nint, htot) integer, intent(in) :: Nint integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2) double precision, intent(out) :: htot - double precision :: hmono,htwoe,hthree + double precision :: hmono, htwoe, hthree integer :: degree + call get_excitation_degree(key_j, key_i, degree, Nint) if(degree.gt.2)then - htot = 0.d0 + htot = 0.d0 else - call htilde_mu_mat_bi_ortho(key_j,key_i, Nint, hmono,htwoe,hthree,htot) + call htilde_mu_mat_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, htot) endif end subroutine htilde_mu_mat_bi_ortho_tot @@ -29,52 +32,63 @@ end subroutine htilde_mu_mat_bi_ortho_tot ! -- subroutine htilde_mu_mat_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, htot) - implicit none - use bitmasks + BEGIN_DOC + ! ! where |key_j> is developed on the LEFT basis and |key_i> is developed on the RIGHT basis !! ! Returns the detail of the matrix element in terms of single, two and three electron contribution. !! WARNING !! ! ! Non hermitian !! + ! END_DOC - integer, intent(in) :: Nint - integer(bit_kind), intent(in) :: key_i(Nint,2),key_j(Nint,2) - double precision, intent(out) :: hmono,htwoe,hthree,htot - integer :: degree - hmono = 0.d0 - htwoe= 0.d0 - htot = 0.d0 + use bitmasks + + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + double precision, intent(out) :: hmono, htwoe, hthree, htot + integer :: degree + + hmono = 0.d0 + htwoe = 0.d0 + htot = 0.d0 hthree = 0.D0 + call get_excitation_degree(key_i, key_j, degree, Nint) - if(degree.gt.2)return + if(degree.gt.2) return + if(degree == 0)then - call diag_htilde_mu_mat_bi_ortho(Nint, key_i, hmono, htwoe, htot) + call diag_htilde_mu_mat_bi_ortho(Nint, key_i, hmono, htwoe, htot) else if (degree == 1)then - call single_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot) + call single_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot) else if(degree == 2)then - call double_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot) - endif - if(three_body_h_tc) then - if(degree == 2) then - if(.not.double_normal_ord) then - call double_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree) - endif - else if(degree == 1)then - call single_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree) - else if(degree == 0)then - call diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree) - endif - endif - htot = hmono + htwoe + hthree - if(degree==0)then - htot += nuclear_repulsion - endif + call double_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot) + endif + + if(three_body_h_tc) then + if(degree == 2) then + if(.not.double_normal_ord) then + call double_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree) + endif + else if(degree == 1) then + call single_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree) + else if(degree == 0) then + call diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree) + endif + endif + + htot = hmono + htwoe + hthree + if(degree==0) then + htot += nuclear_repulsion + endif end +! --- + subroutine diag_htilde_mu_mat_bi_ortho(Nint, key_i, hmono, htwoe, htot) BEGIN_DOC diff --git a/src/tc_bi_ortho/slater_tc_3e.irp.f b/src/tc_bi_ortho/slater_tc_3e.irp.f index 7c2c9c86..f4899a80 100644 --- a/src/tc_bi_ortho/slater_tc_3e.irp.f +++ b/src/tc_bi_ortho/slater_tc_3e.irp.f @@ -207,6 +207,8 @@ subroutine single_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree) end +! --- + subroutine double_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree) BEGIN_DOC @@ -244,7 +246,7 @@ subroutine double_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree) return endif - if(core_tc_op)then + if(core_tc_op) then do i = 1, Nint key_i_core(i,1) = xor(key_i(i,1),core_bitmask(i,1)) key_i_core(i,2) = xor(key_i(i,2),core_bitmask(i,2)) @@ -291,3 +293,6 @@ subroutine double_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree) endif hthree *= phase end + +! --- +