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'
Read (*,*) n, l
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
!----------------------------------------------------------------------
! Initialize constant radial grid
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'
open(7,file='radial.dat',status='replace')
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.