Fortran double precision normalization factor shows single precision error

[Cayo]

Hello. I need help to understand why my code is not giving me double numerical precision.
In this piece of code, I create a function that is analytically normalized, but when I calculate the numerical normalization factor, it seems to be with single precision. Here it is:

program dyna
implicit none
integer i
integer, parameter :: Nq1=61
real(kind=8), parameter :: saw2au=1822.889950851334d0 
real(kind=8), parameter :: nitro=14.0067d0*saw2au 
real(kind=8), parameter :: mredu=nitro**2.0d0/(2.0d0*nitro)
real(kind=8) :: e0,pi,ch,x,x0,stepX,w,expo,c0,c1
complex(kind=8) :: soma,vec0(Nq1)
pi = 3.141592653589793d0
e0=0.005d0
x0=2.09970623d0 
stepX=0.01d0
w=2.d0*e0
c0 = ((mredu*w)/pi)**(1.d0/4.d0)
c1 = (4.d0/pi*(mredu*w)**3.d0)**(1.d0/4.d0)
do i=1,Nq1
  ch=(i-(Nq1-1.d0)/2.d0-1.d0)*stepX
  x=x0+ch 
  expo = dexp(-(x-x0)**2.d0*mredu*w/2.d0)
  vec0(i) = c0 * expo
end do
!----------------------------------------!
!normalizing                             !
soma=0.0d0                               !
do i=1,Nq1                               !
  soma=soma+conjg(vec0(i))*vec0(i)*stepX !
end do                                   !
vec0=vec0/sqrt(soma)                     !
!----------------------------------------!
write(*,'(a24,2(f23.15))')'normalization constant =',soma
end program

I get 0.999998932390816, and it should be 1.0d0.
I use complex vectors because later in the code they will become complex.
I don't understand where is the problem. Can someone help me, please?
Thanks in advance,
Cayo

 

[DrClaude]

kind=8 doesn't guarantee that you will get double precision. The Fortran 90 standard doesn't specify how kind is to be implemented, so it is compiler-dependent.

If you want to insure that all variables are declared as double precision, you can use

integer, parameter :: dp = kind(1.d0)

real(kind=dp) :: x
complex(kind=dp) :: z

 

[Cayo]

Thank you for the reply, you are correct.
But, I changed what you suggested and it didn't changed the result. Any other idea?

 

[DrClaude]

I looked at your code in more details, and I don't understand why you say that

it should be 1.0d0.

You are approximating an integral from $-\infty$ to $\infty$ using a finite grid. Why would you expect that to be an exact? Check the values of vec0(1) and vec0(Nq1).

 

[Cayo]

Hello DrClaude,
You are right. Thats the problem, the size of the grid is too small. Thank you for the help.