실습 · 연습문제 2개

LJ 환산 단위($\varepsilon = \sigma = 1$)에서 계의 퍼텐셜 에너지는 모든 원자 쌍의 합이다.

$$E(\mathbf{R}) = \sum_{i < j} 4\left[\left(\frac{1}{r_{ij}}\right)^{12} - \left(\frac{1}{r_{ij}}\right)^{6}\right], \qquad r_{ij} = \lVert \mathbf{r}_i - \mathbf{r}_j \rVert.$$

원자 $i$에 작용하는 힘은 에너지의 음의 gradient $\mathbf{F}_i = -\partial E / \partial \mathbf{r}_i$다. 함수 lj_forces(R)를 완성해 힘 배열을 반환한다.

인자	shape
R	(n, 3) 원자 좌표

반환값의 shape은 (n, 3)이다.

한 쌍 $(i, j)$가 원자 $i$에 주는 힘은 두 원자를 잇는 방향 $\mathbf{r}_i - \mathbf{r}_j$을 향하고, 크기는 거리로 정해진다. 쌍 퍼텐셜을 미분하면

$$\mathbf{F}_i = \sum_{j \neq i} \frac{24\left(2 r_{ij}^{-12} - r_{ij}^{-6}\right)}{r_{ij}^{2}}\,(\mathbf{r}_i - \mathbf{r}_j)$$

이다. 테스트는 다음을 확인한다.

• 힘이 에너지의 중심 차분 수치 미분과 일치한다(해석적 gradient 검증).
• 전체 힘의 합이 0이다(작용·반작용).
• 평형 거리 $2^{1/6}$에 놓인 두 원자 사이의 힘이 0이다.

힌트: 거리 행렬의 대각을 np.inf로 두면 자기 자신과의 항이 0이 된다. np.einsum("ij,ijk->ik", coeff, diff)로 각 쌍의 계수를 방향 벡터에 곱해 합할 수 있다.

분자동역학의 velocity Verlet 적분기는 위치 $\mathbf{R}$, 속도 $\mathbf{V}$, 힘 $\mathbf{F}$를 시간 간격 $\Delta t$만큼 전진시킨다(질량 1). 한 스텝은 세 단계다.

$$\mathbf{V}_{1/2} = \mathbf{V} + \tfrac{1}{2}\Delta t\, \mathbf{F}(\mathbf{R}), \qquad \mathbf{R}' = \mathbf{R} + \Delta t\, \mathbf{V}_{1/2}, \qquad \mathbf{V}' = \mathbf{V}_{1/2} + \tfrac{1}{2}\Delta t\, \mathbf{F}(\mathbf{R}').$$

함수 verlet_step(R, V, force_fn, dt)를 완성해 갱신된 (R_new, V_new)를 반환한다.

인자	의미
R	(n, 3) 현재 위치
V	(n, 3) 현재 속도
force_fn	위치 배열을 받아 힘 (n, 3)을 반환하는 함수
dt	시간 간격

핵심은 힘을 두 번 쓰는 것이다. 반 스텝 속도 갱신에 현재 위치의 힘을, 나머지 반 스텝에 갱신된 위치의 힘을 쓴다.

테스트는 다음을 확인한다.

• 조화 퍼텐셜에서 한 스텝이 위 식과 정확히 일치한다.
• Lennard-Jones 계에서 여러 스텝을 돌렸을 때 총에너지가 보존된다(보존력에서 velocity Verlet의 핵심 성질).
• 힘이 0이면 등속 운동($\mathbf{R}' = \mathbf{R} + \Delta t\,\mathbf{V}$)이 된다.