分子动力学模拟：涨落偶极子模型的实现

Question

我正在对二氧化硅进行分子动力学模拟。一段时间前，我转向波动偶极子模型，经过很大的努力，我仍然有问题实现它。

简而言之，系统中的所有氧原子都是可极化的，它们的偶极矩取决于它们相对于系统中所有其他原子的位置。更具体地说，我使用TS势(http://digitallibrary.sissa.it/bitstream/handle/1963/2874/tangney.pdf？sequence=2)，其中偶极子在每个时间步迭代找到。

这意味着，在计算作用在原子上的力时，我必须考虑到这个势能对坐标的依赖关系。之前，我使用的是简单的成对势模型，所以我将我的程序设置为使用解析公式来计算力，解析公式是通过微分势能表达式得到的。

那么我如何实现这个模型呢？还有，有没有可能像我以前那样用解析的方法计算力，或者有没有必要用导数的有限差分公式来计算力？

我还没有在文献中找到我的问题的答案，但如果你知道一些文章、网站或书籍，这些材料被突出显示，请告诉我那个来源。

谢谢你抽出时间！

为了分析这种情况，我模拟了一个由两个原子组成的系统：Si和O，它们有相反的电荷，所以它们会振荡。能量时间依赖关系图如下所示：

顶部的曲线代表动能，中间的曲线代表不考虑偶极相互作用的势能，底部的曲线代表系统的势能，其中考虑了偶极相互作用。

您可以从这张图形中清楚地看到，系统正在做它不应该做的事情：沿着潜在的斜坡往上爬。所以我决定这是因为我没有把偶极矩坐标依赖考虑进去。例如，在一个给定的时间点，我们计算力，它们的方向是使两个原子互相移动。但当我们确实把它们朝对方移动（即使是轻微的），偶极矩改变了，我们发现我们实际上结束了比以前更高的势能！在下一个时间步骤中，情况是相同的。

Answer 1

不确定我是否完全理解您的问题，但听起来您可能需要实现一个马尔可夫链类型的解决方案？

更多信息请看这篇文章：http://freakonometrics.hypotheses.org/6803

编辑。我建议这样做的原因是，这听起来像是一个系统，其中每个原子的状态依赖于它的邻居，而邻居的状态又依赖于它们的邻居，以此类推。从概念上讲，这可以被建模为一个巨大的矩阵，您可以根据它的邻居（？？？）迭代更新每个值。这是很难解决的，但是本文展示了如何使用马尔可夫链来解决一个非常大的转移矩阵问题，而不是计算实际的矩阵。