分子动力学模拟:波动偶极子模型实现

Question

我正在进行二氧化硅的分子动力学模拟.前段时间我转向波动的偶极子模型,经过多方努力,我仍然遇到问题.

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

这意味着在评估作用于原子的力时,我必须考虑这种对坐标的潜在能量依赖性.以前,我使用的是简单的成对潜在模型,所以我会设置我的程序来使用通过区分潜在能量表达式获得的分析公式来计算力.

现在我不知所措:如何实现这一新潜力？在我发现的所有文章中,他们只给出了公式,而不是算法.在我看来,当我计算力,作用于某个原子时,我必须考虑到这个原子的偶极子的变化,所有相邻原子的偶极子的变化,然后是更多原子的偶极子的变化等等,因为它们彼此依赖.毕竟,正是由于这种相互依赖性,在每个时间步骤迭代地发现了偶极子.显然,我不能迭代地为每个原子计算力,因为算法的计算复杂度太高了.我应该使用一些简单的函数来解释偶极子的变化吗？这看起来也不是一个好主意,因为偶极子是迭代计算的,具有高精度,然后,它实际上很重要(计算力),我们会使用原油函数吗？

那么我该如何实现这个模型呢？此外,是否可以像以前一样分析计算力,或者是否有必要使用导数的有限差分公式计算它们？

我没有在文献中找到我的问题的答案,但是如果你知道一些文章,网站或书籍,这些材料被突出显示,请指导我.

感谢您的时间!

================================================== ================================

更新:

谢谢您的回答.不幸的是,这不是我的问题.我没有问过如何计算偶极子,但是如果计算这些偶极子随着运动的变化而变化很大.

我试图以直接的方式计算力(不考虑偶极子通过它们的距离相互依赖,只计算每一步的偶极子,然后计算力就像那些偶极是静态的),但我得到的结果在物理上是不正确的.

为了分析这种情况,我建立了一个由两个原子组成的系统的模拟:Si和O.它们具有相反的电荷,因此它们会振荡.能量时间依赖图形看起来像这样:

顶部的曲线代表动能,中间的曲线代表潜在的能量而不考虑偶极相互作用,底部的曲线代表系统的势能,其中偶极相互作用被考虑在内.

从图中可以清楚地看到,系统正在做它不应该做的事情:爬上潜在的斜坡.所以我认为这是因为我没有考虑偶极矩坐标依赖性.例如,在给定的时间点,我们计算力,并且它们被定向以便使两个原子彼此相向移动.但是当我们确实将它们朝向彼此移动时(甚至略微移动),偶极矩会发生变化,我们发现实际上我们的能量比以前更高!在下一个时间步骤中情况是相同的.

所以问题是,如何将这种影响考虑在内,导致我能想到的几种方式无论是计算量太强还是太粗糙.

Answer 1

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

有关更多信息，请参阅这篇好文章：http://freakonometrics.hypotheses.org/6803

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