我正在对二氧化硅进行分子动力学模拟。前段时间我转向波动偶极子模型，经过努力我仍然无法实现它。
简而言之，系统中的所有氧原子都是可极化的，它们的偶极矩取决于它们相对于系统中所有其他原子的位置。更具体地说，我使用 TS 势（http://digitallibrary.sissa.it/bitstream/handle/1963/2874/tangney.pdf?sequence=2），其中在每个时间步迭代地发现偶极子。
这意味着在评估作用在原子上的力时，我必须考虑到这种势能对坐标的依赖性。以前，我使用简单的成对势模型，因此我将程序设置为使用通过对势能表达式求微分获得的解析公式来计算力。
现在我很茫然：如何实现这个新的潜力？在我发现的所有文章中，它们只给你公式，而不是算法。正如我所看到的，当我计算作用在某个原子上的力时，我必须考虑这个原子的偶极子的变化，所有相邻原子的偶极子的变化，然后是更多原子的偶极子的变化，等等，因为它们相互依赖。毕竟，正是由于这种相互依赖性，在每个时间步都迭代地发现了偶极子。显然，我不能为每个原子迭代地计算力，因为算法的计算复杂度太高了。我应该使用一些简单的函数来解释偶极子的变化吗？这看起来也不是一个好主意，因为偶极子是迭代计算的，精度很高，
那么如何实现这个模型呢？此外，是否可以像我之前所做的那样分析计算力，或者是否有必要使用导数的有限差分公式来计算它们？
我还没有在文学中找到我的问题的答案，但是如果你知道一些文章、网站或书籍，其中突出显示了这些材料，请引导我到那个来源。
感谢您的时间！
==================================================== =================================
更新：
谢谢您的回答。不幸的是，这不是我的问题。我没有问如何计算偶极子，而是如何计算力，因为这些偶极子随运动而变化很大。
我试图以直接的方式计算力（不考虑偶极子通过它们的距离相互依赖，只计算每个步骤上的偶极子，然后像这些偶极子是静态的一样计算力），但我得到的结果在物理上并不正确。
为了分析这种情况，我建立了一个仅由两个原子组成的系统的模拟：Si 和 O。它们具有相反的电荷，因此它们会振荡。能量时间依赖性图形如下所示：
顶部曲线表示动能，中间曲线表示未考虑偶极相互作用的势能，底部曲线表示系统的势能，其中考虑了偶极相互作用。
从这张图中可以清楚地看到，系统正在做它不应该做的事情：爬上潜在的斜坡。所以我认为这是因为我没有考虑偶极矩坐标依赖性。例如，在给定的时间点，我们计算力，并且它们被引导以使两个原子相互移动。但是，当我们确实将它们彼此靠近（甚至是稍微移动）时，偶极矩发生了变化，我们发现我们实际上最终获得了比以前更高的势能！在下一个时间步中，情况相同。
所以问题是，如何考虑这种影响，导致我能想到的几种方法要么计算量太大，要么太粗糙。