入力例 1s,2s,2p,3p軌道に固有エネルギーの理論値を代入しています。
縦軸も横軸も無い駄目駄目グラフですが、水素原子の1s,2s,2p,3p軌道の動径波動関数を書いています。
横軸の一目盛りは1ボーア半径、縦軸の単位は任意です。エネルギーの単位はハートリー(2リードベルグ)なので、主量子数nに対してとなります。
s軌道(l=0)は原点で有限の値を持ちますが、p軌道(l=1)は原点で0になります。本来、境界条件から波動関数は遠方では0に収束するはずですが、微分方程式の解法が素朴なオイラー法なので、1s軌道の計算は誤差が溜まって距離が大きいところで発散しています。
MODULE m_plot
USE uhoplot
IMPLICIT NONE
INTEGER, PARAMETER :: nx = 800, ny = 600, noff = 50
CONTAINS
SUBROUTINE axis(imax)
IMPLICIT NONE
INTEGER, INTENT(IN) :: imax
INTEGER :: i, ix, iy
CALL gr_move( noff, ny / 2 )
CALL gr_line( nx - noff, ny / 2 )
CALL gr_move( noff, noff)
CALL gr_line( noff, ny - noff)
CALL gr_show()
DO i = 1, imax
ix = (nx - 2 * noff) * i / imax + noff
iy = ny / 2
CALL gr_move(ix, iy)
CALL gr_line(ix, iy + 10)
END DO
RETURN
END SUBROUTINE axis
SUBROUTINE line(p1, p2, q1, q2, xmin, xmax, ymin, ymax)
REAL, INTENT(IN) :: p1, p2, q1, q2, xmin, xmax, ymin, ymax
INTEGER :: ip1, ip2, iq1, iq2
ip1 = (nx - 2 * noff) * p1 / (xmax - xmin) + noff
ip2 = (ny - 2 * noff) * p2 / (ymax - ymin) + ny / 2
iq1 = (nx - 2 * noff) * q1 / (xmax - xmin) + noff
iq2 = (ny - 2 * noff) * q2 / (ymax - ymin) + ny / 2
CALL gr_move(ip1, ip2)
CALL gr_line(iq1, iq2)
RETURN
END SUBROUTINE line
END MODULE m_plot
PROGRAM hatom
USE m_plot
IMPLICIT NONE
REAL, PARAMETER :: dr = 0.01, rmin = 0.01, rmax = 20.0
INTEGER :: l
REAL :: e, r, x, x1, x2
REAL :: a1, b1, a2, b2
l = 0
CALL gr_on('Wavefunction of Hydrogen', nx, ny)
CALL axis(INT(rmax))
DO
WRITE(*, *) ' Input Angular Momentum L & Energy e (e > 100 to exit) '
READ(*, *) l, e
IF (e > 100.0) EXIT
r = rmin
x = r**(l + 1)
x1 = (l+1) * r**l
a1 = r
b1 = x / r
DO
CALL schroed( x2, x, e, r, l )
r = r + dr
x1 = x1 + dr * x2
x = x + dr * x1
a2 = r
b2 = x / r
CALL line(a1, b1, a2, b2, 0.0, rmax, 2.0, -2.0)
a1 = a2
b1 = b2
IF (r > rmax) EXIT
END DO
CALL gr_show()
END DO
CALL gr_off(0)
STOP
CONTAINS
SUBROUTINE schroed(x2, x, e, r, l)
REAL , INTENT(IN ):: x, e, r
REAL , INTENT(OUT):: x2
INTEGER, INTENT(IN) :: l
REAL :: v, a, b
v = - 1.0 / r
a = 2.0 * (v - e)
b = l * (l + 1) / r**2
x2 = (a + b) * x
RETURN
END SUBROUTINE schroed
END PROGRAM hatom