diff --git a/plugins/local/non_h_ints_mu/jast_1e.irp.f b/plugins/local/non_h_ints_mu/jast_1e.irp.f index 9700c182..b2eef504 100644 --- a/plugins/local/non_h_ints_mu/jast_1e.irp.f +++ b/plugins/local/non_h_ints_mu/jast_1e.irp.f @@ -129,8 +129,7 @@ END_PROVIDER elseif(j1e_type .eq. "Charge_Harmonizer") then - ! The - sign is in the integral over r2 - ! [(N-1)/2N] x \sum_{\mu,\nu} P_{\mu,\nu} \int dr2 [-1 * \grad_r1 J(r1,r2)] \phi_\mu(r2) \phi_nu(r2) + ! -[(N-1)/2N] x \sum_{\mu,\nu} P_{\mu,\nu} \int dr2 [\grad_r1 J_2e(r1,r2)] \phi_\mu(r2) \phi_nu(r2) PROVIDE elec_alpha_num elec_beta_num elec_num PROVIDE mo_coef @@ -151,7 +150,7 @@ END_PROVIDER endif Pt = Pa + Pb - g = 0.5d0 * (dble(elec_num) - 1.d0) / dble(elec_num) + g = -0.5d0 * (dble(elec_num) - 1.d0) / dble(elec_num) call dgemv("T", ao_num*ao_num, n_points_final_grid, g, int2_grad1_u2e_ao(1,1,1,1), ao_num*ao_num, Pt, 1, 0.d0, j1e_gradx, 1) call dgemv("T", ao_num*ao_num, n_points_final_grid, g, int2_grad1_u2e_ao(1,1,1,2), ao_num*ao_num, Pt, 1, 0.d0, j1e_grady, 1) diff --git a/plugins/local/non_h_ints_mu/jast_1e_utils.irp.f b/plugins/local/non_h_ints_mu/jast_1e_utils.irp.f index 80ed8c6e..ba7477cc 100644 --- a/plugins/local/non_h_ints_mu/jast_1e_utils.irp.f +++ b/plugins/local/non_h_ints_mu/jast_1e_utils.irp.f @@ -1,207 +1,6 @@ ! --- -BEGIN_PROVIDER [double precision, int2_u2e_ao, (ao_num, ao_num, n_points_final_grid)] - - BEGIN_DOC - ! - ! int2_u2e_ao(i,j,ipoint,:) = \int dr2 J_2e(r1,r2) \phi_i(r2) \phi_j(r2) - ! - ! where r1 = r(ipoint) - ! - END_DOC - - implicit none - integer :: ipoint, i, j, jpoint - double precision :: time0, time1 - double precision :: x, y, z, r2 - double precision :: dx, dy, dz - double precision :: tmp_ct - double precision :: tmp0, tmp1, tmp2, tmp3 - - PROVIDE j2e_type - PROVIDE Env_type - - call wall_time(time0) - print*, ' providing int2_u2e_ao ...' - - if(tc_integ_type .eq. "semi-analytic") then - - if( (j2e_type .eq. "Mu") .and. & - ( (env_type .eq. "None") .or. (env_type .eq. "Prod_Gauss") .or. (env_type .eq. "Sum_Gauss") ) ) then - - PROVIDE mu_erf - PROVIDE env_type env_val - PROVIDE Ir2_Mu_long_Du_0 Ir2_Mu_long_Du_x Ir2_Mu_long_Du_y Ir2_Mu_long_Du_z Ir2_Mu_long_Du_2 - PROVIDE Ir2_Mu_gauss_Du - - tmp_ct = 0.5d0 / (dsqrt(dacos(-1.d0)) * mu_erf) - - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, i, j, x, y, z, r2, dx, dy, dz, & - !$OMP tmp0, tmp1, tmp2, tmp3) & - !$OMP SHARED (ao_num, n_points_final_grid, final_grid_points, & - !$OMP tmp_ct, env_val, Ir2_Mu_long_Du_0, & - !$OMP Ir2_Mu_long_Du_x, Ir2_Mu_long_Du_y, & - !$OMP Ir2_Mu_long_Du_z, Ir2_Mu_gauss_Du, & - !$OMP Ir2_Mu_long_Du_2, int2_u2e_ao) - !$OMP DO SCHEDULE (static) - 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) - r2 = x*x + y*y + z*z - - dx = x * env_val(ipoint) - dy = y * env_val(ipoint) - dz = z * env_val(ipoint) - - tmp0 = 0.5d0 * env_val(ipoint) * r2 - tmp1 = 0.5d0 * env_val(ipoint) - tmp3 = tmp_ct * env_val(ipoint) - - do j = 1, ao_num - do i = 1, ao_num - - tmp2 = tmp1 * Ir2_Mu_long_Du_2(i,j,ipoint) - dx * Ir2_Mu_long_Du_x(i,j,ipoint) - dy * Ir2_Mu_long_Du_y(i,j,ipoint) - dz * Ir2_Mu_long_Du_z(i,j,ipoint) - - int2_u2e_ao(i,j,ipoint) = tmp0 * Ir2_Mu_long_Du_0(i,j,ipoint) + tmp2 - tmp3 * Ir2_Mu_gauss_Du(i,j,ipoint) - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - - else - - print *, ' Error in int2_u2e_ao: Unknown Jastrow' - stop - - endif ! j2e_type - - else - - write(*, '(A, A, A)') ' Error: The integration type ', trim(tc_integ_type), ' has not been implemented yet' - stop - - endif ! tc_integ_type - - call wall_time(time1) - print*, ' wall time for int2_u2e_ao (min) =', (time1-time0)/60.d0 - call print_memory_usage() - -END_PROVIDER - -! --- - -BEGIN_PROVIDER [double precision, int2_grad1_u2e_ao, (ao_num, ao_num, n_points_final_grid, 3)] - - BEGIN_DOC - ! - ! int2_grad1_u2e_ao(i,j,ipoint,:) = \int dr2 [-1 * \grad_r1 J_2e(r1,r2)] \phi_i(r2) \phi_j(r2) - ! - ! where r1 = r(ipoint) - ! - END_DOC - - implicit none - integer :: ipoint, i, j, m, jpoint - double precision :: time0, time1 - double precision :: x, y, z, r2 - double precision :: dx, dy, dz - double precision :: tmp_ct - double precision :: tmp0, tmp1, tmp2 - double precision :: tmp0_x, tmp0_y, tmp0_z - double precision :: tmp1_x, tmp1_y, tmp1_z - - PROVIDE j2e_type - PROVIDE Env_type - - call wall_time(time0) - print*, ' providing int2_grad1_u2e_ao ...' - - if(tc_integ_type .eq. "semi-analytic") then - - - if( (j2e_type .eq. "Mu") .and. & - ( (env_type .eq. "None") .or. (env_type .eq. "Prod_Gauss") .or. (env_type .eq. "Sum_Gauss") ) ) then - - PROVIDE mu_erf - PROVIDE env_type env_val env_grad - PROVIDE Ir2_Mu_long_Du_0 Ir2_Mu_long_Du_x Ir2_Mu_long_Du_y Ir2_Mu_long_Du_z Ir2_Mu_long_Du_2 - PROVIDE Ir2_Mu_gauss_Du - - tmp_ct = 0.5d0 / (dsqrt(dacos(-1.d0)) * mu_erf) - - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, i, j, x, y, z, r2, dx, dy, dz, tmp1, tmp2, & - !$OMP tmp0_x, tmp0_y, tmp0_z, tmp1_x, tmp1_y, tmp1_z) & - !$OMP SHARED (ao_num, n_points_final_grid, final_grid_points, & - !$OMP tmp_ct, env_val, env_grad, Ir2_Mu_long_Du_0, & - !$OMP Ir2_Mu_long_Du_x, Ir2_Mu_long_Du_y, & - !$OMP Ir2_Mu_long_Du_z, Ir2_Mu_gauss_Du, & - !$OMP Ir2_Mu_long_Du_2, int2_grad1_u2e_ao) - !$OMP DO SCHEDULE (static) - 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) - r2 = x*x + y*y + z*z - - dx = env_grad(1,ipoint) - dy = env_grad(2,ipoint) - dz = env_grad(3,ipoint) - - tmp0_x = 0.5d0 * (env_val(ipoint) * x + r2 * dx) - tmp0_y = 0.5d0 * (env_val(ipoint) * y + r2 * dy) - tmp0_z = 0.5d0 * (env_val(ipoint) * z + r2 * dz) - - tmp1 = 0.5d0 * env_val(ipoint) - - tmp1_x = tmp_ct * dx - tmp1_y = tmp_ct * dy - tmp1_z = tmp_ct * dz - - do j = 1, ao_num - do i = 1, ao_num - - tmp2 = 0.5d0 * Ir2_Mu_long_Du_2(i,j,ipoint) - x * Ir2_Mu_long_Du_x(i,j,ipoint) - y * Ir2_Mu_long_Du_y(i,j,ipoint) - z * Ir2_Mu_long_Du_z(i,j,ipoint) - - int2_grad1_u2e_ao(i,j,ipoint,1) = -Ir2_Mu_long_Du_0(i,j,ipoint) * tmp0_x + tmp1 * Ir2_Mu_long_Du_x(i,j,ipoint) - dx * tmp2 + tmp1_x * Ir2_Mu_gauss_Du(i,j,ipoint) - int2_grad1_u2e_ao(i,j,ipoint,2) = -Ir2_Mu_long_Du_0(i,j,ipoint) * tmp0_y + tmp1 * Ir2_Mu_long_Du_y(i,j,ipoint) - dy * tmp2 + tmp1_y * Ir2_Mu_gauss_Du(i,j,ipoint) - int2_grad1_u2e_ao(i,j,ipoint,3) = -Ir2_Mu_long_Du_0(i,j,ipoint) * tmp0_z + tmp1 * Ir2_Mu_long_Du_z(i,j,ipoint) - dz * tmp2 + tmp1_z * Ir2_Mu_gauss_Du(i,j,ipoint) - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - - else - - print *, ' Error in int2_grad1_u2e_ao: Unknown Jastrow' - stop - - endif ! j2e_type - - else - - write(*, '(A, A, A)') ' Error: The integration type ', trim(tc_integ_type), ' has not been implemented yet' - stop - - endif ! tc_integ_type - - call wall_time(time1) - print*, ' wall time for int2_grad1_u2e_ao (min) =', (time1-time0)/60.d0 - call print_memory_usage() - -END_PROVIDER - -! --- - subroutine get_j1e_coef_fit_ao(dim_fit, coef_fit) implicit none diff --git a/plugins/local/non_h_ints_mu/jast_2e_utils.irp.f b/plugins/local/non_h_ints_mu/jast_2e_utils.irp.f new file mode 100644 index 00000000..8c25b377 --- /dev/null +++ b/plugins/local/non_h_ints_mu/jast_2e_utils.irp.f @@ -0,0 +1,188 @@ + +! --- + +BEGIN_PROVIDER [double precision, int2_u2e_ao, (ao_num, ao_num, n_points_final_grid)] + + BEGIN_DOC + ! + ! int2_u2e_ao(i,j,ipoint,:) = \int dr2 J_2e(r1,r2) \phi_i(r2) \phi_j(r2) + ! + ! where r1 = r(ipoint) + ! + END_DOC + + implicit none + integer :: ipoint, i, j, jpoint + double precision :: time0, time1 + double precision :: x, y, z, r2 + double precision :: dx, dy, dz + double precision :: tmp_ct + double precision :: tmp0, tmp1, tmp2, tmp3 + + PROVIDE j2e_type + PROVIDE Env_type + + call wall_time(time0) + print*, ' providing int2_u2e_ao ...' + + if( (j2e_type .eq. "Mu") .and. & + ( (env_type .eq. "None") .or. (env_type .eq. "Prod_Gauss") .or. (env_type .eq. "Sum_Gauss") ) ) then + + PROVIDE mu_erf + PROVIDE env_type env_val + PROVIDE Ir2_Mu_long_Du_0 Ir2_Mu_long_Du_x Ir2_Mu_long_Du_y Ir2_Mu_long_Du_z Ir2_Mu_long_Du_2 + PROVIDE Ir2_Mu_gauss_Du + + tmp_ct = 0.5d0 / (dsqrt(dacos(-1.d0)) * mu_erf) + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, x, y, z, r2, dx, dy, dz, & + !$OMP tmp0, tmp1, tmp2, tmp3) & + !$OMP SHARED (ao_num, n_points_final_grid, final_grid_points, & + !$OMP tmp_ct, env_val, Ir2_Mu_long_Du_0, & + !$OMP Ir2_Mu_long_Du_x, Ir2_Mu_long_Du_y, & + !$OMP Ir2_Mu_long_Du_z, Ir2_Mu_gauss_Du, & + !$OMP Ir2_Mu_long_Du_2, int2_u2e_ao) + !$OMP DO SCHEDULE (static) + 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) + r2 = x*x + y*y + z*z + + dx = x * env_val(ipoint) + dy = y * env_val(ipoint) + dz = z * env_val(ipoint) + + tmp0 = 0.5d0 * env_val(ipoint) * r2 + tmp1 = 0.5d0 * env_val(ipoint) + tmp3 = tmp_ct * env_val(ipoint) + + do j = 1, ao_num + do i = 1, ao_num + + tmp2 = tmp1 * Ir2_Mu_long_Du_2(i,j,ipoint) - dx * Ir2_Mu_long_Du_x(i,j,ipoint) - dy * Ir2_Mu_long_Du_y(i,j,ipoint) - dz * Ir2_Mu_long_Du_z(i,j,ipoint) + + int2_u2e_ao(i,j,ipoint) = tmp0 * Ir2_Mu_long_Du_0(i,j,ipoint) + tmp2 - tmp3 * Ir2_Mu_gauss_Du(i,j,ipoint) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + else + + print *, ' Error in int2_u2e_ao: Unknown Jastrow' + stop + + endif ! j2e_type + + call wall_time(time1) + print*, ' wall time for int2_u2e_ao (min) =', (time1-time0)/60.d0 + call print_memory_usage() + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, int2_grad1_u2e_ao, (ao_num, ao_num, n_points_final_grid, 3)] + + BEGIN_DOC + ! + ! int2_grad1_u2e_ao(i,j,ipoint,:) = \int dr2 [\grad_r1 J_2e(r1,r2)] \phi_i(r2) \phi_j(r2) + ! + ! where r1 = r(ipoint) + ! + END_DOC + + implicit none + integer :: ipoint, i, j, m, jpoint + double precision :: time0, time1 + double precision :: x, y, z, r2 + double precision :: dx, dy, dz + double precision :: tmp_ct + double precision :: tmp0, tmp1, tmp2 + double precision :: tmp0_x, tmp0_y, tmp0_z + double precision :: tmp1_x, tmp1_y, tmp1_z + + PROVIDE j2e_type + PROVIDE Env_type + + call wall_time(time0) + print*, ' providing int2_grad1_u2e_ao ...' + + if( (j2e_type .eq. "Mu") .and. & + ( (env_type .eq. "None") .or. (env_type .eq. "Prod_Gauss") .or. (env_type .eq. "Sum_Gauss") ) ) then + + PROVIDE mu_erf + PROVIDE env_type env_val env_grad + PROVIDE Ir2_Mu_long_Du_0 Ir2_Mu_long_Du_x Ir2_Mu_long_Du_y Ir2_Mu_long_Du_z Ir2_Mu_long_Du_2 + PROVIDE Ir2_Mu_gauss_Du + + tmp_ct = 0.5d0 / (dsqrt(dacos(-1.d0)) * mu_erf) + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, x, y, z, r2, dx, dy, dz, tmp1, tmp2, & + !$OMP tmp0_x, tmp0_y, tmp0_z, tmp1_x, tmp1_y, tmp1_z) & + !$OMP SHARED (ao_num, n_points_final_grid, final_grid_points, & + !$OMP tmp_ct, env_val, env_grad, Ir2_Mu_long_Du_0, & + !$OMP Ir2_Mu_long_Du_x, Ir2_Mu_long_Du_y, & + !$OMP Ir2_Mu_long_Du_z, Ir2_Mu_gauss_Du, & + !$OMP Ir2_Mu_long_Du_2, int2_grad1_u2e_ao) + !$OMP DO SCHEDULE (static) + 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) + r2 = x*x + y*y + z*z + + dx = env_grad(1,ipoint) + dy = env_grad(2,ipoint) + dz = env_grad(3,ipoint) + + tmp0_x = 0.5d0 * (env_val(ipoint) * x + r2 * dx) + tmp0_y = 0.5d0 * (env_val(ipoint) * y + r2 * dy) + tmp0_z = 0.5d0 * (env_val(ipoint) * z + r2 * dz) + + tmp1 = 0.5d0 * env_val(ipoint) + + tmp1_x = tmp_ct * dx + tmp1_y = tmp_ct * dy + tmp1_z = tmp_ct * dz + + do j = 1, ao_num + do i = 1, ao_num + + tmp2 = 0.5d0 * Ir2_Mu_long_Du_2(i,j,ipoint) - x * Ir2_Mu_long_Du_x(i,j,ipoint) - y * Ir2_Mu_long_Du_y(i,j,ipoint) - z * Ir2_Mu_long_Du_z(i,j,ipoint) + + int2_grad1_u2e_ao(i,j,ipoint,1) = Ir2_Mu_long_Du_0(i,j,ipoint) * tmp0_x - tmp1 * Ir2_Mu_long_Du_x(i,j,ipoint) + dx * tmp2 - tmp1_x * Ir2_Mu_gauss_Du(i,j,ipoint) + int2_grad1_u2e_ao(i,j,ipoint,2) = Ir2_Mu_long_Du_0(i,j,ipoint) * tmp0_y - tmp1 * Ir2_Mu_long_Du_y(i,j,ipoint) + dy * tmp2 - tmp1_y * Ir2_Mu_gauss_Du(i,j,ipoint) + int2_grad1_u2e_ao(i,j,ipoint,3) = Ir2_Mu_long_Du_0(i,j,ipoint) * tmp0_z - tmp1 * Ir2_Mu_long_Du_z(i,j,ipoint) + dz * tmp2 - tmp1_z * Ir2_Mu_gauss_Du(i,j,ipoint) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + FREE Ir2_Mu_long_Du_0 Ir2_Mu_long_Du_x Ir2_Mu_long_Du_y Ir2_Mu_long_Du_z Ir2_Mu_long_Du_2 + FREE Ir2_Mu_gauss_Du + + else + + print *, ' Error in int2_grad1_u2e_ao: Unknown Jastrow' + stop + + endif ! j2e_type + + call wall_time(time1) + print*, ' wall time for int2_grad1_u2e_ao (min) =', (time1-time0)/60.d0 + call print_memory_usage() + +END_PROVIDER + +! --- + diff --git a/plugins/local/non_h_ints_mu/tc_integ.irp.f b/plugins/local/non_h_ints_mu/tc_integ.irp.f index 67ab4c89..2255cb5c 100644 --- a/plugins/local/non_h_ints_mu/tc_integ.irp.f +++ b/plugins/local/non_h_ints_mu/tc_integ.irp.f @@ -3,7 +3,7 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f BEGIN_DOC ! - ! int2_grad1_u12_ao(i,j,ipoint,:) = \int dr2 [-1 * \grad_r1 J(r1,r2)] \phi_i(r2) \phi_j(r2) + ! int2_grad1_u12_ao(i,j,ipoint,:) = \int dr2 [\grad_r1 J(r1,r2)] \phi_i(r2) \phi_j(r2) ! ! where r1 = r(ipoint) ! @@ -63,123 +63,12 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f int2_grad1_u12_ao = 0.d0 - !elseif((j2e_type .eq. "Mu") .and. (env_type .eq. "None")) then - - ! PROVIDE v_ij_erf_rk_cst_mu x_v_ij_erf_rk_cst_mu - - ! int2_grad1_u12_ao = 0.d0 - ! !$OMP PARALLEL & - ! !$OMP DEFAULT (NONE) & - ! !$OMP PRIVATE (ipoint, i, j, x, y, z, tmp1) & - ! !$OMP SHARED ( ao_num, n_points_final_grid, final_grid_points & - ! !$OMP , v_ij_erf_rk_cst_mu, x_v_ij_erf_rk_cst_mu, int2_grad1_u12_ao) - ! !$OMP DO SCHEDULE (static) - ! 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 - ! tmp1 = v_ij_erf_rk_cst_mu(i,j,ipoint) - ! int2_grad1_u12_ao(i,j,ipoint,1) = 0.5d0 * (tmp1 * x - x_v_ij_erf_rk_cst_mu(i,j,ipoint,1)) - ! int2_grad1_u12_ao(i,j,ipoint,2) = 0.5d0 * (tmp1 * y - x_v_ij_erf_rk_cst_mu(i,j,ipoint,2)) - ! int2_grad1_u12_ao(i,j,ipoint,3) = 0.5d0 * (tmp1 * z - x_v_ij_erf_rk_cst_mu(i,j,ipoint,3)) - ! enddo - ! enddo - ! enddo - ! !$OMP END DO - ! !$OMP END PARALLEL - - !elseif((j2e_type .eq. "Mu") .and. (env_type .eq. "Prod_Gauss")) then - - ! PROVIDE env_type env_val env_grad - ! PROVIDE v_ij_erf_rk_cst_mu_env v_ij_u_cst_mu_env_an x_v_ij_erf_rk_cst_mu_env - - ! int2_grad1_u12_ao = 0.d0 - ! !$OMP PARALLEL & - ! !$OMP DEFAULT (NONE) & - ! !$OMP PRIVATE (ipoint, i, j, x, y, z, tmp0, tmp1, tmp2, tmp0_x, tmp0_y, tmp0_z) & - ! !$OMP SHARED (ao_num, n_points_final_grid, final_grid_points, env_val, env_grad, & - ! !$OMP v_ij_erf_rk_cst_mu_env, v_ij_u_cst_mu_env_an, x_v_ij_erf_rk_cst_mu_env, int2_grad1_u12_ao) - ! !$OMP DO SCHEDULE (static) - ! 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) - ! tmp0 = 0.5d0 * env_val(ipoint) - ! tmp0_x = env_grad(1,ipoint) - ! tmp0_y = env_grad(2,ipoint) - ! tmp0_z = env_grad(3,ipoint) - ! do j = 1, ao_num - ! do i = 1, ao_num - ! tmp1 = tmp0 * v_ij_erf_rk_cst_mu_env(i,j,ipoint) - ! tmp2 = v_ij_u_cst_mu_env_an(i,j,ipoint) - ! int2_grad1_u12_ao(i,j,ipoint,1) = tmp1 * x - tmp0 * x_v_ij_erf_rk_cst_mu_env(i,j,ipoint,1) - tmp2 * tmp0_x - ! int2_grad1_u12_ao(i,j,ipoint,2) = tmp1 * y - tmp0 * x_v_ij_erf_rk_cst_mu_env(i,j,ipoint,2) - tmp2 * tmp0_y - ! int2_grad1_u12_ao(i,j,ipoint,3) = tmp1 * z - tmp0 * x_v_ij_erf_rk_cst_mu_env(i,j,ipoint,3) - tmp2 * tmp0_z - ! enddo - ! enddo - ! enddo - ! !$OMP END DO - ! !$OMP END PARALLEL - - !elseif((j2e_type .eq. "Mu") .and. (env_type .eq. "Sum_Gauss")) then - elseif( (j2e_type .eq. "Mu") .and. & ( (env_type .eq. "None") .or. (env_type .eq. "Prod_Gauss") .or. (env_type .eq. "Sum_Gauss") ) ) then - PROVIDE mu_erf - PROVIDE env_type env_val env_grad - PROVIDE Ir2_Mu_long_Du_0 Ir2_Mu_long_Du_x Ir2_Mu_long_Du_y Ir2_Mu_long_Du_z Ir2_Mu_long_Du_2 - PROVIDE Ir2_Mu_gauss_Du - - tmp_ct = 0.5d0 / (dsqrt(dacos(-1.d0)) * mu_erf) - - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, i, j, x, y, z, r2, dx, dy, dz, tmp1, tmp2, & - !$OMP tmp0_x, tmp0_y, tmp0_z, tmp1_x, tmp1_y, tmp1_z) & - !$OMP SHARED (ao_num, n_points_final_grid, final_grid_points, & - !$OMP tmp_ct, env_val, env_grad, Ir2_Mu_long_Du_0, & - !$OMP Ir2_Mu_long_Du_x, Ir2_Mu_long_Du_y, & - !$OMP Ir2_Mu_long_Du_z, Ir2_Mu_gauss_Du, & - !$OMP Ir2_Mu_long_Du_2, int2_grad1_u12_ao) - !$OMP DO SCHEDULE (static) - 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) - r2 = x*x + y*y + z*z - - dx = env_grad(1,ipoint) - dy = env_grad(2,ipoint) - dz = env_grad(3,ipoint) - - tmp0_x = 0.5d0 * (env_val(ipoint) * x + r2 * dx) - tmp0_y = 0.5d0 * (env_val(ipoint) * y + r2 * dy) - tmp0_z = 0.5d0 * (env_val(ipoint) * z + r2 * dz) - - tmp1 = 0.5d0 * env_val(ipoint) - - tmp1_x = tmp_ct * dx - tmp1_y = tmp_ct * dy - tmp1_z = tmp_ct * dz - - do j = 1, ao_num - do i = 1, ao_num - - tmp2 = 0.5d0 * Ir2_Mu_long_Du_2(i,j,ipoint) - x * Ir2_Mu_long_Du_x(i,j,ipoint) - y * Ir2_Mu_long_Du_y(i,j,ipoint) - z * Ir2_Mu_long_Du_z(i,j,ipoint) - - int2_grad1_u12_ao(i,j,ipoint,1) = -Ir2_Mu_long_Du_0(i,j,ipoint) * tmp0_x + tmp1 * Ir2_Mu_long_Du_x(i,j,ipoint) - dx * tmp2 + tmp1_x * Ir2_Mu_gauss_Du(i,j,ipoint) - int2_grad1_u12_ao(i,j,ipoint,2) = -Ir2_Mu_long_Du_0(i,j,ipoint) * tmp0_y + tmp1 * Ir2_Mu_long_Du_y(i,j,ipoint) - dy * tmp2 + tmp1_y * Ir2_Mu_gauss_Du(i,j,ipoint) - int2_grad1_u12_ao(i,j,ipoint,3) = -Ir2_Mu_long_Du_0(i,j,ipoint) * tmp0_z + tmp1 * Ir2_Mu_long_Du_z(i,j,ipoint) - dz * tmp2 + tmp1_z * Ir2_Mu_gauss_Du(i,j,ipoint) - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - + PROVIDE int2_grad1_u2e_ao + int2_grad1_u12_ao = int2_grad1_u2e_ao + else print *, ' Error in int2_grad1_u12_ao: Unknown Jastrow' @@ -195,20 +84,13 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f PROVIDE ao_overlap PROVIDE j1e_gradx j1e_grady j1e_gradz - double precision, allocatable :: int_tmp(:,:,:,:) - - ! minus because we calculate \int [-\grad_1 u(1,2)] - tmp_ct = -1.d0 / (dble(elec_num) - 1.d0) + tmp_ct = 1.d0 / (dble(elec_num) - 1.d0) !$OMP PARALLEL & !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, i, j, tmp0_x, tmp0_y, tmp0_z, int_tmp) & + !$OMP PRIVATE (ipoint, i, j, tmp0_x, tmp0_y, tmp0_z) & !$OMP SHARED (ao_num, n_points_final_grid, tmp_ct, ao_overlap, & !$OMP j1e_gradx, j1e_grady, j1e_gradz, int2_grad1_u12_ao) - - allocate(int_tmp(ao_num,ao_num,n_points_final_grid,3)) - int_tmp = 0.d0 - !$OMP DO do ipoint = 1, n_points_final_grid tmp0_x = tmp_ct * j1e_gradx(ipoint) @@ -216,34 +98,15 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f tmp0_z = tmp_ct * j1e_gradz(ipoint) do j = 1, ao_num do i = 1, ao_num - int_tmp(i,j,ipoint,1) = int_tmp(i,j,ipoint,1) + tmp0_x * ao_overlap(i,j) - int_tmp(i,j,ipoint,2) = int_tmp(i,j,ipoint,2) + tmp0_y * ao_overlap(i,j) - int_tmp(i,j,ipoint,3) = int_tmp(i,j,ipoint,3) + tmp0_z * ao_overlap(i,j) + int2_grad1_u12_ao(i,j,ipoint,1) = int2_grad1_u12_ao(i,j,ipoint,1) + tmp0_x * ao_overlap(i,j) + int2_grad1_u12_ao(i,j,ipoint,2) = int2_grad1_u12_ao(i,j,ipoint,2) + tmp0_y * ao_overlap(i,j) + int2_grad1_u12_ao(i,j,ipoint,3) = int2_grad1_u12_ao(i,j,ipoint,3) + tmp0_z * ao_overlap(i,j) enddo enddo enddo - !$OMP END DO NOWAIT - - !$OMP CRITICAL - int2_grad1_u12_ao = int2_grad1_u12_ao + int_tmp - !$OMP END CRITICAL - - deallocate(int_tmp) + !$OMP END DO !$OMP END PARALLEL - else - - !if((j2e_type .eq. "Mu") .and. (env_type .eq. "None")) then - ! FREE v_ij_erf_rk_cst_mu x_v_ij_erf_rk_cst_mu - !elseif((j2e_type .eq. "Mu") .and. (env_type .eq. "Prod_Gauss")) then - ! FREE v_ij_erf_rk_cst_mu_env v_ij_u_cst_mu_env_an x_v_ij_erf_rk_cst_mu_env - !elseif((j2e_type .eq. "Mu") .and. (env_type .eq. "Sum_Gauss")) then - - if( (j2e_type .eq. "Mu") .and. & - ( (env_type .eq. "None") .or. (env_type .eq. "Prod_Gauss") .or. (env_type .eq. "Sum_Gauss") ) ) then - FREE Ir2_Mu_long_Du_0 Ir2_Mu_long_Du_x Ir2_Mu_long_Du_y Ir2_Mu_long_Du_z Ir2_Mu_gauss_Du Ir2_Mu_long_Du_2 - endif - endif ! j1e_type ! --- @@ -532,7 +395,7 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p else - print *, ' Error in int2_grad1_u12_square_ao: Unknown Jhastrow' + print *, ' Error in int2_grad1_u12_square_ao: Unknown Jastrow' stop endif ! j2e_type @@ -544,75 +407,46 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p PROVIDE elec_num PROVIDE ao_overlap PROVIDE j1e_gradx j1e_grady j1e_gradz + PROVIDE int2_grad1_u2e_ao - double precision, allocatable :: int_tmp(:,:,:) - - tmp_ct1 = 1.d0 / (dsqrt(dacos(-1.d0)) * mu_erf) - tmp_ct2 = 1.d0 / (dble(elec_num) - 1.d0) - - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, i, j, x, y, z, r2, dx1, dy1, dz1, & - !$OMP dx2, dy2, dz2, dr12, tmp0, tmp1, tmp2, tmp3, tmp4, & - !$OMP tmp0_x, tmp0_y, tmp0_z, int_tmp) & - !$OMP SHARED (ao_num, n_points_final_grid, final_grid_points, & - !$OMP tmp_ct1, tmp_ct2, env_val, env_grad, & - !$OMP j1e_gradx, j1e_grady, j1e_gradz, ao_overlap, & - !$OMP Ir2_Mu_long_Du_0, Ir2_Mu_long_Du_2, & - !$OMP Ir2_Mu_long_Du_x, Ir2_Mu_long_Du_y, & - !$OMP Ir2_Mu_long_Du_z, Ir2_Mu_gauss_Du, int2_grad1_u12_square_ao) - - allocate(int_tmp(ao_num,ao_num,n_points_final_grid)) - int_tmp = 0.d0 + tmp_ct1 = 2.d0 / (dble(elec_num) - 1.d0) + tmp_ct2 = 1.d0 / ((dble(elec_num) - 1.d0) * (dble(elec_num) - 1.d0)) + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, dx, dy, dz, r2, & + !$OMP tmp0, tmp0_x, tmp0_y, tmp0_z) & + !$OMP SHARED (ao_num, n_points_final_grid, & + !$OMP tmp_ct1, tmp_ct2, ao_overlap, & + !$OMP j1e_gradx, j1e_grady, j1e_gradz, & + !$OMP int2_grad1_u2e_ao, int2_grad1_u12_square_ao) !$OMP DO 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) - r2 = x*x + y*y + z*z + dx = j1e_gradx(ipoint) + dy = j1e_grady(ipoint) + dz = j1e_gradz(ipoint) + r2 = dx*dx + dy*dy + dz*dz - dx1 = env_grad(1,ipoint) - dy1 = env_grad(2,ipoint) - dz1 = env_grad(3,ipoint) - - dx2 = j1e_gradx(ipoint) - dy2 = j1e_grady(ipoint) - dz2 = j1e_gradz(ipoint) - - dr12 = dx1*dx2 + dy1*dy2 + dz1*dz2 - - tmp0 = tmp_ct2 * (env_val(ipoint) * (dx2*x + dy2*y + dz2*z) + r2*dr12) - tmp1 = tmp_ct2 * dr12 - tmp2 = tmp_ct1 * tmp_ct2 * dr12 - tmp3 = tmp_ct2 * tmp_ct2 * (dx2*dx2 + dy2*dy2 + dz2*dz2) - - tmp0_x = tmp_ct2 * (env_val(ipoint) * dx2 + 2.d0 * dr12 * x) - tmp0_y = tmp_ct2 * (env_val(ipoint) * dy2 + 2.d0 * dr12 * y) - tmp0_z = tmp_ct2 * (env_val(ipoint) * dz2 + 2.d0 * dr12 * z) + tmp0 = tmp_ct2 * r2 + tmp0_x = tmp_ct1 * dx + tmp0_y = tmp_ct1 * dy + tmp0_z = tmp_ct1 * dz do j = 1, ao_num do i = 1, ao_num - - tmp4 = tmp0_x * Ir2_Mu_long_Du_x(i,j,ipoint) + tmp0_y * Ir2_Mu_long_Du_y(i,j,ipoint) + tmp0_z * Ir2_Mu_long_Du_z(i,j,ipoint) - int_tmp(i,j,ipoint) = int_tmp(i,j,ipoint) + tmp0 * Ir2_Mu_long_Du_0(i,j,ipoint) - tmp4 + tmp1 * Ir2_Mu_long_Du_2(i,j,ipoint) & - - tmp2 * Ir2_Mu_gauss_Du(i,j,ipoint) + tmp3 * ao_overlap(i,j) + int2_grad1_u12_square_ao(i,j,ipoint) = int2_grad1_u12_square_ao(i,j,ipoint) & + + tmp0 * ao_overlap(i,j) & + + tmp0_x * int2_grad1_u2e_ao(i,j,ipoint,1) & + + tmp0_y * int2_grad1_u2e_ao(i,j,ipoint,2) & + + tmp0_z * int2_grad1_u2e_ao(i,j,ipoint,3) enddo enddo enddo - !$OMP END DO NOWAIT - - !$OMP CRITICAL - int2_grad1_u12_square_ao = int2_grad1_u12_square_ao + int_tmp - !$OMP END CRITICAL - - deallocate(int_tmp) + !$OMP END DO !$OMP END PARALLEL - FREE Ir2_Mu_long_Du_0 Ir2_Mu_long_Du_x Ir2_Mu_long_Du_y Ir2_Mu_long_Du_z Ir2_Mu_gauss_Du Ir2_Mu_long_Du_2 - endif ! j1e_type ! --- diff --git a/plugins/local/non_h_ints_mu/test_non_h_ints.irp.f b/plugins/local/non_h_ints_mu/test_non_h_ints.irp.f index e349d412..3f88c53f 100644 --- a/plugins/local/non_h_ints_mu/test_non_h_ints.irp.f +++ b/plugins/local/non_h_ints_mu/test_non_h_ints.irp.f @@ -755,7 +755,7 @@ subroutine test_j1e_fit_ao() allocate(x(n_points_final_grid), y(n_points_final_grid), z(n_points_final_grid)) - g = 0.5d0 * (dble(elec_num) - 1.d0) / dble(elec_num) + g = -0.5d0 * (dble(elec_num) - 1.d0) / dble(elec_num) call dgemv("T", ao_num*ao_num, n_points_final_grid, g, int2_grad1_u2e_ao(1,1,1,1), ao_num*ao_num, Pt, 1, 0.d0, x, 1) call dgemv("T", ao_num*ao_num, n_points_final_grid, g, int2_grad1_u2e_ao(1,1,1,2), ao_num*ao_num, Pt, 1, 0.d0, y, 1) diff --git a/plugins/local/non_h_ints_mu/total_tc_int.irp.f b/plugins/local/non_h_ints_mu/total_tc_int.irp.f index 59f5174b..4cedf0e6 100644 --- a/plugins/local/non_h_ints_mu/total_tc_int.irp.f +++ b/plugins/local/non_h_ints_mu/total_tc_int.irp.f @@ -167,12 +167,15 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n !$OMP END PARALLEL do m = 1, 3 - call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 & + call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, -1.d0 & , int2_grad1_u12_ao(1,1,1,m), ao_num*ao_num, b_mat(1,1,1,m), n_points_final_grid & , 1.d0, ao_two_e_tc_tot, ao_num*ao_num) enddo deallocate(b_mat) + FREE int2_grad1_u12_ao + FREE int2_grad1_u2e_ao + endif ! var_tc ! ---