diff --git a/plugins/local/jastrow/EZFIO.cfg b/plugins/local/jastrow/EZFIO.cfg index 2eac6aa2..8728916d 100644 --- a/plugins/local/jastrow/EZFIO.cfg +++ b/plugins/local/jastrow/EZFIO.cfg @@ -89,6 +89,12 @@ doc: linear coef of functions in 1e-Jastrow interface: ezfio size: (jastrow.j1e_size,nuclei.nucl_num) +[j1e_coef_ao] +type: double precision +doc: coefficients of the 1-body Jastrow in AOs +interface: ezfio +size: (nuclei.nucl_num) + [j1e_expo] type: double precision doc: exponenets of functions in 1e-Jastrow diff --git a/plugins/local/jastrow/README.md b/plugins/local/jastrow/README.md index f7ea8e02..22486edd 100644 --- a/plugins/local/jastrow/README.md +++ b/plugins/local/jastrow/README.md @@ -9,7 +9,7 @@ The main keywords are: ## j2e_type Options -1. **none:** No 2e-Jastrow is used. +1. **None:** No 2e-Jastrow is used. 2. **Mu:** 2e-Jastrow inspired by Range Separated Density Functional Theory. It has the following shape:

@@ -28,7 +28,7 @@ The 2-electron Jastrow is multiplied by an envelope \(v\):

-- if `env_type` is **none**: No envelope is used. +- if `env_type` is **None**: No envelope is used. - if `env_type` is **Prod_Gauss**:

@@ -50,7 +50,7 @@ The 1-electron Jastrow used is:

-- if `j1e_type` is **none**: No one-electron Jastrow is used. +- if `j1e_type` is **None**: No one-electron Jastrow is used. - if `j1e_type` is **Gauss**: We use

@@ -65,4 +65,5 @@ are defined by the tables `j1e_coef` and `j1e_expo`, respectively.

+- if `j1e_type` is **Charge_Harmonizer_AO**: The one-electron Jastrow factor **Charge_Harmonizer** is fitted by the atomic orbitals diff --git a/plugins/local/ao_many_one_e_ints/listj1b.irp.f b/plugins/local/jastrow/listj1b.irp.f similarity index 92% rename from plugins/local/ao_many_one_e_ints/listj1b.irp.f rename to plugins/local/jastrow/listj1b.irp.f index 2b049943..49954d47 100644 --- a/plugins/local/ao_many_one_e_ints/listj1b.irp.f +++ b/plugins/local/jastrow/listj1b.irp.f @@ -7,7 +7,11 @@ BEGIN_PROVIDER [integer, List_env1s_size] PROVIDE env_type - if(env_type .eq. "Prod_Gauss") then + if(env_type .eq. "None") then + + List_env1s_size = 1 + + elseif(env_type .eq. "Prod_Gauss") then List_env1s_size = 2**nucl_num @@ -63,11 +67,17 @@ END_PROVIDER provide env_type env_expo env_coef - List_env1s_coef = 0.d0 - List_env1s_expo = 0.d0 - List_env1s_cent = 0.d0 + if(env_type .eq. "None") then - if(env_type .eq. "Prod_Gauss") then + List_env1s_coef( 1) = 1.d0 + List_env1s_expo( 1) = 0.d0 + List_env1s_cent(1:3,1) = 0.d0 + + elseif(env_type .eq. "Prod_Gauss") then + + List_env1s_coef = 0.d0 + List_env1s_expo = 0.d0 + List_env1s_cent = 0.d0 do i = 1, List_env1s_size @@ -150,7 +160,11 @@ BEGIN_PROVIDER [integer, List_env1s_square_size] implicit none double precision :: tmp - if(env_type .eq. "Prod_Gauss") then + if(env_type .eq. "None") then + + List_env1s_square_size = 1 + + elseif(env_type .eq. "Prod_Gauss") then List_env1s_square_size = 3**nucl_num @@ -220,11 +234,17 @@ END_PROVIDER provide env_type env_expo env_coef - List_env1s_square_coef = 0.d0 - List_env1s_square_expo = 0.d0 - List_env1s_square_cent = 0.d0 + if(env_type .eq. "None") then - if(env_type .eq. "Prod_Gauss") then + List_env1s_square_coef( 1) = 1.d0 + List_env1s_square_expo( 1) = 0.d0 + List_env1s_square_cent(1:3,1) = 0.d0 + + elseif(env_type .eq. "Prod_Gauss") then + + List_env1s_square_coef = 0.d0 + List_env1s_square_expo = 0.d0 + List_env1s_square_cent = 0.d0 do i = 1, List_env1s_square_size 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 96275887..c8da0680 100644 --- a/plugins/local/non_h_ints_mu/jast_1e.irp.f +++ b/plugins/local/non_h_ints_mu/jast_1e.irp.f @@ -75,6 +75,7 @@ END_PROVIDER double precision :: a, c, g, tmp_x, tmp_y, tmp_z double precision :: time0, time1 double precision, allocatable :: Pa(:,:), Pb(:,:), Pt(:,:) + double precision, allocatable :: coef_fit(:) PROVIDE j1e_type @@ -133,7 +134,7 @@ END_PROVIDER PROVIDE elec_alpha_num elec_beta_num elec_num PROVIDE mo_coef - PROVIDE int2_grad1_u2b_ao + PROVIDE int2_grad1_u2e_ao allocate(Pa(ao_num,ao_num), Pb(ao_num,ao_num), Pt(ao_num,ao_num)) @@ -152,12 +153,59 @@ END_PROVIDER 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_u2b_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_u2b_ao(1,1,1,2), ao_num*ao_num, Pt, 1, 0.d0, j1e_grady, 1) - call dgemv("T", ao_num*ao_num, n_points_final_grid, g, int2_grad1_u2b_ao(1,1,1,3), ao_num*ao_num, Pt, 1, 0.d0, j1e_gradz, 1) + 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) + call dgemv("T", ao_num*ao_num, n_points_final_grid, g, int2_grad1_u2e_ao(1,1,1,3), ao_num*ao_num, Pt, 1, 0.d0, j1e_gradz, 1) + + FREE int2_grad1_u2e_ao deallocate(Pa, Pb, Pt) + elseif(j1e_type .eq. "Charge_Harmonizer_AO") then + + ! \grad_1 \sum_{\eta} C_{\eta} \chi_{\eta} + ! where + ! \chi_{\eta} are the AOs + ! C_{\eta} are fitted to mimic (j1e_type .eq. "Charge_Harmonizer") + ! + ! The - sign is in the parameters C_{\eta} + + PROVIDE aos_grad_in_r_array + + allocate(coef_fit(ao_num)) + + call get_j1e_coef_fit_ao(ao_num, coef_fit) + call ezfio_set_jastrow_j1e_coef_ao(coef_fit) + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i, ipoint, tmp_x, tmp_y, tmp_z, & + !$OMP c) & + !$OMP SHARED (n_points_final_grid, ao_num, & + !$OMP aos_grad_in_r_array, coef_fit, & + !$OMP j1e_gradx, j1e_grady, j1e_gradz) + !$OMP DO SCHEDULE (static) + do ipoint = 1, n_points_final_grid + + tmp_x = 0.d0 + tmp_y = 0.d0 + tmp_z = 0.d0 + do i = 1, ao_num + c = coef_fit(i) + tmp_x = tmp_x + c * aos_grad_in_r_array(i,ipoint,1) + tmp_y = tmp_y + c * aos_grad_in_r_array(i,ipoint,2) + tmp_z = tmp_z + c * aos_grad_in_r_array(i,ipoint,3) + enddo + + j1e_gradx(ipoint) = tmp_x + j1e_grady(ipoint) = tmp_y + j1e_gradz(ipoint) = tmp_z + enddo + !$OMP END DO + !$OMP END PARALLEL + + deallocate(coef_fit) + else print *, ' Error in j1e_grad: Unknown j1e_type = ', j1e_type 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 1e95f80a..defe8897 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,11 +1,106 @@ ! --- -BEGIN_PROVIDER [double precision, int2_grad1_u2b_ao, (ao_num, ao_num, n_points_final_grid, 3)] +BEGIN_PROVIDER [double precision, int2_u2e_ao, (ao_num, ao_num, n_points_final_grid)] BEGIN_DOC ! - ! int2_grad1_u2b_ao(i,j,ipoint,:) = \int dr2 [-1 * \grad_r1 J_2b(r1,r2)] \phi_i(r2) \phi_j(r2) + ! 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) + + tmp1 = 0.5d0 * env_val(ipoint) + tmp0 = tmp1 * r2 + 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) ! @@ -22,35 +117,23 @@ BEGIN_PROVIDER [double precision, int2_grad1_u2b_ao, (ao_num, ao_num, n_points_f double precision :: tmp1_x, tmp1_y, tmp1_z PROVIDE j2e_type + PROVIDE Env_type call wall_time(time0) + print*, ' providing int2_grad1_u2e_ao ...' - print*, ' providing int2_grad1_u2b_ao ...' - - if(tc_integ_type .eq. "numeric") then - - ! TODO combine 1shot & int2_grad1_u12_ao_num - - PROVIDE int2_grad1_u12_ao_num - int2_grad1_u2b_ao = int2_grad1_u12_ao_num - - !PROVIDE int2_grad1_u12_ao_num_1shot - !int2_grad1_u2b_ao = int2_grad1_u12_ao_num_1shot - - elseif(tc_integ_type .eq. "semi-analytic") then - - ! --- + if(tc_integ_type .eq. "semi-analytic") then if((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_u2b_ao = 0.d0 + int2_grad1_u2e_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_u2b_ao) + !$OMP , v_ij_erf_rk_cst_mu, x_v_ij_erf_rk_cst_mu, int2_grad1_u2e_ao) !$OMP DO SCHEDULE (static) do ipoint = 1, n_points_final_grid x = final_grid_points(1,ipoint) @@ -59,9 +142,9 @@ BEGIN_PROVIDER [double precision, int2_grad1_u2b_ao, (ao_num, ao_num, n_points_f do j = 1, ao_num do i = 1, ao_num tmp1 = v_ij_erf_rk_cst_mu(i,j,ipoint) - int2_grad1_u2b_ao(i,j,ipoint,1) = 0.5d0 * (tmp1 * x - x_v_ij_erf_rk_cst_mu(i,j,ipoint,1)) - int2_grad1_u2b_ao(i,j,ipoint,2) = 0.5d0 * (tmp1 * y - x_v_ij_erf_rk_cst_mu(i,j,ipoint,2)) - int2_grad1_u2b_ao(i,j,ipoint,3) = 0.5d0 * (tmp1 * z - x_v_ij_erf_rk_cst_mu(i,j,ipoint,3)) + int2_grad1_u2e_ao(i,j,ipoint,1) = 0.5d0 * (tmp1 * x - x_v_ij_erf_rk_cst_mu(i,j,ipoint,1)) + int2_grad1_u2e_ao(i,j,ipoint,2) = 0.5d0 * (tmp1 * y - x_v_ij_erf_rk_cst_mu(i,j,ipoint,2)) + int2_grad1_u2e_ao(i,j,ipoint,3) = 0.5d0 * (tmp1 * z - x_v_ij_erf_rk_cst_mu(i,j,ipoint,3)) enddo enddo enddo @@ -73,12 +156,12 @@ BEGIN_PROVIDER [double precision, int2_grad1_u2b_ao, (ao_num, ao_num, n_points_f 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_u2b_ao = 0.d0 + int2_grad1_u2e_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_u2b_ao) + !$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_u2e_ao) !$OMP DO SCHEDULE (static) do ipoint = 1, n_points_final_grid x = final_grid_points(1,ipoint) @@ -92,9 +175,9 @@ BEGIN_PROVIDER [double precision, int2_grad1_u2b_ao, (ao_num, ao_num, n_points_f 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_u2b_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_u2b_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_u2b_ao(i,j,ipoint,3) = tmp1 * z - tmp0 * x_v_ij_erf_rk_cst_mu_env(i,j,ipoint,3) - tmp2 * tmp0_z + int2_grad1_u2e_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_u2e_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_u2e_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 @@ -110,7 +193,7 @@ BEGIN_PROVIDER [double precision, int2_grad1_u2b_ao, (ao_num, ao_num, n_points_f tmp_ct = 0.5d0 / (dsqrt(dacos(-1.d0)) * mu_erf) - int2_grad1_u2b_ao = 0.d0 + int2_grad1_u2e_ao = 0.d0 !$OMP PARALLEL & !$OMP DEFAULT (NONE) & @@ -120,7 +203,7 @@ BEGIN_PROVIDER [double precision, int2_grad1_u2b_ao, (ao_num, ao_num, n_points_f !$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_u2b_ao) + !$OMP Ir2_Mu_long_Du_2, int2_grad1_u2e_ao) !$OMP DO SCHEDULE (static) do ipoint = 1, n_points_final_grid @@ -148,9 +231,9 @@ BEGIN_PROVIDER [double precision, int2_grad1_u2b_ao, (ao_num, ao_num, n_points_f 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_u2b_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_u2b_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_u2b_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) + 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 @@ -159,7 +242,7 @@ BEGIN_PROVIDER [double precision, int2_grad1_u2b_ao, (ao_num, ao_num, n_points_f else - print *, ' Error in int2_grad1_u2b_ao: Unknown Jastrow' + print *, ' Error in int2_grad1_u2e_ao: Unknown Jastrow' stop endif ! j2e_type @@ -172,10 +255,97 @@ BEGIN_PROVIDER [double precision, int2_grad1_u2b_ao, (ao_num, ao_num, n_points_f endif ! tc_integ_type call wall_time(time1) - print*, ' wall time for int2_grad1_u2b_ao (min) =', (time1-time0)/60.d0 + 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 + integer , intent(in) :: dim_fit + double precision, intent(out) :: coef_fit(dim_fit) + + integer :: i, ipoint + double precision :: g + double precision, allocatable :: A(:,:), b(:), A_inv(:,:) + double precision, allocatable :: Pa(:,:), Pb(:,:), Pt(:,:) + double precision, allocatable :: u1e_tmp(:) + + PROVIDE j1e_type + PROVIDE int2_u2e_ao + PROVIDE elec_alpha_num elec_beta_num elec_num + PROVIDE mo_coef + PROVIDE ao_overlap + + ! --- --- --- + ! get u1e(r) + + allocate(Pa(ao_num,ao_num), Pb(ao_num,ao_num), Pt(ao_num,ao_num)) + + call dgemm( 'N', 'T', ao_num, ao_num, elec_alpha_num, 1.d0 & + , mo_coef, size(mo_coef, 1), mo_coef, size(mo_coef, 1) & + , 0.d0, Pa, size(Pa, 1)) + + if(elec_alpha_num .eq. elec_beta_num) then + Pb = Pa + else + call dgemm( 'N', 'T', ao_num, ao_num, elec_beta_num, 1.d0 & + , mo_coef, size(mo_coef, 1), mo_coef, size(mo_coef, 1) & + , 0.d0, Pb, size(Pb, 1)) + endif + Pt = Pa + Pb + + allocate(u1e_tmp(n_points_final_grid)) + + g = 0.5d0 * (dble(elec_num) - 1.d0) / dble(elec_num) + call dgemv("T", ao_num*ao_num, n_points_final_grid, g, int2_u2e_ao(1,1,1), ao_num*ao_num, Pt, 1, 0.d0, u1e_tmp, 1) + + FREE int2_u2e_ao + + deallocate(Pa, Pb, Pt) + + ! --- --- --- + ! get A & b + + allocate(A(ao_num,ao_num), b(ao_num)) + + A(1:ao_num,1:ao_num) = ao_overlap(1:ao_num,1:ao_num) + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i, ipoint) & + !$OMP SHARED (n_points_final_grid, ao_num, & + !$OMP final_weight_at_r_vector, aos_in_r_array_transp, u1e_tmp, b) + !$OMP DO SCHEDULE (static) + do i = 1, ao_num + b(i) = 0.d0 + do ipoint = 1, n_points_final_grid + b(i) = b(i) + final_weight_at_r_vector(ipoint) * aos_in_r_array_transp(ipoint,i) * u1e_tmp(ipoint) + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + deallocate(u1e_tmp) + + ! --- --- --- + ! solve Ax = b + + allocate(A_inv(ao_num,ao_num)) + call get_inverse(A, ao_num, ao_num, A_inv, ao_num) + deallocate(A) + + ! coef_fit = A_inv x b + call dgemv("N", ao_num, ao_num, 1.d0, A_inv, ao_num, b, 1, 0.d0, coef_fit, 1) + deallocate(A_inv, b) + + return +end + +! --- + + + 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 88336485..ed0f8f6b 100644 --- a/plugins/local/non_h_ints_mu/tc_integ.irp.f +++ b/plugins/local/non_h_ints_mu/tc_integ.irp.f @@ -63,67 +63,70 @@ 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 + !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 + ! 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 + ! 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 + !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 + ! 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 + ! 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. "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 @@ -132,8 +135,6 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f tmp_ct = 0.5d0 / (dsqrt(dacos(-1.d0)) * mu_erf) - int2_grad1_u12_ao = 0.d0 - !$OMP PARALLEL & !$OMP DEFAULT (NONE) & !$OMP PRIVATE (ipoint, i, j, x, y, z, r2, dx, dy, dz, tmp1, tmp2, & @@ -220,11 +221,14 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f 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")) 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 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 6a30d909..4ace5d1c 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 @@ -623,7 +623,7 @@ subroutine test_j1e_grad() double precision, allocatable :: pa(:,:), Pb(:,:), Pt(:,:) double precision, allocatable :: x(:), y(:), z(:) - PROVIDE int2_grad1_u2b_ao + PROVIDE int2_grad1_u2e_ao PROVIDE mo_coef allocate(Pa(ao_num,ao_num), Pb(ao_num,ao_num), Pt(ao_num,ao_num)) @@ -652,9 +652,9 @@ subroutine test_j1e_grad() z(ipoint) = 0.d0 do i = 1, ao_num do j = 1, ao_num - x(ipoint) = x(ipoint) + g * Pt(i,j) * int2_grad1_u2b_ao(i,j,ipoint,1) - y(ipoint) = y(ipoint) + g * Pt(i,j) * int2_grad1_u2b_ao(i,j,ipoint,2) - z(ipoint) = z(ipoint) + g * Pt(i,j) * int2_grad1_u2b_ao(i,j,ipoint,3) + x(ipoint) = x(ipoint) + g * Pt(i,j) * int2_grad1_u2e_ao(i,j,ipoint,1) + y(ipoint) = y(ipoint) + g * Pt(i,j) * int2_grad1_u2e_ao(i,j,ipoint,2) + z(ipoint) = z(ipoint) + g * Pt(i,j) * int2_grad1_u2e_ao(i,j,ipoint,3) enddo enddo enddo diff --git a/plugins/local/non_hermit_dav/biorthog.irp.f b/plugins/local/non_hermit_dav/biorthog.irp.f index ab12150f..2229e17d 100644 --- a/plugins/local/non_hermit_dav/biorthog.irp.f +++ b/plugins/local/non_hermit_dav/biorthog.irp.f @@ -142,7 +142,7 @@ subroutine non_hrmt_diag_split_degen(n, A, leigvec, reigvec, n_real_eigv, eigval enddo enddo -end subroutine non_hrmt_diag_split_degen +end ! --- @@ -248,7 +248,7 @@ subroutine non_hrmt_real_diag_new(n, A, leigvec, reigvec, n_real_eigv, eigval) print*,'Your matrix intrinsically contains complex eigenvalues' endif -end subroutine non_hrmt_real_diag_new +end ! --- @@ -519,7 +519,7 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei return -end subroutine non_hrmt_bieig +end ! --- @@ -692,7 +692,7 @@ subroutine non_hrmt_bieig_random_diag(n, A, leigvec, reigvec, n_real_eigv, eigva return -end subroutine non_hrmt_bieig_random_diag +end ! --- @@ -801,7 +801,7 @@ subroutine non_hrmt_real_im(n, A, leigvec, reigvec, n_real_eigv, eigval) deallocate( S ) -end subroutine non_hrmt_real_im +end ! --- @@ -906,7 +906,7 @@ subroutine non_hrmt_generalized_real_im(n, A, B, leigvec, reigvec, n_real_eigv, deallocate( S ) -end subroutine non_hrmt_generalized_real_im +end ! --- @@ -1042,7 +1042,7 @@ subroutine non_hrmt_bieig_fullvect(n, A, leigvec, reigvec, n_real_eigv, eigval) return -end subroutine non_hrmt_bieig_fullvect +end ! --- diff --git a/plugins/local/non_hermit_dav/lapack_diag_non_hermit.irp.f b/plugins/local/non_hermit_dav/lapack_diag_non_hermit.irp.f index 4d51b79e..cb38347e 100644 --- a/plugins/local/non_hermit_dav/lapack_diag_non_hermit.irp.f +++ b/plugins/local/non_hermit_dav/lapack_diag_non_hermit.irp.f @@ -54,7 +54,7 @@ subroutine lapack_diag_non_sym(n, A, WR, WI, VL, VR) deallocate(Atmp, WORK) -end subroutine lapack_diag_non_sym +end subroutine non_sym_diag_inv_right(n,A,leigvec,reigvec,n_real_eigv,eigval) @@ -269,7 +269,7 @@ subroutine lapack_diag_non_sym_new(n, A, WR, WI, VL, VR) deallocate( Atmp ) deallocate( WORK, SCALE_array, RCONDE, RCONDV, IWORK ) -end subroutine lapack_diag_non_sym_new +end ! --- @@ -323,7 +323,7 @@ subroutine lapack_diag_non_sym_right(n, A, WR, WI, VR) ! write(*, '(1000(F16.10,X))') VR(:,i) ! enddo -end subroutine lapack_diag_non_sym_right +end ! --- @@ -437,7 +437,7 @@ subroutine non_hrmt_real_diag(n, A, leigvec, reigvec, n_real_eigv, eigval) print*, ' Notice that if you are interested in ground state it is not a problem :)' endif -end subroutine non_hrmt_real_diag +end ! --- @@ -495,7 +495,7 @@ subroutine lapack_diag_general_non_sym(n, A, B, WR, WI, VL, VR) deallocate( WORK, Atmp ) -end subroutine lapack_diag_general_non_sym +end ! --- @@ -570,7 +570,7 @@ subroutine non_hrmt_general_real_diag(n, A, B, reigvec, leigvec, n_real_eigv, ei enddo enddo -end subroutine non_hrmt_general_real_diag +end ! --- @@ -727,7 +727,7 @@ subroutine impose_biorthog_qr(m, n, thr_d, thr_nd, Vl, Vr) deallocate(tmp) return -end subroutine impose_biorthog_qr +end ! --- @@ -890,7 +890,7 @@ subroutine impose_biorthog_lu(m, n, Vl, Vr, S) !stop return -end subroutine impose_biorthog_lu +end ! --- @@ -996,7 +996,7 @@ subroutine check_EIGVEC(n, m, A, eigval, leigvec, reigvec, thr_diag, thr_norm, s deallocate( Mtmp ) -end subroutine check_EIGVEC +end ! --- @@ -1066,7 +1066,7 @@ subroutine check_degen(n, m, eigval, leigvec, reigvec) stop endif -end subroutine check_degen +end ! --- @@ -1169,7 +1169,7 @@ subroutine impose_weighted_orthog_svd(n, m, W, C) ! --- -end subroutine impose_weighted_orthog_svd +end ! --- @@ -1266,7 +1266,7 @@ subroutine impose_orthog_svd(n, m, C) ! --- -end subroutine impose_orthog_svd +end ! --- @@ -1365,7 +1365,7 @@ subroutine impose_orthog_svd_overlap(n, m, C, overlap) !enddo deallocate(S) -end subroutine impose_orthog_svd_overlap +end ! --- @@ -1442,7 +1442,7 @@ subroutine impose_orthog_GramSchmidt(n, m, C) ! --- -end subroutine impose_orthog_GramSchmidt +end ! --- @@ -1484,7 +1484,7 @@ subroutine impose_orthog_ones(n, deg_num, C) endif enddo -end subroutine impose_orthog_ones +end ! --- @@ -1577,7 +1577,7 @@ subroutine impose_orthog_degen_eigvec(n, e0, C0) endif enddo -end subroutine impose_orthog_degen_eigvec +end ! --- @@ -1661,7 +1661,7 @@ subroutine get_halfinv_svd(n, S) deallocate(S0, Stmp, Stmp2) -end subroutine get_halfinv_svd +end ! --- @@ -1776,7 +1776,7 @@ subroutine check_biorthog_binormalize(n, m, Vl, Vr, thr_d, thr_nd, stop_ifnot) stop endif -end subroutine check_biorthog_binormalize +end ! --- @@ -1840,7 +1840,7 @@ subroutine check_weighted_biorthog(n, m, W, Vl, Vr, thr_d, thr_nd, accu_d, accu_ stop endif -end subroutine check_weighted_biorthog +end ! --- @@ -1907,7 +1907,7 @@ subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, thr_d, thr_nd, stop_ stop endif -end subroutine check_biorthog +end ! --- @@ -1949,7 +1949,7 @@ subroutine check_orthog(n, m, V, accu_d, accu_nd, S) !print*, ' diag acc: ', accu_d !print*, ' nondiag acc: ', accu_nd -end subroutine check_orthog +end ! --- @@ -2067,7 +2067,7 @@ subroutine reorder_degen_eigvec(n, deg_num, e0, L0, R0) ! endif ! enddo ! -end subroutine reorder_degen_eigvec +end ! --- @@ -2188,7 +2188,7 @@ subroutine impose_biorthog_degen_eigvec(n, deg_num, e0, L0, R0) endif enddo -end subroutine impose_biorthog_degen_eigvec +end ! --- @@ -2282,7 +2282,7 @@ subroutine impose_orthog_biorthog_degen_eigvec(n, thr_d, thr_nd, e0, L0, R0) endif enddo -end subroutine impose_orthog_biorthog_degen_eigvec +end ! --- @@ -2420,7 +2420,7 @@ subroutine impose_unique_biorthog_degen_eigvec(n, thr_d, thr_nd, e0, C0, W0, L0, endif enddo -end subroutine impose_unique_biorthog_degen_eigvec +end ! --- @@ -2503,7 +2503,7 @@ subroutine max_overlap_qr(m, n, S0, V) ! --- return -end subroutine max_overlap_qr +end ! --- @@ -2538,7 +2538,7 @@ subroutine max_overlap_invprod(n, m, S, V) deallocate(tmp, invS) return -end subroutine max_overlap_invprod +end ! --- @@ -2623,7 +2623,7 @@ subroutine impose_biorthog_svd(n, m, L, R) deallocate(tmp, U, V, D) -end subroutine impose_biorthog_svd +end ! --- @@ -2668,8 +2668,7 @@ subroutine impose_biorthog_inverse(n, m, L, R) deallocate(S,Lt) -end subroutine impose_biorthog_inverse - +end ! --- @@ -2831,7 +2830,7 @@ subroutine impose_weighted_biorthog_qr(m, n, thr_d, thr_nd, Vl, W, Vr) call check_weighted_biorthog_binormalize(m, n, Vl, W, Vr, thr_d, thr_nd, .false.) return -end subroutine impose_weighted_biorthog_qr +end ! --- @@ -2948,7 +2947,7 @@ subroutine check_weighted_biorthog_binormalize(n, m, Vl, W, Vr, thr_d, thr_nd, s stop endif -end subroutine check_weighted_biorthog_binormalize +end ! --- @@ -3066,7 +3065,7 @@ subroutine impose_weighted_biorthog_svd(n, m, overlap, L, R) deallocate(S) return -end subroutine impose_weighted_biorthog_svd +end ! ---