LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2025-01-03 10:05:57 +01:00
Code Issues Releases Wiki Activity

Beginning logn range integrals

Browse Source
This commit is contained in:
Emmanuel Giner 2017-04-14 19:23:12 +02:00
parent 1df5dced1e
commit 54abf1ddc4
2 changed files with 44 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
src/Integrals_Bielec/EZFIO.cfg
View File

@ -51,3 +51,19 @@ doc: If |<ij|kl>| < ao_integrals_threshold then <pq|rs> is zero
interface: ezfio,provider,ocaml
default: 1.e-15
ezfio_name: threshold_mo
[mu_erf]
type: double precision
doc: cutting of the interaction in the range separated model
interface: ezfio,provider,ocaml
default: 0.5
ezfio_name: mu_erf
[long_range]
type: logical
doc: if true, compute all the integrals using the long range interaction
interface: ezfio,provider,ocaml
default: False
ezfio_name: long_range

37
src/Integrals_Bielec/ao_bi_integrals.irp.f
View File

@ -158,7 +158,7 @@ double precision function ao_bielec_integral_schwartz_accel(i,j,k,l)
q_inv = 1.d0/qq
schwartz_kl(s,r) = general_primitive_integral(dim1, &
Q_new,Q_center,fact_q,qq,q_inv,iorder_q, &
Q_new,Q_center,fact_q,qq,q_inv,iorder_q) &
Q_new,Q_center,fact_q,qq,q_inv,iorder_q) &
* coef2
schwartz_kl(0,r) = max(schwartz_kl(0,r),schwartz_kl(s,r))
enddo
@ -471,9 +471,15 @@ double precision function general_primitive_integral(dim, &
! Gaussian Product
! ----------------
double precision :: p_plus_q
if(long_range)then
p_plus_q = (p+q) * ((p*q)/(p+q) + mu_erf*mu_erf)/(mu_erf*mu_erf)
else
p_plus_q = (p+q)
endif
pq = p_inv*0.5d0*q_inv
pq_inv = 0.5d0/(p+q)
pq_inv = 0.5d0/p_plus_q
p10_1 = q*pq ! 1/(2p)
p01_1 = p*pq ! 1/(2q)
pq_inv_2 = pq_inv+pq_inv
@ -548,7 +554,7 @@ double precision function general_primitive_integral(dim, &
return
endif
rho = p*q *pq_inv_2
rho = p*q *pq_inv_2 ! le rho qui va bien
dist = (P_center(1) - Q_center(1))*(P_center(1) - Q_center(1)) + &
(P_center(2) - Q_center(2))*(P_center(2) - Q_center(2)) + &
(P_center(3) - Q_center(3))*(P_center(3) - Q_center(3))
@ -574,7 +580,8 @@ double precision function general_primitive_integral(dim, &
double precision :: rint_sum
accu = accu + rint_sum(n_pt_out,const,d1)
general_primitive_integral = fact_p * fact_q * accu *pi_5_2*p_inv*q_inv/dsqrt(p+q)
! change p+q in dsqrt
general_primitive_integral = fact_p * fact_q * accu *pi_5_2*p_inv*q_inv/dsqrt(p_plus_q)
end
@ -620,10 +627,16 @@ double precision function ERI(alpha,beta,delta,gama,a_x,b_x,c_x,d_x,a_y,b_y,c_y,
ASSERT (gama >= 0.d0)
p = alpha + beta
q = delta + gama
double precision :: p_plus_q
if(long_range)then
p_plus_q = (p+q) * ((p*q)/(p+q) + mu_erf*mu_erf)/(mu_erf*mu_erf)
else
p_plus_q = (p+q)
endif
ASSERT (p+q >= 0.d0)
n_pt = ishft( nx+ny+nz,1 )
coeff = pi_5_2 / (p * q * dsqrt(p+q))
coeff = pi_5_2 / (p * q * dsqrt(p_plus_q))
if (n_pt == 0) then
ERI = coeff
return
@ -650,15 +663,21 @@ subroutine integrale_new(I_f,a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z,p,q
integer :: i, n_iter, n_pt, j
double precision :: I_f, pq_inv, p10_1, p10_2, p01_1, p01_2,rho,pq_inv_2
integer :: ix,iy,iz, jx,jy,jz, sx,sy,sz
double precision :: p_plus_q
j = ishft(n_pt,-1)
ASSERT (n_pt > 1)
pq_inv = 0.5d0/(p+q)
if(long_range)then
p_plus_q = (p+q) * ((p*q)/(p+q) + mu_erf*mu_erf)/(mu_erf*mu_erf)
else
p_plus_q = (p+q)
endif
pq_inv = 0.5d0/(p_plus_q)
pq_inv_2 = pq_inv + pq_inv
p10_1 = 0.5d0/p
p01_1 = 0.5d0/q
p10_2 = 0.5d0 * q /(p * q + p * p)
p01_2 = 0.5d0 * p /(q * q + q * p)
p10_2 = 0.5d0 * q /(p * p_plus_q)
p01_2 = 0.5d0 * p /(q * p_plus_q)
double precision :: B00(n_pt_max_integrals)
double precision :: B10(n_pt_max_integrals), B01(n_pt_max_integrals)
double precision :: t1(n_pt_max_integrals), t2(n_pt_max_integrals)