From d32b170fd33d99e6a642cc540a5cba7e2e0abc54 Mon Sep 17 00:00:00 2001 From: eginer Date: Sat, 22 Oct 2022 16:21:37 +0200 Subject: [PATCH] documentation --- src/ao_many_one_e_ints/ao_erf_gauss.irp.f | 2 +- .../prim_int_erf_gauss.irp.f | 2 +- src/non_h_ints_mu/README.rst | 4 + src/non_h_ints_mu/fit_j.irp.f | 2 +- src/non_h_ints_mu/grad_squared.irp.f | 97 +++++++++++++++++++ src/non_h_ints_mu/new_grad_tc.irp.f | 2 +- src/non_h_ints_mu/test_non_h_ints.irp.f | 43 +++++++- 7 files changed, 143 insertions(+), 9 deletions(-) diff --git a/src/ao_many_one_e_ints/ao_erf_gauss.irp.f b/src/ao_many_one_e_ints/ao_erf_gauss.irp.f index 39be352f..0c5430d5 100644 --- a/src/ao_many_one_e_ints/ao_erf_gauss.irp.f +++ b/src/ao_many_one_e_ints/ao_erf_gauss.irp.f @@ -90,7 +90,7 @@ end subroutine erfc_mu_gauss_xyz_ij_ao(i,j,mu, C_center, delta,gauss_ints) implicit none BEGIN_DOC - ! gauss_ints(m) = \int dr exp(-delta (r - C)^2 ) x/y/z * ( 1 - erf(mu |r-r'|))/ |r-r'| * AO_i(r') * AO_j(r') + ! gauss_ints(m) = \int dr exp(-delta (r - C)^2 ) x/y/z * ( 1 - erf(mu |r-C|))/ |r-C| * AO_i(r) * AO_j(r) ! ! with m = 1 ==> x, m = 2, m = 3 ==> z ! diff --git a/src/ao_many_one_e_ints/prim_int_erf_gauss.irp.f b/src/ao_many_one_e_ints/prim_int_erf_gauss.irp.f index 641d25fe..4c2b65a6 100644 --- a/src/ao_many_one_e_ints/prim_int_erf_gauss.irp.f +++ b/src/ao_many_one_e_ints/prim_int_erf_gauss.irp.f @@ -142,7 +142,7 @@ double precision function erf_mu_gauss(D_center,delta,mu,A_center,B_center,power ! ! .. math:: ! - ! \int dr exp(-delta (r - D)^2 ) erf(mu*|r-r'|)/ |r-r'| * (x-A_x)^a (x-B_x)^b \exp(-\alpha (x-A_x)^2 - \beta (x-B_x)^2 ) + ! \int dr exp(-delta (r - D)^2 ) erf(mu*|r-D|)/ |r-D| * (x-A_x)^a (x-B_x)^b \exp(-\alpha (x-A_x)^2 - \beta (x-B_x)^2 ) ! END_DOC diff --git a/src/non_h_ints_mu/README.rst b/src/non_h_ints_mu/README.rst index 6a36bb98..fb1c25b1 100644 --- a/src/non_h_ints_mu/README.rst +++ b/src/non_h_ints_mu/README.rst @@ -9,3 +9,7 @@ The two providers are : +) ao_non_hermit_term_chemist which returns the non hermitian part of the two-electron TC Hamiltonian on the MO basis. +) mo_non_hermit_term_chemist which returns the non hermitian part of the two-electron TC Hamiltonian on the BI-ORTHO MO basis. + +!\sum_mm = 1,3 \sum_R phi_i(R) \phi_k(R) grad_1_u_ij_mu(j,l,R,mm) grad_1_u_ij_mu(m,n,R,mm) +!\sum_mm+= 1,3 \sum_R phi_j(R) \phi_l(R) grad_1_u_ij_mu(i,k,R,mm) grad_1_u_ij_mu(m,n,R,mm) +!\sum_mm+= 1,3 \sum_R phi_m(R) \phi_n(R) grad_1_u_ij_mu(i,k,R,mm) grad_1_u_ij_mu(j,l,R,mm) diff --git a/src/non_h_ints_mu/fit_j.irp.f b/src/non_h_ints_mu/fit_j.irp.f index 695ead7f..6624459f 100644 --- a/src/non_h_ints_mu/fit_j.irp.f +++ b/src/non_h_ints_mu/fit_j.irp.f @@ -20,7 +20,7 @@ END_PROVIDER ! ! 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) +! F(x) is fitted by - 1/sqrt(pi) * exp(-alpha * x) exp(-beta * x^2) (see expo_j_xmu) ! ! The slater function exp(-alpha * x) is fitted with n_max_fit_slat gaussians ! diff --git a/src/non_h_ints_mu/grad_squared.irp.f b/src/non_h_ints_mu/grad_squared.irp.f index a88521a1..01b5d8d6 100644 --- a/src/non_h_ints_mu/grad_squared.irp.f +++ b/src/non_h_ints_mu/grad_squared.irp.f @@ -32,6 +32,103 @@ print*,'Wall time for grad_1_squared_u_ij_mu = ',time1 - time0 END_PROVIDER + BEGIN_PROVIDER [ double precision, grad_1_squared_u_ij_mu_new, (n_points_final_grid, ao_num, ao_num)] + 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 + include 'constants.include.F' + double precision :: r(3),delta,coef + double precision :: overlap_gauss_r12_ao,time0,time1 + integer :: num_a,num_b,power_A(3), power_B(3),l,k + double precision :: A_center(3), B_center(3),overlap_gauss_r12,alpha,beta,analytical_j + 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, coef_i, coef_j, k_ab,center_new(3),p_new,c_tmp,coef_last ! constant factor + double precision :: coefxy, coefx, coefy, coefz,coefxyz + integer :: d(3),lx,ly,lz,iorder_tmp(3),dim1 + double precision :: overlap,overlap_x,overlap_y,overlap_z,thr + dim1=100 + thr = 0.d0 + print*,'providing grad_1_squared_u_ij_mu_new ...' + grad_1_squared_u_ij_mu_new = 0.d0 + 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 + d = 0 + do i = 1, ao_num + do j = 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) + if(ao_overlap_abs(j,i).lt.1.d-12)then + cycle + endif + num_A = ao_nucl(i) + power_A(1:3)= ao_power(i,1:3) + A_center(1:3) = nucl_coord(num_A,1:3) + num_B = ao_nucl(j) + power_B(1:3)= ao_power(j,1:3) + B_center(1:3) = nucl_coord(num_B,1:3) + do l=1,ao_prim_num(i) + coef_i = ao_coef_normalized_ordered_transp(l,i) + alpha = ao_expo_ordered_transp(l,i) + do k=1,ao_prim_num(j) + beta = ao_expo_ordered_transp(k,j) + coef_j = ao_coef_normalized_ordered_transp(k,j) + + ! New gaussian/polynom defined by :: new pol new center new expo cst fact new order + ! from gaussian_A * gaussian_B + call give_explicit_poly_and_gaussian(A_new , A_center_new , alpha_new, fact_a_new , iorder_a_new , & + beta,alpha,power_B,power_A,B_center,A_center,n_pt_max_integrals) + c_tmp = coef_i*coef_j*fact_a_new + if(dabs(c_tmp).lt.thr)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) + do igauss = 1, n_max_fit_slat + delta = expo_gauss_1_erf_x_2(igauss) + coef = coef_gauss_1_erf_x_2(igauss) + coef_last = c_tmp * coef + if(dabs(coef_last).lt.thr)cycle + do lx = 0, iorder_a_new(1) + coefx = A_new(lx,1) + coefx *= coef_last +! 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 = coefz * coefxy +! if(dabs(coefxyz).lt.thr)cycle + iorder_tmp(3) = lz +! call gaussian_product(alpha_new,A_center_new,delta,r,k_ab,p_new,center_new) +! if(dabs(coef_last*k_ab).lt.thr)cycle + call overlap_gaussian_xyz(A_center_new,r,alpha_new,delta,iorder_tmp,d,overlap_x,overlap_y,overlap_z,overlap,dim1) + grad_1_squared_u_ij_mu_new(ipoint,j,i) += -0.25 * coefxyz * overlap + enddo ! igauss + enddo ! ipoint + enddo ! lz + enddo ! ly + enddo ! lx + enddo ! k + enddo ! l + enddo ! j + enddo ! i + call wall_time(time1) + print*,'Wall time for grad_1_squared_u_ij_mu_new = ',time1 - time0 + END_PROVIDER + + BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao_num)] implicit none BEGIN_DOC diff --git a/src/non_h_ints_mu/new_grad_tc.irp.f b/src/non_h_ints_mu/new_grad_tc.irp.f index 068381b4..f205b781 100644 --- a/src/non_h_ints_mu/new_grad_tc.irp.f +++ b/src/non_h_ints_mu/new_grad_tc.irp.f @@ -28,7 +28,7 @@ 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) = + ! 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) ! 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 index c535d0c5..c098b0f5 100644 --- a/src/non_h_ints_mu/test_non_h_ints.irp.f +++ b/src/non_h_ints_mu/test_non_h_ints.irp.f @@ -1,13 +1,14 @@ 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 +! 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 +! call routine_fit + call routine_grad_squared_new end subroutine routine_lapl_grad @@ -85,6 +86,38 @@ subroutine routine_grad_squared end +subroutine routine_grad_squared_new + implicit none + integer :: i,j,k,l,ipoint + double precision :: grad_squared, get_ao_tc_sym_two_e_pot,new,accu,contrib + double precision :: count_n,accu_relat + accu = 0.d0 + accu_relat = 0.d0 + count_n = 0.d0 + do i = 1, ao_num + do j = 1, ao_num + do ipoint = 1, n_points_final_grid + grad_squared = grad_1_squared_u_ij_mu(j,i,ipoint) + new = grad_1_squared_u_ij_mu_new(ipoint,j,i) + 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,ipoint + print*,grad_squared,new,contrib + print*,2.0d0*contrib/dabs(grad_squared+new+1.d-12) + endif + accu += contrib + enddo + enddo + enddo + print*,'accu = ',accu/count_n + print*,'accu/rel = ',accu_relat/count_n + +end + subroutine routine_fit implicit none integer :: i,nx