diff --git a/src/ao_many_one_e_ints/grad2_jmu_manu.irp.f b/src/ao_many_one_e_ints/grad2_jmu_manu.irp.f index 65966c81..9bb16475 100644 --- a/src/ao_many_one_e_ints/grad2_jmu_manu.irp.f +++ b/src/ao_many_one_e_ints/grad2_jmu_manu.irp.f @@ -7,121 +7,6 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p ! END_DOC - implicit none - integer :: i, j, ipoint, i_1s, i_fit - 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, external :: NAI_pol_mult_erf_ao_with1s - double precision :: j12_mu_r12 - double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2 - dsqpi_3_2 = (dacos(-1.d0))**(3/2) - - provide mu_erf final_grid_points j1b_pen ao_overlap_abs - call wall_time(wall0) - - - int2_u_grad1u_j1b2_test = 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, dist, & - !$OMP sigma_ij, beta_ij, factor_ij_1s,center_ij_1s, dist_ij_ipoint, & - !$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 ao_prod_dist_grid, ao_prod_sigma, ao_overlap_abs_grid,ao_prod_center,dsqpi_3_2, & - !$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, & - !$OMP List_all_comb_b3_cent, int2_u_grad1u_j1b2_test) - !$OMP DO - do ipoint = 1, n_points_final_grid - do i = 1, ao_num - do j = i, ao_num - if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-10)cycle - dist_ij_ipoint = ao_prod_dist_grid(j,i,ipoint) ! distance to the grid point for the distribution |chi_i(r)chi_j(r)| - sigma_ij = ao_prod_sigma(j,i) ! typical spatial extension of the distribution |chi_i(r)chi_j(r)| - 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) - beta = List_all_comb_b3_expo (i_1s) -! if(beta.gt.1.d3)cycle - if(dabs(coef).lt.1.d-10)cycle - 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)) - sigma_ij = 1.d0/sigma_ij - sigma_ij *= sigma_ij - sigma_ij *= 0.5d0 - double precision :: beta_ij, factor_ij_1s, center_ij_1s(3) -! call gaussian_product(sigma_ij,ao_prod_center(1:3,j,i),beta,B_center,factor_ij_1s,beta_ij,center_ij_1s) -! if(factor_ij_1s*ao_overlap_abs_grid(j,i).lt.1.d-15)cycle -! if(factor_ij_1s*dsqpi_3_2*(beta_ij)**(-3/2)*ao_overlap_abs_grid(j,i).lt.1.d-20)cycle - - do i_fit = 1, n_max_fit_slat - - expo_fit = expo_gauss_j_mu_1_erf(i_fit) - call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s) - if(factor_ij_1s*ao_overlap_abs_grid(j,i).lt.1.d-15)cycle - coef_fit = coef_gauss_j_mu_1_erf(i_fit) - - 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 * dist - if(expo_coef_1s .gt. 20.d0) cycle - coef_tmp = coef * coef_fit * dexp(-expo_coef_1s) - if(dabs(coef_tmp) .lt. 1d-08) cycle - - int_fit = NAI_pol_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r) - - tmp += coef_tmp * int_fit - enddo - enddo - - int2_u_grad1u_j1b2_test(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_u_grad1u_j1b2_test(j,i,ipoint) = int2_u_grad1u_j1b2_test(i,j,ipoint) - enddo - enddo - enddo - - call wall_time(wall1) - print*, ' wall time for int2_u_grad1u_j1b2_test', wall1 - wall0 - -END_PROVIDER - -! --- - - -BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test_2, (ao_num, ao_num, n_points_final_grid)] - - BEGIN_DOC - ! - ! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 u_12^mu [\grad_1 u_12^mu] - ! - END_DOC - implicit none integer :: i, j, ipoint, i_1s, i_fit double precision :: r(3), int_fit, expo_fit, coef_fit, coef_tmp @@ -138,7 +23,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test_2, (ao_num, ao_num, n call wall_time(wall0) - int2_u_grad1u_j1b2_test_2 = 0.d0 + int2_u_grad1u_j1b2_test = 0.d0 !$OMP PARALLEL DEFAULT (NONE) & !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & @@ -150,7 +35,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test_2, (ao_num, ao_num, n !$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, & !$OMP ao_prod_dist_grid, ao_prod_sigma, ao_overlap_abs_grid,ao_prod_center,dsqpi_3_2, & !$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, ao_abs_comb_b3_j1b, & - !$OMP List_comb_thr_b3_cent, int2_u_grad1u_j1b2_test_2) + !$OMP List_comb_thr_b3_cent, int2_u_grad1u_j1b2_test) !$OMP DO do ipoint = 1, n_points_final_grid do i = 1, ao_num @@ -198,7 +83,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test_2, (ao_num, ao_num, n enddo enddo - int2_u_grad1u_j1b2_test_2(j,i,ipoint) = tmp + int2_u_grad1u_j1b2_test(j,i,ipoint) = tmp enddo enddo enddo @@ -208,15 +93,228 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test_2, (ao_num, ao_num, n do ipoint = 1, n_points_final_grid do i = 2, ao_num do j = 1, i-1 - int2_u_grad1u_j1b2_test_2(j,i,ipoint) = int2_u_grad1u_j1b2_test_2(i,j,ipoint) + int2_u_grad1u_j1b2_test(j,i,ipoint) = int2_u_grad1u_j1b2_test(i,j,ipoint) enddo enddo enddo call wall_time(wall1) - print*, ' wall time for int2_u_grad1u_j1b2_test_2', wall1 - wall0 + print*, ' wall time for int2_u_grad1u_j1b2_test', wall1 - wall0 END_PROVIDER ! --- +BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_no_v, (ao_num, ao_num, n_points_final_grid)] + + BEGIN_DOC + ! + ! -\frac{1}{4} x int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 [1 - erf(mu r12)]^2 + ! + END_DOC + + implicit none + integer :: i, j, ipoint, i_1s, i_fit + double precision :: r(3), expo_fit, coef_fit + double precision :: coef, beta, B_center(3) + double precision :: tmp + double precision :: wall0, wall1 + + double precision, allocatable :: int_fit_v(:) + double precision, external :: overlap_gauss_r12_ao_with1s + double precision :: factor_ij_1s,beta_ij,center_ij_1s(3),int_j1b,int_gauss,dsqpi_3_2 + dsqpi_3_2 = (dacos(-1.d0))**(3/2) + + provide mu_erf final_grid_points_transp j1b_pen List_comb_thr_b3_coef + call wall_time(wall0) + + int2_grad1u2_grad2u2_j1b2_test_no_v(:,:,:) = 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_v, tmp,int_gauss,int_j1b,factor_ij_1s,beta_ij,center_ij_1s) & +! !$OMP SHARED (n_points_final_grid, ao_num, final_grid_points,List_comb_b3_size_thr,& +! !$OMP final_grid_points_transp, n_max_fit_slat, & +! !$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, & +! !$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, & +! !$OMP List_comb_thr_b3_cent, int2_grad1u2_grad2u2_j1b2_test_no_v, ao_abs_comb_b3_j1b,& +! !$OMP ao_overlap_abs,dsqpi_3_2) +! !$OMP DO SCHEDULE(dynamic) +! do i = 1, ao_num +! do j = 1, ao_num + do i = 14,14 + do j = 17,17 + if(ao_overlap_abs(j,i) .lt. 1.d-12) then + cycle + endif + +! if(ipoint==1)then +! if(i+j.lt.10)then +! print*,j,i +! endif +! endif +! do i_1s = 1, List_comb_b3_size_thr(j,i) + do i_1s = 1, 1 + + coef = List_comb_thr_b3_coef (i_1s,j,i) + beta = List_comb_thr_b3_expo (i_1s,j,i) + int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i) +! if(dabs(coef)*dabs(int_j1b).lt.1.d-15)cycle + B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i) + B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i) + B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i) +! if(ipoint==1)then +! if(i+j.lt.10)then +! print*,coef,beta +! print*,B_center +! endif +! endif + +! do i_fit = 1, n_max_fit_slat + do i_fit = 15,15 + if(j==17.and.i==14)then + print*,i_fit,i_1s + endif +! do ipoint = 1, n_points_final_grid + do ipoint = 4,4 + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + + expo_fit = expo_gauss_1_erf_x_2(i_fit) +! call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s) + coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit) * coef +! if(dabs(coef_fit)*factor_ij_1s*dabs(int_j1b).lt.1.d-15)cycle + +! call overlap_gauss_r12_ao_with1s_v(B_center, beta, final_grid_points_transp, & +! expo_fit, i, j, int_fit_v, n_points_final_grid) +! if(ipoint == 4)then +! print*,'ipoint == 4 !!' +! endif + int_gauss = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) + + int2_grad1u2_grad2u2_j1b2_test_no_v(j,i,ipoint) += coef_fit * int_gauss + + enddo + enddo + enddo + enddo + enddo + +! !$OMP END DO +! !$OMP END PARALLEL + +! do ipoint = 1, n_points_final_grid +! do i = 1, ao_num +! do j = 1, i-1 +! int2_grad1u2_grad2u2_j1b2_test_no_v(j,i,ipoint) = int2_grad1u2_grad2u2_j1b2_test_no_v(i,j,ipoint) +! enddo +! enddo +! enddo + + call wall_time(wall1) + print*, ' wall time for int2_grad1u2_grad2u2_j1b2_test_no_v', wall1 - wall0 + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_num, n_points_final_grid)] + + BEGIN_DOC + ! + ! -\frac{1}{4} x int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 [1 - erf(mu r12)]^2 + ! + END_DOC + + implicit none + integer :: i, j, ipoint, i_1s, i_fit + double precision :: r(3), expo_fit, coef_fit + double precision :: coef, beta, B_center(3) + double precision :: tmp + double precision :: wall0, wall1 + + double precision, allocatable :: int_fit_v(:) + double precision, external :: overlap_gauss_r12_ao_with1s + + provide mu_erf final_grid_points_transp j1b_pen + call wall_time(wall0) + + double precision :: int_j1b + int2_grad1u2_grad2u2_j1b2_test(:,:,:) = 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_v, tmp,int_j1b) & +! !$OMP SHARED (n_points_final_grid, ao_num, List_comb_b3_size_thr,& +! !$OMP final_grid_points_transp, n_max_fit_slat, & +! !$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, & +! !$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, & +! !$OMP List_comb_thr_b3_cent, int2_grad1u2_grad2u2_j1b2_test,& +! !$OMP ao_abs_comb_b3_j1b,ao_overlap_abs) + + allocate(int_fit_v(n_points_final_grid)) + !$OMP DO SCHEDULE(dynamic) +! do i = 1, ao_num +! do j = i, ao_num + do i = 14,14 + do j = 17,17 + + if(ao_overlap_abs(j,i) .lt. 1.d-12) then + cycle + endif +! if(i+j.lt.10)then +! print*,j,i +! endif + +! do i_1s = 1, List_comb_b3_size_thr(j,i) + do i_1s = 1, 1 + + coef = List_comb_thr_b3_coef (i_1s,j,i) + beta = List_comb_thr_b3_expo (i_1s,j,i) + int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i) +! if(dabs(coef)*dabs(int_j1b).lt.1.d-15)cycle + B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i) + B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i) + B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i) +! if(i+j.lt.10)then +! print*,coef,beta +! print*,B_center +! endif + +! do i_fit = 1, n_max_fit_slat + do i_fit = 15,15 + + expo_fit = expo_gauss_1_erf_x_2(i_fit) + coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit) * coef + + if(j==17.and.i==14)then + print*,i_fit,i_1s + endif + call overlap_gauss_r12_ao_with1s_v(B_center, beta, final_grid_points_transp, & + expo_fit, i, j, int_fit_v, n_points_final_grid) + + do ipoint = 1, n_points_final_grid + int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) += coef_fit * int_fit_v(ipoint) + enddo + + enddo + + enddo + enddo + enddo +! !$OMP END DO + deallocate(int_fit_v) +! !$OMP END PARALLEL + + do ipoint = 1, n_points_final_grid + do i = 2, ao_num + do j = 1, i-1 + int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) = int2_grad1u2_grad2u2_j1b2_test(i,j,ipoint) + enddo + enddo + enddo + + call wall_time(wall1) + print*, ' wall time for int2_grad1u2_grad2u2_j1b2_test', wall1 - wall0 + +END_PROVIDER + diff --git a/src/ao_many_one_e_ints/listj1b_sorted.irp.f b/src/ao_many_one_e_ints/listj1b_sorted.irp.f index 58c77f5c..4ef5db99 100644 --- a/src/ao_many_one_e_ints/listj1b_sorted.irp.f +++ b/src/ao_many_one_e_ints/listj1b_sorted.irp.f @@ -13,6 +13,7 @@ coef = List_all_comb_b2_coef (i_1s) if(dabs(coef).lt.1.d-10)cycle beta = List_all_comb_b2_expo (i_1s) + beta = max(beta,1.d-10) center(1:3) = List_all_comb_b2_cent(1:3,i_1s) int_j1b = 0.d0 do ipoint = 1, n_points_final_grid @@ -92,46 +93,48 @@ END_PROVIDER END_PROVIDER + BEGIN_PROVIDER [ integer, List_comb_b3_size_thr, (ao_num, ao_num)] &BEGIN_PROVIDER [ integer, max_List_comb_b3_size_thr] implicit none integer :: i_1s,i,j,ipoint double precision :: coef,beta,center(3),int_j1b,thr double precision :: r(3),weight,dist - thr = 1.d-14 + thr = 1.d-10 List_comb_b3_size_thr = 0 do i = 1, ao_num - do j = i, ao_num + do j = 1, ao_num do i_1s = 1, List_all_comb_b3_size coef = List_all_comb_b3_coef (i_1s) - if(dabs(coef).lt.thr)cycle beta = List_all_comb_b3_expo (i_1s) center(1:3) = List_all_comb_b3_cent(1:3,i_1s) - int_j1b = 0.d0 - do ipoint = 1, n_points_final_grid - r(1:3) = final_grid_points(1:3,ipoint) - weight = final_weight_at_r_vector(ipoint) - dist = ( center(1) - r(1) )*( center(1) - r(1) ) - dist += ( center(2) - r(2) )*( center(2) - r(2) ) - dist += ( center(3) - r(3) )*( center(3) - r(3) ) - int_j1b += dabs(aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,j))*dexp(-beta*dist) * weight - enddo - if(dabs(coef)*dabs(int_j1b).gt.thr)then + if(dabs(coef).lt.thr)cycle +! int_j1b = 0.d0 +! do ipoint = 1, n_points_final_grid +! r(1:3) = final_grid_points(1:3,ipoint) +! weight = final_weight_at_r_vector(ipoint) +! dist = ( center(1) - r(1) )*( center(1) - r(1) ) +! dist += ( center(2) - r(2) )*( center(2) - r(2) ) +! dist += ( center(3) - r(3) )*( center(3) - r(3) ) +! int_j1b += dabs(aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,j))*dexp(-beta*dist) * weight +! enddo +! if(dabs(coef)*dabs(int_j1b).gt.thr)then List_comb_b3_size_thr(j,i) += 1 - endif +! endif enddo enddo enddo - do i = 1, ao_num - do j = 1, i-1 - List_comb_b3_size_thr(j,i) = List_comb_b3_size_thr(i,j) - enddo - enddo +! do i = 1, ao_num +! do j = 1, i-1 +! List_comb_b3_size_thr(j,i) = List_comb_b3_size_thr(i,j) +! enddo +! enddo integer :: list(ao_num) do i = 1, ao_num list(i) = maxval(List_comb_b3_size_thr(:,i)) enddo max_List_comb_b3_size_thr = maxval(list) + print*,'max_List_comb_b3_size_thr = ',max_List_comb_b3_size_thr END_PROVIDER @@ -143,46 +146,46 @@ END_PROVIDER integer :: i_1s,i,j,ipoint,icount double precision :: coef,beta,center(3),int_j1b,thr double precision :: r(3),weight,dist - thr = 1.d-14 + thr = 1.d-10 ao_abs_comb_b3_j1b = 10000000.d0 do i = 1, ao_num - do j = i, ao_num + do j = 1, ao_num icount = 0 do i_1s = 1, List_all_comb_b3_size coef = List_all_comb_b3_coef (i_1s) - if(dabs(coef).lt.thr)cycle beta = List_all_comb_b3_expo (i_1s) + beta = max(beta,1.d-10) center(1:3) = List_all_comb_b3_cent(1:3,i_1s) - int_j1b = 0.d0 - do ipoint = 1, n_points_final_grid - r(1:3) = final_grid_points(1:3,ipoint) - weight = final_weight_at_r_vector(ipoint) - dist = ( center(1) - r(1) )*( center(1) - r(1) ) - dist += ( center(2) - r(2) )*( center(2) - r(2) ) - dist += ( center(3) - r(3) )*( center(3) - r(3) ) - int_j1b += dabs(aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,j))*dexp(-beta*dist) * weight - enddo - if(dabs(coef)*dabs(int_j1b).gt.thr)then + if(dabs(coef).lt.thr)cycle +! int_j1b = 0.d0 +! do ipoint = 1, n_points_final_grid +! r(1:3) = final_grid_points(1:3,ipoint) +! weight = final_weight_at_r_vector(ipoint) +! dist = ( center(1) - r(1) )*( center(1) - r(1) ) +! dist += ( center(2) - r(2) )*( center(2) - r(2) ) +! dist += ( center(3) - r(3) )*( center(3) - r(3) ) +! int_j1b += dabs(aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,j))*dexp(-beta*dist) * weight +! enddo +! if(dabs(coef)*dabs(int_j1b).gt.thr)then icount += 1 List_comb_thr_b3_coef(icount,j,i) = coef List_comb_thr_b3_expo(icount,j,i) = beta List_comb_thr_b3_cent(1:3,icount,j,i) = center(1:3) - ao_abs_comb_b3_j1b(icount,j,i) = int_j1b - endif +! ao_abs_comb_b3_j1b(icount,j,i) = int_j1b +! endif enddo enddo enddo - do i = 1, ao_num - do j = 1, i-1 - do icount = 1, List_comb_b3_size_thr(j,i) - List_comb_thr_b3_coef(icount,j,i) = List_comb_thr_b3_coef(icount,i,j) - List_comb_thr_b3_expo(icount,j,i) = List_comb_thr_b3_expo(icount,i,j) - List_comb_thr_b3_cent(1,icount,j,i) = List_comb_thr_b3_cent(1,icount,i,j) - List_comb_thr_b3_cent(2,icount,j,i) = List_comb_thr_b3_cent(2,icount,i,j) - List_comb_thr_b3_cent(3,icount,j,i) = List_comb_thr_b3_cent(3,icount,i,j) - enddo - enddo - enddo - + +! do i = 1, ao_num +! do j = 1, i-1 +! do icount = 1, List_comb_b3_size_thr(j,i) +! List_comb_thr_b3_coef(icount,j,i) = List_comb_thr_b3_coef(icount,i,j) +! List_comb_thr_b3_expo(icount,j,i) = List_comb_thr_b3_expo(icount,i,j) +! List_comb_thr_b3_cent(1:3,icount,j,i) = List_comb_thr_b3_cent(1:3,icount,i,j) +! enddo +! enddo +! enddo END_PROVIDER + diff --git a/src/tc_scf/test_int.irp.f b/src/tc_scf/test_int.irp.f index 69953f02..c470d1b4 100644 --- a/src/tc_scf/test_int.irp.f +++ b/src/tc_scf/test_int.irp.f @@ -11,39 +11,47 @@ program test_ints my_grid_becke = .True. ! my_n_pt_r_grid = 30 ! my_n_pt_a_grid = 50 - my_n_pt_r_grid = 10 ! small grid for quick debug - my_n_pt_a_grid = 26 ! small grid for quick debug + my_n_pt_r_grid = 3 ! small grid for quick debug + my_n_pt_a_grid = 6 ! small grid for quick debug touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid ! call routine_int2_u_grad1u_j1b2 ! call routine_v_ij_erf_rk_cst_mu_j1b ! call routine_x_v_ij_erf_rk_cst_mu_tmp_j1b - call routine_v_ij_u_cst_mu_j1b +! call routine_v_ij_u_cst_mu_j1b + +! ! call routine_test_j1b + call routine_int2_grad1u2_grad2u2_j1b2 end subroutine routine_test_j1b implicit none integer :: i,icount,j icount = 0 -! do i = 1, List_all_comb_b2_size -! if(dabs(List_all_comb_b2_coef(i)).gt.1.d-10)then -! icount += 1 -! endif -! print*,i,List_all_comb_b2_expo(i),List_all_comb_b2_coef(i) -! enddo -! print*,'List_all_comb_b2_coef,icount = ',List_all_comb_b2_size + do i = 1, List_all_comb_b3_size + if(dabs(List_all_comb_b3_coef(i)).gt.1.d-10)then + print*,'' + print*,List_all_comb_b3_expo(i),List_all_comb_b3_coef(i) + print*,List_all_comb_b3_cent(1:3,i) + print*,'' + icount += 1 + endif + + enddo + print*,'List_all_comb_b3_coef,icount = ',List_all_comb_b3_size,icount do i = 1, ao_num do j = 1, ao_num do icount = 1, List_comb_b3_size_thr(j,i) - print*,List_comb_thr_b3_cent(1:3,icount,j,i) -! print*,'',j,i -! print*,List_comb_b2_size_thr(j,i),List_comb_b3_size_thr(j,i),ao_overlap_abs_grid(j,i) + print*,'',j,i + print*,List_comb_thr_b3_expo(icount,j,i),List_comb_thr_b3_coef(icount,j,i) + print*,List_comb_thr_b3_cent(1:3,icount,j,i) + print*,'' enddo +! enddo enddo enddo - print*,'max_List_comb_b2_size_thr = ',max_List_comb_b2_size_thr,List_all_comb_b2_size - print*,'max_List_comb_b2_size_thr = ',max_List_comb_b3_size_thr,List_all_comb_b3_size + print*,'max_List_comb_b3_size_thr = ',max_List_comb_b3_size_thr,List_all_comb_b3_size end @@ -52,19 +60,6 @@ subroutine routine_int2_u_grad1u_j1b2 integer :: i,j,ipoint,k,l double precision :: weight,accu_relat, accu_abs, contrib double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) -! print*,'ao_overlap_abs = ' -! do i = 1, ao_num -! write(*,'(100(F10.5,X))')ao_overlap_abs(i,:) -! enddo -! print*,'center = ' -! do i = 1, ao_num -! write(*,'(100(F10.5,X))')ao_prod_center(2,i,:) -! enddo -! print*,'sigma = ' -! do i = 1, ao_num -! write(*,'(100(F10.5,X))')ao_prod_sigma(i,:) -! enddo - allocate(array(ao_num, ao_num, ao_num, ao_num)) array = 0.d0 @@ -76,15 +71,7 @@ subroutine routine_int2_u_grad1u_j1b2 do l = 1, ao_num do i = 1, ao_num do j = 1, ao_num - array(j,i,l,k) += int2_u_grad1u_j1b2_test_2(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight -! if(dabs(int2_u_grad1u_j1b2(j,i,ipoint)).gt.1.d-6)then -! if(dabs(int2_u_grad1u_j1b2_test_2(j,i,ipoint)-int2_u_grad1u_j1b2(j,i,ipoint)).gt.1.d-6)then -! print*,int2_u_grad1u_j1b2(j,i,ipoint), int2_u_grad1u_j1b2_test_2(j,i,ipoint),dabs(int2_u_grad1u_j1b2_test_2(j,i,ipoint)-int2_u_grad1u_j1b2(j,i,ipoint)) -! print*,i,j -! print*,final_grid_points(:,i) -! stop -! endif -! endif + array(j,i,l,k) += int2_u_grad1u_j1b2_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight array_ref(j,i,l,k) += int2_u_grad1u_j1b2(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight enddo enddo @@ -290,4 +277,70 @@ subroutine routine_v_ij_u_cst_mu_j1b +end + +subroutine routine_int2_grad1u2_grad2u2_j1b2 + implicit none + integer :: i,j,ipoint,k,l + double precision :: weight,accu_relat, accu_abs, contrib + double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) + double precision, allocatable :: ints(:,:,:) + allocate(ints(ao_num, ao_num, n_points_final_grid)) + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + do j = 1, ao_num + read(33,*)ints(j,i,ipoint) + enddo + enddo + enddo + + allocate(array(ao_num, ao_num, ao_num, ao_num)) + array = 0.d0 + allocate(array_ref(ao_num, ao_num, ao_num, ao_num)) + array_ref = 0.d0 + do ipoint = 1, n_points_final_grid + weight = final_weight_at_r_vector(ipoint) + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + array(j,i,l,k) += int2_grad1u2_grad2u2_j1b2_test_no_v(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight +! array(j,i,l,k) += int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight + array_ref(j,i,l,k) += int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight +! array(j,i,l,k) += ints(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight +! array_ref(j,i,l,k) += int2_grad1u2_grad2u2_j1b2(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight +! array_ref(j,i,l,k) += ints(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight + if(dabs(int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint)).gt.1.d-6)then + if(dabs(int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) - int2_grad1u2_grad2u2_j1b2_test_no_v(j,i,ipoint)).gt.1.d-6)then + print*,j,i,ipoint + print*,int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) , int2_grad1u2_grad2u2_j1b2_test_no_v(j,i,ipoint), dabs(int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) - int2_grad1u2_grad2u2_j1b2_test_no_v(j,i,ipoint)) + print*,int2_grad1u2_grad2u2_j1b2_test(i,j,ipoint) , int2_grad1u2_grad2u2_j1b2_test_no_v(i,j,ipoint), dabs(int2_grad1u2_grad2u2_j1b2_test(i,j,ipoint) - int2_grad1u2_grad2u2_j1b2_test_no_v(i,j,ipoint)) + stop + endif + endif + enddo + enddo + enddo + enddo + enddo + accu_relat = 0.d0 + accu_abs = 0.d0 + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k)) + accu_abs += contrib + if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then + accu_relat += contrib/dabs(array_ref(j,i,l,k)) + endif + enddo + enddo + enddo + enddo + print*,'accu_abs = ',accu_abs/dble(ao_num)**4 + print*,'accu_relat = ',accu_relat/dble(ao_num)**4 + + + end