PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-11-03 12:43:48 +01:00
Code Releases Activity

Fixed division by zero in RSDFT

Browse Source
This commit is contained in:
Anthony Scemama 2020-05-20 11:45:41 +02:00
parent 764aed423e
commit 512525508b
3 changed files with 34 additions and 25 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
src/dft_utils_func/routines_exc_sr_lda.irp.f
View File

@ -138,6 +138,8 @@ subroutine ex_lda_sr(mu,rho_a,rho_b,ex,vx_a,vx_b)
!Density and kF
rho_a_2=rho_a*2.D0
akf = ckf*(rho_a_2**f13)
! Avoid division by zero
if (akf == 0.d0) akf = 1.d-20
a = mu/(z2*akf)
a2 = a*a
a3 = a2*a
@ -169,6 +171,7 @@ subroutine ex_lda_sr(mu,rho_a,rho_b,ex,vx_a,vx_b)
!Density and kF
rho_b_2= rho_b * 2.d0
akf = ckf*(rho_b_2**f13)
if (akf == 0.d0) akf = 1.d-20
a = mu/(z2*akf)
a2 = a*a
a3 = a2*a

51
src/utils/need.irp.f
View File

@ -92,42 +92,47 @@
end
double precision function erf0(x)
implicit double precision (a-h,o-z)
if(x.lt.0.d0)then
erf0=-gammp(0.5d0,x**2)
implicit none
double precision, intent(in) :: x
double precision, external :: gammp
if(x < 0.d0)then
erf0=-gammp(0.5d0,x*x)
else
erf0=gammp(0.5d0,x**2)
erf0=gammp(0.5d0,x*x)
endif
end
double precision function gammp(a,x)
implicit double precision (a-h,o-z)
if(x.lt.0..or.a.le.0.)stop 'error in gammp'
if(x.lt.a+1.)then
implicit none
double precision, intent(in) :: a, x
double precision :: gln, gammcf
if(x<0.d0.or.a<=0.d0) then
stop 'error in gammp'
endif
if(x < a+1.d0)then
call gser(gammp,a,x,gln)
else
call gcf(gammcf,a,x,gln)
gammp=1.-gammcf
gammp=1.d0-gammcf
endif
return
end
subroutine gser(gamser,a,x,gln)
implicit double precision (a-h,o-z)
parameter (itmax=100,eps=3.e-7)
parameter (itmax=100,eps=3.d-7)
gln=gammln(a)
if(x.le.0.)then
if(x.lt.0.) stop 'error in gser'
gamser=0.
if(x.le.0.d0)then
if(x.lt.0.d0) stop 'error in gser'
gamser=0.d0
return
endif
ap=a
sum=1./a
sum=1.d0/a
del=sum
do 11 n=1,itmax
ap=ap+1.
ap=ap+1.d0
del=del*x/ap
sum=sum+del
if(abs(del).lt.abs(sum)*eps)go to 1
@ -139,14 +144,14 @@
subroutine gcf(gammcf,a,x,gln)
implicit double precision (a-h,o-z)
parameter (itmax=100,eps=3.e-7)
parameter (itmax=100,eps=3.d-7)
gln=gammln(a)
gold=0.
a0=1.
gold=0.d0
a0=1.d0
a1=x
b0=0.
b1=1.
fac=1.
b0=0.d0
b1=1.d0
fac=1.d0
do 11 n=1,itmax
an=float(n)
ana=an-a
@ -155,8 +160,8 @@
anf=an*fac
a1=x*a0+anf*a1
b1=x*b0+anf*b1
if(a1.ne.0.)then
fac=1./a1
if(a1.ne.0.d0)then
fac=1.d0/a1
g=b1*fac
if(abs((g-gold)/g).lt.eps)go to 1
gold=g

5
src/utils/prim_in_r.irp.f
View File

@ -24,8 +24,9 @@ double precision function primitive_value_explicit(power_prim,center_prim,alpha,
end
double precision function give_pol_in_r(r,pol,center, alpha,iorder, max_dim)
double precision :: r(3), center(3), alpha,pol(0:max_dim,3)
implicit none
integer, intent(in) :: iorder(3), max_dim
double precision :: r(3), center(3), alpha,pol(0:max_dim,3)
integer :: i,m
double precision :: gauss(3), x
gauss = 0.d0
@ -33,7 +34,7 @@ double precision function give_pol_in_r(r,pol,center, alpha,iorder, max_dim)
do m = 1, 3
x = r(m) - center(m)
do i = 0, iorder(m)
gauss(m) += pol(i,m) * dexp(-alpha *x**2 ) * x**i
gauss(m) += pol(i,m) * dexp(-alpha *x*x ) * x**i
enddo
enddo
give_pol_in_r = gauss(1) * gauss(2) * gauss(3)