シリカの分子動力学シミュレーションを行っています。しばらく前に、変動双極子モデルに目を向けましたが、多大な努力を払ったにもかかわらず、まだ実装に問題があります。
つまり、システム内のすべての酸素原子は分極可能であり、それらの双極子モーメントは、システム内の他のすべての原子に対する位置に依存します。より具体的には、TS ポテンシャル ( http://digitallibrary.sissa.it/bitstream/handle/1963/2874/tangney.pdf?sequence=2 ) を使用します。双極子は各時間ステップで繰り返し検出されます。
これは、原子に作用する力を評価する際に、座標に対するこのポテンシャル エネルギーの依存性を考慮に入れる必要があることを意味します。以前は単純なペアワイズ ポテンシャル モデルを使用していたので、ポテンシャル エネルギーの式を微分して得られる解析式を使用して力を計算するようにプログラムを設定していました。
今、私は途方に暮れています: この新しい可能性をどのように実装するのですか? 私が見つけたすべての記事では、数式のみが提供され、アルゴリズムは提供されていません。私が見ているように、特定の原子に作用する力を計算するとき、この原子の双極子の変化、すべての隣接原子の双極子の変化、さらに多くの原子の双極子の変化を考慮に入れる必要がありますなど、互いに依存しているためです。結局のところ、各時間ステップで双極子が繰り返し検出されるのは、この相互依存性のためです。明らかに、アルゴリズムの計算が非常に複雑になるため、原子ごとに力を繰り返し計算することはできません。双極子の変化を説明するためにいくつかの単純な関数を使用する必要がありますか? 双極子は高精度で反復的に計算されるため、これも良い考えのようには見えません。
では、このモデルをどのように実装すればよいでしょうか。また、以前のように解析的に力を計算することは可能ですか、それとも導関数の有限差分式を使用して計算する必要がありますか?
私の質問に対する答えは文献で見つかりませんでしたが、この資料が強調されている記事、サイト、または本を知っている場合は、その情報源を教えてください.
お時間をいただきありがとうございます！
================================================== ================================
アップデート：
ご回答ありがとうございます。残念ながら、これは私の質問ではありませんでした。双極子を計算する方法は尋ねませんでしたが、これらの双極子が動きによって大きく変化することを考えると、力を計算する方法を尋ねました。
私は単純な方法で力を計算しようとしました (距離による双極子の相互依存性を考慮せず、各ステップで双極子を計算し、それらの双極子が静的であるかのように力を計算します)、得られた結果は物理的に正しくありませんでした。
この状況を分析するために、Si と O の 2 つの原子だけで構成される系のシミュレーションを設定しました。これらは反対の電荷を持っているため、振動します。エネルギー時間依存のグラフは次のようになります。
上の曲線は運動エネルギーを表し、真ん中の曲線は双極子相互作用を考慮しない場合のポテンシャルエネルギーを表し、下の曲線は双極子相互作用が考慮されたシステムのポテンシャルエネルギーを表します。
この図から、システムがすべきでないこと、つまり潜在的な勾配を上っていることがはっきりとわかります。これは、双極子モーメントの座標依存性を考慮していなかったことが原因であると判断しました。たとえば、特定の時点で力を計算すると、両方の原子が互いに近づくように力が向けられます。しかし、それらを互いに近づけると (少しでも)、双極子モーメントが変化し、実際には以前よりも高い位置エネルギーになっていることがわかります! 次の時間ステップの間、状況は同じです。
問題は、この効果をどのように考慮に入れるかということです。私が考えることができるいくつかの方法は、計算量が多すぎるか、粗雑すぎるためです。