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

improved the building of tc_grad_and_lapl_ao_test and #tc_grad_square_ao_test

Browse Source
This commit is contained in:
eginer 2022-12-13 18:40:37 +01:00
parent 23ec3ba18a
commit a34653b5d1
6 changed files with 138 additions and 96 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -1,5 +1,5 @@
BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_no_v, (ao_num, ao_num, n_points_final_grid)] BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC BEGIN_DOC
! !
@ -16,13 +16,14 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_no_v, (ao_num,
double precision, allocatable :: int_fit_v(:) double precision, allocatable :: int_fit_v(:)
double precision, external :: overlap_gauss_r12_ao_with1s 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 double precision :: int_gauss,dsqpi_3_2,int_j1b
double precision :: factor_ij_1s,beta_ij,center_ij_1s(3)
dsqpi_3_2 = (dacos(-1.d0))**(3/2) dsqpi_3_2 = (dacos(-1.d0))**(3/2)
provide mu_erf final_grid_points_transp j1b_pen List_comb_thr_b3_coef provide mu_erf final_grid_points_transp j1b_pen List_comb_thr_b3_coef
call wall_time(wall0) call wall_time(wall0)
int2_grad1u2_grad2u2_j1b2_test_no_v(:,:,:) = 0.d0 int2_grad1u2_grad2u2_j1b2_test(:,:,:) = 0.d0
!$OMP PARALLEL DEFAULT (NONE) & !$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center,& !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center,&
@ -31,7 +32,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_no_v, (ao_num,
!$OMP final_grid_points_transp, ng_fit_jast, & !$OMP final_grid_points_transp, ng_fit_jast, &
!$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, & !$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_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 List_comb_thr_b3_cent, int2_grad1u2_grad2u2_j1b2_test, ao_abs_comb_b3_j1b,&
!$OMP ao_overlap_abs,dsqpi_3_2) !$OMP ao_overlap_abs,dsqpi_3_2)
!$OMP DO SCHEDULE(dynamic) !$OMP DO SCHEDULE(dynamic)
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
@ -49,7 +50,6 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_no_v, (ao_num,
coef = List_comb_thr_b3_coef (i_1s,j,i) coef = List_comb_thr_b3_coef (i_1s,j,i)
beta = List_comb_thr_b3_expo (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) 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(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(2) = List_comb_thr_b3_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i) B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
@ -57,15 +57,16 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_no_v, (ao_num,
do i_fit = 1, ng_fit_jast do i_fit = 1, ng_fit_jast
expo_fit = expo_gauss_1_erf_x_2(i_fit) 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) !DIR$ FORCEINLINE
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 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 if(dabs(coef_fit*factor_ij_1s*int_j1b).lt.1.d-10)cycle
! call overlap_gauss_r12_ao_with1s_v(B_center, beta, final_grid_points_transp, & ! 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) ! expo_fit, i, j, int_fit_v, n_points_final_grid)
int_gauss = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) 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 int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) += coef_fit * int_gauss
enddo enddo
enddo enddo
@ -79,17 +80,17 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_no_v, (ao_num,
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
do i = 1, ao_num do i = 1, ao_num
do j = 1, i-1 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) int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) = int2_grad1u2_grad2u2_j1b2_test(i,j,ipoint)
enddo enddo
enddo enddo
enddo enddo
call wall_time(wall1) call wall_time(wall1)
print*, ' wall time for int2_grad1u2_grad2u2_j1b2_test_no_v', wall1 - wall0 print*, ' wall time for int2_grad1u2_grad2u2_j1b2_test', wall1 - wall0
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_num, n_points_final_grid)] BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_v, (ao_num, ao_num, n_points_final_grid)]
! !
! BEGIN_DOC ! BEGIN_DOC
! ! ! !
@ -104,15 +105,15 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
double precision :: tmp double precision :: tmp
double precision :: wall0, wall1 double precision :: wall0, wall1
double precision, allocatable :: int_fit_v(:) double precision, allocatable :: int_fit_v(:),big_array(:,:,:)
double precision, external :: overlap_gauss_r12_ao_with1s double precision, external :: overlap_gauss_r12_ao_with1s
provide mu_erf final_grid_points_transp j1b_pen provide mu_erf final_grid_points_transp j1b_pen
call wall_time(wall0) call wall_time(wall0)
double precision :: int_j1b double precision :: int_j1b
int2_grad1u2_grad2u2_j1b2_test(:,:,:) = 0.d0 big_array(:,:,:) = 0.d0
! allocate(big_array(n_points_final_grid,ao_num, ao_num))
!$OMP PARALLEL DEFAULT (NONE) & !$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center,& !$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 coef_fit, expo_fit, int_fit_v, tmp,int_j1b) &
@ -120,7 +121,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
!$OMP final_grid_points_transp, ng_fit_jast, & !$OMP final_grid_points_transp, ng_fit_jast, &
!$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, & !$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_coef, List_comb_thr_b3_expo, &
!$OMP List_comb_thr_b3_cent, int2_grad1u2_grad2u2_j1b2_test,& !$OMP List_comb_thr_b3_cent, big_array,&
!$OMP ao_abs_comb_b3_j1b,ao_overlap_abs) !$OMP ao_abs_comb_b3_j1b,ao_overlap_abs)
! !
allocate(int_fit_v(n_points_final_grid)) allocate(int_fit_v(n_points_final_grid))
@ -151,7 +152,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
expo_fit, i, j, int_fit_v, size(int_fit_v,1),n_points_final_grid) expo_fit, i, j, int_fit_v, size(int_fit_v,1),n_points_final_grid)
do ipoint = 1, 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) big_array(ipoint,j,i) += coef_fit * int_fit_v(ipoint)
enddo enddo
enddo enddo
@ -162,17 +163,24 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
!$OMP END DO !$OMP END DO
deallocate(int_fit_v) deallocate(int_fit_v)
!$OMP END PARALLEL !$OMP END PARALLEL
do i = 1, ao_num
do j = i, ao_num
do ipoint = 1, n_points_final_grid
int2_grad1u2_grad2u2_j1b2_test_v(j,i,ipoint) = big_array(ipoint,j,i)
enddo
enddo
enddo
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
do i = 2, ao_num do i = 2, ao_num
do j = 1, i-1 do j = 1, i-1
int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) = int2_grad1u2_grad2u2_j1b2_test(i,j,ipoint) int2_grad1u2_grad2u2_j1b2_test_v(j,i,ipoint) = big_array(ipoint,i,j)
enddo enddo
enddo enddo
enddo enddo
call wall_time(wall1) call wall_time(wall1)
print*, ' wall time for int2_grad1u2_grad2u2_j1b2_test', wall1 - wall0 print*, ' wall time for int2_grad1u2_grad2u2_j1b2_test_v', wall1 - wall0
END_PROVIDER END_PROVIDER
@ -192,6 +200,7 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
double precision, external :: overlap_gauss_r12_ao double precision, external :: overlap_gauss_r12_ao
double precision, external :: overlap_gauss_r12_ao_with1s double precision, external :: overlap_gauss_r12_ao_with1s
double precision :: factor_ij_1s,beta_ij,center_ij_1s(3)
provide mu_erf final_grid_points j1b_pen provide mu_erf final_grid_points j1b_pen
call wall_time(wall0) call wall_time(wall0)
@ -200,7 +209,7 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
!$OMP PARALLEL DEFAULT (NONE) & !$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & !$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 coef_fit, expo_fit, int_fit, tmp, int_j1b,factor_ij_1s,beta_ij,center_ij_1s) &
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b3_size, & !$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b3_size, &
!$OMP final_grid_points, ng_fit_jast, & !$OMP final_grid_points, ng_fit_jast, &
!$OMP expo_gauss_j_mu_x_2, coef_gauss_j_mu_x_2, & !$OMP expo_gauss_j_mu_x_2, coef_gauss_j_mu_x_2, &
@ -231,6 +240,9 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
expo_fit = expo_gauss_j_mu_x_2(i_fit) expo_fit = expo_gauss_j_mu_x_2(i_fit)
coef_fit = coef_gauss_j_mu_x_2(i_fit) coef_fit = coef_gauss_j_mu_x_2(i_fit)
!DIR$ FORCEINLINE
call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s)
if(dabs(coef_fit*coef*factor_ij_1s*int_j1b).lt.1.d-10)cycle
! --- ! ---
@ -323,17 +335,6 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (3, ao_num, ao_num
expo_fit = expo_gauss_j_mu_1_erf(i_fit) expo_fit = expo_gauss_j_mu_1_erf(i_fit)
coef_fit = coef_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)) & dist = (B_center(1) - r(1)) * (B_center(1) - r(1)) &
+ (B_center(2) - r(2)) * (B_center(2) - r(2)) & + (B_center(2) - r(2)) * (B_center(2) - r(2)) &
+ (B_center(3) - r(3)) * (B_center(3) - r(3)) + (B_center(3) - r(3)) * (B_center(3) - r(3))
@ -347,7 +348,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (3, ao_num, ao_num
expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist
coef_tmp = coef * coef_fit * dexp(-expo_coef_1s) coef_tmp = coef * coef_fit * dexp(-expo_coef_1s)
if(dabs(coef_tmp) .lt. 1d-10) cycle if(dabs(coef_tmp*int_j1b) .lt. 1d-10) cycle
call NAI_pol_x_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r, int_fit) call NAI_pol_x_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r, int_fit)

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

@ -119,7 +119,7 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp_j1b_test, (3, ao_num
double precision :: coef, beta, B_center(3), r(3), ints(3), ints_coulomb(3) double precision :: coef, beta, B_center(3), r(3), ints(3), ints_coulomb(3)
double precision :: tmp_x, tmp_y, tmp_z double precision :: tmp_x, tmp_y, tmp_z
double precision :: wall0, wall1 double precision :: wall0, wall1
double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2,int_j1b double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2,int_j1b,factor_ij_1s,beta_ij,center_ij_1s
dsqpi_3_2 = (dacos(-1.d0))**(3/2) dsqpi_3_2 = (dacos(-1.d0))**(3/2)
call wall_time(wall0) call wall_time(wall0)
@ -128,7 +128,7 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp_j1b_test, (3, ao_num
!$OMP PARALLEL DEFAULT (NONE) & !$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, ints, ints_coulomb, & !$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 int_j1b, tmp_x, tmp_y, tmp_z,factor_ij_1s,beta_ij,center_ij_1s) &
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b2_size, 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 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 x_v_ij_erf_rk_cst_mu_tmp_j1b_test, mu_erf,ao_abs_comb_b2_j1b, &
@ -157,6 +157,10 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp_j1b_test, (3, ao_num
B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i) B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i) B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
! approximate 1 - erf(mu r12) = exp(-2 mu r12^2)
! !DIR$ FORCEINLINE
! call gaussian_product(expo_good_j_mu_1gauss,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s)
! if(dabs(coef * factor_ij_1s*int_j1b).lt.1.d-10)cycle
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, mu_erf, r, ints )
call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, ints_coulomb) call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, ints_coulomb)

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

@ -5,13 +5,13 @@
integer :: i_1s,i,j,ipoint integer :: i_1s,i,j,ipoint
double precision :: coef,beta,center(3),int_j1b,thr double precision :: coef,beta,center(3),int_j1b,thr
double precision :: r(3),weight,dist double precision :: r(3),weight,dist
thr = 1.d-12 thr = 1.d-15
List_comb_thr_b2_size = 0 List_comb_thr_b2_size = 0
do i = 1, ao_num do i = 1, ao_num
do j = i, ao_num do j = i, ao_num
do i_1s = 1, List_all_comb_b2_size do i_1s = 1, List_all_comb_b2_size
coef = List_all_comb_b2_coef (i_1s) coef = List_all_comb_b2_coef (i_1s)
if(dabs(coef).lt.1.d-12)cycle if(dabs(coef).lt.1.d-15)cycle
beta = List_all_comb_b2_expo (i_1s) beta = List_all_comb_b2_expo (i_1s)
beta = max(beta,1.d-12) beta = max(beta,1.d-12)
center(1:3) = List_all_comb_b2_cent(1:3,i_1s) center(1:3) = List_all_comb_b2_cent(1:3,i_1s)
@ -51,7 +51,7 @@ END_PROVIDER
integer :: i_1s,i,j,ipoint,icount integer :: i_1s,i,j,ipoint,icount
double precision :: coef,beta,center(3),int_j1b,thr double precision :: coef,beta,center(3),int_j1b,thr
double precision :: r(3),weight,dist double precision :: r(3),weight,dist
thr = 1.d-12 thr = 1.d-15
ao_abs_comb_b2_j1b = 10000000.d0 ao_abs_comb_b2_j1b = 10000000.d0
do i = 1, ao_num do i = 1, ao_num
do j = i, ao_num do j = i, ao_num
@ -100,7 +100,7 @@ END_PROVIDER
integer :: i_1s,i,j,ipoint integer :: i_1s,i,j,ipoint
double precision :: coef,beta,center(3),int_j1b,thr double precision :: coef,beta,center(3),int_j1b,thr
double precision :: r(3),weight,dist double precision :: r(3),weight,dist
thr = 1.d-12 thr = 1.d-15
List_comb_thr_b3_size = 0 List_comb_thr_b3_size = 0
do i = 1, ao_num do i = 1, ao_num
do j = 1, ao_num do j = 1, ao_num
@ -146,7 +146,7 @@ END_PROVIDER
integer :: i_1s,i,j,ipoint,icount integer :: i_1s,i,j,ipoint,icount
double precision :: coef,beta,center(3),int_j1b,thr double precision :: coef,beta,center(3),int_j1b,thr
double precision :: r(3),weight,dist double precision :: r(3),weight,dist
thr = 1.d-12 thr = 1.d-15
ao_abs_comb_b3_j1b = 10000000.d0 ao_abs_comb_b3_j1b = 10000000.d0
do i = 1, ao_num do i = 1, ao_num
do j = 1, ao_num do j = 1, ao_num
@ -164,7 +164,7 @@ END_PROVIDER
dist = ( center(1) - r(1) )*( center(1) - r(1) ) dist = ( center(1) - r(1) )*( center(1) - r(1) )
dist += ( center(2) - r(2) )*( center(2) - r(2) ) dist += ( center(2) - r(2) )*( center(2) - r(2) )
dist += ( center(3) - r(3) )*( center(3) - r(3) ) 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 int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight
enddo enddo
if(dabs(coef)*dabs(int_j1b).gt.thr)then if(dabs(coef)*dabs(int_j1b).gt.thr)then
icount += 1 icount += 1

27
src/non_h_ints_mu/grad_squared_manu.irp.f
View File

@ -9,8 +9,11 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test, (ao_num, ao_num, ao_nu
implicit none implicit none
integer :: ipoint, i, j, k, l integer :: ipoint, i, j, k, l
double precision :: weight1, ao_ik_r, ao_i_r double precision :: weight1, ao_ik_r, ao_i_r,contrib,contrib2
double precision, allocatable :: ac_mat(:,:,:,:), bc_mat(:,:,:,:) double precision, allocatable :: ac_mat(:,:,:,:), bc_mat(:,:,:,:)
double precision :: wall1, wall0
provide u12sq_j1bsq_test u12_grad1_u12_j1b_grad1_j1b_test grad12_j12_test
call wall_time(wall0)
allocate(ac_mat(ao_num,ao_num,ao_num,ao_num)) allocate(ac_mat(ao_num,ao_num,ao_num,ao_num))
ac_mat = 0.d0 ac_mat = 0.d0
@ -20,16 +23,18 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test, (ao_num, ao_num, ao_nu
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
weight1 = final_weight_at_r_vector(ipoint) weight1 = final_weight_at_r_vector(ipoint)
do i = 1, ao_num do j = 1, ao_num
ao_i_r = weight1 * aos_in_r_array_transp(ipoint,i) do l = 1, ao_num
contrib = u12sq_j1bsq_test(l,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b_test(l,j,ipoint)
contrib2=grad12_j12_test(l,j,ipoint)
do i = 1, ao_num
ao_i_r = weight1 * aos_in_r_array(i,ipoint)
do k = 1, ao_num
ao_ik_r = ao_i_r * aos_in_r_array(k,ipoint)
do k = 1, ao_num ac_mat(k,i,l,j) += ao_ik_r * contrib
ao_ik_r = ao_i_r * aos_in_r_array_transp(ipoint,k) bc_mat(k,i,l,j) += ao_ik_r * contrib2
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_test(l,j,ipoint)
enddo enddo
enddo enddo
enddo enddo
@ -45,6 +50,8 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test, (ao_num, ao_num, ao_nu
enddo enddo
enddo enddo
enddo enddo
call wall_time(wall1)
print*,'wall time for tc_grad_square_ao_test',wall1 - wall0
deallocate(ac_mat) deallocate(ac_mat)
deallocate(bc_mat) deallocate(bc_mat)

46
src/non_h_ints_mu/new_grad_tc_manu.irp.f
View File

@ -95,34 +95,38 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao_test, (ao_num, ao_num, ao_
double precision :: weight1, contrib_x, contrib_y, contrib_z, tmp_x, tmp_y, tmp_z 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 :: ao_k_r, ao_i_r, ao_i_dx, ao_i_dy, ao_i_dz
double precision, allocatable :: ac_mat(:,:,:,:) double precision, allocatable :: ac_mat(:,:,:,:)
double precision :: wall0, wall1
provide int2_grad1_u12_ao_test
call wall_time(wall0)
allocate(ac_mat(ao_num,ao_num,ao_num,ao_num)) allocate(ac_mat(ao_num,ao_num,ao_num,ao_num))
ac_mat = 0.d0 ac_mat = 0.d0
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
weight1 = 0.5d0 * final_weight_at_r_vector(ipoint) weight1 = 0.5d0 * final_weight_at_r_vector(ipoint)
do j = 1, ao_num
do l = 1, ao_num
contrib_x = int2_grad1_u12_ao_test(1,l,j,ipoint)
contrib_y = int2_grad1_u12_ao_test(2,l,j,ipoint)
contrib_z = int2_grad1_u12_ao_test(3,l,j,ipoint)
do i = 1, ao_num
ao_i_r = weight1 * aos_in_r_array (i,ipoint)
ao_i_dx = weight1 * aos_grad_in_r_array_transp(1,i,ipoint)
ao_i_dy = weight1 * aos_grad_in_r_array_transp(2,i,ipoint)
ao_i_dz = weight1 * aos_grad_in_r_array_transp(3,i,ipoint)
do k = 1, ao_num
ao_k_r = aos_in_r_array(k,ipoint)
tmp_x = ao_k_r * ao_i_dx - ao_i_r * aos_grad_in_r_array_transp(1,k,ipoint)
tmp_y = ao_k_r * ao_i_dy - ao_i_r * aos_grad_in_r_array_transp(2,k,ipoint)
tmp_z = ao_k_r * ao_i_dz - ao_i_r * aos_grad_in_r_array_transp(3,k,ipoint)
do i = 1, ao_num tmp_x *= contrib_x
ao_i_r = weight1 * aos_in_r_array_transp (ipoint,i) tmp_y *= contrib_y
ao_i_dx = weight1 * aos_grad_in_r_array_transp_bis(ipoint,i,1) tmp_z *= contrib_z
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 ac_mat(k,i,l,j) += tmp_x + tmp_y + tmp_z
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 enddo
@ -139,6 +143,8 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao_test, (ao_num, ao_num, ao_
enddo enddo
enddo enddo
call wall_time(wall1)
print*,'wall time for tc_grad_and_lapl_ao_test',wall1 - wall0
deallocate(ac_mat) deallocate(ac_mat)
END_PROVIDER END_PROVIDER

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

@ -11,8 +11,8 @@ program test_ints
my_grid_becke = .True. my_grid_becke = .True.
my_n_pt_r_grid = 30 my_n_pt_r_grid = 30
my_n_pt_a_grid = 50 my_n_pt_a_grid = 50
! my_n_pt_r_grid = 10 ! small grid for quick debug ! my_n_pt_r_grid = 15 ! small grid for quick debug
! my_n_pt_a_grid = 14 ! 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 touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
my_extra_grid_becke = .True. my_extra_grid_becke = .True.
my_n_pt_r_extra_grid = 30 my_n_pt_r_extra_grid = 30
@ -40,10 +40,21 @@ program test_ints
! call test_total_grad_lapl ! call test_total_grad_lapl
! call test_total_grad_square ! call test_total_grad_square
! call test_ao_tc_int_chemist ! call test_ao_tc_int_chemist
call test_grid_points_ao ! call test_grid_points_ao
call test_tc_scf
end end
subroutine test_tc_scf
implicit none
! provide tc_grad_square_ao_test
provide tc_grad_and_lapl_ao_test
! provide int2_u_grad1u_x_j1b2_test
! provide x_v_ij_erf_rk_cst_mu_tmp_j1b_test
! print*,'TC_HF_energy = ',TC_HF_energy
! print*,'grad_non_hermit = ',grad_non_hermit
end
subroutine test_ao_tc_int_chemist subroutine test_ao_tc_int_chemist
implicit none implicit none
provide ao_tc_int_chemist provide ao_tc_int_chemist
@ -309,17 +320,18 @@ end
subroutine routine_int2_grad1u2_grad2u2_j1b2 subroutine routine_int2_grad1u2_grad2u2_j1b2
implicit none implicit none
integer :: i,j,ipoint,k,l integer :: i,j,ipoint,k,l
integer :: ii , jj
double precision :: weight,accu_relat, accu_abs, contrib double precision :: weight,accu_relat, accu_abs, contrib
double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
double precision, allocatable :: ints(:,:,:) double precision, allocatable :: ints(:,:,:)
allocate(ints(ao_num, ao_num, n_points_final_grid)) allocate(ints(ao_num, ao_num, n_points_final_grid))
! do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
! do i = 1, ao_num do i = 1, ao_num
! do j = 1, ao_num do j = 1, ao_num
! read(33,*)ints(j,i,ipoint) read(33,*)ints(j,i,ipoint)
! enddo enddo
! enddo enddo
! enddo enddo
allocate(array(ao_num, ao_num, ao_num, ao_num)) allocate(array(ao_num, ao_num, ao_num, ao_num))
array = 0.d0 array = 0.d0
@ -331,17 +343,17 @@ subroutine routine_int2_grad1u2_grad2u2_j1b2
do l = 1, ao_num do l = 1, ao_num
do i = 1, ao_num do i = 1, ao_num
do j = 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(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) += 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(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 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)).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 ! if(dabs(int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) - int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint)).gt.1.d-6)then
! print*,j,i,ipoint ! 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(j,i,ipoint) , int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint), dabs(int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) - int2_grad1u2_grad2u2_j1b2_test(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)) ! print*,int2_grad1u2_grad2u2_j1b2_test(i,j,ipoint) , int2_grad1u2_grad2u2_j1b2_test(i,j,ipoint), dabs(int2_grad1u2_grad2u2_j1b2_test(i,j,ipoint) - int2_grad1u2_grad2u2_j1b2_test(i,j,ipoint))
! stop ! stop
! endif ! endif
! endif ! endif
@ -350,23 +362,35 @@ subroutine routine_int2_grad1u2_grad2u2_j1b2
enddo enddo
enddo enddo
enddo enddo
double precision :: e_ref, e_new
accu_relat = 0.d0 accu_relat = 0.d0
accu_abs = 0.d0 accu_abs = 0.d0
do k = 1, ao_num e_ref = 0.d0
do l = 1, ao_num e_new = 0.d0
do i = 1, ao_num do ii = 1, elec_alpha_num
do j = 1, ao_num do jj = ii, elec_alpha_num
contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k)) do k = 1, ao_num
accu_abs += contrib do l = 1, ao_num
if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then do i = 1, ao_num
accu_relat += contrib/dabs(array_ref(j,i,l,k)) do j = 1, ao_num
endif e_ref += mo_coef(j,ii) * mo_coef(i,ii) * array_ref(j,i,l,k) * mo_coef(l,jj) * mo_coef(k,jj)
e_new += mo_coef(j,ii) * mo_coef(i,ii) * array(j,i,l,k) * mo_coef(l,jj) * mo_coef(k,jj)
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 enddo
enddo enddo
enddo enddo
print*,'accu_abs = ',accu_abs/dble(ao_num)**4 enddo
print*,'accu_relat = ',accu_relat/dble(ao_num)**4 print*,'e_ref = ',e_ref
print*,'e_new = ',e_new
! print*,'accu_abs = ',accu_abs/dble(ao_num)**4
! print*,'accu_relat = ',accu_relat/dble(ao_num)**4