Fortran Help with fortran numerical precision

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:

Fortran:
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?
Cayo

Related Programming and Computer Science News on Phys.org

DrClaude

Mentor
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
Fortran:
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

Mentor
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.

"Help with fortran numerical precision"

Physics Forums Values

We Value Quality
• Topics based on mainstream science
• Proper English grammar and spelling
We Value Civility
• Positive and compassionate attitudes
• Patience while debating
We Value Productivity
• Disciplined to remain on-topic
• Recognition of own weaknesses
• Solo and co-op problem solving