電解質水溶液中のイオンの溶媒和自由エネルギーの計算について

[Masaru Matsugami]

松林先生，ユーザの皆様

お世話になります。松上＠熊本高専です。

お忙しいところ，すみません。

今回は，電解質溶液中のイオンの溶媒和自由エネルギー計算について質問させて頂きます。

現在，溶媒和イオン液体中へ各種イオンを挿入し，溶媒和自由エネルギーを計算しております(ERmod 0.3.5)。

例えば，挿入イオン: Li+, 溶媒: LiG3TFSAなどを計算を行なっています。

その計算結果に違和感がありました。

・ヒストグラム（エネルギー分布関数）のふるまいがおかしい？(Li+ [1st com.], G3[2nd com.], TFSA- [3rd com])

※Li+-Li+の2体ポテンシャルは，rの全領域で正となるが，エネルギーが負の領域が大きい。

イオンーイオン間のヒストグラムは，E=0 kcal/molが最大とならず，-5 kcal/mol, 5kcal/mol付近にピークが現れる。

(添付ファイル：engsln_Li-IL.01, engref_Li-IL.01)

・溶媒和自由エネルギーを溶媒種別に見ると，1st component (Li+)の結果が負となる。

(添付ファイル：slvfe_Li-IL.txt)

 cumulative average & 95% error for solvation free energy

              total             1st component         2nd component         3rd component

  1  -145.4704              -99.7811              -49.6419               15.6197

  2  -141.4065     8.1277   -95.0011     9.5601   -50.2959     1.3081    15.5575     0.1242

  3  -141.5180     4.6978   -92.3599     7.6399   -51.4740     2.4743    13.9830     3.1499

  4  -136.4334    10.6980   -90.3083     6.7839   -51.6518     1.7853    17.1936     6.7966

  5  -135.5069     8.4913   -89.7562     5.3695   -51.8888     1.4619    17.8051     5.4048

  6  -135.8892     6.9751   -89.5283     4.4078   -51.9506     1.2000    17.2568     4.5472

  7  -136.8060     6.1736   -90.5196     4.2200   -51.7917     1.0628    17.1723     3.8468

  8  -139.3073     7.3220   -91.3685     4.0297   -51.9370     0.9652    15.6651     4.4928

  9  -142.1017     8.5401   -92.9897     4.8108   -51.9185     0.8520    14.4735     4.6238

 10  -144.9586     9.5390   -95.0245     5.9225   -51.8814     0.7657    13.6143     4.4784

もう少し単純な系でテスト計算を行なってみました。

NaCl溶液（約20w%, 298K）中へNa+を挿入し，溶媒和自由エネルギーを計算してみました。

solute: Na+ 1個, solvent: Na+(1st) 50個, Cl-(2nd) 50個, SPC/E Water(3rd) 650個

GROMACS 5.0.4, PME, Cut-off=1.35 nm

こちらでも，イオン間のヒストグラムが同様のふるまいをしました。

(添付ファイル：engsln_Na-NaClaq.01, engref_Na-NaClaq.01)

また，溶媒和自由エネルギーの値もおかしいような気がします。

(添付ファイル：slvfe_Na-NaClaq.out)

そこで質問ですが，

・添付ファイルのような，ヒストグラム（エネルギー分布関数）のふるまいは正しいのでしょうか。

・MD計算の計算法に問題があるのでしょうか。

・あるいは，ERmodの計算のパラメーター設定に問題があるのでしょうか。

もし，電解質溶液中のイオンの溶媒和自由エネルギーを計算する際の注意点/計算法がある場合は，教えて頂けたら助かります。

お忙しい中，大変ご迷惑をおかけしますが，アドバイスよろしくお願い申し上げます。

松上

[石塚良介]

松上さん:

お久しぶりです。石塚です。

イオン液体も電解液もエネルギー分布関数の0 kcal/molで

ピークが出てこないのは妙ですね。

ただ、情報が少ないので、何が問題かはわかりません。

こういう時は、次の点を確認します。

１．MDの平衡化ができているか？

２．系のサイズは十分とれているか？

→　系のサイズを２倍にしても溶媒和自由エネルギーに

変化がないことを確認する。

３．溶質ー溶媒動径分布関数に異常は見られないか？

→　ちゃんと1に収束しているかの確認。

４．４に異常がなければ、分子間相互作用と動径分布関数から

solvation energyを計算してERmodの計算したsolvation energyと比較する。

もし、両者が一致しなければ、エネルギー分布関数の計算が

おかしい可能性がある。

５．ERmodの溶媒和自由エネルギーが収束しているか？

→　サンプリングが十分とれているかの確認。特に、有効電荷や

分極を考慮していないイオン液体はドロドロなので、

長時間シミュレーションが必要になる（50nsでも足りないことがある）。

engdiv、numref、numslnを大きめにとり、収束チェック。

こんなところでしょうか？

石塚

2017年7月1日 8:43 Masaru Matsugami <masar...@gmail.com>:

--
You received this message because you are subscribed to the Google Groups "ermod-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ermod-users+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Masaru Matsugami]

石塚様

お忙しい中，コメント頂き，ありがとうございました。

アドバイス頂いた点を1つずつ確認してみます。

溶媒和イオン液体の系では，系のサイズを2倍にしてもだいたい同様の結果が得られています。

NaClの系は，少しサイズが小さいのでもっと大きくして再度計算をしてみます。

また，何か分かりましたらご報告させて頂きます。

よろしくお願い申し上げます。

松上

2017年7月2日日曜日 0時49分53秒 UTC+9 石塚良介:

To unsubscribe from this group and stop receiving emails from it, send an email to ermod-users...@googlegroups.com.

[Masaru Matsugami]

石塚様、松林先生

お世話になります。松上です。

先日の石塚様のコメントを参考にして再度20w%NaClaq系でテストをしてみましたので報告・ご相談させて頂きます。

大変お忙しいと思いますが、コメントを頂けたら幸いです。

問題点：電解質溶液中のイオン種の溶媒和エネルギーを計算する際に、

　　　　エネルギー分布関数の0 kcal/molでピークが出ない。

前回は、Na+(solute)1個, Na+Cl- 50個, SPC/E water 650個と系が小さく、計算時間も2 nsと短かった。

今回は、系のサイズを大きくした（2倍, 4倍, 8倍）。

　　　　50 ns (50000 snapshot)計算を行った。

その結果、エネルギー分布関数の振る舞いがどう変化するか確認してみました。

※とりあえずは、時間がなかったのでsoln系のみエネルギー分布関数の計算を行っています。

結果は、系のサイズを大きくすると、エネルギー分布関数の振る舞いに改善が見られた（添付ファイル：Fig. 1）。

しかしながら、8倍してもエネルギー分布関数の0 kcal/molでピークが出なかった。

溶質-溶媒のg(r)も1に収束しているように見えますが、まだ計算時間が足りないでしょうか（添付ファイル：Fig. 2）。

何か改善策がありましたら教えて頂けると助かります。

ご迷惑をおかけしますが、よろしくお願い申し上げます。

松上

2017年7月2日日曜日 16時07分56秒 UTC+9 Masaru Matsugami:

[Nobuyuki MATUBAYASI]

松上さん

イオンーイオンの対エネルギーが０にピークを持たず、

しかし、系を大きくすれば、ピークが０に近づいて行くことは正しい挙動です。

Ewaldは、対エネルギーを全MD cellで積分したら０になるように定式化されています。

ですので、距離が近いときにマイナスのエネルギーであれば、遠距離でプラスになり、

距離が近いときにプラスであれば、遠距離でマイナスです。

つまり、同種イオンの場合に、小さいマイナスエネルギーにピークが出て、

異種イオンで、小さいプラスエネルギーにピークが出ます

上記に対応して、溶媒和計算をすると、

カチオンーカチオンの寄与が負で、カチオンーアニオンの寄与が正となることが結構あります

そこで、系のサイズを変えた場合に、

20w%NaClaq系でのsolvation energy および solvation free energy はどうなりましたか？

松林

[Masaru Matsugami]

松林先生

ご連絡が遅くなり、申し訳ありません。

詳しくご解説頂き、ありがとうございました。

> イオンーイオンの対エネルギーが０にピークを持たず、

> しかし、系を大きくすれば、ピークが０に近づいて行くことは正しい挙動です。

> Ewaldは、対エネルギーを全MD cellで積分したら０になるように定式化されています。

> ですので、距離が近いときにマイナスのエネルギーであれば、遠距離でプラスになり、

> 距離が近いときにプラスであれば、遠距離でマイナスです。

> つまり、同種イオンの場合に、小さいマイナスエネルギーにピークが出て、

> 異種イオンで、小さいプラスエネルギーにピークが出ます

イオンーイオンのエネルギー分布関数の振る舞いとしては正しい挙動とのこと了解いたしました。

> そこで、系のサイズを変えた場合に、

> 20w%NaClaq系でのsolvation energy および solvation free energy はどうなりましたか？

こちらは、全系列で計算がまだ終わっておりません（マシンの数とパワーがないので・・・）。

終わり次第、ご報告させて頂きます。

松上

2017年7月13日木曜日 15時51分02秒 UTC+9 Nobuyuki MATUBAYASI: