LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-11-12 16:33:37 +01:00
Code Releases Activity

working on the grad squared

Browse Source
This commit is contained in:
eginer 2022-12-10 15:23:50 +01:00
parent ee987554e7
commit 5358bd9a61
10 changed files with 921 additions and 199 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -1,110 +1,4 @@
BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (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
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,int_j1b
double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2
double precision :: beta_ij,center_ij_1s(3),factor_ij_1s
dsqpi_3_2 = (dacos(-1.d0))**(3/2)
provide mu_erf final_grid_points j1b_pen ao_overlap_abs List_comb_thr_b3_cent
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 beta_ij,center_ij_1s,factor_ij_1s, &
!$OMP int_j1b,alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp) &
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_b3_size_thr, &
!$OMP final_grid_points, ng_fit_jast, &
!$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)
!$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
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_comb_b3_size_thr(j,i)
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-10)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)
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))
do i_fit = 1, ng_fit_jast
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*dabs(coef*int_j1b)*dsqpi_3_2*beta_ij**(-3/2).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_grad1u2_grad2u2_j1b2_test_no_v, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
@ -133,7 +27,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_no_v, (ao_num,
!$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 SHARED (n_points_final_grid, ao_num, final_grid_points,List_comb_thr_b3_size,&
!$OMP final_grid_points_transp, ng_fit_jast, &
!$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, &
!$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, &
@ -150,7 +44,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_no_v, (ao_num,
cycle
endif
do i_1s = 1, List_comb_b3_size_thr(j,i)
do i_1s = 1, List_comb_thr_b3_size(j,i)
coef = List_comb_thr_b3_coef (i_1s,j,i)
beta = List_comb_thr_b3_expo (i_1s,j,i)
@ -222,7 +116,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
!$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 SHARED (n_points_final_grid, ao_num, List_comb_thr_b3_size,&
!$OMP final_grid_points_transp, ng_fit_jast, &
!$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, &
!$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, &
@ -238,7 +132,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
cycle
endif
do i_1s = 1, List_comb_b3_size_thr(j,i)
do i_1s = 1, List_comb_thr_b3_size(j,i)
coef = List_comb_thr_b3_coef (i_1s,j,i)
beta = List_comb_thr_b3_expo (i_1s,j,i)
@ -282,3 +176,317 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
END_PROVIDER
BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (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]^2
!
END_DOC
implicit none
integer :: i, j, ipoint, i_1s, i_fit
double precision :: r(3), int_fit, expo_fit, coef_fit
double precision :: coef, beta, B_center(3), tmp
double precision :: wall0, wall1,int_j1b
double precision, external :: overlap_gauss_r12_ao
double precision, external :: overlap_gauss_r12_ao_with1s
provide mu_erf final_grid_points j1b_pen
call wall_time(wall0)
int2_u2_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, int_j1b) &
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b3_size, &
!$OMP final_grid_points, ng_fit_jast, &
!$OMP expo_gauss_j_mu_x_2, coef_gauss_j_mu_x_2, &
!$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, &
!$OMP List_comb_thr_b3_cent, int2_u2_j1b2_test,ao_abs_comb_b3_j1b)
!$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_1s = 1, List_comb_thr_b3_size(j,i)
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-10)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)
do i_fit = 1, ng_fit_jast
expo_fit = expo_gauss_j_mu_x_2(i_fit)
coef_fit = coef_gauss_j_mu_x_2(i_fit)
! ---
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
tmp += coef * coef_fit * int_fit
enddo
! ---
enddo
int2_u2_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_u2_j1b2_test(j,i,ipoint) = int2_u2_j1b2_test(i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for int2_u2_j1b2_test', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (3, 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] r2
!
END_DOC
implicit none
integer :: i, j, ipoint, i_1s, i_fit
double precision :: r(3), int_fit(3), expo_fit, coef_fit
double precision :: coef, beta, B_center(3), dist
double precision :: alpha_1s, alpha_1s_inv, centr_1s(3), expo_coef_1s, coef_tmp
double precision :: tmp_x, tmp_y, tmp_z, int_j1b
double precision :: wall0, wall1
provide mu_erf final_grid_points j1b_pen
call wall_time(wall0)
int2_u_grad1u_x_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, alpha_1s, dist, &
!$OMP alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp, &
!$OMP tmp_x, tmp_y, tmp_z,int_j1b) &
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b3_size, &
!$OMP final_grid_points, ng_fit_jast, &
!$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, &
!$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, &
!$OMP List_comb_thr_b3_cent, int2_u_grad1u_x_j1b2_test,ao_abs_comb_b3_j1b)
!$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_x = 0.d0
tmp_y = 0.d0
tmp_z = 0.d0
do i_1s = 1, List_comb_thr_b3_size(j,i)
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-10)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)
do i_fit = 1, ng_fit_jast
expo_fit = expo_gauss_j_mu_1_erf(i_fit)
coef_fit = coef_gauss_j_mu_1_erf(i_fit)
! ---
! call NAI_pol_x_mult_erf_ao_with1s(i, j, expo_fit, r, 1.d+9, r, int_fit)
! tmp_x += coef_fit * int_fit(1)
! tmp_y += coef_fit * int_fit(2)
! tmp_z += coef_fit * int_fit(3)
! if( (dabs(int_fit(1)) + dabs(int_fit(2)) + dabs(int_fit(3))) .lt. 3d-10 ) cycle
! ---
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))
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
coef_tmp = coef * coef_fit * dexp(-expo_coef_1s)
if(dabs(coef_tmp) .lt. 1d-10) cycle
call NAI_pol_x_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r, int_fit)
tmp_x += coef_tmp * int_fit(1)
tmp_y += coef_tmp * int_fit(2)
tmp_z += coef_tmp * int_fit(3)
enddo
! ---
enddo
int2_u_grad1u_x_j1b2_test(1,j,i,ipoint) = tmp_x
int2_u_grad1u_x_j1b2_test(2,j,i,ipoint) = tmp_y
int2_u_grad1u_x_j1b2_test(3,j,i,ipoint) = tmp_z
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_x_j1b2_test(1,j,i,ipoint) = int2_u_grad1u_x_j1b2_test(1,i,j,ipoint)
int2_u_grad1u_x_j1b2_test(2,j,i,ipoint) = int2_u_grad1u_x_j1b2_test(2,i,j,ipoint)
int2_u_grad1u_x_j1b2_test(3,j,i,ipoint) = int2_u_grad1u_x_j1b2_test(3,i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for int2_u_grad1u_x_j1b2_test', wall1 - wall0
END_PROVIDER
BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (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
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,int_j1b
double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2
double precision :: beta_ij,center_ij_1s(3),factor_ij_1s
dsqpi_3_2 = (dacos(-1.d0))**(3/2)
provide mu_erf final_grid_points j1b_pen ao_overlap_abs List_comb_thr_b3_cent
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 beta_ij,center_ij_1s,factor_ij_1s, &
!$OMP int_j1b,alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp) &
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b3_size, &
!$OMP final_grid_points, ng_fit_jast, &
!$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)
!$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
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_comb_thr_b3_size(j,i)
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-10)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)
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))
do i_fit = 1, ng_fit_jast
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*dabs(coef*int_j1b)*dsqpi_3_2*beta_ij**(-3/2).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
! ---

14
src/ao_many_one_e_ints/grad2_jmu_modif.irp.f
View File

@ -51,7 +51,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n
int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j)
tmp += -0.25d0 * coef_fit * int_fit
if(dabs(int_fit) .lt. 1d-10) cycle
! if(dabs(coef_fit*int_fit) .lt. 1d-12) cycle
! ---
@ -143,7 +143,7 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final
int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j)
tmp += coef_fit * int_fit
if(dabs(int_fit) .lt. 1d-10) cycle
! if(dabs(coef_fit*int_fit) .lt. 1d-12) cycle
! ---
@ -241,7 +241,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2, (3, ao_num, ao_num, n_p
tmp_x += coef_fit * int_fit(1)
tmp_y += coef_fit * int_fit(2)
tmp_z += coef_fit * int_fit(3)
if( (dabs(int_fit(1)) + dabs(int_fit(2)) + dabs(int_fit(3))) .lt. 3d-10 ) cycle
! if( dabs(coef_fit)*(dabs(int_fit(1)) + dabs(int_fit(2)) + dabs(int_fit(3))) .lt. 3d-10 ) cycle
! ---
@ -265,7 +265,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2, (3, ao_num, ao_num, n_p
expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist
coef_tmp = coef * coef_fit * dexp(-expo_coef_1s)
if(dabs(coef_tmp) .lt. 1d-10) cycle
! if(dabs(coef_tmp) .lt. 1d-12) cycle
call NAI_pol_x_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r, int_fit)
@ -351,7 +351,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points
! ---
int_fit = NAI_pol_mult_erf_ao_with1s(i, j, expo_fit, r, 1.d+9, r)
! if(dabs(int_fit) .lt. 1d-10) cycle
! if(dabs(coef_fit)*dabs(int_fit) .lt. 1d-12) cycle
tmp += coef_fit * int_fit
@ -375,9 +375,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-12) cycle
int_fit = NAI_pol_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r)

12
src/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f
View File

@ -25,7 +25,7 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num,
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, int_mu, int_coulomb, tmp, int_j1b)&
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_b2_size_thr, final_grid_points, &
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b2_size, final_grid_points, &
!$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo, List_comb_thr_b2_cent,ao_abs_comb_b2_j1b, &
!$OMP v_ij_erf_rk_cst_mu_j1b_test, mu_erf, &
!$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2)
@ -41,7 +41,7 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num,
if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle
tmp = 0.d0
do i_1s = 1, List_comb_b2_size_thr(j,i)
do i_1s = 1, List_comb_thr_b2_size(j,i)
coef = List_comb_thr_b2_coef (i_1s,j,i)
beta = List_comb_thr_b2_expo (i_1s,j,i)
@ -129,7 +129,7 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp_j1b_test, (3, ao_num
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, ints, ints_coulomb, &
!$OMP int_j1b, tmp_x, tmp_y, tmp_z) &
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_b2_size_thr, final_grid_points,&
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b2_size, final_grid_points,&
!$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo, List_comb_thr_b2_cent, &
!$OMP x_v_ij_erf_rk_cst_mu_tmp_j1b_test, mu_erf,ao_abs_comb_b2_j1b, &
!$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2)
@ -147,7 +147,7 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp_j1b_test, (3, ao_num
tmp_x = 0.d0
tmp_y = 0.d0
tmp_z = 0.d0
do i_1s = 1, List_comb_b2_size_thr(j,i)
do i_1s = 1, List_comb_thr_b2_size(j,i)
coef = List_comb_thr_b2_coef (i_1s,j,i)
beta = List_comb_thr_b2_expo (i_1s,j,i)
@ -223,7 +223,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
!$OMP SHARED (n_points_final_grid, ao_num, &
!$OMP final_grid_points, ng_fit_jast, &
!$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, &
!$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_b2_size_thr, &
!$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, &
!$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_test,ao_abs_comb_b2_j1b, &
!$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2)
!$OMP DO
@ -238,7 +238,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle
tmp = 0.d0
do i_1s = 1, List_comb_b2_size_thr(j,i)
do i_1s = 1, List_comb_thr_b2_size(j,i)
coef = List_comb_thr_b2_coef (i_1s,j,i)
beta = List_comb_thr_b2_expo (i_1s,j,i)

6
src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f
View File

@ -49,7 +49,7 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po
int_mu = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r)
int_coulomb = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r)
if(dabs(int_mu - int_coulomb) .lt. 1d-10) cycle
! if(dabs(coef)*dabs(int_mu - int_coulomb) .lt. 1d-12) cycle
tmp += coef * (int_mu - int_coulomb)
@ -169,7 +169,7 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp_j1b, (3, ao_num, ao_
call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, ints )
call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, ints_coulomb)
if( (dabs(ints(1)-ints_coulomb(1)) + dabs(ints(2)-ints_coulomb(2)) + dabs(ints(3)-ints_coulomb(3))) .lt. 3d-10) cycle
! if( dabs(coef)*(dabs(ints(1)-ints_coulomb(1)) + dabs(ints(2)-ints_coulomb(2)) + dabs(ints(3)-ints_coulomb(3))) .lt. 3d-10) cycle
tmp_x += coef * (ints(1) - ints_coulomb(1))
tmp_y += coef * (ints(2) - ints_coulomb(2))
@ -277,7 +277,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_
B_center(3) = List_all_comb_b2_cent(3,1)
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
if(dabs(int_fit) .lt. 1d-10) cycle
! if(dabs(int_fit*coef) .lt. 1d-12) cycle
tmp += coef * coef_fit * int_fit

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

@ -1,70 +1,70 @@
BEGIN_PROVIDER [ integer, List_comb_b2_size_thr, (ao_num, ao_num)]
&BEGIN_PROVIDER [ integer, max_List_comb_b2_size_thr]
BEGIN_PROVIDER [ integer, List_comb_thr_b2_size, (ao_num, ao_num)]
&BEGIN_PROVIDER [ integer, max_List_comb_thr_b2_size]
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-10
List_comb_b2_size_thr = 0
thr = 1.d-12
List_comb_thr_b2_size = 0
do i = 1, ao_num
do j = i, ao_num
do i_1s = 1, List_all_comb_b2_size
coef = List_all_comb_b2_coef (i_1s)
if(dabs(coef).lt.1.d-10)cycle
if(dabs(coef).lt.1.d-12)cycle
beta = List_all_comb_b2_expo (i_1s)
beta = max(beta,1.d-10)
beta = max(beta,1.d-12)
center(1:3) = List_all_comb_b2_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)
do ipoint = 1, n_points_extra_final_grid
r(1:3) = final_grid_points_extra(1:3,ipoint)
weight = final_weight_at_r_vector_extra(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_b2_size_thr(j,i) += 1
List_comb_thr_b2_size(j,i) += 1
endif
enddo
enddo
enddo
do i = 1, ao_num
do j = 1, i-1
List_comb_b2_size_thr(j,i) = List_comb_b2_size_thr(i,j)
List_comb_thr_b2_size(j,i) = List_comb_thr_b2_size(i,j)
enddo
enddo
integer :: list(ao_num)
do i = 1, ao_num
list(i) = maxval(List_comb_b2_size_thr(:,i))
list(i) = maxval(List_comb_thr_b2_size(:,i))
enddo
max_List_comb_b2_size_thr = maxval(list)
max_List_comb_thr_b2_size = maxval(list)
END_PROVIDER
BEGIN_PROVIDER [ double precision, List_comb_thr_b2_coef, ( max_List_comb_b2_size_thr,ao_num, ao_num )]
&BEGIN_PROVIDER [ double precision, List_comb_thr_b2_expo, ( max_List_comb_b2_size_thr,ao_num, ao_num )]
&BEGIN_PROVIDER [ double precision, List_comb_thr_b2_cent, (3, max_List_comb_b2_size_thr,ao_num, ao_num )]
&BEGIN_PROVIDER [ double precision, ao_abs_comb_b2_j1b, ( max_List_comb_b2_size_thr ,ao_num, ao_num)]
BEGIN_PROVIDER [ double precision, List_comb_thr_b2_coef, ( max_List_comb_thr_b2_size,ao_num, ao_num )]
&BEGIN_PROVIDER [ double precision, List_comb_thr_b2_expo, ( max_List_comb_thr_b2_size,ao_num, ao_num )]
&BEGIN_PROVIDER [ double precision, List_comb_thr_b2_cent, (3, max_List_comb_thr_b2_size,ao_num, ao_num )]
&BEGIN_PROVIDER [ double precision, ao_abs_comb_b2_j1b, ( max_List_comb_thr_b2_size ,ao_num, ao_num)]
implicit none
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-10
thr = 1.d-12
ao_abs_comb_b2_j1b = 10000000.d0
do i = 1, ao_num
do j = i, ao_num
icount = 0
do i_1s = 1, List_all_comb_b2_size
coef = List_all_comb_b2_coef (i_1s)
if(dabs(coef).lt.1.d-10)cycle
if(dabs(coef).lt.1.d-12)cycle
beta = List_all_comb_b2_expo (i_1s)
center(1:3) = List_all_comb_b2_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)
do ipoint = 1, n_points_extra_final_grid
r(1:3) = final_grid_points_extra(1:3,ipoint)
weight = final_weight_at_r_vector_extra(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) )
@ -83,7 +83,7 @@ END_PROVIDER
do i = 1, ao_num
do j = 1, i-1
do icount = 1, List_comb_b2_size_thr(j,i)
do icount = 1, List_comb_thr_b2_size(j,i)
List_comb_thr_b2_coef(icount,j,i) = List_comb_thr_b2_coef(icount,i,j)
List_comb_thr_b2_expo(icount,j,i) = List_comb_thr_b2_expo(icount,i,j)
List_comb_thr_b2_cent(1:3,icount,j,i) = List_comb_thr_b2_cent(1:3,icount,i,j)
@ -94,14 +94,14 @@ 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]
BEGIN_PROVIDER [ integer, List_comb_thr_b3_size, (ao_num, ao_num)]
&BEGIN_PROVIDER [ integer, max_List_comb_thr_b3_size]
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-10
List_comb_b3_size_thr = 0
thr = 1.d-12
List_comb_thr_b3_size = 0
do i = 1, ao_num
do j = 1, ao_num
do i_1s = 1, List_all_comb_b3_size
@ -110,43 +110,43 @@ END_PROVIDER
center(1:3) = List_all_comb_b3_cent(1:3,i_1s)
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)
do ipoint = 1, n_points_extra_final_grid
r(1:3) = final_grid_points_extra(1:3,ipoint)
weight = final_weight_at_r_vector_extra(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
List_comb_thr_b3_size(j,i) += 1
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)
! List_comb_thr_b3_size(j,i) = List_comb_thr_b3_size(i,j)
! enddo
! enddo
integer :: list(ao_num)
do i = 1, ao_num
list(i) = maxval(List_comb_b3_size_thr(:,i))
list(i) = maxval(List_comb_thr_b3_size(:,i))
enddo
max_List_comb_b3_size_thr = maxval(list)
print*,'max_List_comb_b3_size_thr = ',max_List_comb_b3_size_thr
max_List_comb_thr_b3_size = maxval(list)
print*,'max_List_comb_thr_b3_size = ',max_List_comb_thr_b3_size
END_PROVIDER
BEGIN_PROVIDER [ double precision, List_comb_thr_b3_coef, ( max_List_comb_b3_size_thr,ao_num, ao_num )]
&BEGIN_PROVIDER [ double precision, List_comb_thr_b3_expo, ( max_List_comb_b3_size_thr,ao_num, ao_num )]
&BEGIN_PROVIDER [ double precision, List_comb_thr_b3_cent, (3, max_List_comb_b3_size_thr,ao_num, ao_num )]
&BEGIN_PROVIDER [ double precision, ao_abs_comb_b3_j1b, ( max_List_comb_b3_size_thr ,ao_num, ao_num)]
BEGIN_PROVIDER [ double precision, List_comb_thr_b3_coef, ( max_List_comb_thr_b3_size,ao_num, ao_num )]
&BEGIN_PROVIDER [ double precision, List_comb_thr_b3_expo, ( max_List_comb_thr_b3_size,ao_num, ao_num )]
&BEGIN_PROVIDER [ double precision, List_comb_thr_b3_cent, (3, max_List_comb_thr_b3_size,ao_num, ao_num )]
&BEGIN_PROVIDER [ double precision, ao_abs_comb_b3_j1b, ( max_List_comb_thr_b3_size ,ao_num, ao_num)]
implicit none
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-10
thr = 1.d-12
ao_abs_comb_b3_j1b = 10000000.d0
do i = 1, ao_num
do j = 1, ao_num
@ -154,13 +154,13 @@ END_PROVIDER
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)
beta = max(beta,1.d-10)
beta = max(beta,1.d-12)
center(1:3) = List_all_comb_b3_cent(1:3,i_1s)
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)
do ipoint = 1, n_points_extra_final_grid
r(1:3) = final_grid_points_extra(1:3,ipoint)
weight = final_weight_at_r_vector_extra(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) )
@ -179,7 +179,7 @@ END_PROVIDER
! do i = 1, ao_num
! do j = 1, i-1
! do icount = 1, List_comb_b3_size_thr(j,i)
! do icount = 1, List_comb_thr_b3_size(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)

1
src/non_h_ints_mu/grad_squared.irp.f
View File

@ -290,6 +290,7 @@ BEGIN_PROVIDER [ double precision, u12sq_j1bsq, (ao_num, ao_num, n_points_final_
END_PROVIDER
! ---
! ---
BEGIN_PROVIDER [ double precision, u12_grad1_u12_j1b_grad1_j1b, (ao_num, ao_num, n_points_final_grid) ]

134
src/non_h_ints_mu/grad_squared_manu.irp.f Normal file
View File

@ -0,0 +1,134 @@
BEGIN_PROVIDER [double precision, tc_grad_square_ao_test, (ao_num, ao_num, ao_num, ao_num)]
BEGIN_DOC
!
! tc_grad_square_ao_test(k,i,l,j) = -1/2 <kl | |\grad_1 u(r1,r2)|^2 + |\grad_1 u(r1,r2)|^2 | ij>
!
END_DOC
implicit none
integer :: ipoint, i, j, k, l
double precision :: weight1, ao_ik_r, ao_i_r
double precision, allocatable :: ac_mat(:,:,:,:), bc_mat(:,:,:,:)
allocate(ac_mat(ao_num,ao_num,ao_num,ao_num))
ac_mat = 0.d0
allocate(bc_mat(ao_num,ao_num,ao_num,ao_num))
bc_mat = 0.d0
do ipoint = 1, n_points_final_grid
weight1 = final_weight_at_r_vector(ipoint)
do i = 1, ao_num
ao_i_r = weight1 * aos_in_r_array_transp(ipoint,i)
do k = 1, ao_num
ao_ik_r = ao_i_r * aos_in_r_array_transp(ipoint,k)
do j = 1, ao_num
do l = 1, ao_num
ac_mat(k,i,l,j) += ao_ik_r * ( u12sq_j1bsq_test(l,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b_test(l,j,ipoint) )
bc_mat(k,i,l,j) += ao_ik_r * grad12_j12(l,j,ipoint)
enddo
enddo
enddo
enddo
enddo
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
tc_grad_square_ao_test(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i) + bc_mat(k,i,l,j)
enddo
enddo
enddo
enddo
deallocate(ac_mat)
deallocate(bc_mat)
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, u12sq_j1bsq_test, (ao_num, ao_num, n_points_final_grid) ]
implicit none
integer :: ipoint, i, j
double precision :: tmp_x, tmp_y, tmp_z
double precision :: tmp1
double precision :: time0, time1
print*, ' providing u12sq_j1bsq_test ...'
call wall_time(time0)
do ipoint = 1, n_points_final_grid
tmp_x = v_1b_grad(1,ipoint)
tmp_y = v_1b_grad(2,ipoint)
tmp_z = v_1b_grad(3,ipoint)
tmp1 = -0.5d0 * (tmp_x * tmp_x + tmp_y * tmp_y + tmp_z * tmp_z)
do j = 1, ao_num
do i = 1, ao_num
u12sq_j1bsq_test(i,j,ipoint) = tmp1 * int2_u2_j1b2_test(i,j,ipoint)
enddo
enddo
enddo
call wall_time(time1)
print*, ' Wall time for u12sq_j1bsq_test = ', time1 - time0
END_PROVIDER
BEGIN_PROVIDER [ double precision, u12_grad1_u12_j1b_grad1_j1b_test, (ao_num, ao_num, n_points_final_grid) ]
implicit none
integer :: ipoint, i, j, m, igauss
double precision :: x, y, z
double precision :: tmp_v, tmp_x, tmp_y, tmp_z
double precision :: tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9
double precision :: time0, time1
double precision, external :: overlap_gauss_r12_ao
print*, ' providing u12_grad1_u12_j1b_grad1_j1b_test ...'
call wall_time(time0)
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)
tmp_v = v_1b (ipoint)
tmp_x = v_1b_grad(1,ipoint)
tmp_y = v_1b_grad(2,ipoint)
tmp_z = v_1b_grad(3,ipoint)
tmp3 = tmp_v * tmp_x
tmp4 = tmp_v * tmp_y
tmp5 = tmp_v * tmp_z
tmp6 = -x * tmp3
tmp7 = -y * tmp4
tmp8 = -z * tmp5
do j = 1, ao_num
do i = 1, ao_num
tmp9 = int2_u_grad1u_j1b2(i,j,ipoint)
u12_grad1_u12_j1b_grad1_j1b_test(i,j,ipoint) = tmp6 * tmp9 + tmp3 * int2_u_grad1u_x_j1b2_test(1,i,j,ipoint) &
+ tmp7 * tmp9 + tmp4 * int2_u_grad1u_x_j1b2_test(2,i,j,ipoint) &
+ tmp8 * tmp9 + tmp5 * int2_u_grad1u_x_j1b2_test(3,i,j,ipoint)
enddo
enddo
enddo
call wall_time(time1)
print*, ' Wall time for u12_grad1_u12_j1b_grad1_j1b_test = ', time1 - time0
END_PROVIDER
! ---

2
src/non_h_ints_mu/new_grad_tc.irp.f
View File

@ -82,6 +82,7 @@ END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num, ao_num)]
BEGIN_DOC
@ -149,3 +150,4 @@ END_PROVIDER
! ---

146
src/non_h_ints_mu/new_grad_tc_manu.irp.f Normal file
View File

@ -0,0 +1,146 @@
BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_test, (3, ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! int2_grad1_u12_ao_test(:,i,j,ipoint) = \int dr2 [-1 * \grad_r1 J(r1,r2)] \phi_i(r2) \phi_j(r2)
!
! where r1 = r(ipoint)
!
! if J(r1,r2) = u12:
!
! int2_grad1_u12_ao_test(:,i,j,ipoint) = 0.5 x \int dr2 [(r1 - r2) (erf(mu * r12)-1)r_12] \phi_i(r2) \phi_j(r2)
! = 0.5 * [ v_ij_erf_rk_cst_mu(i,j,ipoint) * r(:) - x_v_ij_erf_rk_cst_mu(i,j,ipoint,:) ]
!
! if J(r1,r2) = u12 x v1 x v2
!
! int2_grad1_u12_ao_test(:,i,j,ipoint) = v1 x [ 0.5 x \int dr2 [(r1 - r2) (erf(mu * r12)-1)r_12] v2 \phi_i(r2) \phi_j(r2) ]
! - \grad_1 v1 x [ \int dr2 u12 v2 \phi_i(r2) \phi_j(r2) ]
! = 0.5 v_1b(ipoint) * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) * r(:)
! - 0.5 v_1b(ipoint) * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,:)
! - v_1b_grad[:,ipoint] * v_ij_u_cst_mu_j1b(i,j,ipoint)
!
!
END_DOC
implicit none
integer :: ipoint, i, j
double precision :: x, y, z, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2
PROVIDE j1b_type
if(j1b_type .eq. 3) then
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 * v_1b(ipoint)
tmp_x = v_1b_grad(1,ipoint)
tmp_y = v_1b_grad(2,ipoint)
tmp_z = v_1b_grad(3,ipoint)
do j = 1, ao_num
do i = 1, ao_num
tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint)
tmp2 = v_ij_u_cst_mu_j1b_test(i,j,ipoint)
int2_grad1_u12_ao_test(1,i,j,ipoint) = tmp1 * x - tmp0 * x_v_ij_erf_rk_cst_mu_tmp_j1b_test(1,i,j,ipoint) - tmp2 * tmp_x
int2_grad1_u12_ao_test(2,i,j,ipoint) = tmp1 * y - tmp0 * x_v_ij_erf_rk_cst_mu_tmp_j1b_test(2,i,j,ipoint) - tmp2 * tmp_y
int2_grad1_u12_ao_test(3,i,j,ipoint) = tmp1 * z - tmp0 * x_v_ij_erf_rk_cst_mu_tmp_j1b_test(3,i,j,ipoint) - tmp2 * tmp_z
enddo
enddo
enddo
else
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_test(1,i,j,ipoint) = tmp1 * x - x_v_ij_erf_rk_cst_mu_tmp(1,i,j,ipoint)
int2_grad1_u12_ao_test(2,i,j,ipoint) = tmp1 * y - x_v_ij_erf_rk_cst_mu_tmp(2,i,j,ipoint)
int2_grad1_u12_ao_test(3,i,j,ipoint) = tmp1 * z - x_v_ij_erf_rk_cst_mu_tmp(3,i,j,ipoint)
enddo
enddo
enddo
int2_grad1_u12_ao_test *= 0.5d0
endif
END_PROVIDER
BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao_test, (ao_num, ao_num, ao_num, ao_num)]
BEGIN_DOC
!
! tc_grad_and_lapl_ao_test(k,i,l,j) = < k l | -1/2 \Delta_1 u(r1,r2) - \grad_1 u(r1,r2) | ij >
!
! = 1/2 \int dr1 (phi_k(r1) \grad_r1 phi_i(r1) - phi_i(r1) \grad_r1 phi_k(r1)) . \int dr2 \grad_r1 u(r1,r2) \phi_l(r2) \phi_j(r2)
!
! This is obtained by integration by parts.
!
END_DOC
implicit none
integer :: ipoint, i, j, k, l
double precision :: weight1, contrib_x, contrib_y, contrib_z, tmp_x, tmp_y, tmp_z
double precision :: ao_k_r, ao_i_r, ao_i_dx, ao_i_dy, ao_i_dz
double precision, allocatable :: ac_mat(:,:,:,:)
allocate(ac_mat(ao_num,ao_num,ao_num,ao_num))
ac_mat = 0.d0
do ipoint = 1, n_points_final_grid
weight1 = 0.5d0 * final_weight_at_r_vector(ipoint)
do i = 1, ao_num
ao_i_r = weight1 * aos_in_r_array_transp (ipoint,i)
ao_i_dx = weight1 * aos_grad_in_r_array_transp_bis(ipoint,i,1)
ao_i_dy = weight1 * aos_grad_in_r_array_transp_bis(ipoint,i,2)
ao_i_dz = weight1 * aos_grad_in_r_array_transp_bis(ipoint,i,3)
do k = 1, ao_num
ao_k_r = aos_in_r_array_transp(ipoint,k)
tmp_x = ao_k_r * ao_i_dx - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,1)
tmp_y = ao_k_r * ao_i_dy - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,2)
tmp_z = ao_k_r * ao_i_dz - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,3)
do j = 1, ao_num
do l = 1, ao_num
contrib_x = int2_grad1_u12_ao_test(1,l,j,ipoint) * tmp_x
contrib_y = int2_grad1_u12_ao_test(2,l,j,ipoint) * tmp_y
contrib_z = int2_grad1_u12_ao_test(3,l,j,ipoint) * tmp_z
ac_mat(k,i,l,j) += contrib_x + contrib_y + contrib_z
enddo
enddo
enddo
enddo
enddo
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
tc_grad_and_lapl_ao_test(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i)
enddo
enddo
enddo
enddo
deallocate(ac_mat)
END_PROVIDER
! ---

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

@ -9,20 +9,37 @@ program test_ints
print *, 'starting ...'
my_grid_becke = .True.
! my_n_pt_r_grid = 30
my_n_pt_r_grid = 10
! 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 = 10 ! small grid for quick debug
my_n_pt_a_grid = 14 ! 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_test_j1b
my_n_pt_r_extra_grid = 30
my_n_pt_a_extra_grid = 50 ! small extra_grid for quick debug
touch my_extra_grid_becke my_n_pt_r_extra_grid my_n_pt_a_extra_grid
!! OK
!call routine_int2_u_grad1u_j1b2
!! OK
!call routine_v_ij_erf_rk_cst_mu_j1b
!! OK
! call routine_x_v_ij_erf_rk_cst_mu_tmp_j1b
!! OK
! call routine_v_ij_u_cst_mu_j1b
!! OK
!call routine_int2_u2_j1b2
!! OK
!call routine_int2_u_grad1u_x_j1b2
!! OK
! call routine_int2_grad1u2_grad2u2_j1b2
! call routine_int2_u_grad1u_j1b2
! call test_total_grad_lapl
call test_total_grad_square
end
subroutine routine_test_j1b
@ -42,7 +59,7 @@ subroutine routine_test_j1b
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)
do icount = 1, List_comb_thr_b3_size(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)
@ -51,7 +68,7 @@ subroutine routine_test_j1b
! enddo
enddo
enddo
print*,'max_List_comb_b3_size_thr = ',max_List_comb_b3_size_thr,List_all_comb_b3_size
print*,'max_List_comb_thr_b3_size = ',max_List_comb_thr_b3_size,List_all_comb_b3_size
end
@ -221,7 +238,7 @@ end
subroutine routine_v_ij_u_cst_mu_j1b
subroutine routine_v_ij_u_cst_mu_j1b_test
implicit none
integer :: i,j,ipoint,k,l
double precision :: weight,accu_relat, accu_abs, contrib
@ -286,13 +303,13 @@ subroutine routine_int2_grad1u2_grad2u2_j1b2
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
! 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
@ -306,18 +323,18 @@ subroutine routine_int2_grad1u2_grad2u2_j1b2
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_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) += 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))
! 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
! stop
! endif
! endif
enddo
enddo
enddo
@ -343,4 +360,218 @@ subroutine routine_int2_grad1u2_grad2u2_j1b2
end
subroutine routine_int2_u2_j1b2
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_u2_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_u2_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
subroutine routine_int2_u_grad1u_x_j1b2
implicit none
integer :: i,j,ipoint,k,l,m
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
do m = 1, 3
array(j,i,l,k) += int2_u_grad1u_x_j1b2_test(m,j,i,ipoint) * aos_grad_in_r_array_transp(m,k,ipoint) * aos_in_r_array(l,ipoint) * weight
array_ref(j,i,l,k) += int2_u_grad1u_x_j1b2(m,j,i,ipoint) * aos_grad_in_r_array_transp(m,k,ipoint) * aos_in_r_array(l,ipoint) * weight
enddo
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
subroutine routine_v_ij_u_cst_mu_j1b
implicit none
integer :: i,j,ipoint,k,l
double precision :: weight,accu_relat, accu_abs, contrib
double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
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) += v_ij_u_cst_mu_j1b_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
array_ref(j,i,l,k) += v_ij_u_cst_mu_j1b(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
subroutine test_total_grad_lapl
implicit none
integer :: i,j,ipoint,k,l
double precision :: weight,accu_relat, accu_abs, contrib
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(tc_grad_and_lapl_ao_test(j,i,l,k) - tc_grad_and_lapl_ao(j,i,l,k))
accu_abs += contrib
if(dabs(tc_grad_and_lapl_ao(j,i,l,k)).gt.1.d-10)then
accu_relat += contrib/dabs(tc_grad_and_lapl_ao(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
subroutine test_total_grad_square
implicit none
integer :: i,j,ipoint,k,l
double precision :: weight,accu_relat, accu_abs, contrib
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(tc_grad_square_ao_test(j,i,l,k) - tc_grad_square_ao(j,i,l,k))
accu_abs += contrib
if(dabs(tc_grad_square_ao(j,i,l,k)).gt.1.d-10)then
accu_relat += contrib/dabs(tc_grad_square_ao(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