PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-08-24 21:21:47 +02:00
Code Releases Activity

minor modifs

Browse Source
This commit is contained in:
eginer 2022-12-08 17:44:40 +01:00
parent 5bd7b7ca6b
commit 7ce04cdd6a
3 changed files with 357 additions and 203 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

338
src/ao_many_one_e_ints/grad2_jmu_manu.irp.f
View File

@ -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

95
src/ao_many_one_e_ints/listj1b_sorted.irp.f
View File

@ -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

127
src/tc_scf/test_int.irp.f
View File

@ -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