분자동역학(molecular dynamics)에 분할 알고리즘 적용하기

D.P. 란다우 그룹에서 나온 논문 하나를 훑어봤습니다. 제목은 "Molecular and spin dynamics simulations using modern integration methods(현대 적분 방법을 이용한 분자와 스핀동역학 시늉내기)"이고, <아메리칸 저널 오브 피직스>에 2005년에 실렸습니다. 내용 중 일부를 소개합니다.

질량이 각각 m_i인 N개의 입자를 생각합시다. 이들의 위치와 속도는 r_i, v_i이고 두 입자 사이의 거리의 함수인 포텐셜 u(r_ij)로 상호작용합니다. 해밀토니안은 다음처럼 씁니다.

\(H=\sum_i \frac{1}{2}m_iv_i^2+\sum_{i<j}u(r_{ij})\)

j에 의해 i에 가해지는 힘은

\(f_{ij}=-\nabla_{r_i} u(r_{ij})\)

이고, 이 힘들의 합에 의해 입자의 위치가 변하겠죠.

\(m_i\frac{d^2r_i}{dt^2}=\sum_{j\neq i}f_{ij}\equiv f_i\)

모든 입자의 위치와 속도의 집합을 배열(configuration) y(t)={r_i, v_i}라고 하면, 운동방정식은 뿌아송 괄호(Poisson bracket)를 이용해 다음처럼 간단히 씌어집니다.

\(\frac{dy(t)}{dt}=[y(t),H]\equiv \sum_i\frac{1}{m_i}\left(\frac{\partial y(t)}{\partial r_i} \frac{\partial H}{\partial v_i} - \frac{\partial y(t)}{\partial v_i} \frac{\partial H}{\partial r_i} \right)\)

리우빌 연산자(Liouville operator)로 다시 쓰면,

\(\frac{dy(t)}{dt}=\hat L y(t),\ \hat L\equiv \sum_i\left(v_i \frac{\partial}{\partial r_i} + \frac{f_i}{m_i} \frac{\partial}{\partial v_i}\right)\)

입니다. 위 식의 해는 다음과 같습니다.

\(y(t+\tau)=e^{\hat L\tau}y(t)\)

L을 지수 위로 올리니 L 위의 모자(hat)가 보기에 좋지 않네요;;; 여튼 이 문제를 수치적으로 푸는 여러 방법 중에 분할 알고리즘(decomposition algorithm)을 소개합니다. 위와 같은 지수 연산자를 분할한다는 말인데요, 다음처럼 1차 또는 2차까지 분할할 수 있습니다.

\(e^{(A+B)\tau}=e^{A\tau}e^{B\tau}+\mathcal{O}(\tau^2),\ e^{(A+B)\tau}=e^{B\tau/2}e^{A\tau}e^{B\tau/2}+\mathcal{O}(\tau^3)\)

여기서 A와 B는 숫자가 아니라 연산자입니다. 위의 리우빌 연산자를 다음처럼 두 부분으로 나눕니다.

\(\hat L_1=\sum_i v_i\frac{\partial}{\partial r_i},\ \hat L_2=\sum_i \frac{f_i}{m_i}\frac{\partial}{\partial v_i}\)

2차 분할(decomposition to 2nd order)을 적용합니다.

\(e^{(\hat L_1+\hat L_2)\tau}=e^{\hat L_2\tau/2}e^{\hat L_1\tau}e^{\hat L_2\tau/2}+\mathcal{O}(\tau^3)\)

맨 오른쪽 연산자를 우선 적용해봅니다.

\(e^{\hat L_2\tau/2}y=\exp\left(\frac{\tau}{2} \sum_j \frac{f_j}{m_j}\frac{\partial}{\partial v_j}\right)\{r_i(t),v_i(t)\}= \left\{r_i(t),v_i(t)+\frac{f_i(t)}{m_i}\frac{\tau}{2}\right\}\)

즉 시간 τ/2 동안 위치는 그대로인데 속도만 변했습니다. 이렇게 얻은 t+τ/2에서의 속도를 이용해 위치를 변화시킵니다.

\(e^{\hat L_1\tau}e^{\hat L_2\tau/2}y=\left\{r_i(t)+v_i(t)\tau+\frac{f_i(t)}{m_i}\frac{\tau^2}{2}, v_i(t)+\frac{f_i(t)}{m_i}\frac{\tau}{2} \right\}\)

이제 시간 τ 동안 속도는 그대로인데 위치만 변했죠. 마지막으로 맨 왼쪽 연산자를 적용합니다(수식 쳐넣기 귀찮네요). 그 결과 시간 τ/2 동안 위치는 그대로인데 속도만 변했습니다. 다만 이때 속도를 변화시킨 힘은 t+τ에서의 위치의 함수이므로 역시 t+τ에서의 값을 이용해야 합니다. 최종 결과는 다음과 같습니다.

\(r_i(t+\tau)&=&r_i(t)+v_i(t)\tau+\frac{f_i(t)}{m_i}\frac{\tau^2}{2},\\ v_i(t+\tau)&=&v_i(t)+\frac{f_i(t)}{m_i}\frac{\tau}{2}+\frac{f_i(t+\tau)}{m_i}\frac{\tau}{2}\)

L₁만으로 이루어진 연산자는 속도는 그대로 놔두고 이 속도에 의존하여 위치만 바꿉니다. 이것만 그냥 다시 써보겠습니다.

\(e^{\hat L_1\tau}y &=& \exp\left(\tau\sum_j v_j\frac{\partial}{\partial r_j}\right)\{r_i(t),v_i(t)\} \\ &=& \left\{r_i(t)+v_i(t)\tau,v_i(t)\right\} \equiv \{r_i(t+\tau),v_i(t+\tau)\}\)

야코비안을 구합니다.

\(J_i = \left|\begin{array}{cc} \frac{\partial r_i(t+\tau)}{\partial r_i(t)} & \frac{\partial r_i(t+\tau)}{\partial v_i(t)} \\ \frac{\partial v_i(t+\tau)}{\partial r_i(t)} & \frac{\partial v_i(t+\tau)}{\partial v_i(t)} \end{array}\right|=1\)

속도가 위치와 무관하게 그대로 있었기 때문에 야코비안은 1이 되었고, 이는 곧 리우빌 연산자에 의해 상태공간의 부피가 변하지 않았음을 뜻합니다. L₂에 대해서도 마찬가지로 야코비안은 1이 됩니다.

다음으로 시간대칭이 성립하는 걸 보이겠습니다.

\(U(\tau)\equiv e^{B\tau/2}e^{A\tau}e^{B\tau/2},\ U(\tau)U(-\tau)=1\)

위 왼쪽처럼 정의하면 위 오른쪽과 같은 결과를 얻는데 이는 곧 시간대칭을 뜻합니다.

결론을 말하면, 상태공간 부피 보존과 시간대칭은 특히 해밀토니안 시스템에서 중요한 특성인데 컴퓨터로 시늉내기를 할 때에도 이러한 특성들이 잘 지켜지도록 하기에 분할 알고리즘이 적절하다고 합니다.