diff --git a/src/ao_two_e_ints/two_e_integrals.irp.f b/src/ao_two_e_ints/two_e_integrals.irp.f index e60e6eeb..b4b21f5c 100644 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ b/src/ao_two_e_ints/two_e_integrals.irp.f @@ -444,14 +444,17 @@ BEGIN_PROVIDER [ logical, ao_two_e_integrals_in_map ] END_PROVIDER -BEGIN_PROVIDER [ double precision, ao_two_e_integral_schwartz,(ao_num,ao_num) ] - implicit none +! --- + +BEGIN_PROVIDER [ double precision, ao_two_e_integral_schwartz, (ao_num, ao_num) ] + BEGIN_DOC ! Needed to compute Schwartz inequalities END_DOC - integer :: i,k - double precision :: ao_two_e_integral,cpu_1,cpu_2, wall_1, wall_2 + implicit none + integer :: i, k + double precision :: ao_two_e_integral,cpu_1,cpu_2, wall_1, wall_2 ao_two_e_integral_schwartz(1,1) = ao_two_e_integral(1,1,1,1) !$OMP PARALLEL DO PRIVATE(i,k) & @@ -468,6 +471,7 @@ BEGIN_PROVIDER [ double precision, ao_two_e_integral_schwartz,(ao_num,ao_num) ] END_PROVIDER +! --- double precision function general_primitive_integral(dim, & P_new,P_center,fact_p,p,p_inv,iorder_p, & diff --git a/src/utils/one_e_integration.irp.f b/src/utils/one_e_integration.irp.f index cacc3bf7..9c1d2445 100644 --- a/src/utils/one_e_integration.irp.f +++ b/src/utils/one_e_integration.irp.f @@ -92,41 +92,48 @@ subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,& overlap = overlap_x * overlap_y * overlap_z end + +! --- + +subroutine overlap_x_abs(A_center, B_center, alpha, beta, power_A, power_B, overlap_x, lower_exp_val, dx, nx) - -subroutine overlap_x_abs(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx) - implicit none BEGIN_DOC ! .. math :: ! ! \int_{-infty}^{+infty} (x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) dx ! END_DOC - integer :: i,j,k,l - integer,intent(in) :: power_A,power_B - double precision, intent(in) :: lower_exp_val - double precision,intent(in) :: A_center, B_center,alpha,beta - double precision, intent(out) :: overlap_x,dx - integer, intent(in) :: nx - double precision :: x_min,x_max,domain,x,factor,dist,p,p_inv,rho - double precision :: P_center - if(power_A.lt.0.or.power_B.lt.0)then + + implicit none + + integer, intent(in) :: power_A, power_B, nx + double precision, intent(in) :: lower_exp_val, A_center, B_center, alpha, beta + double precision, intent(out) :: overlap_x, dx + + integer :: i, j, k, l + double precision :: x_min, x_max, domain, x, factor, dist, p, p_inv, rho + double precision :: P_center + double precision :: tmp + + if(power_A.lt.0 .or. power_B.lt.0) then overlap_x = 0.d0 dx = 0.d0 return endif - p = alpha + beta - p_inv= 1.d0/p - rho = alpha * beta * p_inv - dist = (A_center - B_center)*(A_center - B_center) + + p = alpha + beta + p_inv = 1.d0/p + rho = alpha * beta * p_inv + dist = (A_center - B_center)*(A_center - B_center) P_center = (alpha * A_center + beta * B_center) * p_inv - if(rho*dist.gt.80.d0)then + + if(rho*dist.gt.80.d0) then overlap_x= 0.d0 return endif + factor = dexp(-rho * dist) - double precision :: tmp tmp = dsqrt(lower_exp_val/p) x_min = P_center - tmp