분자동역학(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)라고 하면, 운동방정식은 뿌아송 괄호를 이용해 다음처럼 간단히 씌어집니다.

\(\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 H}{\partial r_i} \frac{\partial y(t)}{\partial v_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)\)

수치적으로 푸는 여러 방법 중에 분할 알고리즘(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차 분할