# Forthcoming book: Quantum Mechanics: Detailed Historical, Mathematical and Computational Approaches

Quantum Mechanics: Detailed Historical, Mathematical and Computational Approaches
By Caio Lima Firme
eBook Published16 June 2022
Imprint: CRC Press
Pages: 462

vi) Fundamentals of Fortran and numerical calculation along with the source codes for numerical solutions of several mathematical and quantum problems. All source codes are in the author’s site: (https://www.fortrancodes.com/);
vii) Chapters devoted to linear algebra and differential equations applied to quantum mechanics and their numerical solutions;
viii) Complete solution for the one-electron and two-electron problems using Schrödinger’s time independent equation along with their source codes.

The code site above is not active yet, but an example of the author’s codes is here. Some comments:

(1) It uses both real*8 and double precision declarations, the first being non-standard and the second archaic.
(2) The `allocate` statement can allocate several variables in the same line, and the `deallocate` statement can deallocate several variables in the same line. `Deallocate` statements just before the end of a program are unnecessary.
(3) Some code has unnecessary parentheses, for example

``````Do i=0,mesh
oldr = rmin + dr*i !oldr is in Bohr unit
r(i)= (exp(oldr))
r2(i)= (r(i) * r(i))
End do
``````
full code -- (I added the declaration of oldr to make it compile)
``````!Program name: HYDROGEN-RADIAL

!Radial equation of the hydrogen-like atom

!Objective: plot the radial wave function for each n,l pair

! The equation is atomic units

!---------------------------------------------------------------

!The Coulomb potential energy equation is:

! V(r) = -2Z/r

! The ceentrifugal potential equation is:

! VL=l(l+1)/r^2

! The energy is:

! E=-2*Z/n^2

!---------------------------------------------------------------

! Use of Numerov method to find the wave function

! Function fn from Numerov equation is

! Fn = 1 + K2*h^2/12 , where K2=-2(Veff - E)

! In the code, h is changed into dr

!---------------------------------------------------------------

Program hydrogen

Implicit none

integer :: mesh, n, l, i, j, zeta, m, maxiter=100

Real*8 rmin, dr, zmesh, y_out_m, oldr

double precision :: e, fh12, norm, eps=1.0E-6

double precision, allocatable :: r(:), r2(:), y(:), vc(:), vcent(:), &

veff(:), k2(:), fn(:), f(:)

!----------------------------------------------------

zeta = 1

!------------------------------------------------------

! Initial data

zmesh= zeta !zmesh is zeta as a real number

rmin= 1.E-5

dr= 0.005d0

!----------------------------------------------------------------

! Number of points for the calculation (mesh)=1000 for Z=1

mesh = 1000

Allocate (r(mesh))

Allocate (r2(mesh))

Allocate (y(mesh))

Allocate (vc(mesh))

Allocate (f(mesh))

Allocate (fn(mesh))

Allocate (vcent(mesh))

Allocate (veff(mesh))

Allocate (k2(mesh))

!----------------------------------------------------------------

!User entering data

Write (*,*) 'Enter n and l. n > 0 and l=n-1'

If (n<1) then

Write (*,*) 'n < 1. Unphysical.'

Stop

Else if (n < l+1) then

Write (*,*) 'n < l+1. It has to be n=l+1!'

Stop

Else if (l < 0) then

Write (*,*) 'l < 0. Unphysical.'

Stop

End if

!-----------------------------------------------------------------

! Initialize exponential radial grid only for the potential

! This grid does not work for 1s and 2p states

Do i=0,mesh

oldr = rmin + dr*i !oldr is in Bohr unit

r(i)= (exp(oldr))

r2(i)= (r(i) * r(i))

End do

!-----------------------------------------------------

! Initialize the potential energy grid for the y(i)

Do i=0,mesh

vc(i)= -2*zmesh/r(i)

vcent(i)=(l*(l+1.d0))/(r2(i))

veff(i)= vc(i)+vcent(i)

End do

!----------------------------------------------------------------------

dr= 0.02d0

Do i=0,mesh

r(i) = rmin + dr*i !oldr is in Bohr unit

r2(i)= (r(i) * r(i))

End do

!----------------------------------------------------------------------

!Set the  eigenvalue

e= -zmesh/n**2

!-----------------------------------------------------------------------

!calculate part of Numerov function f, fh12 and K2

fh12=dr*dr/12.0d0

k2= vc-e+vcent

f(0)= fh12*k2(0)

!-----------------------------------------------------------------

!Start Numerov integration

!----------------------------------------------------------------

!set boundary conditions

y(0)=r(0)**(l+1)

y(1)=r(1)**(l+1)

y(mesh)=dr

m=0

!---------------------------------------------------------------

!find the matching point m checking the change of sign of f

do i=1,mesh

f(i)= fh12*k2(i)

if (f(i) == 0.d0) f(i)=1.E-20

if (f(i) /= sign(f(i),f(i-1))) m=i

end do

!--------------------------------------------------------------

! With the values of f, obtain fn

!obtain the y(mesh-1) for inward integration

fn=1.d0-f

y(mesh-1)=y(mesh)*(12.d0-10.d0*fn(mesh))/fn(mesh-1)

!--------------------------------------------------------------

!Start Numerov outward integration

do i=1,m-1

y(i+1)=y(i)*(12.d0-10.d0*fn(i))-(fn(i-1)*y(i-1))/fn(i+1)

if (y(i) == 0.0d0) y(i)=1.d-20

end do

y_out_m = y(m)

!print *, y(m)

!------------------------------------------------------------------

! Start Numerov inward integration

do i = mesh-1,m+1,-1

y(i-1)=y(i)*(12.d0-10.d0*fn(i))-(fn(i+1)*y(i+1))/fn(i-1)

if (y(i-1) > 1.E10) then

do j=mesh,i-1,-1

y(j)=y(j)/y(i-1)

end do

end if

end do

!------------------------------------------------------------------

! rescale function to match at the turning point

y(:m-1) = y(:m-1)/y_out_m

y(m:) = y(m:)/y(m)

!----------------------------------------------------------------

!Normalization process

norm = sqrt(dot_product(y,y))

y=y/norm

!----------------------------------------------------------------

!Print the result of the wave function

do i=0,mesh

print *, r(i),y(i)

End do

!-------------------------------------------------------------

!Printing results r(i) and y(i) in the 'radial'

do i=0,mesh

write (7,'(3e16.8,f12.6,3e16.8,f12.6)') r(i),y(i)

End do

write (7,'(/)')

Close(7)

!--------------------------------------------------------------------------

Deallocate (r)

Deallocate (r2)

Deallocate (vc)

Deallocate (vcent)

Deallocate (veff)

Deallocate (k2)

Deallocate (y)

Deallocate (f)

Deallocate (fn)

End program
``````

This text will be hidden

I am not an academic and do not have an account at ResearchGate. Someone who does could contact the author here and suggest that the codes be put on GitHub for public comment before the book is published.

4 Likes