Range of DVR parameters for the Calculation of Vibrational Energy of Anharmonic Oscillators

비조화 진동자 진동에너지 계산에 적합한 DVR 계산 변수 결정

Abstract

We summarize the discrete variable representation method which is a simple numerical method enabling us to calculate the vibrational energies and wave functions of anharmonic oscillators. The ranges of its parameters well-performing for the calculation of fundamental and overtone transition energies are predicted by analyzing the model of Morse oscillator.

비조화진동자의 진동에너지와 파동함수를 수치적으로 간단히 계산할 수 있는 Discrete Variable Representation 방법을 요약하고, Morse 진동자를 모형으로 활용하여 실제 분자계의 기본전이 및 들뜬 전이에너지를 구하는데 적용할 수 있는 최적의 변수 범위를 예측하였다.

서 론

양자화학적 Schrödinger 방정식을 수치적으로 푸는 여러가지 방법들 중에서 Discrete Variable Representation (DVR) 방법은 간단한 원리에 기반하면서 여러가지 장점을 가지고 있어서 널리 활용되고 있다.1 특히, 응축상 동역학 연구에 최근 그 중요성이 커지고 있는 진동 분광학 실험에 대한 이론적 해석 연구를 수행할 때 매우 유용한 방법이 양자/고전 혼합방법인데, 이 이론적 방법을 개발하기 위해서는 분자들의 여러가지 위치 상태에 대해 진동상태를 양자화학적으로 풀어야 한다.2,3 이러한 연구를 위해서 효과적이고 정확한 진동상태를 구할 수 있는 DVR 방법의 유용성은 매우 크다고 할 수 있다.

본 논문에서는 보통의 분자계에서 보이는 크기의 비조화성을 가지는 진동자들의 기본(0-1) 전이와 들뜬(1-2) 전이에너지를 효과적으로 계산할 수 있는 DVR 방법에 대해 요약하고, Morse 진동자로 모형화된 비조화진동자들의 전이에너지를 DVR 방법으로 계산할 때 필요한 파라미터의 범위를 예측하고 최적화된 값들을 찾아본다.

조화진동자

해밀토니안 작용자가 아래와 같이 주어지는 임의의 일차원 조화진동자를 고려해 보자.

여기에서 이고 는 진동자의 평형점으로부터의 벗어난 변위를 나타내며 평형점에서는 그 값이 0 이다. 는 그 운동량 작용자, μ는 질량이다. k는 힘 상수로서, 조화포텐셜 에너지의 평형점에서의 이차 미분 값이다. 이제, 다음과 같이 단위가 없는 작용자들을 도입하자.

ω는 조화진동자의 각진동수로서 힘 상수 k와 의 관계를 가진다. 위치와 운동량 작용자들의 양자역학적 교환자 특성, , 에 의해 이 작용자들의 교환자는 를 만족함을 쉽게 확인할 수 있다. 이 작용자들을 식 (1)에 대입하면, 단위가 없는 조화진동자 해밀토니안(ĥ ≡ Ĥ0/hω)은

와 같이 나타낼 수 있으며 이다. 조화진동자에 대한 보존(Boson) 작용자는

와 같이 정의되며 이 작용자들의 교환자는 을 만족함을 알 수 있다. 변위와 운동량 작용자를 이 보존 작용자들로 나타내면

이 되고, 이를 식 (3)에 대입하면 조화진동자의 해밀토니안을 다음과 같이 나타낼 수 있다.

이차 양자화 이론 방법을 이용하여 이 해밀토니안의 eigenstate의 성질을 다음과 같이 구할 수 있다.4

여기서 이다. 따라서, 인접 상태간 전이에너지는

이 되며, 인접 상태간 전이에너지는 상태의 에너지가 증가해도 일정하다. 조화진동자의 eigenstate 들은 Boson 작용자들에 대해 다음 성질을 만족함을 보일 수 있다.4

여기서 n = 0, 1, 2, ⋯이다. 이 성질들과 eigenstate 들의 직교 성질 ( = δmn )을 이용하여, , , 의 조화진동자 상태 표현을 아래와 같이 얻을 수 있다.

운동량과 변위 작용자의 이차항 합은 대각화되어 있음을 확인할 수 있다.

Morse 진동자

Morse 진동자의 해밀토니안은 다음과 같이 정의된다.

여기서

은 Morse 포텐셜 에너지, D는 포텐셜 에너지 깊이, Qe는 평형점을 나타낸다. 평형점에서의 힘 상수는 k = 2a2D로 얻어진다. 따라서, 평형점에서의 조화진동수는 로 얻어진다.

단위가 없는 Morse 해밀토니안(ĥMO ≡ ĤMO/ħω)의 eigenvalue와 eigenstate를 고려하면

속박된 상태에 대한 해밀토니안의 eigenvalue는 해석적으로 다음과 같이 잘 알려져 있다.5

여기서 이며 비조화성 상수라 한다. n = 0,1, 2,..., [δ −1−1/2]이며 [x]는 x 보다 작은 최대 정수를 의미한다. 따라서, 인접상태간 전이에너지는

이 되며, 인접 상태간 전이에너지는 상태의 에너지가 증가함에 따라 비조화성 만큼 감소하게 된다. 따라서 실험적으로 측정이 편리한 기본 전이진동수 ω10 ≡ ω(ε1– ε0) 와 비조화 진동수 Δ ≡ ω {ε1– ε0 – (ε2 – ε1)} 의 항으로써 Morse 포텐셜 에너지의 두 매개변수인 ω와 δ를 나타낼 수 있다.

Morse 진동자의 전이 변위는 상태간 전이와 관련된 여러가지 물리량을 결정하는 중요한 인자이다. 이에 대한 해석적 해가 m > n인 경우에 다음과 같이 주어진다.6

Discrete Variable Representation (DVR)

일반적 비조화 진동자를 고려해 보자. 이 진동자의 해밀토니안은 아래와 같다.

식 (2)에 의해 주어진 조화 진동자의 단위가 없는 운동량과 변위 작용자를 이용하여 식 (20)을 단위가 없는 해밀토니안(ĥ ≡ Ĥ/ħω)의 형태로 쓸 수 있다.

예를 들어 Morse 포텐셜 에너지의 경우

로 나타낼 수 있다. 이제, 해밀토니안 (21)을 조화진동자의 eigenstate representation 으로 표현할 수 있다. 조화진동자 eigenstate 들로 구성된 무한 차원의 Hilbert 공간 대신, 근사적으로 유한 차원의 공간을 고려하자. 보통 바닥상태 근처의 비조화 진동자의 상태를 기술할 것이므로, 조화진동자의 바닥상태로부터 N개의 상태{|0〉, |1〉 ,..., |N – 1〉}가 basis를 구성하는 Hilbert 공간을 고려하면, 다음의 completeness가 근사적으로 성립한다. . 이 completeness를 이용하면, 임의의 작용자들을 조화진동자의 유한한 Hilbert 공간에서 행렬로 표현하는 것이 가능하다. 이제, 식 (12)를 이용하여 변위 작용자를 그러한 행렬로 표현한 후, 그 행렬의 eigenvalue 방정식을 수치적으로 풀어서, 를 만족하는 작용자 의 eigenstate N개를 구한다. Eigenstate 들 의 complete를 이용하여 를 아래와 같이 나타낸다.

이제, 포텐셜 작용자를 Taylor 전개 형태로 표현하면,

이 되고, 여기에 식 (23)를 대입하면,

가 된다. 이 식을 이용하면, N개의 eigenvalue에 해당하는 변위에서 포텐셜 에너지(v(qn))를 계산하여 포텐셜에너지의 변위작용자 표현을 얻게 되고, 이를 다시 조화진동자 표현으로 변환하면,

을 얻는다. 식 (26)을 식 (10)과 결합하여 비조화 진동자 해밀토니안의 조화진동자 표현 행렬 을 얻게 된다. 여기서 이다. 이제 행렬의 eigenvalue를 구하는 표준 수치 방법을 활용하여,

식을 만족하는 비조화 진동자의 에너지와 진동상태를 얻게 된다. 수치적으로 구한 비조화 진동자의 eigenstate를 아래와 같이 표현하고,

식 (12)를 이용하면, 비조화 진동자 전이 변위를 다음과 같이 얻을 수 있다.

 

결과 및 고찰

DVR 방법을 이용하여 진동상태를 계산하기 위해 필요한 첫 단계는 선택된 N에 대해 변위작용자 의 eigenvalue들 N개를 구하는 것이다. N = 2~14에 대해 수치적으로 구한 eigenvalue 들을 Fig. 1에 나타내었다.

Figure 1.N 값(좌측 세로축)에 따른 변위 작용자 의 eigenvalue (수평축, Solid circle). 곡선은 변위 q에 따른 조화 포텐셜 에너지(V (q) = q2/2, 우측 세로축)를 나타낸다.

변위 작용자의 eigenvalue 들은 q = 0를 중심으로 대칭적으로 분포하며, N이 증가함에 따라 그 분포 범위도 늘어남을 볼 수 있다. 그리고 주어진 차원 N에 해당하는 변위 eigenvalue 들은 N 개의 양자진동에너지에 해당하는 포텐셜 장벽에 갇혀 있음을 알 수 있다. 즉, N 차원의 Hilbert 공간에서 변위 eigenvalue 의 최댓값(q*)은 조화 진동자가 N개의 양자에너지를 가질 때 진동할 수 있는 최대 변위와 비례함을 볼 수 있으며 몇 개의 N에 대해 그 값들을 Table 1에 나타내었다.

Table 1.N 차원의 DVR 계산에서 변위 작용자의 최대 eigenvalue 의 크기

DVR 계산 방법에서 N 의 값이 증가할수록 유한 공간의 completeness의 타당성은 늘어나게 되어 좀 더 정확한 계산이 가능할 것으로 기대된다. 그러나, Fig. 1과 Table 1에서 볼 수 있는 바와 같이 N이 증가할수록 변위 작용자의 eigenvalue 들의 값의 범위도 증가한다. 화학 결합을 비롯한 일반적인 비조화포텐셜에너지를 가지는 진동자의 계산을 위해서는 식 (25)에서 볼 수 있듯이 eigenvalue 에서의 포텐셜에너지를 계산할 필요가 있으나, 이를 위해 양자화학 계산을 수행할 경우 비현실적인 변위에서 양자화학 계산에 오류가 발생할 수 있다. 특히, 평형점보다 짧은 변위값의 경우 원자간 반발력이 급격히 증가하여 정확한 포텐셜에너지를 얻을 수 없는 문제가 발생한다. 뿐만 아니라, 포텐셜에너지를 구하기 위해서는 일반적으로 많은 양자화학 계산 시간이 소요되기 때문에, 그러한 양자화학 계산을 수반한 실제 분자계를 위한 계산에서는 임의로 큰 N을 설정하는 것이 바람직하지 않으며 계산의 정확성을 담보할 수 있는 적절한 수준에서 결정되어야 한다.

일반적인 비조화 진동자의 진동 특성은 Morse 진동자로서 잘 대변될 수 있으므로, 식 (14)에 정의된 Morse 포텐셜에너지 를 모델로 이용하여 진동자의 양자화학적 기본 진동에너지와 비조화성의 특성에 따라 N 값의 범위가 어떻게 결정되는지 살펴보기로 한다. 위에서 언급한 바와 같이, 포텐셜에너지의 양자화학적 계산이 가능한 최대의 수축 변위(qα)는 그 위치에서의 포텐셜에너지에 의해 결정된다고 가정하자. 그러므로 qα를 VMO(qα) = αħω 인 변위로 정의한다(여기서 α는 양자화학적 계산의 가능 범위를 한정하는 파라미터이다. 이 변위는 반발력 영역의 높은 에너지에 해당하므로, 2δα > 1을 만족해야 한다). 즉, qα는 Morse 진동에너지(ħω)의 α배에 해당하는 포텐셜 에너지를 가지는 반발력 영역의 최대 수축 변위로 정의되며 를 만족한다. 따라서, 주어진 α에 대해 최대 수축변위는 진동자의 비조화성 δ에 의해 아래와 같이 결정된다.

한편, DVR 계산의 타당성을 위해서, 관심의 진동상태(예를 들어 3 차 비선형 분광학에서는 바닥상태로부터 이차 들뜬 상태까지)의 파동함수가 펼쳐진 변위 영역의 포텐셜 에너지는 정확히 기술되어야 할 필요가 있다. 그러한 영역들의 최대 확장 변위를 qβ라 표시하고, VMO(qβ) = βħω 인 변위로 정의하자(여기서 β는 관심대상의 진동상태가 정의된 진동자 변위의 최대 범위를 한정하는 파라미터이며, 최대 확장 변위에서의 에너지는 진동자의 해리 에너지보다 클 수 없으므로 2δβ < 1을 만족해야만 한다). 따라서, 주어진 β에 대해 최대 확장 변위는 진동자의 비조화성 δ에 의해 아래와 같이 결정된다.

여기서 이다. 이 두 변위값을 비조화성의 함수로서 Fig. 2에 나타내었다.

Figure 2.Morse 진동자의 비조화성 δ의 변화에 따른 경계변위 (붉은 선)과 (검은 선). qe= 0.

이 두 변위 사이의 변위들에 대한 양자화학적 상태를 정확히 계산하기 위해서는 이 변위 구간의 포텐셜에너지 정보가 DVR 계산에 정확히 반영되어야 한다. 따라서, 변위 작용자의 eigevnvalue 들의 분포는 이 변위 구간을 모두 포함할 수 있도록 설정되어야 할 것이다. 그런데, 변위 작용자의 eigenvalue 들은 Fig. 1에서 보듯이 q = 0를 중심으로 대칭적으로 분포하며 이 값의 분포 범위를 -q*~q* 라고 하면,

의 두 부등식을 만족하는 양수의 q*가 존재 해야 함을 의미한다. 이 두 부등식을 q*에 대한 부등식으로 바꾸면

이 된다. max[A, B]는 A, B 중 더 큰 수를 의미한다.

Fig. 3은 α = 20, qe = 0 일 때, 세 개의 진동변위영역 파라미터(β = 2,3,4)에 대해 위에서 도입된 세 변위 경계값 들을 진동자 비조화성(δ)의 함수로서 나타낸 그림이다. 이 그림이 보여주는 바와 같이, α = 20, qe = 0, β = 4의 파라미터 경우에는 진동자의 비조화성(δ)이 0 에서부터 약 0.04 보다 작을 때 q*≈4.1 이 식 (33)를 만족하나, 비조화성이 0.04 보다 커지면 식 (33)을 만족하는 q*가 존재하지 않음을 알 수 있다.

Figure 3.Morse 진동자의 비조화성(δ)에 따른 (dashed lines), (dash-dot lines)과 −qα (solid line). α = 20, qe = 0.

따라서, δ = 0~0.1의 비조화성을 가지는 모든 경우에 식 (33)을 만족하는 q* 값이 존재하기 위해서는, 그림에서 보이듯이 좀 더 작은 β의 값이 설정되면 되는데(예를 들어, β = 2이면 식 (33)을 만족하는 q*는 약 3.3 정도 이다), 너무 작은 β의 경우 1-2 들뜬 전이에너지를 계산하는데 필요한 포텐셜 에너지 영역을 충분히 포함하지 못해서 계산의 정확성을 낮추게 될 것이다. 좀 더 바람직한 방법은 적절한 크기의 β에 대하여 을 좀 더 낮추기 위해 음수의 qe를 도입하는 것이다.

예를 들어 다음 Fig. 4는 α와 β는 Fig. 3에서와 동일하지만 qe = −1.0로 변화하여 얻은 그림이다. 이 그림은 Morse 포텐셜 에너지의 평형점과 다른 조화진동자의 eigenstate를 활용한 DVR 계산법의 경우 훨씬 넓은 범위의 비조화성 (δ = 0~0.1)에 대해서 식 (33)을 만족하는 q*가 존재함을 보여준다. 즉, β = 3의 경우, q*≈3.8~4.5 범위에서 이 그림에 나타난 모든 비조화성에 대해서 식 (33)을 만족하는 값임을 알 수 있다. 양자화학 계산을 수행할 eigenvalue 들의 갯수를 최소화해야 할 필요성이 있으므로 이들 중 최솟값인 q*≈3.8이 선호된다.

Figure 4.Morse 진동자의 비조화성(δ)에 따른 (dashed lines), (dash-dot lines), −qα (solid lines). 각 선이 나타내는 바는 Fig. 3 에서와 동일하다. α = 20, qe = −1.0.

Fig. 5는 α와 β는 Fig. 4에서와 동일하지만 qe = −1.5로 변화하여 얻은 그림이다. Fig. 4로부터 −qα와 이 다소 상승하고, 은 다소 하락하여 결국 식 (33)을 만족하는 q*의 범위는 β = 3에 대해 q*≈3.8~5.0 범위에서 존재하여 앞에서와 같은 이유로 q*≈3.8이 선호된다. qe를 계속 낮추면 −qα와 의 차이는 늘어나지만, 이 상승하여 식 (33)을 만족하는 q*의 최솟값은 점점 늘어나며 그에 해당하는 변위작용자 eigenvalue의 갯수는 증가하여 양자화학 계산의 소요 시간은 늘어나므로 바람직하지 않다. 따라서, 최적의 q*의 값은 3.8로 추정되며, qe= −1.0~−1.5 범위에서 결정될 것이다. 따라서 Table 1을 참고하면, δ = 0~0.1 범위의 비조화성을 가지는 진동자의 기본 및 1-2 전이에너지 계산에 최적인 DVR의 차원 N은 12로 추정된다.

Figure 5.Morse 진동자의 비조화성(δ)에 따른 (dashed lines), (dash-dot lines), −qα (solid lines). 각 선이 나타내는 바는 Fig. 4에서와 동일하다. α = 20, qe = −1.5.

 

결 론

본 논문을 통해, 비조화진동자의 양자화학적 진동상태를 수치적으로 간단히 구할수 있는 DVR 방법에 대해 요약하고, 이를 실제 분자계에 적용하기 위해 필요한 파라미터 범위를 고려하여 최적의 파라미터를 찾았다. 진동상태의 기본전이와 1-2 전이와 관련된 실제 화학 진동현상의 비조화성을 Morse 진동자가 잘 대변한다고 가정하고 이를 모델로 하여 분자계에 일반적으로 나타나는 크기의 다양한 비조화성에 적용될 수 있도록 파라미터들이 결정되었다.