我只会开发对给定原子的力的计算（BalazsToth 答案的第 1 部分），因为这一点实际上有点棘手，我认为它应该得到更多解释。您可以在任何一本分子模拟书籍中找到有用的评论，我的个人参考是 Allen 和 Tildesley 的“液体计算模拟”以及 W. Hoover 的“分子动力学”物理讲义。

一对粒子相互作用的力沿着由粒子坐标差定义的向量指向。我们称这个向量$\mathbf{r}_{ij}$和具有相同方向的酉向量：

$$\mathbf{r}_{ij} = \mathbf{r}_i - \mathbf{r}_j\ ;\ \mathbf{n}_{ij}=\frac{\mathbf{r}_{ij}}{|\mathbf{r}_{ij}|}=\frac{\mathbf{r}_{ij}}{r}$$

和$|\mathbf{r}_{ij}|=\sqrt{r_{ij,x}^2+r_{ij,y}^2+r_{ij,z}^2}=r$

如果你计算力$f$关于$r$：

$$f_{ij} = -\frac{dV}{dr}$$ 你会注意到：

1. 这个值是正的，当$r<r_0$和消极的时候$r>r_0$.
2. 在第一种情况下，力将原子彼此推开，在第二种情况下使它们更靠近。

因此作用在原子上的力$i$与矢量方向相反$\mathbf{n}_{ij}$正如这里定义的那样。一些作者用相反的值定义它，或者例如，BalazsToth 使用$\mathbf{n}_{ji}$.

然后，您可以将力添加到原子上的总力$i$：

$$\mathbf{F}_{i} = \sum_{j=1\ ;\ j \ne i}^N -f_{ij}\mathbf{n}_{ij}$$

关于你的算法的一句话。鉴于$f_{ij}=-f_{ji}$由于牛顿第三定律，您可以在计算力的函数之外创建循环。您实际上可以使您的代码更像：

void ComputeForces{
for(int i=0;i<N_particles;i++){
  for(int j=i+1;j<N_particles,j++){
    Computeforce(atm[i], atm[j]);
  }
}

在函数 ComputeForce 中，采用两个参数可以让您在不断进行时缩短循环，因为您可以更新原子上的力$j$也是。

我有以下评论：

1.力计算

运动方程$i^\text{th}$有质量的原子$m_i$读作

$\textbf{a}_i=\frac{d\textbf{v}_i}{dt}=\frac{1}{m_i}\cdot f(\vert\textbf{r}_j-\textbf{r}_i\vert)\textbf{n}_{ji}$,

其中粒子间力$f(\vert\textbf{r}_j-\textbf{r}_i\vert)$只取决于粒子之间的距离$i$和$j$. 在您的实现中，您使用每个方向上的距离来计算力，这是不正确的。相反，您只需计算每对粒子的力一次（使用它们之间的标量距离）并将其与归一化方向向量相乘$\textbf{n}_{ji}=\frac{\textbf{r}_j-\textbf{r}_i}{\vert\textbf{r}_j-\textbf{r}_i\vert}$以获得力矢量的分量。

2. 邻居

在大量原子的情况下效率不高，因为$\mathcal{O}(n^2)$计算复杂度。

由于原子之间的相互作用随着距离的增加而变得越来越小，因此应用某个阈值或截止距离以仅考虑每个原子的相邻原子是有益的。根据阈值，您可以将域划分为阈值大小的单元并执行邻居搜索。这是一种无需实际计算任何一对原子之间的距离即可找到潜在邻居的技术。根据要考虑的潜在邻居列表，您可以计算距离并检查它们是否真的在影响半径内。有几种方法可以找到潜在的邻居，但总的来说$\mathcal{O}(n)$可以实现复杂度。在此处阅读有关邻居搜索算法的更多信息。

当然，如果您的相互作用定律具有无限半径，则忽略影响半径之外的粒子可能会导致模拟错误。因此，阈值的应用是准确性和计算性能要求之间的权衡。

3. 粒子容器

您当前的实现将粒子存储为 AoS（结构数组）。从编程的角度来看，这是一种方便的方法，但效率低于 SoA（数组结构），尤其是对于并行计算。在后者中，您将创建一个大型结构，其中包含存储所有粒子属性的数组。它最终得到内存中对齐良好的数据，这对于这种类型的计算来说是更可取的。