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 cdc86456..c9c3b259 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 @@ -344,9 +344,9 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points 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. 80.d0) cycle + 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 + 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) diff --git a/src/ao_many_one_e_ints/listj1b.irp.f b/src/ao_many_one_e_ints/listj1b.irp.f index 1178cc31..32308f59 100644 --- a/src/ao_many_one_e_ints/listj1b.irp.f +++ b/src/ao_many_one_e_ints/listj1b.irp.f @@ -168,7 +168,7 @@ END_PROVIDER do j = 1, nucl_num tmp_alphaj = dble(List_all_comb_b3(j,i)) * j1b_pen(j) - print*,List_all_comb_b3(j,i),j1b_pen(j) +! print*,List_all_comb_b3(j,i),j1b_pen(j) List_all_comb_b3_expo(i) += tmp_alphaj List_all_comb_b3_cent(1,i) += tmp_alphaj * nucl_coord(j,1) List_all_comb_b3_cent(2,i) += tmp_alphaj * nucl_coord(j,2) diff --git a/src/dft_utils_in_r/ao_prod_mlti_pl.irp.f b/src/dft_utils_in_r/ao_prod_mlti_pl.irp.f new file mode 100644 index 00000000..22e2fe3d --- /dev/null +++ b/src/dft_utils_in_r/ao_prod_mlti_pl.irp.f @@ -0,0 +1,111 @@ +BEGIN_PROVIDER [ double precision, ao_overlap_abs_grid, (ao_num, ao_num)] + implicit none + integer :: i,j,ipoint + double precision :: contrib, weight,r(3) + ao_overlap_abs_grid = 0.D0 + do ipoint = 1,n_points_final_grid + r(:) = final_grid_points(:,ipoint) + weight = final_weight_at_r_vector(ipoint) + do i = 1, ao_num + do j = 1, ao_num + contrib = dabs(aos_in_r_array(j,ipoint) * aos_in_r_array(i,ipoint)) * weight + ao_overlap_abs_grid(j,i) += contrib + enddo + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_prod_center, (3, ao_num, ao_num)] + implicit none + BEGIN_DOC +! ao_prod_center(1:3,j,i) = \int dr |phi_i(r) phi_j(r)| x/y/z / \int |phi_i(r) phi_j(r)| +! +! if \int |phi_i(r) phi_j(r)| < 1.d-15 then ao_prod_center = 0. + END_DOC + integer :: i,j,m,ipoint + double precision :: contrib, weight,r(3) + ao_prod_center = 0.D0 + do ipoint = 1,n_points_final_grid + r(:) = final_grid_points(:,ipoint) + weight = final_weight_at_r_vector(ipoint) + do i = 1, ao_num + do j = 1, ao_num + contrib = dabs(aos_in_r_array(j,ipoint) * aos_in_r_array(i,ipoint)) * weight + do m = 1, 3 + ao_prod_center(m,j,i) += contrib * r(m) + enddo + enddo + enddo + enddo + do i = 1, ao_num + do j = 1, ao_num + if(dabs(ao_overlap_abs_grid(j,i)).gt.1.d-10)then + do m = 1, 3 + ao_prod_center(m,j,i) *= 1.d0/ao_overlap_abs_grid(j,i) + enddo + endif + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_prod_sigma, (ao_num, ao_num)] + implicit none + BEGIN_DOC +! ao_prod_sigma(i,j) = \int |phi_i(r) phi_j(r)| dsqrt((x - <|i|x|j|>)^2 + (y - <|i|y|j|>)^2 +(z - <|i|z|j|>)^2) / \int |phi_i(r) phi_j(r)| +! +! gives you a precise idea of the spatial extension of the distribution phi_i(r) phi_j(r) + END_DOC + ao_prod_sigma = 0.d0 + integer :: i,j,m,ipoint + double precision :: contrib, weight,r(3),contrib_x2 + do ipoint = 1,n_points_final_grid + r(:) = final_grid_points(:,ipoint) + weight = final_weight_at_r_vector(ipoint) + do i = 1, ao_num + do j = 1, ao_num + contrib = dabs(aos_in_r_array(j,ipoint) * aos_in_r_array(i,ipoint)) * weight + contrib_x2 = 0.d0 + do m = 1, 3 + contrib_x2 += (r(m) - ao_prod_center(m,j,i)) * (r(m) - ao_prod_center(m,j,i)) + enddo + contrib_x2 = dsqrt(contrib_x2) + ao_prod_sigma(j,i) += contrib * contrib_x2 + enddo + enddo + enddo + + do i = 1, ao_num + do j = 1, ao_num + if(dabs(ao_overlap_abs_grid(j,i)).gt.1.d-10)then + ao_prod_sigma(j,i) *= 1.d0/ao_overlap_abs_grid(j,i) + endif + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_prod_dist_grid, (ao_num, ao_num, n_points_final_grid)] + implicit none + BEGIN_DOC + ! ao_prod_dist_grid(j,i,ipoint) = distance between the center of |phi_i(r) phi_j(r)| and the grid point r(ipoint) + END_DOC + integer :: i,j,m,ipoint + double precision :: distance,r(3) + do ipoint = 1, n_points_final_grid + r(:) = final_grid_points(:,ipoint) + do i = 1, ao_num + do j = 1, ao_num + distance = 0.d0 + do m = 1, 3 + distance += (ao_prod_center(m,j,i) - r(m))*(ao_prod_center(m,j,i) - r(m)) + enddo + distance = dsqrt(distance) + ao_prod_dist_grid(j,i,ipoint) = distance + 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 f3b68f43..a515e0b8 100644 --- a/src/non_h_ints_mu/j12_nucl_utils.irp.f +++ b/src/non_h_ints_mu/j12_nucl_utils.irp.f @@ -237,6 +237,23 @@ end function j12_mu ! --- +double precision function j12_mu_r12(r12) + + include 'constants.include.F' + + implicit none + double precision, intent(in) :: r12 + double precision :: mu_r12 + + mu_r12 = mu_erf * r12 + + j12_mu_r12 = 0.5d0 * r12 * (1.d0 - derf(mu_r12)) - inv_sq_pi_2 * dexp(-mu_r12*mu_r12) / mu_erf + + return +end function j12_mu_r12 + +! --- + double precision function j12_mu_gauss(r1, r2) implicit none diff --git a/src/tc_scf/test_int.irp.f b/src/tc_scf/test_int.irp.f new file mode 100644 index 00000000..29fcb2a9 --- /dev/null +++ b/src/tc_scf/test_int.irp.f @@ -0,0 +1,77 @@ +program test_ints + + BEGIN_DOC +! TODO : Put the documentation of the program here + END_DOC + + implicit none + + print *, 'starting ...' + + 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 + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + call routine + +end + +subroutine routine + implicit none + 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 + 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_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 + 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