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 8196614f..722ff2ff 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 @@ -1,4 +1,72 @@ + +! --- + +BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2, (ao_num, ao_num, n_points_final_grid)] + + BEGIN_DOC + ! + ! -\frac{1}{4} x int dr2 phi_i(r2) phi_j(r2) [1 - erf(mu r12)]^2 + ! + END_DOC + + implicit none + integer :: i, j, ipoint, i_fit + double precision :: r(3), expo_fit, coef_fit + double precision :: tmp + double precision :: wall0, wall1 + + double precision, external :: overlap_gauss_r12_ao + + print*, ' providing int2_grad1u2_grad2u2 ...' + call wall_time(wall0) + + provide mu_erf final_grid_points j1b_pen + + int2_grad1u2_grad2u2 = 0.d0 + + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, i_fit, r, coef_fit, expo_fit, tmp) & + !$OMP SHARED (n_points_final_grid, ao_num, final_grid_points, ng_fit_jast, & + !$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2,int2_grad1u2_grad2u2) + !$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 + + tmp = 0.d0 + do i_fit = 1, ng_fit_jast + + expo_fit = expo_gauss_1_erf_x_2(i_fit) + coef_fit = coef_gauss_1_erf_x_2(i_fit) + + tmp += -0.25d0 * coef_fit * overlap_gauss_r12_ao(r, expo_fit, i, j) + enddo + + int2_grad1u2_grad2u2(j,i,ipoint) = tmp + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + do ipoint = 1, n_points_final_grid + do i = 2, ao_num + do j = 1, i-1 + int2_grad1u2_grad2u2(j,i,ipoint) = int2_grad1u2_grad2u2(i,j,ipoint) + enddo + enddo + enddo + + call wall_time(wall1) + print*, ' wall time for int2_grad1u2_grad2u2 =', wall1 - wall0 + +END_PROVIDER + ! --- BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n_points_final_grid)] 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 7a567979..963a49a6 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 @@ -53,13 +53,13 @@ subroutine compute_ao_tc_sym_two_e_pot_jl(j, l, n_integrals, buffer_i, buffer_va integral_erf = ao_two_e_integral_erf(i, k, j, l) integral = integral_erf + integral_pot - if( j1b_type .eq. 1 ) then - !print *, ' j1b type 1 is added' - integral = integral + j1b_gauss_2e_j1(i, k, j, l) - elseif( j1b_type .eq. 2 ) then - !print *, ' j1b type 2 is added' - integral = integral + j1b_gauss_2e_j2(i, k, j, l) - endif + !if( j1b_type .eq. 1 ) then + ! !print *, ' j1b type 1 is added' + ! integral = integral + j1b_gauss_2e_j1(i, k, j, l) + !elseif( j1b_type .eq. 2 ) then + ! !print *, ' j1b type 2 is added' + ! integral = integral + j1b_gauss_2e_j2(i, k, j, l) + !endif if(abs(integral) < thr) then cycle 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 b0455570..5f2795f1 100644 --- a/src/bi_ort_ints/one_e_bi_ort.irp.f +++ b/src/bi_ort_ints/one_e_bi_ort.irp.f @@ -8,22 +8,22 @@ BEGIN_PROVIDER [double precision, ao_one_e_integrals_tc_tot, (ao_num,ao_num)] ao_one_e_integrals_tc_tot = ao_one_e_integrals - provide j1b_type + !provide j1b_type - if( (j1b_type .eq. 1) .or. (j1b_type .eq. 2) ) then - - print *, ' do things properly !' - stop + !if( (j1b_type .eq. 1) .or. (j1b_type .eq. 2) ) then + ! + ! print *, ' do things properly !' + ! stop - !do i = 1, ao_num - ! do j = 1, ao_num - ! ao_one_e_integrals_tc_tot(j,i) += ( j1b_gauss_hermI (j,i) & - ! + j1b_gauss_hermII (j,i) & - ! + j1b_gauss_nonherm(j,i) ) - ! enddo - !enddo + ! !do i = 1, ao_num + ! ! do j = 1, ao_num + ! ! ao_one_e_integrals_tc_tot(j,i) += ( j1b_gauss_hermI (j,i) & + ! ! + j1b_gauss_hermII (j,i) & + ! ! + j1b_gauss_nonherm(j,i) ) + ! ! enddo + ! !enddo - endif + !endif END_PROVIDER diff --git a/src/non_h_ints_mu/jast_deriv.irp.f b/src/non_h_ints_mu/jast_deriv.irp.f index d1a044cb..4cfd13d2 100644 --- a/src/non_h_ints_mu/jast_deriv.irp.f +++ b/src/non_h_ints_mu/jast_deriv.irp.f @@ -19,8 +19,12 @@ END_DOC implicit none - integer :: ipoint, jpoint - double precision :: r1(3), r2(3) + integer :: ipoint, jpoint + double precision :: r1(3), r2(3) + double precision :: v1b_r1, v1b_r2, u2b_r12 + double precision :: grad1_v1b(3), grad1_u2b(3) + double precision :: dx, dy, dz + double precision, external :: j12_mu, j1b_nucl PROVIDE j1b_type PROVIDE final_grid_points_extra @@ -28,12 +32,43 @@ grad1_u12_num = 0.d0 grad1_u12_squared_num = 0.d0 - if((j1b_type .ge. 100) .and. (j1b_type .lt. 200)) then + if(j1b_type .eq. 100) then - double precision :: v1b_r1, v1b_r2, u2b_r12 - double precision :: grad1_v1b(3), grad1_u2b(3) - double precision :: dx, dy, dz - double precision, external :: j12_mu, j1b_nucl + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, jpoint, r1, r2, v1b_r1, v1b_r2, u2b_r12, grad1_v1b, grad1_u2b, dx, dy, dz) & + !$OMP SHARED (n_points_final_grid, n_points_extra_final_grid, final_grid_points, & + !$OMP final_grid_points_extra, grad1_u12_num, grad1_u12_squared_num) + !$OMP DO SCHEDULE (static) + do ipoint = 1, n_points_final_grid ! r1 + + r1(1) = final_grid_points(1,ipoint) + r1(2) = final_grid_points(2,ipoint) + r1(3) = final_grid_points(3,ipoint) + + do jpoint = 1, n_points_extra_final_grid ! r2 + + r2(1) = final_grid_points_extra(1,jpoint) + r2(2) = final_grid_points_extra(2,jpoint) + r2(3) = final_grid_points_extra(3,jpoint) + + call grad1_j12_mu(r1, r2, grad1_u2b) + + dx = grad1_u2b(1) + dy = grad1_u2b(2) + dz = grad1_u2b(3) + + grad1_u12_num(jpoint,ipoint,1) = dx + grad1_u12_num(jpoint,ipoint,2) = dy + grad1_u12_num(jpoint,ipoint,3) = dz + + grad1_u12_squared_num(jpoint,ipoint) = dx*dx + dy*dy + dz*dz + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + elseif((j1b_type .gt. 100) .and. (j1b_type .lt. 200)) then !$OMP PARALLEL & !$OMP DEFAULT (NONE) & diff --git a/src/non_h_ints_mu/tc_integ.irp.f b/src/non_h_ints_mu/tc_integ.irp.f index 784947b4..1a86a2e7 100644 --- a/src/non_h_ints_mu/tc_integ.irp.f +++ b/src/non_h_ints_mu/tc_integ.irp.f @@ -41,7 +41,34 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f else - if(j1b_type .eq. 3) then + if(j1b_type .eq. 0) 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(j1b_type .eq. 3) then PROVIDE v_1b_grad v_ij_erf_rk_cst_mu_j1b v_ij_u_cst_mu_j1b x_v_ij_erf_rk_cst_mu_j1b @@ -172,7 +199,27 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p PROVIDE j1b_type - if(j1b_type .eq. 3) then + if(j1b_type .eq. 0) then + + PROVIDE int2_grad1u2_grad2u2 + + int2_grad1_u12_square_ao = 0.d0 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i, j, ipoint) & + !$OMP SHARED (int2_grad1_u12_square_ao, ao_num, n_points_final_grid, int2_grad1u2_grad2u2) + !$OMP DO SCHEDULE (static) + do ipoint = 1, n_points_final_grid + do j = 1, ao_num + do i = 1, ao_num + int2_grad1_u12_square_ao(i,j,ipoint) = int2_grad1u2_grad2u2(i,j,ipoint) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + elseif(j1b_type .eq. 3) then PROVIDE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12 diff --git a/src/non_h_ints_mu/total_tc_int.irp.f b/src/non_h_ints_mu/total_tc_int.irp.f index a60c99da..450bbef0 100644 --- a/src/non_h_ints_mu/total_tc_int.irp.f +++ b/src/non_h_ints_mu/total_tc_int.irp.f @@ -11,6 +11,13 @@ BEGIN_PROVIDER [double precision, ao_vartc_int_chemist, (ao_num, ao_num, ao_num, call wall_time(wall0) if(test_cycle_tc) then + + PROVIDE j1b_type + if(j1b_type .ne. 3) then + print*, ' TC integrals with cycle can not be used for j1b_type =', j1b_type + stop + endif + do j = 1, ao_num do l = 1, ao_num do i = 1, ao_num @@ -20,7 +27,9 @@ BEGIN_PROVIDER [double precision, ao_vartc_int_chemist, (ao_num, ao_num, ao_num, enddo enddo enddo + else + do j = 1, ao_num do l = 1, ao_num do i = 1, ao_num @@ -30,6 +39,7 @@ BEGIN_PROVIDER [double precision, ao_vartc_int_chemist, (ao_num, ao_num, ao_num, enddo enddo enddo + endif call wall_time(wall1) @@ -50,6 +60,12 @@ BEGIN_PROVIDER [double precision, ao_tc_int_chemist, (ao_num, ao_num, ao_num, ao if(test_cycle_tc) then + PROVIDE j1b_type + if(j1b_type .ne. 3) then + print*, ' TC integrals with cycle can not be used for j1b_type =', j1b_type + stop + endif + ao_tc_int_chemist = ao_tc_int_chemist_test else diff --git a/src/tc_keywords/EZFIO.cfg b/src/tc_keywords/EZFIO.cfg index fad6f1c2..41c07d0b 100644 --- a/src/tc_keywords/EZFIO.cfg +++ b/src/tc_keywords/EZFIO.cfg @@ -176,7 +176,7 @@ default: 1.e-7 type: logical doc: If |true|, the integrals of the three-body jastrow are computed with cycles interface: ezfio,provider,ocaml -default: Flase +default: False [thresh_biorthog_diag] type: Threshold