간단한 양자문제: 1 Dimensional infinite potential well 문제를 생각해보자.
Potential은 아래와 같이 주어진다.
Schrodinger equation을 이용하면,
$$ - \frac{\hbar^2}{2m} \frac{d^2 \psi}{dx^2} = E\psi $$
이고, 이를 풀면 아래와 같은 규격화된 solution을 갖는다.
$$\psi_n (x) = \sqrt{\frac{2}{a_0}} sin(\frac{n \pi}{a_0} x)$$
이 때 에너지는 다음과 같이 양자화된다.
$$ E_n = \frac{\hbar^2 {k_n}^2}{2m} = \frac{n^2 \pi^2 \hbar^2}{2ma_0^2}, \ \ \ n=1,2,3,...$$
수학적인 관점에서, 방금 해석적으로 푼 1차원 우물 문제는 경계값 문제 (Boundary value problem) 이다. 이런 경계값 문제를 수치적으로는 어떻게 풀어낼 수 있을까? 경계값 문제를 수치적으로 풀어내는 방법 중 하나가 바로 Bisection method (혹은 Shooting method)라고 부르는 방법이다. 원리는 굉장히 간단하다.
우리가 infinite potential 문제에서 풀어내려고 하는 것은, 파동함수에 대해 주어진 경계조건을 만족시키는 양자화 된 에너지 값을 구해내는 것이다. 정답인 에너지 값 이외에는, 파동함수가 발산하기 때문에 경계조건에서 $\psi(x)$ 가 0이 되지 못한다. (이에 대한 자세한 설명은 내 지난 글 [예제 연습장/양자역학] - 에너지가 양자화 되지 않는다면, Eigenfunction은 어떻게 될까? 을 참조하라)
계산 순서는 아래와 같다.
- 임의의 두 에너지 값 $(E_1, E_2)$ 를 설정한다.
- 그리고 두 에너지 값의 midpoint $E' (E'=\frac{ E_1+E_2}{2})$ 를 계산하고, 그 때의 $\psi(a_0)$를 계산한다. (슈뢰딩거 방정식 $H\psi = E\psi$ , 즉 2차 미분 방정식을 RK4와 같은 수치적 방법으로 계산하면 된다.)
- 만약 $\psi(a_0)=0$의 값이 0 이라면, 우리는 정답을 구했다!
- 만약 0이 아니고 같은 부호를 가지고 있다면, $E_1$을 $E'$으로 설정하고, 반대의 경우에는 $E_2$를 $E'$으로 설정한 뒤에, $\psi(a_0)=0$이 되거나 $|E_1 - E_2|$이 미리 정해둔 target accuracy보다 작을 때 까지 2번을 반복한다.
수치적으로 계산하기 위해서 다음 조건을 이용했다.
$$ m = 9.1094e-31 \ (kg) \\ \hbar = 1.0546e-34 \ (J\cdot s) \\ e=1.6022e-19 \ (C) \\ a_0 = 5.292e-11 \ (m) $$
(무차원 변수를 사용해야 하지만, 직관성을 위해 차원을 포함한 변수로 계산을 진행했다. e값은 에너지를 eV단위로 나타내기 위해서 사용했다.)
실제 파이썬으로 계산한 결과는 다음 그래프와 같다. (파동함수는 해답을 구한 뒤 규격화 해주었다.)
구현한 코드를 간략하게 소개해본다.
필요한 패키지를 불러와주고,
import numpy as np아래 함수들을 먼저 정의해준다.
def f(r,kfactor):
psi = r[0]
fpsi = r[1]
fphi = -kfactor*psi
answer = np.array([fpsi,fphi], float)
return answer
def RK4(r,kfactor):
#Container
RK4v_psilist = np.empty(N,float) #x location
RK4v_philist = np.empty(N,float) #x velocity
RK4v_xlist = np.empty(N,float)
#initial
RK4v_psilist[0] = r[0]
RK4v_philist[0] = r[1]
RK4v_xlist[0] = 0
for i in range(N-1):
#RK4 calculation
k1 = dx * f(r,kfactor)
k2 = dx * f(r+k1/2,kfactor)
k3 = dx * f(r+k2/2,kfactor)
k4 = dx * f(r+k3,kfactor)
r += (k1+(2*k2)+(2*k3)+k4)/6
RK4v_psilist[i+1] = r[0]
RK4v_philist[i+1] = r[1]
RK4v_xlist[i+1] = RK4v_xlist[i] + dx
ans = [RK4v_psilist, RK4v_philist, RK4v_xlist]
return ans
def final_psi(Energy):
ksquare = (2*m*Energy)/(hbar**2) ##k^2 = 2mE/hbar^2
r = np.array([0,1],float) #r = [psi,phi]
return RK4(r,ksquare)RK4는 미분방정식을 수치적으로 푸는 방법으로써, 이번 게시글에서는 자세한 작동원리는 소개하지 않는다.
위의 f(r,kfactor)와 RK4(r,kfactor)는 RK4를 구현하기 위한 함수이다. RK4 함수에 초기값 r과 에너지 값(계산상 편의를 위해 $k=\frac{2mE}{\hbar^2}$ 값을 사용했다.)을 넣어주면 경계조건 범위 $[0,a_0]$안의 N개의 좌표에 해당하는 위치의 값을 계산해내고, 리스트로 $\psi$ 값 (psilist)과 $\frac{d\psi}{dx}$ (philist), $x$ 값 (xlist)를 제공해준다.
사실 에너지만 구하려면 마지막 값만 살리면 되는데, 굳이 모든 값을 리스트로 저장해서 리턴하는 이유는, 그래프로 그려주기 위함이다.
마지막 final_psi에 에너지값을 넣어주면, 자동으로 k^2과 초기조건이 들어가서 RK4를 돌린 결과를 리턴해준다.
다음, 초기조건을 설정해주자
#initial value
m = 9.1094e-31 #kg
hbar = 1.0546e-34 #Js
e = 1.6022e-19 #C
a0 = 5.292e-11 #m
dx = a0/1000 #time_step
N = int(a0/dx)+1
#Arbitrary initial energy
#1
E1 = 100*e
E2 = 200*e
target_acc =abs(E1-E2)*0.000001 #target_accuracy이제 Bisection method로 적절한 값이 나올 때 까지 반복한다!
#first check
h1 = final_psi(E1)[0][-1]
h2 = final_psi(E2)[0][-1]
#perform binary search
while abs(E2-E1) > target_acc:
temp_E = (E1+E2)*0.5
checker = final_psi(temp_E)[0][-1] #Last psi value
if h1 * checker > 0:
E1, h1 = temp_E, checker
else:
E2, h2 = temp_E, checker이 While문을 거치고 살아남은 E_1과 E_2의 중간값이 바로 우리가 구한 정답이고, 그 때의 psi값도 구해줄 수 있다.
#Final value
final_E = (E1+E2)*0.5
final_h = final_psi(final_E)
print('E=',final_E/e,"ev")이러면 다음과 같이 최종 에너지가 출력된다.
범위를 100~200eV에서만 찾았는데, 마찬가지로 400~600, 1000~1400, 1800~2200eV 등등등에서 찾아보면 위의 내 결과와 같이 여러 양자화된 에너지 값을 찾을 수 있을 것이다.
그래프까지 그리는 방법은 나중에 기회가 되면 또 작성해보겠다.
