diff --git a/src/ao_many_one_e_ints/ao_gaus_gauss.irp.f b/src/ao_many_one_e_ints/ao_gaus_gauss.irp.f index cd9a486d..681d1e6f 100644 --- a/src/ao_many_one_e_ints/ao_gaus_gauss.irp.f +++ b/src/ao_many_one_e_ints/ao_gaus_gauss.irp.f @@ -116,6 +116,7 @@ double precision function overlap_gauss_r12_ao(D_center,delta,i,j) if(ao_overlap_abs(j,i).lt.1.d-12)then return endif + ! TODO :: PUT CYCLES IN LOOPS num_A = ao_nucl(i) power_A(1:3)= ao_power(i,1:3) A_center(1:3) = nucl_coord(num_A,1:3) diff --git a/src/ao_tc_eff_map/map_integrals_eff_pot.irp.f b/src/ao_tc_eff_map/map_integrals_eff_pot.irp.f index aea4644f..4a6128b9 100644 --- a/src/ao_tc_eff_map/map_integrals_eff_pot.irp.f +++ b/src/ao_tc_eff_map/map_integrals_eff_pot.irp.f @@ -78,7 +78,7 @@ double precision function get_ao_tc_sym_two_e_pot(i,j,k,l,map) result(result) use map_module implicit none BEGIN_DOC - ! Gets one |AO| two-electron integral from the |AO| map + ! Gets one |AO| two-electron integral from the |AO| map in PHYSICIST NOTATION END_DOC integer, intent(in) :: i,j,k,l integer(key_kind) :: idx diff --git a/src/ao_two_e_ints/map_integrals.irp.f b/src/ao_two_e_ints/map_integrals.irp.f index 55b2d5e2..de4195ba 100644 --- a/src/ao_two_e_ints/map_integrals.irp.f +++ b/src/ao_two_e_ints/map_integrals.irp.f @@ -326,7 +326,9 @@ double precision function get_ao_two_e_integral(i,j,k,l,map) result(result) use map_module implicit none BEGIN_DOC - ! Gets one AO bi-electronic integral from the AO map + ! Gets one AO bi-electronic integral from the AO map in PHYSICIST NOTATION + ! + ! <1:k, 2:l |1:i, 2:j> END_DOC integer, intent(in) :: i,j,k,l integer(key_kind) :: idx diff --git a/src/non_h_ints_mu/fit_j.irp.f b/src/non_h_ints_mu/fit_j.irp.f new file mode 100644 index 00000000..695ead7f --- /dev/null +++ b/src/non_h_ints_mu/fit_j.irp.f @@ -0,0 +1,91 @@ +BEGIN_PROVIDER [ double precision, expo_j_xmu, (n_fit_1_erf_x) ] + implicit none + BEGIN_DOC + ! F(x) = x * (1 - erf(x)) - 1/sqrt(pi) * exp(-x**2) is fitted with a gaussian and a Slater + ! + ! \approx - 1/sqrt(pi) * exp(-alpha * x ) exp(-beta * x**2) + ! + ! where alpha = expo_j_xmu(1) and beta = expo_j_xmu(2) + END_DOC + expo_j_xmu(1) = 1.7477d0 + expo_j_xmu(2) = 0.668662d0 + +END_PROVIDER + + BEGIN_PROVIDER [double precision, expo_gauss_j_mu_x, (n_max_fit_slat)] +&BEGIN_PROVIDER [double precision, coef_gauss_j_mu_x, (n_max_fit_slat)] + implicit none + BEGIN_DOC +! J(mu,r12) = 1/2 r12 * (1 - erf(mu*r12)) - 1/(2 sqrt(pi)*mu) exp(-(mu*r12)^2) is expressed as +! +! J(mu,r12) = 0.5/mu * F(r12*mu) where F(x) = x * (1 - erf(x)) - 1/sqrt(pi) * exp(-x**2) +! +! F(x) is fitted by - 1/sqrt(pi) * exp(-alpha * x) exp(-beta*mu^2x^2) (see expo_j_xmu) +! +! The slater function exp(-alpha * x) is fitted with n_max_fit_slat gaussians +! +! See Appendix 2 of JCP 154, 084119 (2021) +! + END_DOC + integer :: i + double precision :: expos(n_max_fit_slat),alpha,beta + alpha = expo_j_xmu(1) * mu_erf + call expo_fit_slater_gam(alpha,expos) + beta = expo_j_xmu(2) * mu_erf**2.d0 + + do i = 1, n_max_fit_slat + expo_gauss_j_mu_x(i) = expos(i) + beta + coef_gauss_j_mu_x(i) = coef_fit_slat_gauss(i) / (2.d0 * mu_erf) * (- 1/dsqrt(dacos(-1.d0))) + enddo +END_PROVIDER + +double precision function F_x_j(x) + implicit none + BEGIN_DOC + ! F_x_j(x) = dimension-less correlation factor = x (1 - erf(x)) - 1/sqrt(pi) exp(-x^2) + END_DOC + double precision, intent(in) :: x + F_x_j = x * (1.d0 - derf(x)) - 1/dsqrt(dacos(-1.d0)) * dexp(-x**2) + +end + +double precision function j_mu_F_x_j(x) + implicit none + BEGIN_DOC + ! j_mu_F_x_j(x) = correlation factor = 1/2 r12 * (1 - erf(mu*r12)) - 1/(2 sqrt(pi)*mu) exp(-(mu*r12)^2) + ! + ! = 1/(2*mu) * F_x_j(mu*x) + END_DOC + double precision :: F_x_j + double precision, intent(in) :: x + j_mu_F_x_j = 0.5d0/mu_erf * F_x_j(x*mu_erf) +end + +double precision function j_mu(x) + implicit none + double precision, intent(in) :: x + BEGIN_DOC + ! j_mu(x) = correlation factor = 1/2 r12 * (1 - erf(mu*r12)) - 1/(2 sqrt(pi)*mu) exp(-(mu*r12)^2) + END_DOC + j_mu = 0.5d0* x * (1.d0 - derf(mu_erf*x)) - 0.5d0/( dsqrt(dacos(-1.d0))*mu_erf) * dexp(-(mu_erf*x)*(mu_erf*x)) + +end + +double precision function j_mu_fit_gauss(x) + implicit none + BEGIN_DOC + ! j_mu_fit_gauss(x) = correlation factor = 1/2 r12 * (1 - erf(mu*r12)) - 1/(2 sqrt(pi)*mu) exp(-(mu*r12)^2) + ! + ! but fitted with gaussians + END_DOC + double precision, intent(in) :: x + integer :: i + double precision :: alpha,coef + j_mu_fit_gauss = 0.d0 + do i = 1, n_max_fit_slat + alpha = expo_gauss_j_mu_x(i) + coef = coef_gauss_j_mu_x(i) + j_mu_fit_gauss += coef_gauss_j_mu_x(i) * dexp(-expo_gauss_j_mu_x(i)*x*x) + enddo + +end diff --git a/src/non_h_ints_mu/grad_squared.irp.f b/src/non_h_ints_mu/grad_squared.irp.f new file mode 100644 index 00000000..a88521a1 --- /dev/null +++ b/src/non_h_ints_mu/grad_squared.irp.f @@ -0,0 +1,72 @@ + BEGIN_PROVIDER [ double precision, grad_1_squared_u_ij_mu, ( ao_num, ao_num,n_points_final_grid)] + implicit none + integer :: ipoint,i,j,m,igauss + BEGIN_DOC + ! grad_1_squared_u_ij_mu(j,i,ipoint) = -1/2 \int dr2 phi_j(r2) phi_i(r2) |\grad_r1 u(r1,r2,\mu)|^2 + ! |\grad_r1 u(r1,r2,\mu)|^2 = 1/4 * (1 - erf(mu*r12))^2 + ! ! (1 - erf(mu*r12))^2 = \sum_i coef_gauss_1_erf_x_2(i) * exp(-expo_gauss_1_erf_x_2(i) * r12^2) + END_DOC + double precision :: r(3),delta,coef + double precision :: overlap_gauss_r12_ao,time0,time1 + print*,'providing grad_1_squared_u_ij_mu ...' + call wall_time(time0) + !TODO : strong optmization : write the loops in a different way + ! : for each couple of AO, the gaussian product are done once for all + 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 + ! \int dr2 phi_j(r2) phi_i(r2) (1 - erf(mu*r12))^2 + ! = \sum_i coef_gauss_1_erf_x_2(i) \int dr2 phi_j(r2) phi_i(r2) exp(-expo_gauss_1_erf_x_2(i) * (r_1 - r_2)^2) + do igauss = 1, n_max_fit_slat + delta = expo_gauss_1_erf_x_2(igauss) + coef = coef_gauss_1_erf_x_2(igauss) + grad_1_squared_u_ij_mu(j,i,ipoint) += -0.25 * coef * overlap_gauss_r12_ao(r,delta,i,j) + enddo + enddo + enddo + enddo + call wall_time(time1) + print*,'Wall time for grad_1_squared_u_ij_mu = ',time1 - time0 + END_PROVIDER + +BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao_num)] + implicit none + BEGIN_DOC + ! tc_grad_square_ao(k,i,l,j) = -1/2 + ! + END_DOC + integer :: ipoint,i,j,k,l + double precision :: contrib,weight1 + double precision, allocatable :: ac_mat(:,:,:,:) + allocate(ac_mat(ao_num, ao_num, ao_num, ao_num)) + ac_mat = 0.d0 + do ipoint = 1, n_points_final_grid + weight1 = final_weight_at_r_vector(ipoint) + do j = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do k = 1, ao_num + contrib = weight1 *0.5D0* (aos_in_r_array_transp(ipoint,k) * aos_in_r_array_transp(ipoint,i)) + ! \int dr1 phi_k(r1) phi_i(r1) . \int dr2 |\grad_1 u(r1,r2)|^2 \phi_l(r2) \phi_j(r2) + ac_mat(k,i,l,j) += grad_1_squared_u_ij_mu(l,j,ipoint) * contrib + 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(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i) + enddo + enddo + enddo + enddo + +END_PROVIDER + diff --git a/src/non_h_ints_mu/grad_tc_int.irp.f b/src/non_h_ints_mu/grad_tc_int.irp.f index dd60e724..40600335 100644 --- a/src/non_h_ints_mu/grad_tc_int.irp.f +++ b/src/non_h_ints_mu/grad_tc_int.irp.f @@ -11,6 +11,8 @@ END_DOC double precision, allocatable :: b_mat(:,:,:,:),ac_mat(:,:,:,:) ! provide v_ij_erf_rk_cst_mu provide v_ij_erf_rk_cst_mu x_v_ij_erf_rk_cst_mu +! ao_non_hermit_term_chemist = non_h_ints +! return call wall_time(wall0) allocate(b_mat(n_points_final_grid,ao_num,ao_num,3),ac_mat(ao_num, ao_num, ao_num, ao_num)) !$OMP PARALLEL & @@ -35,6 +37,9 @@ END_DOC !$OMP END DO !$OMP END PARALLEL + + ! (A) b_mat(ipoint,k,i,m) X v_ij_erf_rk_cst_mu(j,l,r1) + ! 1/2 \int dr1 x1 phi_k(1) d/dx1 phi_i(1) \int dr2 (1 - erf(mu_r12))/r12 phi_j(2) phi_l(2) ac_mat = 0.d0 do m = 1, 3 ! A B^T dim(A,1) dim(B,2) dim(A,2) alpha * A LDA @@ -60,6 +65,8 @@ END_DOC !$OMP END DO !$OMP END PARALLEL + ! (B) b_mat(ipoint,k,i,m) X x_v_ij_erf_rk_cst_mu(j,l,r1,m) + ! 1/2 \int dr1 phi_k(1) d/dx1 phi_i(1) \int dr2 x2(1 - erf(mu_r12))/r12 phi_j(2) phi_l(2) do m = 1, 3 ! A B^T dim(A,1) dim(B,2) dim(A,2) alpha * A LDA call dgemm("N","N",ao_num*ao_num,ao_num*ao_num,n_points_final_grid,-1.d0,x_v_ij_erf_rk_cst_mu(1,1,1,m),ao_num*ao_num & @@ -75,6 +82,7 @@ END_DOC do l = 1, ao_num do i = 1, ao_num do k = 1, ao_num + ! (ki|lj) (ki|lj) (lj|ki) ao_non_hermit_term_chemist(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i) enddo enddo diff --git a/src/non_h_ints_mu/new_grad_tc.irp.f b/src/non_h_ints_mu/new_grad_tc.irp.f new file mode 100644 index 00000000..068381b4 --- /dev/null +++ b/src/non_h_ints_mu/new_grad_tc.irp.f @@ -0,0 +1,70 @@ +BEGIN_PROVIDER [ double precision, grad_1_u_ij_mu, ( ao_num, ao_num,n_points_final_grid,3)] + implicit none + BEGIN_DOC + ! grad_1_u_ij_mu(i,j,ipoint) = -1 * \int dr2 \grad_r1 u(r1,r2) \phi_i(r2) \phi_j(r2) + ! + ! where r1 = r(ipoint) + ! + ! grad_1_u_ij_mu(i,j,ipoint) = \int dr2 (r1 - r2) (erf(mu * r12)-1)/2 r_12 \phi_i(r2) \phi_j(r2) + END_DOC + integer :: ipoint,i,j,m + double precision :: r(3) + do m = 1, 3 + 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 + grad_1_u_ij_mu(i,j,ipoint,m) = v_ij_erf_rk_cst_mu(i,j,ipoint) * r(m) - x_v_ij_erf_rk_cst_mu(i,j,ipoint,m) + enddo + enddo + enddo + enddo + grad_1_u_ij_mu *= 0.5d0 + +END_PROVIDER + +BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num, ao_num)] + implicit none + BEGIN_DOC + ! tc_grad_and_lapl_ao(k,i,l,j) = + ! + ! = 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 + integer :: ipoint,i,j,k,l,m + double precision :: contrib,weight1 + double precision, allocatable :: ac_mat(:,:,:,:) + allocate(ac_mat(ao_num, ao_num, ao_num, ao_num)) + ac_mat = 0.d0 + do m = 1, 3 + do ipoint = 1, n_points_final_grid + weight1 = final_weight_at_r_vector(ipoint) + do j = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do k = 1, ao_num + contrib = weight1 *0.5D0* (aos_in_r_array_transp(ipoint,k) * aos_grad_in_r_array_transp_bis(ipoint,i,m) & + -aos_in_r_array_transp(ipoint,i) * aos_grad_in_r_array_transp_bis(ipoint,k,m) ) + ! \int dr1 phi_k(r1) \grad_r1 phi_i(r1) . \int dr2 \grad_r1 u(r1,r2) \phi_l(r2) \phi_j(r2) + ac_mat(k,i,l,j) += grad_1_u_ij_mu(l,j,ipoint,m) * contrib + enddo + 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(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i) + enddo + enddo + enddo + enddo + +END_PROVIDER diff --git a/src/non_h_ints_mu/test_non_h_ints.irp.f b/src/non_h_ints_mu/test_non_h_ints.irp.f new file mode 100644 index 00000000..c535d0c5 --- /dev/null +++ b/src/non_h_ints_mu/test_non_h_ints.irp.f @@ -0,0 +1,102 @@ +program test_non_h + implicit none + my_grid_becke = .True. + my_n_pt_r_grid = 50 + my_n_pt_a_grid = 74 +! my_n_pt_r_grid = 10 ! small grid for quick debug +! my_n_pt_a_grid = 26 ! small grid for quick debug + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid +!call routine_grad_squared + call routine_fit +end + +subroutine routine_lapl_grad + implicit none + integer :: i,j,k,l + double precision :: grad_lapl, get_ao_tc_sym_two_e_pot,new,accu,contrib + double precision :: ao_two_e_integral_erf,get_ao_two_e_integral,count_n,accu_relat +! !!!!!!!!!!!!!!!!!!!!! WARNING +! THIS ROUTINE MAKES SENSE ONLY IF HAND MODIFIED coef_gauss_eff_pot(1:n_max_fit_slat) = 0. to cancel (1-erf(mu*r12))^2 + accu = 0.d0 + accu_relat = 0.d0 + count_n = 0.d0 + do i = 1, ao_num + do j = 1, ao_num + do k = 1, ao_num + do l = 1, ao_num + grad_lapl = get_ao_tc_sym_two_e_pot(i,j,k,l,ao_tc_sym_two_e_pot_map) ! pure gaussian part : comes from Lapl + grad_lapl += ao_two_e_integral_erf(i, k, j, l) ! erf(mu r12)/r12 : comes from Lapl + grad_lapl += ao_non_hermit_term_chemist(k,i,l,j) ! \grad u(r12) . grad + new = tc_grad_and_lapl_ao(k,i,l,j) + new += get_ao_two_e_integral(i,j,k,l,ao_integrals_map) + contrib = dabs(new - grad_lapl) + if(dabs(grad_lapl).gt.1.d-12)then + count_n += 1.d0 + accu_relat += 2.0d0 * contrib/dabs(grad_lapl+new) + endif + if(contrib.gt.1.d-10)then + print*,i,j,k,l + print*,grad_lapl,new,contrib + print*,2.0d0*contrib/dabs(grad_lapl+new+1.d-12) + endif + accu += contrib + enddo + enddo + enddo + enddo + print*,'accu = ',accu/count_n + print*,'accu/rel = ',accu_relat/count_n + +end + +subroutine routine_grad_squared + implicit none + integer :: i,j,k,l + double precision :: grad_squared, get_ao_tc_sym_two_e_pot,new,accu,contrib + double precision :: count_n,accu_relat +! !!!!!!!!!!!!!!!!!!!!! WARNING +! THIS ROUTINE MAKES SENSE ONLY IF HAND MODIFIED coef_gauss_eff_pot(n_max_fit_slat:n_max_fit_slat+1) = 0. to cancel exp(-'mu*r12)^2) + accu = 0.d0 + accu_relat = 0.d0 + count_n = 0.d0 + do i = 1, ao_num + do j = 1, ao_num + do k = 1, ao_num + do l = 1, ao_num + grad_squared = get_ao_tc_sym_two_e_pot(i,j,k,l,ao_tc_sym_two_e_pot_map) ! pure gaussian part : comes from Lapl + new = tc_grad_square_ao(k,i,l,j) + contrib = dabs(new - grad_squared) + if(dabs(grad_squared).gt.1.d-12)then + count_n += 1.d0 + accu_relat += 2.0d0 * contrib/dabs(grad_squared+new) + endif + if(contrib.gt.1.d-10)then + print*,i,j,k,l + print*,grad_squared,new,contrib + print*,2.0d0*contrib/dabs(grad_squared+new+1.d-12) + endif + accu += contrib + enddo + enddo + enddo + enddo + print*,'accu = ',accu/count_n + print*,'accu/rel = ',accu_relat/count_n + +end + +subroutine routine_fit + implicit none + integer :: i,nx + double precision :: dx,xmax,x,j_mu,j_mu_F_x_j,j_mu_fit_gauss + nx = 500 + xmax = 5.d0 + dx = xmax/dble(nx) + x = 0.d0 + print*,'coucou',mu_erf + do i = 1, nx + write(33,'(100(F16.10,X))') x,j_mu(x),j_mu_F_x_j(x),j_mu_fit_gauss(x) + x += dx + enddo + +end