ONE_PART_QHO
!---------------------------------------------------------------
!Program
name: ONE_PART_QHO
!Original
code: Giannozzi,P., Ercolessi, F. and Gironcoli, S.
!Modified
by Caio Lima Firme
!One
particle quantum harmonic oscillator
!PARABOLIC
POTENTIAL WELL
!Objective:
to plot the wave function for each node
!---------------------------------------------------------------
!The
dimensionless physical quantities:
!Potential
energy: V(x)=0.5X^2
!
Energy: E=1/2+n
!---------------------------------------------------------------
!
Use of Numerov method to find the wave function
!
Function fn from Numerov equation is
! Fn = 1 + K2*h^2/12 , where
K2=2m/((H/2pi)^2)*(E-V)
! H=Planck constant
!---------------------------------------------------------------
integer ::
grid, i, icl
integer ::
n, hn, ncross, j, iter
double
precision :: xmax, h, fh12, norm, djump, fac
double precision ::
Vupper, Vlower, e, k
double precision :: x(0:300), y(0:300),
V(0:300), fn(0:300)
character (len=20) :: fileout
character (len=1) :: A
!------------------------------------------------------
!
Boundary conditions:
y(0) = 0.0d0
x(0) = 0.0d0
xmax=10
grid=300
!
The number 300 refers to the grid everywhere in the code
h =
xmax/grid
! Make explicit one term of the fn function
from the Numerov method: fh12
fh12=h*h/12.0d0
!------------------------------------------------------
!Input data file
print *, 'output file
name is (include .dat at the end): '
read *, fileout
!-------------------------------------------------------------
! Record in the computer memory each value of
x(i) and V(i)
do i = 0, grid
x(i) = i
* h
V(i) =
0.5d0 * x(i)*x(i)
write
(*,*) i, 'X(i)=', x(i),'V(i)=', V(i)
end do
!---------------------------------------------------------------
!
Beginning of searching process
10 search: do
!Input
the number of nodes, n
20 print *,
'give the number of nodes, n (n=0,1,2,3,): '
read *, n
if (n <
0) then
go to 20
end if
!---------------------------------------------------------------
! set initial lower and upper bounds of the
potential V(x)
! Set
trial energy E
Vupper=maxval (V(0:300))
Vlower=minval
(V(0:300))
e = 0.5d0
* (Vlower + Vupper)
iter =
100
Write (*,*) "Vupper=", Vupper,
"e(initial)=", e
!------------------------------------------------------------------
!
BEGINNING OF THE RECURSIVE PROCESS
(STEP=j)
iterate: do j = 1, iter
! set up the fn-function used by the Numerov
method
! fn=1+[2(V-E)]h^2/12 without m (mass) and
H/2pi (Planck constant)
! determine the position of its last crossing,
i.e. change of sign
! fn
< 0 means classically allowed region
! fn
> 0 means classically forbidden region
fn(0)=fh12*(2.0d0*(V(0)-e))
icl=-1
do i=1,grid
fn(i)=fh12*2.0d0*(V(i)-e)
if ( fn(i) == 0.0d0) fn(i)=1.d-15
! store the index 'icl' where the last change
of sign has been found
if (
fn(i) /= sign(fn(i),fn(i-1)) ) icl=i
end do
! fn as
required by the Numerov method
fn =
1.0d0 - fn
! determination of the wave-function in the
first two points
!--------------------------------------------------------------------
! Routine to determine whether n (number of
nodes) is even or odd
! Check
the program EVEN-ODD
hn = n/2
! if n is even, the wave function is even
! if
n is odd, the wave function is odd
if (2*hn
== n) then
! even
number of n: wavefunction is even
y(0) =
1.0d0
! assume f(-1) = f(1)
y(1) =
0.5d0*(12.0d0-10.0d0*fn(0))*y(0)/fn(1)
else
! odd number of n: wavefunction is odd
y(0) =
0.0d0
y(1) =
h
end if
!-------------------------------------------------
!
outward integration and count number of crossings
ncross=0
do i
=1,icl-1
y(i+1)=((12.0d0-10.0d0*fn(i))*y(i)-fn(i-1)*y(i-1))/fn(i+1)
if ( y(i)
/= sign(y(i),y(i+1)) ) ncross=ncross+1
end do
fac =
y(icl)
!-------------------------------------------------
! Routine to determine whether n (number of
nodes) is even or odd
! Check
the program EVEN-ODD
if (2*hn
== n) then
! even
number of n: no node in x=0
ncross
= 2*ncross
else
! odd
number of n: node in x=0
ncross
= 2*ncross+1
end if
!--------------------------------------------
!
check number of crossings
!
FIRST ITERATION PROCESS
!THE NUMBER OF CROSSINGS HAS TO BE EQUAL TO THE
NUMBER OF NODES
!Bisection
method (see BISECTION3 code in chapter two)
if ( iter
> 1 ) then
if
(ncross /= n) then
!
Incorrect number of crossings: adjust energy
if ( j == 1) & print
'("step Energy n Discontinuity Vupper
Vlower icl")'
write (*,5)
j, e, ncross, Vupper, Vlower, icl
5 format (i5,f25.15,i5,f35.15,f35.15,i40)
if
(ncross > n) then
! Too
many crossings: current energy is too high,
!
lower the upper bound
Vupper = e
else
! Too few crossings: current energy is too
low,
! raise the lower bound
Vlower = e
end
if
! New trial value:
e =
0.5d0 * (Vupper+Vlower)
! go to beginning of do loop, don't perform
inward integration
cycle
end
if
end if
!The equality ncross = n was reached
!END OF THE FIRST ITERATION PROCESS
!--------------------------------------------------------------------------
!
Correct number of crossings: proceed to inward integration
!
Determination of the wave-function in the last two points
!
assuming y(grid+1) = 0
y(grid) =
h
y(grid-1) =
(12.0d0-10.0d0*fn(grid))*y(grid)/fn(grid-1)
do i = grid-1,icl+1,-1
y(i-1)=((12.0d0-10.0d0*fn(i))*y(i)-fn(i+1)*y(i+1))/fn(i-1)
end do
! -------------------------------------------------------------------------
!
rescale function to match at the classical turning point (icl)
fac =
fac/y(icl)
Write
(*,*) "fac=", fac
y(icl:) =
y(icl:)*fac
!-------------------------------------------------------------------------
!Normalization of the wave function
!DOT_PRODUCT : scalar product of vectors
!See
the code: NORMALIZATION
norm =
dot_product (y, y)
y = y /
sqrt(norm)
!----------------------------------------------------------------------
! SECOND ITERATION PROCESS
! Bisection method (See BISECTION3 code in
chapter one)
if ( iter
> 1 ) then
!
calculate the discontinuity in the first derivative
! y'(i;RIGHT) - y'(i;LEFT)
djump = ( y(icl+1) + y(icl-1) -
(14.0d0-12.0d0*fn(icl))*y(icl) ) / h
write
(*,30) j, e, n, djump, Vupper, Vlower,icl
30 format
(i5,f25.15,i5,f14.8,f22.15,f34.15,i40)
if
(djump*y(icl) > 0.0d0) then
! Energy
is too high --> choose lower energy range
Vupper = e
else
! Energy is too low --> choose upper energy range
Vlower = e
endif
e =
0.5d0 * (Vupper+Vlower)
! -------------------- convergence test -----------------------------
if (
Vupper-Vlower < 1.d-10) exit iterate
endif
!END OF THE SECOND ITERATION PROCESS
end do
iterate
!END OF
THE RECURSIVE PROCESS
!-------------Convergence
achieved-----------------------
print *, 'Do you want to do another calc
(Y/N)?'
read *, A
Select case (A)
Case ('Y','y')
go to 10
case default
continue
end select
!------------------------------------------------------------------
!Printing
results x(i) and y(i) in the 'fileout'
! x<0 region:
open(7,file=fileout,status='replace')
do i=grid,1,-1
write
(7,'(f7.3,3e16.8,f12.6)') -x(i),
((-1)**n)*y(i)
enddo
! x>0
region:
do i=0,grid
write
(7,'(f7.3,3e16.8,f12.6)') x(i), y(i)
enddo
write (7,'(/)')
close(7)
stop
!--------------------------------------------------------------------------
end do search
stop
Comentários
Postar um comentário