Numerov-Cooley Method on a Potential of NO Molecule

산화질소 분자 퍼텐셜에 적용한 Numerov-Cooley 방법

Abstract

In applying Numerov-Cooley method, Excel tool ‘Solver' is used to match those two wave functions propagated inward and outward, respectively. It is numerically confirmed that the same eigenvalue is obtained by using the average of two energy values of each inward and outward wave functions. This method is applied to a NO molecule potential, and we calculated the variations of the average bond distance and tunneling for a given vibrational energy. It is found that the average bond lengths increase proportionately to the vibrational energy, while the tunneling is not so sensitive to the energy changes. Rather substantial amount of tunnel effect is found for every vibrational state.

Numerov-Cooley 방법을 적용할 때 양쪽에서 전파된 파동함수를 연결시키는 것이 관건인데, 잘 알려진 스프레드시트인 엑셀에서 해찾기 도구를 활용하여 Numerov 방법의 파동함수 연결을 해결하였다. 한쪽방향으로 전파하여 얻은 에너지 값과 반대쪽으로 전파하여 얻은 에너지 값들의 평균을 이용해 상태함수의 고유값을 정하는 방법도 같은 결과를 보임을 확인하였다. 산화질소 분자 퍼텐셜 위에서 이 방법을 적용하였고, 진동에너지에 따른 원자간 평균거리와 터널링 변화를 계산하였다. 진동에너지가 증가하면서 분자 결합 길이는 비례하여 늘어나지만, 터널 효과는 모든 진동 상태가 어느 정도의 확률을 가지며 에너지 증가에 둔감하다는 것을 확인하였다.

서 론

일반적인 퍼텐셜 에너지 함수 V(x) 문제에서 묶인 상태(bound states) Schrödinger 방정식의 수치해석 근을 구하는 방법 중 하나가 Numerov 방법1,2이다. 그 유도과정은 Ref. 2에 상세히 소개되어 있으며, 식으로 나타내면 다음과 같다.

여기서 s는 스텝 사이즈이고,

이다. 이 식은 러시아의 천문학자 Numerov에 의해서 고안되었으며, Ref. 2에서는 스프레드시트에서 어떻게 Numerov 방법을 적용하는가를 잘 보여주고 있다. 하지만 이 방법의 단점은 공간에서 한쪽 방향으로만 쏘아 올린다는 것이다. 한쪽으로만 전파되며 파동함수를 계산하다보니 중간에 발생하는 오차는 누적되어가고 거의 마지막 단계에서는 극에 달하게 되어, 상태함수의 모양도 이상해질 뿐만 아니라 그렇게 얻은 고유값이 정확하지 않게 된다는 약점이 있다. 파동함수는 고유값의 미세한 차이에도 민감하게 변화를 보이는 특성이 있다는 것은 잘 알려져 있다. 실제로, 한 스텝 당 오차는 O(s6) 이지만, 각 스텝의 오차가 누적되면 총 오차는 O(s5)이 된다. 이를 개선시키고자 Cooley3는 (i) 각 양쪽 끝에서 쏘아 두개의 ψ를 얻고, (ii) 이 두개 ψ와 각 ψ의 1차 미분 값들이 중간지점 m에서 부드럽게 연결되도록 Newton-Raphson 방법을 사용하여 고유값 E를 조정하는 방법을 제안 하였다. 이렇게 개량된 방법들이 여러 사람들에 의해 만들어졌고, 이것들을 묶어서 간단히 Numerov-Cooley 방법2-8 이라고 부른다. 특히 Johnson7의 ‘renormalized Numerov method’는 프로그래밍하기가 쉽다는 장점과 효율성 때문에 널리 쓰이고 있는 개량된 Numerov 방법 이다. 하지만, 직접 프로그래밍을 하려면 프로그래밍 언어를 알아야하고 또한 작성된 프로그램이 원하는 결과를 제대로 생산하는 것이 결코 쉽지 않다는 것이다. 다행히 스프레드시트 중 쉽게 사용할 수 있는 엑셀에는 우리가 Numerov-Cooley 방법을 적용하는데 필요한 기능들이 들어있다. Ref. 2에서 사용한 Numerov 방법에서 더 나아가, 여기서 우리는 Numerov-Cooley 방법의 아이디어를 살려 스프레드시트에서 쉽게 이용할 수 있는 개량된 Numerov 방법을 제안 한다. 여기서 제안하는 방법은 파동함수 전파에 따른 누적오차가 확실히 줄어들고, 또한 프로그래밍 언어를 모르는 사람도 쉽게 활용할 수 있다는 것이 장점이다.

 

방법 및 결과

식 (1)을 보면, 지수 k가 작은 값에서 큰 값으로(즉, 왼쪽에서 오른쪽으로) 전파되어 얻어지는 ψL과 큰 값에서 작은 값으로(오른쪽에서 왼쪽으로) 전파되는 파동함수 ψR가 같은 함수이어야 하고, 그것이 제대로 된 상태함수라면 각각 얻어진 고유값(즉, 에너지 값) EL과 ER이 같은 값을 가져야한다. 그러나 실제로 전파시켜보면 파동함수 ψL과 ψR가 숫자상으로 정확히 일치하지 않는 값을 갖거니와 각각 얻어진 에너지 값 EL과 ER도 소수점 이하 몇 자리부터 일치하지 않는다는 것을 쉽게 알 수 있다. 이 문제를 해결하는 방법은 두 가지로 접근할 수 있다. 첫째 방법은, Cooley의 아이디어처럼 ψL과 ψR이 중간지점 m에서 부드럽게 연결되도록 하는 것이다. 그러면서 자연스럽게 상태에너지에 접근하는 고유값 E를 얻게 된다. 엑셀에서는 이 두개 파동함수를 연결시키는 작업을 해찾기 기능을 이용해 해결할 수 있다. ψL과 ψR의 값이 표시된 열에서 적절한 중간지점 m을 선택하여 파동함수 값의 차 ψL(xm)-ψR(xm)를 어느 한 셀에 표시한 후, 해찾기로 이 값이 0이 되는 조건을 찾으면서 에너지 값을 바꾸어 주도록 실행하는 것이다. 해찾기에는 Leon Lasdon(University of Texas at Austin)과 Allan Waren(Cleveland State University)이 개발한 비선형 최적 코드 Generalized Reduced Gradient (GRG2)9가 사용된다. 둘째 방법은, 조금 더 단순하게 두 개의 시트(sheet)에 EL과 ER을 따로 구한 후 세 번째 시트에E = (1/2)(EL+ER)을 해찾기의 추정 값으로 넣고 파동함수를 구하는 것이다. 우리는 이 두 가지 방법을 모두 시도해보았는데 엑셀에서 이 두 방법은 숫자상으로 동일한 결과를 가져온다는 것이 확인되었다. 두 가지 방법 중 사용자가 편리한 것을 택해 사용하면 된다. 또 하나 유용한 현상은, 첫째 방법을 사용할 때 중간지점 m을 선택하는 것이 Cooley가 원래 제안했던 것처럼 파동함수의 진폭이 가장 큰 값인 곳일 필요 없이 중간에 적절한 어느 곳을 택해도 같은 결과를 가져온다는 것이 숫자로 확인되었다.

제안된 방법의 적용 예로 산화질소 분자의 Morse 퍼텐셜을 택하여 분자진동 상태함수를 구했다.

식 (3)에 나타나는 De는 분자의 해리에너지, Re은 평형분자결합거리이며, βe는 평형진동수 ve와 관련지어 결정되는데, 우리는 이미 발표한 결과10를 이용하여 De=6.27eV, βe=2.845Å-1, Re=1.159Å으로 선택하였다. NO 분자의 환산질량 μNO= 7.4664 amu=1.2398×10-26 kg을 이용하면, =5.7641×1013s-1이다. 익숙한 파수 단위로 나타내면, =1,922.7 cm-1, = 50,570.9 cm-1이 된다. Morse 퍼텐셜의 해석적 해(analytic solution)는

으로, 우리가 수치 해(numerical solution)를 구했을 때 그 정확도를 쉽게 비교할 수 있는 장점이 있다. 물론, Numerov-Cooley 방법은 묶인 상태를 구하는 문제라면 어떤 일반적인 퍼텐셜에도 적용될 수 있다.

Schrödinger 방정식을 컴퓨터상에서 풀기위해 물리적 차원 없는 수로 환산하여 나타내면

이 되고, 여기서

이다. 환산된 값들은 xr= βe(R-Re), 이고, De,r=1,383.593이 된다. 물리적으로 가능한 원자들 사이 최단거리는 0 이므로 xr,min= -3.297이다. xr,max는 구하고자하는 상태 양자수 n과 Morse 퍼텐셜의 특성을 고려해 |xr,min|과 ∞ 사이에서 적절히 택하면 된다.

Morse 퍼텐셜은 한정된 개수의 묶인 상태들을 갖는데, 그 양자수는 n = 0, 1, 2, ..., nmax이고

인 자연수가 된다. 따라서 nmax= 52이다. 양자수에 따라 증가하는 비조화성을 생각한다면, 마지막 묶인 상태 nmax를 수치적으로 구해 그 오차를 비교하고, 응용 예를 보이는 것이 제안된 방법의 유용함을 보이는 한 가지 방편이 될 수 있을 것이다. n = 0와 n = 52 상태를 위에서 제시한 방법으로 구했다. 계산 첫 단계에서는 스텝사이즈를 0.02로 사용하고 마지막 단계에서는 스텝사이즈를 0.01로 사용하여 효율적으로 고유값과 상태함수를 구할 수 있었다. 얻은 고유값은 =956.781 cm-1, = 50,570.70 cm-1이다. 식 (4)를 이용하여 계산된 값은, =956.781cm-1, =50,570.87cm-1이다. n = 0 상태의 고유값은 두개의 값이 일치한다. 얻어진 n = 52 상태의 고유값이 참 고유값 보다 작은 것은 이 방법이 변분법이 아니기 때문에 가능하고, 고유값에 약간의 오차(0.0003%)가 나타나는 것은 Morse 퍼텐셜의 최대 값이 xmax=∞이기 때문이다. 원한다면 xr을 조금 더 길게 잡아 오차를 줄일 수 있지만, 누적오차와 계산 시간을 고려하여 선택할 문제이다. 성공적으로 구한 n = 52 상태함수의 제곱을 Fig. 1에 그려 나타내었다.

Fig. 1에서 볼 수 있는 것처럼, Morse 진동자는 진동에너지가 증가하여 해리에너지에 가까이 접근할 때 바깥 반환점에서 많은 시간을 보낸다. 응용 예로, 구한 진동 상태함수를 이용하여, 진동에너지에 따른 원자간 평균거리는 n=0일 때 0= 0.01436(즉, 0= 1.164Å)에서 n = 52일 때 52= 13.56(즉, 52= 5.924Å)으로 분자 결합 길이가 5배 이상 늘어나는 것을 확인하였다. 또한, 고전적 반환점들 밖의 면적을 적분하여, 진동에너지에 따른 양자 효과인 터널링 확률도 계산했는데, n = 0일 때 0.18609에서 n = 52일 때 0.30385로 증가한다. 터널링은 바닥상태(n = 0)의 터널 효과도 무시할 수 없을 만큼 크며, 진동에너지 증가에 그다지 민감하지 않다는 것을 확인할 수 있다. 또 다른 응용으로는 상태함수를 알면 계산 가능한 여러 가지 문제, 예를 들면, Franck-Condon 인자 등의 계산에 이 방법이 활용될 수 있다.

Fig. 1.Normalized |ψ52(xr)|2.

 

결 론

요약하면, Numerov-Cooley 방법을 적용할 때 양쪽에서 전파된 파동함수를 연결시키는 것이 관건인데, 잘 알려진 스프레드시트인 엑셀에서 해찾기 도구를 활용하여 Numerov 방법의 파동함수 연결을 해결하였다. 또한, 한쪽방향으로 전파하여 얻은 에너지 값과 반대쪽으로 전파하여 얻은 에너지 값들의 평균을 구해 상태함수의 고유값을 정하는 방법도 같은 결과를 보임을 확인하였다. 산화질소 분자 퍼텐셜 위에서 이 방법을 적용하였고, 진동에너지에 따른 분자결합 길이 변화와 터널링 변화를 효과적으로 계산하였다. NO 분자 경우, 분자 결합 길이는 진동에너지가 증가하면서 무려 5배 이상 늘어나지만, 터널 효과는 모든 진동 상태가 어느 정도의 확률을 가지며 에너지 증가에 둔감하다는 것을 확인할 수 있었다. 여기서 제시한 방법은 프로그램 언어를 모르는 사람도 쉽게 활용할 수 있다는 장점이 있어, 화학교육적인 목적이나 조금 더 복잡한 화학연구문제의 상태함수 계산 등에 쉽게 이용될 것으로 기대된다.