Molecular dynamics (MD) simulations numerically solve the equations of motion for a multi-particle system described by Newton's second law:
$$
m_i \frac{d^2 \vec{r}_i}{dt^2} = \vec{F}_i(\{\vec{r}_j\})
$$
where $m_i$ and $\vec{r}_i$ are the mass and position of particle $i$, respectively, and the force $\vec{F}_i$ is typically given by the negative gradient of the potential energy function $U(\{\vec{r}_j\})$, i.e., $\vec{F}_i = -\nabla_{\vec{r}_i} U$.

Below are four numerical integration schemes for solving these equations of motion (where acceleration $\vec{a}_n = \vec{F}_n / m$).

**1.**
$$
\begin{cases}
\vec{r}_{n+1} = \vec{r}_n + \vec{v}_n \Delta t \\
\vec{v}_{n+1} = \vec{v}_n+\vec{a}_n \Delta t
\end{cases}
$$

**2.**
$$
\vec{r}_{n+1} = 2\vec{r}_n - \vec{r}_{n-1} + \vec{a}_n(\Delta t)^2
$$

**3.**
$$
\begin{cases}
\vec{r}_{n+1} = \vec{r}_n + \vec{v}_n \Delta t + \frac{1}{2}\vec{a}_n(\Delta t)^2 \\
\vec{v}_{n+1} = \vec{v}_n + \frac{\vec{a}_n + \vec{a}_{n+1}}{2}\Delta t
\end{cases}
$$

**4.**
$$
\begin{cases}
\vec{v}_{n+1/2} = \vec{v}_n + \frac{\Delta t}{2} \vec{a}_n \\
\vec{r}_{n+1/2} = \vec{r}_n + \frac{\Delta t}{2} \vec{v}_n \\
\vec{a}_{n+1/2} = \vec{a}(\vec{r}_{n+1/2}) \\
\vec{v}_{n+1} = \vec{v}_n + \Delta t \cdot \vec{a}_{n+1/2} \\
\vec{r}_{n+1} = \vec{r}_n + \Delta t \cdot \vec{v}_{n+1/2}
\end{cases}
$$

Under ideal numerical precision and using identical initial conditions (e.g., $\vec{r}_0, \vec{v}_0$), **which of the following schemes** would generate exactly the same **discretized** particle trajectories $\{\vec{r}_n\}$?
分子动力学（MD）模拟通过数值方法求解由牛顿第二定律描述的多粒子体系的运动方程：
$$
m_i \frac{d^2 \vec{r}_i}{dt^2} = \vec{F}_i(\{\vec{r}_j\})
$$
其中 $m_i$ 和 $\vec{r}_i$ 分别是粒子 $i$ 的质量和位置，力 $\vec{F}_i$ 通常由势能函数 $U(\{\vec{r}_j\})$ 的负梯度给出，即 $\vec{F}_i = -\nabla_{\vec{r}_i} U$

下面列出了四种用于求解该运动方程的数值积分方案（其中加速度 $\vec{a}_n = \vec{F}_n / m$）。

**1.**
$$
\begin{cases}
\vec{r}_{n+1} = \vec{r}_n + \vec{v}_n \Delta t \\
\vec{v}_{n+1} = \vec{v}_n+\vec{a}_n \Delta t
\end{cases}
$$

**2.**
$$
\vec{r}_{n+1} = 2\vec{r}_n - \vec{r}_{n-1} + \vec{a}_n(\Delta t)^2
$$

**3.**
$$
\begin{cases}
\vec{r}_{n+1} = \vec{r}_n + \vec{v}_n \Delta t + \frac{1}{2}\vec{a}_n(\Delta t)^2 \\
\vec{v}_{n+1} = \vec{v}_n + \frac{\vec{a}_n + \vec{a}_{n+1}}{2}\Delta t
\end{cases}
$$

**4.**
$$
\begin{cases}
\vec{v}_{n+1/2} = \vec{v}_n + \frac{\Delta t}{2} \vec{a}_n \\
\vec{r}_{n+1/2} = \vec{r}_n + \frac{\Delta t}{2} \vec{v}_n \\
\vec{a}_{n+1/2} = \vec{a}(\vec{r}_{n+1/2}) \\
\vec{v}_{n+1} = \vec{v}_n + \Delta t \cdot \vec{a}_{n+1/2} \\
\vec{r}_{n+1} = \vec{r}_n + \Delta t \cdot \vec{v}_{n+1/2}
\end{cases}
$$

在理想的数值精度下，并使用完全相同的初始条件（如 $\vec{r}_0, \vec{v}_0$），**下列哪些方案**会生成完全相同的**离散化**粒子运动轨迹 $\{\vec{r}_n\}$？
1 and 2
1 and 3
1 and 4
2 and 3
2 and 4
3 and 4
Except for 1
Except for 2
Except for 3
Except for 4
All of the above
does not exist
1和2
1和3
1和4
2和3
2和4
3和4
除1之外
除2之外
除3之外
除4之外
以上全部
不存在
We can prove through algebraic transformation that the position Verlet algorithm (2) can be derived from the velocity Verlet algorithm (3).

<Checkpoint>[The position Verlet algorithm (2) can be derived from the velocity Verlet algorithm (3)](1)</Checkpoint>

**Scheme 3 (Velocity Verlet) equations:**
$$
\begin{cases}
\vec{r}_{n+1} = \vec{r}_n + \vec{v}_n \Delta t + \frac{1}{2}\vec{a}_n(\Delta t)^2 & (3a) \\
\vec{v}_{n+1} = \vec{v}_n + \frac{\vec{a}_n + \vec{a}_{n+1}}{2}\Delta t & (3b)
\end{cases}
$$

**Scheme 2 (Verlet) equation:**
$$
\vec{r}_{n+1} = 2\vec{r}_n - \vec{r}_{n-1} + \vec{a}_n(\Delta t)^2 \quad (2)
$$

The most concise proof method utilizes the time-reversibility of these two algorithms. The position update formula (3a) in Scheme 3 is essentially a forward Taylor expansion. We can similarly write a backward Taylor expansion to represent the previous position $\vec{r}_{n-1}$.

1. **Write the forward position step:**  
   This is equation (3a) itself in Scheme 3:  
   $$ \vec{r}_{n+1} = \vec{r}_n + \vec{v}_n \Delta t + \frac{1}{2}\vec{a}_n(\Delta t)^2 \quad (\text{Eq. A}) $$

2. **Write the backward position step:**  
   By replacing the time step $\Delta t$ with $-\Delta t$, we can derive an expression for the previous state $\vec{r}_{n-1}$ from the current state $(\vec{r}_n, \vec{v}_n)$:  
   $$ \vec{r}_{n-1} = \vec{r}_n + \vec{v}_n(-\Delta t) + \frac{1}{2}\vec{a}_n(-\Delta t)^2 $$  
   Simplifying yields:  
   $$ \vec{r}_{n-1} = \vec{r}_n - \vec{v}_n \Delta t + \frac{1}{2}\vec{a}_n(\Delta t)^2 \quad (\text{Eq. B}) $$

3. **Add the forward and backward equations:**  
   We sum (Eq. A) and (Eq. B):  
   $$ \vec{r}_{n+1} + \vec{r}_{n-1} = (\vec{r}_n + \vec{v}_n \Delta t + \frac{1}{2}\vec{a}_n(\Delta t)^2) + (\vec{r}_n - \vec{v}_n \Delta t + \frac{1}{2}\vec{a}_n(\Delta t)^2) $$

4. **Simplify:**  
   The velocity terms $\vec{v}_n \Delta t$ and $-\vec{v}_n \Delta t$ cancel out:  
   $$ \vec{r}_{n+1} + \vec{r}_{n-1} = 2\vec{r}_n + 2 \cdot \frac{1}{2}\vec{a}_n(\Delta t)^2 $$  
   $$ \vec{r}_{n+1} + \vec{r}_{n-1} = 2\vec{r}_n + \vec{a}_n(\Delta t)^2 $$

5. **Rearrange to obtain Scheme 2:**  
   Rearrange the equation to isolate $\vec{r}_{n+1}$ on the left-hand side:  
   $$ \vec{r}_{n+1} = 2\vec{r}_n - \vec{r}_{n-1} + \vec{a}_n(\Delta t)^2 $$  
   This is precisely the expression for Scheme 2. Q.E.D.

<Checkpoint>[2 and 3 are equivalent](2)</Checkpoint>

1 is a first-order method and is clearly not equivalent to the others.

<Checkpoint>[1 is a first-order method and is not equivalent to the others](1)</Checkpoint>

For 4, the so-called improved Euler or midpoint method, its position update is  
$$
\vec{r}_{n+1} = \vec{r}_n + \vec{v}_n \Delta t + \frac{1}{2}\vec{a}_n(\Delta t)^2
$$  
and its velocity update is  
$$
v_{n+1}^{ \text{midpoint}} = \vec{v}_n + \Delta t \cdot \vec{a}\left(\vec{r}_n + \frac{\Delta t}{2}\vec{v}_n\right)
$$  
which differs from that in 3. Thus, from the second step onward, it will produce a different trajectory.

<Checkpoint>[4 will produce a different trajectory from the second step onward and is not equivalent to 2 or 3](1)</Checkpoint>
我们可以通过代数变换证明，从速度Verlet算法（3）可以推导出位置Verlet算法（2）。

<Checkpoint>[从速度Verlet算法（3）可以推导出位置Verlet算法（2）](1)</Checkpoint>

**方案3 (速度Verlet) 的方程组:**
$$
\begin{cases}
\vec{r}_{n+1} = \vec{r}_n + \vec{v}_n \Delta t + \frac{1}{2}\vec{a}_n(\Delta t)^2 & (3a) \\
\vec{v}_{n+1} = \vec{v}_n + \frac{\vec{a}_n + \vec{a}_{n+1}}{2}\Delta t & (3b)
\end{cases}
$$

**方案2 (Verlet) 的方程:**
$$
\vec{r}_{n+1} = 2\vec{r}_n - \vec{r}_{n-1} + \vec{a}_n(\Delta t)^2 \quad (2)
$$

最简洁的证明方法是利用这两种算法的时间可逆性（time-reversibility）。方案3中的位置更新公式(3a)本质上是向前的一个泰勒展开。我们同样可以写出向后的泰勒展开来表示上一步的位置 $\vec{r}_{n-1}$。

1.  **写出向前一步的位置：**
    这是方案3中的方程(3a)本身：
    $$ \vec{r}_{n+1} = \vec{r}_n + \vec{v}_n \Delta t + \frac{1}{2}\vec{a}_n(\Delta t)^2 \quad (\text{Eq. A}) $$

2.  **写出向后一步的位置：**
    通过将时间步 $\Delta t$ 替换为 $-\Delta t$，我们可以得到从当前状态 $(\vec{r}_n, \vec{v}_n)$ 推算上一步状态 $\vec{r}_{n-1}$ 的表达式：
    $$ \vec{r}_{n-1} = \vec{r}_n + \vec{v}_n(-\Delta t) + \frac{1}{2}\vec{a}_n(-\Delta t)^2 $$
    简化后得到：
    $$ \vec{r}_{n-1} = \vec{r}_n - \vec{v}_n \Delta t + \frac{1}{2}\vec{a}_n(\Delta t)^2 \quad (\text{Eq. B}) $$

3.  **将向前和向后的方程相加：**
    我们将 (Eq. A) 和 (Eq. B) 两式相加：
    $$ \vec{r}_{n+1} + \vec{r}_{n-1} = (\vec{r}_n + \vec{v}_n \Delta t + \frac{1}{2}\vec{a}_n(\Delta t)^2) + (\vec{r}_n - \vec{v}_n \Delta t + \frac{1}{2}\vec{a}_n(\Delta t)^2) $$

4.  **化简：**
    等式右边的速度项 $\vec{v}_n \Delta t$ 和 $-\vec{v}_n \Delta t$ 相互抵消：
    $$ \vec{r}_{n+1} + \vec{r}_{n-1} = 2\vec{r}_n + 2 \cdot \frac{1}{2}\vec{a}_n(\Delta t)^2 $$
    $$ \vec{r}_{n+1} + \vec{r}_{n-1} = 2\vec{r}_n + \vec{a}_n(\Delta t)^2 $$

5.  **整理得到方案2：**
    将上式重新整理，把 $\vec{r}_{n+1}$ 单独放在等式左边：
    $$ \vec{r}_{n+1} = 2\vec{r}_n - \vec{r}_{n-1} + \vec{a}_n(\Delta t)^2 $$
    这正是方案2的表达式。证明完毕。

<Checkpoint>[2和3是等价的](2)</Checkpoint>

1是一阶方法，与其他方法显然不等价。

<Checkpoint>[1是一阶方法，与其他方法不等价](1)</Checkpoint>

对于4，即所谓的改进欧拉法或中点法，其位置更新为
$$
\vec{r}_{n+1} = \vec{r}_n + \vec{v}_n \Delta t + \frac{1}{2}\vec{a}_n(\Delta t)^2
$$
速度更新为
$$
v_{n+1}^{ \text{midpoint}} = \vec{v}_n + \Delta t \cdot \vec{a}\left(\vec{r}_n + \frac{\Delta t}{2}\vec{v}_n\right)
$$
与3中不同。因此从第二步开始将给出不相同的轨迹。

<Checkpoint>[4从第二步开始将给出不相同的轨迹，与2，3不等价](1)</Checkpoint>