LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-11-09 06:53:38 +01:00
Code Releases Activity

Merge branch 'eginer-dev' into good-dev-tc

Browse Source
This commit is contained in:
AbdAmmar 2022-12-22 14:53:15 +01:00
commit 6c0fb305b9
25 changed files with 2023 additions and 324 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

47
src/ao_many_one_e_ints/ao_gaus_gauss.irp.f
View File

@ -156,6 +156,53 @@ end function overlap_gauss_r12_ao
! --
double precision function overlap_abs_gauss_r12_ao(D_center, delta, i, j)
BEGIN_DOC
! \int dr AO_i(r) AO_j(r) e^{-delta |r-D_center|^2}
END_DOC
implicit none
integer, intent(in) :: i, j
double precision, intent(in) :: D_center(3), delta
integer :: power_A(3), power_B(3), l, k
double precision :: A_center(3), B_center(3), alpha, beta, coef, coef1, analytical_j
double precision, external :: overlap_abs_gauss_r12
overlap_abs_gauss_r12_ao = 0.d0
if(ao_overlap_abs(j,i).lt.1.d-12) then
return
endif
power_A(1:3) = ao_power(i,1:3)
power_B(1:3) = ao_power(j,1:3)
A_center(1:3) = nucl_coord(ao_nucl(i),1:3)
B_center(1:3) = nucl_coord(ao_nucl(j),1:3)
do l = 1, ao_prim_num(i)
alpha = ao_expo_ordered_transp (l,i)
coef1 = ao_coef_normalized_ordered_transp(l,i)
do k = 1, ao_prim_num(j)
beta = ao_expo_ordered_transp(k,j)
coef = coef1 * ao_coef_normalized_ordered_transp(k,j)
if(dabs(coef) .lt. 1d-12) cycle
analytical_j = overlap_abs_gauss_r12(D_center, delta, A_center, B_center, power_A, power_B, alpha, beta)
overlap_abs_gauss_r12_ao += dabs(coef * analytical_j)
enddo
enddo
end function overlap_gauss_r12_ao
! --
subroutine overlap_gauss_r12_ao_v(D_center, LD_D, delta, i, j, resv, LD_resv, n_points)
BEGIN_DOC

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

@ -1,4 +1,396 @@
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
double precision :: int_gauss,dsqpi_3_2,int_j1b
double precision :: factor_ij_1s,beta_ij,center_ij_1s(3),sq_pi_3_2
sq_pi_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(:,:,:) = 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_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, &
!$OMP List_comb_thr_b3_cent, int2_grad1u2_grad2u2_j1b2_test, ao_abs_comb_b3_j1b,&
!$OMP ao_overlap_abs,sq_pi_3_2)
!$OMP DO SCHEDULE(dynamic)
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
if(ao_overlap_abs(j,i) .lt. 1.d-12) then
cycle
endif
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)
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_1_erf_x_2(i_fit)
!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
! if(dabs(coef_fit*factor_ij_1s*int_j1b).lt.1.d-10)cycle ! old version
if(dabs(coef_fit*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-3/2)).lt.1.d-10)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)
int_gauss = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
int2_grad1u2_grad2u2_j1b2_test(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(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
BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_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(:),big_array(:,:,:)
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
big_array(:,:,:) = 0.d0
allocate(big_array(n_points_final_grid,ao_num, 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_j1b) &
!$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, &
!$OMP List_comb_thr_b3_cent, big_array,&
!$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
if(ao_overlap_abs(j,i) .lt. 1.d-12) then
cycle
endif
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-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)
do i_fit = 1, ng_fit_jast
expo_fit = expo_gauss_1_erf_x_2(i_fit)
coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit) * coef
call overlap_gauss_r12_ao_with1s_v(B_center, beta, final_grid_points_transp, size(final_grid_points_transp,1),&
expo_fit, i, j, int_fit_v, size(int_fit_v,1),n_points_final_grid)
do ipoint = 1, n_points_final_grid
big_array(ipoint,j,i) += coef_fit * int_fit_v(ipoint)
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
deallocate(int_fit_v)
!$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 i = 2, ao_num
do j = 1, i-1
int2_grad1u2_grad2u2_j1b2_test_v(j,i,ipoint) = big_array(ipoint,i,j)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for int2_grad1u2_grad2u2_j1b2_test_v', wall1 - wall0
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
double precision :: factor_ij_1s,beta_ij,center_ij_1s(3),sq_pi_3_2
sq_pi_3_2 = (dacos(-1.d0))**(3/2)
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,factor_ij_1s,beta_ij,center_ij_1s) &
!$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,sq_pi_3_2, &
!$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)
!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 ! old version
if(dabs(coef_fit*coef*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-3/2)).lt.1.d-10)cycle
! ---
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, sq_pi_3_2,sq_alpha
sq_pi_3_2 = dacos(-1.D0)**(3/2)
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,sq_alpha) &
!$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,sq_pi_3_2)
!$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)
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)
sq_alpha = alpha_1s_inv * dsqrt(alpha_1s_inv)
! if(dabs(coef_tmp*int_j1b) .lt. 1d-10) cycle ! old version
if(dabs(coef_tmp*int_j1b*sq_pi_3_2*sq_alpha) .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
@ -30,7 +422,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
!$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 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, &
@ -46,7 +438,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
r(3) = final_grid_points(3,ipoint)
tmp = 0.d0
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)
@ -63,7 +455,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
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
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
@ -104,181 +496,3 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
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, 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, &
!$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 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
if(ao_overlap_abs(j,i) .lt. 1.d-12) then
cycle
endif
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-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)
do i_fit = 1, ng_fit_jast
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)
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, 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, &
!$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
if(ao_overlap_abs(j,i) .lt. 1.d-12) then
cycle
endif
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-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)
do i_fit = 1, ng_fit_jast
expo_fit = expo_gauss_1_erf_x_2(i_fit)
coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit) * coef
call overlap_gauss_r12_ao_with1s_v(B_center, beta, final_grid_points_transp, size(final_grid_points_transp,1),&
expo_fit, i, j, int_fit_v, size(int_fit_v,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)
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

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)

132
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)
@ -50,7 +50,7 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num,
B_center(1) = List_comb_thr_b2_cent(1,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)
! TODO :: cycle on the 1 - erf(mur12)
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)
@ -94,9 +94,9 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_nu
do ipoint = 1, n_points_final_grid
do i = 1, ao_num
do j = 1, ao_num
x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,1) = x_v_ij_erf_rk_cst_mu_tmp_j1b(1,j,i,ipoint)
x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,2) = x_v_ij_erf_rk_cst_mu_tmp_j1b(2,j,i,ipoint)
x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,3) = x_v_ij_erf_rk_cst_mu_tmp_j1b(3,j,i,ipoint)
x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,1) = x_v_ij_erf_rk_cst_mu_tmp_j1b_test(1,j,i,ipoint)
x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,2) = x_v_ij_erf_rk_cst_mu_tmp_j1b_test(2,j,i,ipoint)
x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,3) = x_v_ij_erf_rk_cst_mu_tmp_j1b_test(3,j,i,ipoint)
enddo
enddo
enddo
@ -119,22 +119,23 @@ 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 :: tmp_x, tmp_y, tmp_z
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)
provide expo_erfc_mu_gauss ao_prod_sigma ao_prod_center
call wall_time(wall0)
x_v_ij_erf_rk_cst_mu_tmp_j1b_test = 0.d0
!$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 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 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)
!$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma)
! !$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,expo_erfc_mu_gauss)
!$OMP DO
!do ipoint = 1, 10
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
r(2) = final_grid_points(2,ipoint)
@ -142,12 +143,12 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp_j1b_test, (3, ao_num
do i = 1, ao_num
do j = i, ao_num
if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle
if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-10)cycle
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)
@ -157,6 +158,14 @@ 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(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
! if(ao_prod_center(1,j,i).ne.10000.d0)then
! ! approximate 1 - erf(mu r12) by a gaussian * 10
! !DIR$ FORCEINLINE
! call gaussian_product(expo_erfc_mu_gauss,r, &
! ao_prod_sigma(j,i),ao_prod_center(1,j,i), &
! factor_ij_1s,beta_ij,center_ij_1s)
! if(dabs(coef * factor_ij_1s*int_j1b*10.d0 * dsqpi_3_2 * beta_ij**(-3/2)).lt.1.d-10)cycle
! endif
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)
@ -223,7 +232,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 +247,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)
@ -285,3 +294,96 @@ END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2) u(mu, r12) with u(mu,r12) \approx 1/2 mu e^{-2.5 * mu (r12)^2}
!
END_DOC
implicit none
integer :: i, j, ipoint, i_1s
double precision :: r(3), int_fit, expo_fit, coef_fit
double precision :: coef, beta, B_center(3)
double precision :: tmp
double precision :: wall0, wall1
double precision, external :: overlap_gauss_r12_ao_with1s
double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2,int_j1b
dsqpi_3_2 = (dacos(-1.d0))**(3/2)
provide mu_erf final_grid_points j1b_pen
call wall_time(wall0)
v_ij_u_cst_mu_j1b_ng_1_test = 0.d0
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, &
!$OMP beta_ij_u, factor_ij_1s_u, center_ij_1s_u, &
!$OMP coef_fit, expo_fit, int_fit, tmp,coeftot,int_j1b) &
!$OMP SHARED (n_points_final_grid, ao_num, &
!$OMP final_grid_points, expo_good_j_mu_1gauss,coef_good_j_mu_1gauss, &
!$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, &
!$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_ng_1_test,ao_abs_comb_b2_j1b, &
!$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2)
!$OMP DO
!do ipoint = 1, 10
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
if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle
tmp = 0.d0
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)
int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i)
if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle
B_center(1) = List_comb_thr_b2_cent(1,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)
! do i_fit = 1, ng_fit_jast
expo_fit = expo_good_j_mu_1gauss
coef_fit = 1.d0
coeftot = coef * coef_fit
if(dabs(coeftot).lt.1.d-15)cycle
double precision :: beta_ij_u, factor_ij_1s_u, center_ij_1s_u(3),coeftot
call gaussian_product(beta,B_center,expo_fit,r,factor_ij_1s_u,beta_ij_u,center_ij_1s_u)
if(factor_ij_1s_u*ao_overlap_abs_grid(j,i).lt.1.d-15)cycle
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
tmp += coef * coef_fit * int_fit
! enddo
enddo
v_ij_u_cst_mu_j1b_ng_1_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
v_ij_u_cst_mu_j1b_ng_1_test(j,i,ipoint) = v_ij_u_cst_mu_j1b_ng_1_test(i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for v_ij_u_cst_mu_j1b_ng_1_test', wall1 - wall0
END_PROVIDER
! ---

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

59
src/ao_many_one_e_ints/list_grid.irp.f Normal file
View File

@ -0,0 +1,59 @@
BEGIN_PROVIDER [ integer, n_pts_grid_ao_prod, (ao_num, ao_num)]
&BEGIN_PROVIDER [ integer, max_n_pts_grid_ao_prod]
implicit none
integer :: i,j,ipoint
double precision :: overlap, r(3),thr, overlap_abs_gauss_r12_ao,overlap_gauss_r12_ao
double precision :: sigma,dist,center_ij(3),fact_gauss, alpha, center(3)
n_pts_grid_ao_prod = 0
thr = 1.d-11
print*,' expo_good_j_mu_1gauss = ',expo_good_j_mu_1gauss
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, r, overlap, thr,fact_gauss, alpha, center,dist,sigma,center_ij) &
!$OMP SHARED (n_points_final_grid, ao_num, ao_overlap_abs_grid,n_pts_grid_ao_prod,expo_good_j_mu_1gauss,&
!$OMP final_grid_points,ao_prod_center,ao_prod_sigma,ao_nucl)
!$OMP DO
do i = 1, ao_num
! do i = 3,3
do j = 1, ao_num
! do i = 22,22
! do j = 9,9
center_ij(1:3) = ao_prod_center(1:3,j,i)
sigma = ao_prod_sigma(j,i)
sigma *= sigma
sigma = 0.5d0 /sigma
! if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-10)cycle
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)
dist = (center_ij(1) - r(1))*(center_ij(1) - r(1))
dist += (center_ij(2) - r(2))*(center_ij(2) - r(2))
dist += (center_ij(3) - r(3))*(center_ij(3) - r(3))
dist = dsqrt(dist)
call gaussian_product(sigma, center_ij, expo_good_j_mu_1gauss, r, fact_gauss, alpha, center)
! print*,''
! print*,j,i,ao_overlap_abs_grid(j,i),ao_overlap_abs(j,i)
! print*,r
! print*,dist,sigma
! print*,fact_gauss
if( fact_gauss*ao_overlap_abs_grid(j,i).lt.1.d-11)cycle
if(ao_nucl(i) == ao_nucl(j))then
overlap = overlap_abs_gauss_r12_ao(r, expo_good_j_mu_1gauss, i, j)
else
overlap = overlap_gauss_r12_ao(r, expo_good_j_mu_1gauss, i, j)
endif
! print*,overlap
if(dabs(overlap).lt.thr)cycle
n_pts_grid_ao_prod(j,i) += 1
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
integer :: list(ao_num)
do i = 1, ao_num
list(i) = maxval(n_pts_grid_ao_prod(:,i))
enddo
max_n_pts_grid_ao_prod = maxval(list)
END_PROVIDER

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

@ -1,74 +1,74 @@
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-15
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-15)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
int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_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-15
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) )
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
if(dabs(coef)*dabs(int_j1b).gt.thr)then
icount += 1
@ -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-15
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
int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_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-15
ao_abs_comb_b3_j1b = 10000000.d0
do i = 1, ao_num
do j = 1, ao_num
@ -154,17 +154,17 @@ 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) )
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
if(dabs(coef)*dabs(int_j1b).gt.thr)then
icount += 1
@ -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)

74
src/ao_many_one_e_ints/prim_int_gauss_gauss.irp.f
View File

@ -32,16 +32,17 @@ double precision function overlap_gauss_r12(D_center, delta, A_center, B_center,
dim1 = 100
thr = 1.d-10
d(:) = 0 ! order of the polynom for the gaussian exp(-delta (r - D)^2 ) == 0
overlap_gauss_r12 = 0.d0
! New gaussian/polynom defined by :: new pol new center new expo cst fact new order
call give_explicit_poly_and_gaussian( A_new, A_center_new , alpha_new, fact_a_new, iorder_a_new &
, delta, alpha, d, power_A, D_center, A_center, n_pt_max_integrals)
call give_explicit_poly_and_gaussian(A_new , A_center_new , alpha_new, fact_a_new , iorder_a_new ,&
delta,alpha,d,power_A,D_center,A_center,n_pt_max_integrals)
if(fact_a_new.lt.thr)return
! The new gaussian exp(-delta (r - D)^2 ) (x-A_x)^a \exp(-\alpha (x-A_x)^2
accu = 0.d0
do lx = 0, iorder_a_new(1)
coefx = A_new(lx,1)
if(dabs(coefx) .lt. thr) cycle
coefx = A_new(lx,1)*fact_a_new
if(dabs(coefx).lt.thr)cycle
iorder_tmp(1) = lx
do ly = 0, iorder_a_new(2)
@ -63,9 +64,70 @@ double precision function overlap_gauss_r12(D_center, delta, A_center, B_center,
enddo
enddo
enddo
overlap_gauss_r12 = accu
end
overlap_gauss_r12 = fact_a_new * accu
!---
double precision function overlap_abs_gauss_r12(D_center,delta,A_center,B_center,power_A,power_B,alpha,beta)
BEGIN_DOC
! Computes the following integral :
!
! .. math ::
!
! \int dr exp(-delta (r - D)^2 ) |(x-A_x)^a (x-B_x)^b \exp(-\alpha (x-A_x)^2 - \beta (x-B_x)^2 )|
!
END_DOC
implicit none
include 'constants.include.F'
double precision, intent(in) :: D_center(3), delta ! pure gaussian "D"
double precision, intent(in) :: A_center(3),B_center(3),alpha,beta ! gaussian/polynoms "A" and "B"
integer, intent(in) :: power_A(3),power_B(3)
double precision :: overlap_x,overlap_y,overlap_z,overlap
! First you multiply the usual gaussian "A" with the gaussian exp(-delta (r - D)^2 )
double precision :: A_new(0:max_dim,3)! new polynom
double precision :: A_center_new(3) ! new center
integer :: iorder_a_new(3) ! i_order(i) = order of the new polynom ==> should be equal to power_A
double precision :: alpha_new ! new exponent
double precision :: fact_a_new ! constant factor
double precision :: accu,coefx,coefy,coefz,coefxy,coefxyz,thr,dx,lower_exp_val
integer :: d(3),i,lx,ly,lz,iorder_tmp(3),dim1
dim1=50
lower_exp_val = 40.d0
thr = 1.d-12
d(:) = 0 ! order of the polynom for the gaussian exp(-delta (r - D)^2 ) == 0
overlap_abs_gauss_r12 = 0.d0
! New gaussian/polynom defined by :: new pol new center new expo cst fact new order
call give_explicit_poly_and_gaussian(A_new , A_center_new , alpha_new, fact_a_new , iorder_a_new ,&
delta,alpha,d,power_A,D_center,A_center,n_pt_max_integrals)
if(fact_a_new.lt.thr)return
! The new gaussian exp(-delta (r - D)^2 ) (x-A_x)^a \exp(-\alpha (x-A_x)^2
accu = 0.d0
do lx = 0, iorder_a_new(1)
coefx = A_new(lx,1)*fact_a_new
! if(dabs(coefx).lt.thr)cycle
iorder_tmp(1) = lx
do ly = 0, iorder_a_new(2)
coefy = A_new(ly,2)
coefxy = coefx * coefy
if(dabs(coefxy).lt.thr)cycle
iorder_tmp(2) = ly
do lz = 0, iorder_a_new(3)
coefz = A_new(lz,3)
coefxyz = coefxy * coefz
if(dabs(coefxyz).lt.thr)cycle
iorder_tmp(3) = lz
call overlap_x_abs(A_center_new(1),B_center(1),alpha_new,beta,iorder_tmp(1),power_B(1),overlap_x,lower_exp_val,dx,dim1)
call overlap_x_abs(A_center_new(2),B_center(2),alpha_new,beta,iorder_tmp(2),power_B(2),overlap_y,lower_exp_val,dx,dim1)
call overlap_x_abs(A_center_new(3),B_center(3),alpha_new,beta,iorder_tmp(3),power_B(3),overlap_z,lower_exp_val,dx,dim1)
accu += dabs(coefxyz * overlap_x * overlap_y * overlap_z)
enddo
enddo
enddo
overlap_abs_gauss_r12= accu
>>>>>>> f7e58e4a636af0ab066aa644a74ab56cb4de6041
end
!---

35
src/ao_tc_eff_map/fit_j.irp.f
View File

@ -1,5 +1,40 @@
BEGIN_PROVIDER [ double precision, expo_j_xmu_1gauss ]
&BEGIN_PROVIDER [ double precision, coef_j_xmu_1gauss ]
implicit none
BEGIN_DOC
! Upper bound long range fit of F(x) = x * (1 - erf(x)) - 1/sqrt(pi) * exp(-x**2)
!
! with a single gaussian.
!
! Such a function can be used to screen integrals with F(x).
END_DOC
expo_j_xmu_1gauss = 0.5d0
coef_j_xmu_1gauss = 1.d0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, expo_erfc_gauss ]
implicit none
expo_erfc_gauss = 1.41211d0
END_PROVIDER
BEGIN_PROVIDER [ double precision, expo_erfc_mu_gauss ]
implicit none
expo_erfc_mu_gauss = expo_erfc_gauss * mu_erf * mu_erf
END_PROVIDER
BEGIN_PROVIDER [ double precision, expo_good_j_mu_1gauss ]
&BEGIN_PROVIDER [ double precision, coef_good_j_mu_1gauss ]
implicit none
BEGIN_DOC
! exponent of Gaussian in order to obtain an upper bound of J(r12,mu)
!
! Can be used to scree integrals with J(r12,mu)
END_DOC
expo_good_j_mu_1gauss = 2.D0 * mu_erf * expo_j_xmu_1gauss
coef_good_j_mu_1gauss = 0.5d0/mu_erf * coef_j_xmu_1gauss
END_PROVIDER
BEGIN_PROVIDER [ double precision, expo_j_xmu, (n_fit_1_erf_x) ]
implicit none
BEGIN_DOC

34
src/bi_ort_ints/semi_num_ints_mo.irp.f
View File

@ -108,15 +108,27 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_transp, (ao_num, ao_num, 3,
double precision :: wall0, wall1
call wall_time(wall0)
do ipoint = 1, n_points_final_grid
do i = 1, ao_num
do j = 1, ao_num
int2_grad1_u12_ao_transp(j,i,1,ipoint) = int2_grad1_u12_ao(1,j,i,ipoint)
int2_grad1_u12_ao_transp(j,i,2,ipoint) = int2_grad1_u12_ao(2,j,i,ipoint)
int2_grad1_u12_ao_transp(j,i,3,ipoint) = int2_grad1_u12_ao(3,j,i,ipoint)
enddo
enddo
enddo
if(test_cycle_tc)then
do ipoint = 1, n_points_final_grid
do i = 1, ao_num
do j = 1, ao_num
int2_grad1_u12_ao_transp(j,i,1,ipoint) = int2_grad1_u12_ao_test(1,j,i,ipoint)
int2_grad1_u12_ao_transp(j,i,2,ipoint) = int2_grad1_u12_ao_test(2,j,i,ipoint)
int2_grad1_u12_ao_transp(j,i,3,ipoint) = int2_grad1_u12_ao_test(3,j,i,ipoint)
enddo
enddo
enddo
else
do ipoint = 1, n_points_final_grid
do i = 1, ao_num
do j = 1, ao_num
int2_grad1_u12_ao_transp(j,i,1,ipoint) = int2_grad1_u12_ao(1,j,i,ipoint)
int2_grad1_u12_ao_transp(j,i,2,ipoint) = int2_grad1_u12_ao(2,j,i,ipoint)
int2_grad1_u12_ao_transp(j,i,3,ipoint) = int2_grad1_u12_ao(3,j,i,ipoint)
enddo
enddo
enddo
endif
call wall_time(wall1)
print *, ' wall time for int2_grad1_u12_ao_transp ', wall1 - wall0
@ -178,8 +190,8 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_t, (n_points_final_grid, 3,
integer :: i, j, ipoint
do ipoint = 1, n_points_final_grid
do i = 1, mo_num
do j = 1, mo_num
do i = 1, ao_num
do j = 1, ao_num
int2_grad1_u12_ao_t(ipoint,1,j,i) = int2_grad1_u12_ao(1,j,i,ipoint)
int2_grad1_u12_ao_t(ipoint,2,j,i) = int2_grad1_u12_ao(2,j,i,ipoint)
int2_grad1_u12_ao_t(ipoint,3,j,i) = int2_grad1_u12_ao(3,j,i,ipoint)

1
src/determinants/filter_connected.irp.f
View File

@ -96,7 +96,6 @@ subroutine filter_not_connected(key1,key2,Nint,sze,idx)
idx(0) = l-1
end
subroutine filter_connected(key1,key2,Nint,sze,idx)
use bitmasks
implicit none

164
src/determinants/sparse_mat.irp.f Normal file
View File

@ -0,0 +1,164 @@
use bitmasks
subroutine filter_connected_array(key1,key2,ld,Nint,sze,idx)
use bitmasks
implicit none
BEGIN_DOC
! Filters out the determinants that are not connected by H
!
! returns the array idx which contains the index of the
!
! determinants in the array key1 that interact
!
! via the H operator with key2.
!
! idx(0) is the number of determinants that interact with key1
END_DOC
integer, intent(in) :: Nint, ld,sze
integer(bit_kind), intent(in) :: key1(Nint,2,ld)
integer(bit_kind), intent(in) :: key2(Nint,2)
integer, intent(out) :: idx(0:sze)
integer :: i,j,l
integer :: degree_x2
ASSERT (Nint > 0)
ASSERT (sze >= 0)
l=1
if (Nint==1) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt( xor( key1(1,1,i), key2(1,1))) &
+ popcnt( xor( key1(1,2,i), key2(1,2)))
! print*,degree_x2
if (degree_x2 > 4) then
cycle
else
idx(l) = i
l = l+1
endif
enddo
else if (Nint==2) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(2,1,i), key2(2,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))) + &
popcnt(xor( key1(2,2,i), key2(2,2)))
if (degree_x2 > 4) then
cycle
else
idx(l) = i
l = l+1
endif
enddo
else if (Nint==3) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))) + &
popcnt(xor( key1(2,1,i), key2(2,1))) + &
popcnt(xor( key1(2,2,i), key2(2,2))) + &
popcnt(xor( key1(3,1,i), key2(3,1))) + &
popcnt(xor( key1(3,2,i), key2(3,2)))
if (degree_x2 > 4) then
cycle
else
idx(l) = i
l = l+1
endif
enddo
else
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = 0
!DIR$ LOOP COUNT MIN(4)
do j=1,Nint
degree_x2 = degree_x2+ popcnt(xor( key1(j,1,i), key2(j,1))) +&
popcnt(xor( key1(j,2,i), key2(j,2)))
if (degree_x2 > 4) then
exit
endif
enddo
if (degree_x2 <= 5) then
idx(l) = i
l = l+1
endif
enddo
endif
idx(0) = l-1
! print*,'idx(0) = ',idx(0)
end
BEGIN_PROVIDER [ integer, n_sparse_mat]
&BEGIN_PROVIDER [ integer, n_connected_per_det, (N_det)]
&BEGIN_PROVIDER [ integer, n_max_connected_per_det]
implicit none
BEGIN_DOC
! n_sparse_mat = total number of connections in the CI matrix
!
! n_connected_per_det(i) = number of connected determinants to the determinant psi_det(1,1,i)
!
! n_max_connected_per_det = maximum number of connected determinants
END_DOC
integer, allocatable :: idx(:)
allocate(idx(0:N_det))
integer :: i
n_sparse_mat = 0
do i = 1, N_det
call filter_connected_array(psi_det_sorted,psi_det_sorted(1,1,i),psi_det_size,N_int,N_det,idx)
n_connected_per_det(i) = idx(0)
n_sparse_mat += idx(0)
enddo
n_max_connected_per_det = maxval(n_connected_per_det)
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), connected_det_per_det, (N_int,2,n_max_connected_per_det,N_det)]
&BEGIN_PROVIDER [ integer(bit_kind), list_connected_det_per_det, (n_max_connected_per_det,N_det)]
implicit none
BEGIN_DOC
! connected_det_per_det(:,:,j,i) = jth connected determinant to the determinant psi_det(:,:,i)
!
! list_connected_det_per_det(j,i) = index of jth determinant in psi_det which is connected to psi_det(:,:,i)
END_DOC
integer, allocatable :: idx(:)
allocate(idx(0:N_det))
integer :: i,j
do i = 1, N_det
call filter_connected_array(psi_det_sorted,psi_det_sorted(1,1,i),psi_det_size,N_int,N_det,idx)
do j = 1, idx(0)
connected_det_per_det(1:N_int,1:2,j,i) = psi_det_sorted(1:N_int,1:2,idx(j))
list_connected_det_per_det(j,i) = idx(j)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, sparse_h_mat, (n_max_connected_per_det, N_det)]
implicit none
BEGIN_DOC
! sparse matrix format
!
! sparse_h_mat(j,i) = matrix element between the jth connected determinant and psi_det(:,:,i)
END_DOC
integer :: i,j
double precision :: hij
do i = 1, N_det
do j = 1, n_connected_per_det(i)
call i_H_j(psi_det(1,1,i),connected_det_per_det(1,1,j,i),N_int,hij)
sparse_h_mat(j,i) = hij
enddo
enddo
END_PROVIDER

41
src/dft_utils_in_r/ao_in_r.irp.f
View File

@ -40,6 +40,47 @@
END_PROVIDER
BEGIN_PROVIDER[double precision, aos_in_r_array_extra, (ao_num,n_points_extra_final_grid)]
implicit none
BEGIN_DOC
! aos_in_r_array_extra(i,j) = value of the ith ao on the jth grid point
END_DOC
integer :: i,j
double precision :: aos_array(ao_num), r(3)
!$OMP PARALLEL DO &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,r,aos_array,j) &
!$OMP SHARED(aos_in_r_array_extra,n_points_extra_final_grid,ao_num,final_grid_points_extra)
do i = 1, n_points_extra_final_grid
r(1) = final_grid_points_extra(1,i)
r(2) = final_grid_points_extra(2,i)
r(3) = final_grid_points_extra(3,i)
call give_all_aos_at_r(r,aos_array)
do j = 1, ao_num
aos_in_r_array_extra(j,i) = aos_array(j)
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
BEGIN_PROVIDER[double precision, aos_in_r_array_extra_transp, (n_points_extra_final_grid,ao_num)]
implicit none
BEGIN_DOC
! aos_in_r_array_extra_transp(i,j) = value of the jth ao on the ith grid point
END_DOC
integer :: i,j
double precision :: aos_array(ao_num), r(3)
do i = 1, n_points_extra_final_grid
do j = 1, ao_num
aos_in_r_array_extra_transp(i,j) = aos_in_r_array_extra(j,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, aos_grad_in_r_array, (ao_num,n_points_final_grid,3)]
implicit none

65
src/dft_utils_in_r/ao_prod_mlti_pl.irp.f
View File

@ -1,5 +1,28 @@
BEGIN_PROVIDER [ double precision, ao_abs_int_grid, (ao_num)]
implicit none
BEGIN_DOC
! ao_abs_int_grid(i) = \int dr |phi_i(r) |
END_DOC
integer :: i,j,ipoint
double precision :: contrib, weight,r(3)
ao_abs_int_grid = 0.D0
do ipoint = 1,n_points_final_grid
r(:) = final_grid_points(:,ipoint)
weight = final_weight_at_r_vector(ipoint)
do i = 1, ao_num
contrib = dabs(aos_in_r_array(i,ipoint)) * weight
ao_abs_int_grid(i) += contrib
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_overlap_abs_grid, (ao_num, ao_num)]
implicit none
BEGIN_DOC
! ao_overlap_abs_grid(j,i) = \int dr |phi_i(r) phi_j(r)|
END_DOC
integer :: i,j,ipoint
double precision :: contrib, weight,r(3)
ao_overlap_abs_grid = 0.D0
@ -21,7 +44,7 @@ BEGIN_PROVIDER [ double precision, ao_prod_center, (3, ao_num, ao_num)]
BEGIN_DOC
! ao_prod_center(1:3,j,i) = \int dr |phi_i(r) phi_j(r)| x/y/z / \int |phi_i(r) phi_j(r)|
!
! if \int |phi_i(r) phi_j(r)| < 1.d-15 then ao_prod_center = 0.
! if \int |phi_i(r) phi_j(r)| < 1.d-10 then ao_prod_center = 10000.
END_DOC
integer :: i,j,m,ipoint
double precision :: contrib, weight,r(3)
@ -44,20 +67,23 @@ BEGIN_PROVIDER [ double precision, ao_prod_center, (3, ao_num, ao_num)]
do m = 1, 3
ao_prod_center(m,j,i) *= 1.d0/ao_overlap_abs_grid(j,i)
enddo
else
do m = 1, 3
ao_prod_center(m,j,i) = 10000.d0
enddo
endif
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_prod_sigma, (ao_num, ao_num)]
BEGIN_PROVIDER [ double precision, ao_prod_abs_r, (ao_num, ao_num)]
implicit none
BEGIN_DOC
! ao_prod_sigma(i,j) = \int |phi_i(r) phi_j(r)| dsqrt((x - <|i|x|j|>)^2 + (y - <|i|y|j|>)^2 +(z - <|i|z|j|>)^2) / \int |phi_i(r) phi_j(r)|
! ao_prod_abs_r(i,j) = \int |phi_i(r) phi_j(r)| dsqrt((x - <|i|x|j|>)^2 + (y - <|i|y|j|>)^2 +(z - <|i|z|j|>)^2) / \int |phi_i(r) phi_j(r)|
!
! gives you a precise idea of the spatial extension of the distribution phi_i(r) phi_j(r)
END_DOC
ao_prod_sigma = 0.d0
ao_prod_abs_r = 0.d0
integer :: i,j,m,ipoint
double precision :: contrib, weight,r(3),contrib_x2
do ipoint = 1,n_points_final_grid
@ -71,21 +97,34 @@ BEGIN_PROVIDER [ double precision, ao_prod_sigma, (ao_num, ao_num)]
contrib_x2 += (r(m) - ao_prod_center(m,j,i)) * (r(m) - ao_prod_center(m,j,i))
enddo
contrib_x2 = dsqrt(contrib_x2)
ao_prod_sigma(j,i) += contrib * contrib_x2
ao_prod_abs_r(j,i) += contrib * contrib_x2
enddo
enddo
enddo
do i = 1, ao_num
do j = 1, ao_num
if(dabs(ao_overlap_abs_grid(j,i)).gt.1.d-10)then
ao_prod_sigma(j,i) *= 1.d0/ao_overlap_abs_grid(j,i)
endif
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, ao_prod_sigma, (ao_num, ao_num)]
implicit none
BEGIN_DOC
! Gaussian exponent reproducing the product |chi_i(r) chi_j(r)|
!
! Therefore |chi_i(r) chi_j(r)| \approx e^{-ao_prod_sigma(j,i) (r - ao_prod_center(1:3,j,i))**2}
END_DOC
integer :: i,j
double precision :: pi,alpha
pi = dacos(-1.d0)
do i = 1, ao_num
do j = 1, ao_num
! if(dabs(ao_overlap_abs_grid(j,i)).gt.1.d-5)then
alpha = 1.d0/pi * (2.d0*ao_overlap_abs_grid(j,i)/ao_prod_abs_r(j,i))**2
ao_prod_sigma(j,i) = alpha
! endif
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_prod_dist_grid, (ao_num, ao_num, n_points_final_grid)]
implicit none
BEGIN_DOC

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) ]

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

@ -0,0 +1,194 @@
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,contrib,contrib2
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))
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 j = 1, ao_num
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)
ac_mat(k,i,l,j) += ao_ik_r * contrib
bc_mat(k,i,l,j) += ao_ik_r * contrib2
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
call wall_time(wall1)
print*,'wall time for tc_grad_square_ao_test',wall1 - wall0
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
provide int2_u_grad1u_x_j1b2_test
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_test(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
! ---
BEGIN_PROVIDER [ double precision, grad12_j12_test, (ao_num, ao_num, n_points_final_grid) ]
implicit none
integer :: ipoint, i, j, m, igauss
double precision :: r(3), delta, coef
double precision :: tmp1
double precision :: time0, time1
double precision, external :: overlap_gauss_r12_ao
provide int2_grad1u2_grad2u2_j1b2_test
print*, ' providing grad12_j12_test ...'
call wall_time(time0)
PROVIDE j1b_type
if(j1b_type .eq. 3) then
do ipoint = 1, n_points_final_grid
tmp1 = v_1b(ipoint)
tmp1 = tmp1 * tmp1
do j = 1, ao_num
do i = 1, ao_num
grad12_j12_test(i,j,ipoint) = tmp1 * int2_grad1u2_grad2u2_j1b2_test(i,j,ipoint)
enddo
enddo
enddo
else
grad12_j12_test = 0.d0
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 j = 1, ao_num
do i = 1, ao_num
do igauss = 1, n_max_fit_slat
delta = expo_gauss_1_erf_x_2(igauss)
coef = coef_gauss_1_erf_x_2(igauss)
grad12_j12_test(i,j,ipoint) += -0.25d0 * coef * overlap_gauss_r12_ao(r, delta, i, j)
enddo
enddo
enddo
enddo
endif
call wall_time(time1)
print*, ' Wall time for grad12_j12_test = ', time1 - time0
END_PROVIDER
! ---

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

@ -253,3 +253,4 @@ END_PROVIDER
! ---

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

@ -0,0 +1,152 @@
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(:,:,:,:)
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))
ac_mat = 0.d0
do ipoint = 1, n_points_final_grid
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)
tmp_x *= contrib_x
tmp_y *= contrib_y
tmp_z *= contrib_z
ac_mat(k,i,l,j) += tmp_x + tmp_y + tmp_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
call wall_time(wall1)
print*,'wall time for tc_grad_and_lapl_ao_test',wall1 - wall0
deallocate(ac_mat)
END_PROVIDER
! ---

45
src/non_h_ints_mu/total_tc_int.irp.f
View File

@ -9,15 +9,19 @@ BEGIN_PROVIDER [double precision, ao_tc_int_chemist, (ao_num, ao_num, ao_num, ao
call wall_time(wall0)
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
ao_tc_int_chemist(k,i,l,j) = tc_grad_square_ao(k,i,l,j) + tc_grad_and_lapl_ao(k,i,l,j) + ao_two_e_coul(k,i,l,j)
enddo
enddo
enddo
enddo
if(test_cycle_tc)then
ao_tc_int_chemist = ao_tc_int_chemist_test
else
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
ao_tc_int_chemist(k,i,l,j) = tc_grad_square_ao(k,i,l,j) + tc_grad_and_lapl_ao(k,i,l,j) + ao_two_e_coul(k,i,l,j)
enddo
enddo
enddo
enddo
endif
call wall_time(wall1)
print *, ' wall time for ao_tc_int_chemist ', wall1 - wall0
@ -26,6 +30,29 @@ END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, ao_tc_int_chemist_test, (ao_num, ao_num, ao_num, ao_num)]
implicit none
integer :: i, j, k, l
double precision :: wall1, wall0
call wall_time(wall0)
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
ao_tc_int_chemist_test(k,i,l,j) = tc_grad_square_ao_test(k,i,l,j) + tc_grad_and_lapl_ao_test(k,i,l,j) + ao_two_e_coul(k,i,l,j)
enddo
enddo
enddo
enddo
call wall_time(wall1)
print *, ' wall time for ao_tc_int_chemist_test ', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, ao_two_e_coul, (ao_num, ao_num, ao_num, ao_num) ]
BEGIN_DOC

82
src/tc_bi_ortho/test_tc_fock.irp.f
View File

@ -15,7 +15,8 @@ program test_tc_fock
!call routine_2
! call routine_3()
call test_3e
! call test_3e
call routine_tot
end
! ---
@ -84,8 +85,8 @@ subroutine routine_3()
print*, i, a
stop
endif
!print*, ' excited det'
!call debug_det(det_i, N_int)
print*, ' excited det'
call debug_det(det_i, N_int)
call htilde_mu_mat_bi_ortho(det_i, ref_bitmask, N_int, hmono, htwoe, hthree, htilde_ij)
if(dabs(hthree).lt.1.d-10)cycle
@ -116,3 +117,78 @@ subroutine routine_3()
end subroutine routine_3
! ---
subroutine routine_tot()
use bitmasks ! you need to include the bitmasks_module.f90 features
implicit none
integer :: i, a, i_ok, s1,other_spin(2)
double precision :: hmono, htwoe, hthree, htilde_ij
double precision :: err_ai, err_tot, ref, new
integer(bit_kind), allocatable :: det_i(:,:)
allocate(det_i(N_int,2))
other_spin(1) = 2
other_spin(2) = 1
err_tot = 0.d0
! do s1 = 1, 2
s1 = 2
det_i = ref_bitmask
call debug_det(det_i, N_int)
print*, ' HF det'
call debug_det(det_i, N_int)
! do i = 1, elec_num_tab(s1)
! do a = elec_num_tab(s1)+1, mo_num ! virtual
do i = 1, elec_beta_num
do a = elec_beta_num+1, elec_alpha_num! virtual
! do i = elec_beta_num+1, elec_alpha_num
! do a = elec_alpha_num+1, mo_num! virtual
print*,i,a
det_i = ref_bitmask
call do_single_excitation(det_i, i, a, s1, i_ok)
if(i_ok == -1) then
print*, 'PB !!'
print*, i, a
stop
endif
call htilde_mu_mat_bi_ortho(det_i, ref_bitmask, N_int, hmono, htwoe, hthree, htilde_ij)
print*,htilde_ij
if(dabs(htilde_ij).lt.1.d-10)cycle
print*, ' excited det'
call debug_det(det_i, N_int)
if(s1 == 1)then
new = Fock_matrix_tc_mo_alpha(a,i)
else
new = Fock_matrix_tc_mo_beta(a,i)
endif
ref = htilde_ij
! if(s1 == 1)then
! new = fock_a_tot_3e_bi_orth(a,i)
! else if(s1 == 2)then
! new = fock_b_tot_3e_bi_orth(a,i)
! endif
err_ai = dabs(dabs(ref) - dabs(new))
if(err_ai .gt. 1d-7) then
print*,'s1 = ',s1
print*, ' warning on', i, a
print*, ref,new,err_ai
endif
print*, ref,new,err_ai
err_tot += err_ai
write(22, *) htilde_ij
enddo
enddo
! enddo
print *, ' err_tot = ', err_tot
deallocate(det_i)
end subroutine routine_3

6
src/tc_keywords/EZFIO.cfg
View File

@ -166,3 +166,9 @@ doc: Thresholds on the Imag part of energy
interface: ezfio,provider,ocaml
default: 1.e-7
[test_cycle_tc]
type: logical
doc: If |true|, the integrals of the three-body jastrow are computed with cycles
interface: ezfio,provider,ocaml
default: False

23
src/tc_scf/fock_tc_mo_tot.irp.f
View File

@ -73,6 +73,29 @@
+ (Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
enddo
enddo
if(three_body_h_tc)then
! C-O
do j = 1, elec_beta_num
do i = elec_beta_num+1, elec_alpha_num
Fock_matrix_tc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j))
Fock_matrix_tc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i))
enddo
enddo
! C-V
do j = 1, elec_beta_num
do i = elec_alpha_num+1, mo_num
Fock_matrix_tc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j))
Fock_matrix_tc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i))
enddo
enddo
! O-V
do j = elec_beta_num+1, elec_alpha_num
do i = elec_alpha_num+1, mo_num
Fock_matrix_tc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j))
Fock_matrix_tc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i))
enddo
enddo
endif
endif

2
src/tc_scf/fock_three.irp.f
View File

@ -128,6 +128,8 @@ BEGIN_PROVIDER [double precision, diag_three_elem_hf]
call give_abb_contrib(integral_abb)
call give_bbb_contrib(integral_bbb)
diag_three_elem_hf = integral_aaa + integral_aab + integral_abb + integral_bbb
! print*,'integral_aaa + integral_aab + integral_abb + integral_bbb'
! print*,integral_aaa , integral_aab , integral_abb , integral_bbb
endif

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

@ -6,24 +6,44 @@ program test_ints
implicit none
print *, 'starting ...'
print *, ' starting 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 = 15 ! small grid for quick debug
! my_n_pt_a_grid = 26 ! small grid for quick debug
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
!call routine_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_extra_grid_becke = .True.
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
!call routine_int2_grad1u2_grad2u2_j1b2
!! 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
! call test_ao_tc_int_chemist
! call test_grid_points_ao
! call test_tc_scf
call test_int_gauss
!call test_fock_3e_uhf_ao()
call test_fock_3e_uhf_mo()
@ -32,6 +52,32 @@ end
! ---
subroutine test_tc_scf
implicit none
integer :: i
! provide int2_u_grad1u_x_j1b2_test
provide x_v_ij_erf_rk_cst_mu_tmp_j1b_test
! do i = 1, ng_fit_jast
! print*,expo_gauss_1_erf_x_2(i),coef_gauss_1_erf_x_2(i)
! enddo
! 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
implicit none
provide ao_tc_int_chemist
! provide ao_tc_int_chemist_test
! provide tc_grad_square_ao_test
! provide tc_grad_and_lapl_ao_test
end
! ---
subroutine routine_test_j1b
implicit none
integer :: i,icount,j
@ -49,7 +95,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)
@ -58,7 +104,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
@ -228,7 +274,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
@ -289,6 +335,7 @@ end
subroutine routine_int2_grad1u2_grad2u2_j1b2
implicit none
integer :: i,j,ipoint,k,l
integer :: ii , jj
double precision :: weight,accu_relat, accu_abs, contrib
double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
double precision, allocatable :: ints(:,:,:)
@ -311,20 +358,90 @@ subroutine routine_int2_grad1u2_grad2u2_j1b2
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(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) += 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
! 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(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(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(i,j,ipoint), dabs(int2_grad1u2_grad2u2_j1b2_test(i,j,ipoint) - int2_grad1u2_grad2u2_j1b2_test(i,j,ipoint))
! stop
! endif
! endif
enddo
enddo
enddo
enddo
enddo
double precision :: e_ref, e_new
accu_relat = 0.d0
accu_abs = 0.d0
e_ref = 0.d0
e_new = 0.d0
do ii = 1, elec_alpha_num
do jj = ii, elec_alpha_num
do k = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do j = 1, ao_num
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
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
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
@ -350,6 +467,110 @@ subroutine routine_int2_grad1u2_grad2u2_j1b2
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
! ---
@ -491,7 +712,139 @@ end subroutine test_fock_3e_uhf_mo()
! ---
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
subroutine test_grid_points_ao
implicit none
integer :: i,j,ipoint,icount,icount_good, icount_bad,icount_full
double precision :: thr
thr = 1.d-10
! print*,'max_n_pts_grid_ao_prod = ',max_n_pts_grid_ao_prod
! print*,'n_pts_grid_ao_prod'
do i = 1, ao_num
do j = i, ao_num
icount = 0
icount_good = 0
icount_bad = 0
icount_full = 0
do ipoint = 1, n_points_final_grid
! if(dabs(int2_u_grad1u_x_j1b2_test(1,j,i,ipoint)) &
! + dabs(int2_u_grad1u_x_j1b2_test(2,j,i,ipoint)) &
! + dabs(int2_u_grad1u_x_j1b2_test(2,j,i,ipoint)) )
! if(dabs(int2_u2_j1b2_test(j,i,ipoint)).gt.thr)then
! icount += 1
! endif
if(dabs(v_ij_u_cst_mu_j1b_ng_1_test(j,i,ipoint)).gt.thr*0.1d0)then
icount_full += 1
endif
if(dabs(v_ij_u_cst_mu_j1b_test(j,i,ipoint)).gt.thr)then
icount += 1
if(dabs(v_ij_u_cst_mu_j1b_ng_1_test(j,i,ipoint)).gt.thr*0.1d0)then
icount_good += 1
else
print*,j,i,ipoint
print*,dabs(v_ij_u_cst_mu_j1b_test(j,i,ipoint)),dabs(v_ij_u_cst_mu_j1b_ng_1_test(j,i,ipoint)),dabs(v_ij_u_cst_mu_j1b_ng_1_test(j,i,ipoint))/dabs(v_ij_u_cst_mu_j1b_test(j,i,ipoint))
icount_bad += 1
endif
endif
! if(dabs(v_ij_u_cst_mu_j1b_ng_1_test(j,i,ipoint)).gt.thr)then
! endif
enddo
print*,''
print*,j,i
print*,icount,icount_full, icount_bad!,n_pts_grid_ao_prod(j,i)
print*,dble(icount)/dble(n_points_final_grid),dble(icount_full)/dble(n_points_final_grid)
! dble(n_pts_grid_ao_prod(j,i))/dble(n_points_final_grid)
! if(icount.gt.n_pts_grid_ao_prod(j,i))then
! print*,'pb !!'
! endif
enddo
enddo
end
subroutine test_int_gauss
implicit none
integer :: i,j
print*,'center'
do i = 1, ao_num
do j = i, ao_num
print*,j,i
print*,ao_prod_sigma(j,i),ao_overlap_abs_grid(j,i)
print*,ao_prod_center(1:3,j,i)
enddo
enddo
print*,''
double precision :: weight, r(3),integral_1,pi,center(3),f_r,alpha,distance,integral_2
center = 0.d0
pi = dacos(-1.d0)
integral_1 = 0.d0
integral_2 = 0.d0
alpha = 0.75d0
do i = 1, n_points_final_grid
! you get x, y and z of the ith grid point
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
weight = final_weight_at_r_vector(i)
distance = dsqrt( (r(1) - center(1))**2 + (r(2) - center(2))**2 + (r(3) - center(3))**2 )
f_r = dexp(-alpha * distance*distance)
! you add the contribution of the grid point to the integral
integral_1 += f_r * weight
integral_2 += f_r * distance * weight
enddo
print*,'integral_1 =',integral_1
print*,'(pi/alpha)**1.5 =',(pi / alpha)**1.5
print*,'integral_2 =',integral_2
print*,'(pi/alpha)**1.5 =',2.d0*pi / (alpha)**2
end

90
src/tools/print_hmat.irp.f Normal file
View File

@ -0,0 +1,90 @@
program print_h_mat
implicit none
BEGIN_DOC
! program that prints out the CI matrix in sparse form
END_DOC
read_wf = .True.
touch read_wf
call print_wf_dets
call print_wf_coef
call sparse_mat
call full_mat
call test_sparse_mat
end
subroutine print_wf_dets
implicit none
integer :: i,j
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
output=trim(ezfio_filename)//'.wf_det'
i_unit_output = getUnitAndOpen(output,'w')
write(i_unit_output,*)N_det,N_int
do i = 1, N_det
write(i_unit_output,*)psi_det_sorted(1:N_int,1,i)
write(i_unit_output,*)psi_det_sorted(1:N_int,2,i)
enddo
end
subroutine print_wf_coef
implicit none
integer :: i,j
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
output=trim(ezfio_filename)//'.wf_coef'
i_unit_output = getUnitAndOpen(output,'w')
write(i_unit_output,*)N_det,N_states
do i = 1, N_det
write(i_unit_output,*)psi_coef_sorted(i,1:N_states)
enddo
end
subroutine sparse_mat
implicit none
integer :: i,j
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
output=trim(ezfio_filename)//'.hmat_sparse'
i_unit_output = getUnitAndOpen(output,'w')
do i = 1, N_det
write(i_unit_output,*)i,n_connected_per_det(i)
do j =1, n_connected_per_det(i)
write(i_unit_output,*)list_connected_det_per_det(j,i),sparse_h_mat(j,i)
enddo
enddo
end
subroutine full_mat
implicit none
integer :: i,j
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
output=trim(ezfio_filename)//'.hmat_full'
i_unit_output = getUnitAndOpen(output,'w')
do i = 1, N_det
do j = i, N_det
write(i_unit_output,*)i,j,H_matrix_all_dets(j,i)
enddo
enddo
end
subroutine test_sparse_mat
implicit none
integer :: i,j
double precision, allocatable :: eigvec(:,:), eigval(:), hmat(:,:)
allocate(eigval(N_det), eigvec(N_det,N_det),hmat(N_det,N_det))
hmat = 0.d0
do i = 1, N_det
do j =1, n_connected_per_det(i)
hmat(list_connected_det_per_det(j,i),i) = sparse_h_mat(j,i)
enddo
enddo
call lapack_diag(eigval,eigvec,hmat,N_det,N_det)
print*,'The two energies should be the same '
print*,'eigval(1) = ',eigval(1)
print*,'psi_energy= ',CI_electronic_energy(1)
end