BioSupercomputing Newsletter Vol.7

Home -> Newsletter -> Vol.7

研究報告
水の誘電率計算から得られる
古くて新しい問題

中村 春木

大阪大学蛋白質研究所
中村 春木
(分子スケールWG)

 「水は真空に比べて80倍ほどの大きな誘電率を持ち、大きな分極効果が観測される」と教科書には書かれています。この純水の誘電率の分子シミュレーションによる計算には長い歴史があり、最近になって問題点がクリアーになったこともあります。この古くて新しい問題を振り返ることにより、分子だけでなく、シミュレーション計算全般において我々が注意すべきことがあるようにも思い、以下に紹介をいたします。
 まず、ある純水の系の誘電率εは、その系にある各水分子{i} の電気双極子モーメント {} の総和の統計平均により、
式3
として算出されます。ここで、μ0 は一つの水分子がもつ電気双極子モーメントでありGK はKirkwood因子と呼ばれ、
式1
で定義されるスカラ量です。一方、Distance dependent Kirkwood因子と呼ばれる
式2
は、動径分布関数に比べ極めて敏感に水分子の配向と構造を反映する重要な指針を与えます(図1〈表紙〉)。ある水分子を取り巻く第1層目、第2層目の水分子の電気双極子モーメントとの相関、さらに離れた層との相関が見て取れます。式[1]第2項は十分大きな統計をとるとゼロとなる部分ですので、 r を無限大にした時の式[2]のGK(r) が、誘電率に対応するGK に一致します。すなわち、図1の右端の十分遠方の部分での値が誘電率に対応し、標準的な計算法である周期境界条件を使ったPME(particle-mesh Ewald)法での1気圧、300Kでの誘電率の値は96程度になります。 80からのずれは、用いたTIP3Pの水分子モデル1)に依存すると考えられます。
 ところで、液体の純水では水素結合ネットワークができており、そのダイナミクスとしての緩和は遅いため、誘電率の値が収束するには図2(a)〈表紙〉に見られるように1~2ns程度の短い計算では全く足らず、少なくとも6ns程度以上の長いシミュレーションが必要です。しかし、この現象は2011年のGereben & Pusztai論文2)により初めて系統的に指摘されたもので、それ以前の他の論文における短い計算時間での値は信用なりません。私どもも自ら計算していて、数nsのトラジェクトリでは同一の計算手法であっても誘電率の値が異なることが多く、悩んでいた問題でもあったため、この論文の指摘によって極めて明快に問題がクリアーできました。さらに、式[1]のGK の第2項についても図2(b)〈表紙〉のように1~2ns程度では時間平均としての値が無視できるほど小さくなっていないことがわかりました。以上は、計算資源が豊富になったために明らかにされた点です。
 一方、もう一つの問題として、純水における遠距離的な静電力の取り扱いがあげられます。今ではPME法が標準的な手法となっていますが、以前は計算資源の問題から、あるカットオフ距離 dc 内の水分子同士の相互作用のみを考慮する方法がよく使われていました。その際、単純に考えるとなるべく大きな dc をとれば良いと思われますが、Yonetaniの論文3)により、たとえ dc =18Åと大きくしてもGK(r) は正負に一桁ほども大きく振れてずれてしまい、単に定量的に誘電率が異なるだけでなく、定性的な水の構造も正しく再現できないことが示されています。
 歴史的には、Neumann4)によって、 dc を半径とする球の外側に誘電率 εRF の誘電体があるとし、そこからFröhlichによる反作用場(Reaction field)5)を受けるとしたReaction field法によって、上記のカットオフによるartefactを取り去る手法が提案されています。この手法はその後多くの研究者によって試みられましたが、蛋白質水溶液のような均一でない系では、パラメータとして与えるべき εRF の見積もりが困難なこともあり、近年は盛んではありません。一方、最近、理研の福田育夫博士らによりNon-Ewald法の一つとして提案されたZero-dipole summation(ZD) 法6-8)は、Wolfによる電荷の中性条件9)だけでなく電気双極子モーメントの中性条件も課すことによって遠距離力の効果を繰り込み、簡潔なアルゴリズムでありながら高い精度の計算を実現する優れた手法です。この手法を用いると、 dc =12~14Åという常識的には近距離的な相互作用しか考えていないような短い距離で相互作用をカットオフしても、図1および図2に示されるように、PME法とほとんど同一の誘電率やGK(r)が得られます。興味深いことに、このZD法は、ある条件ではReaction field法における εRF→∞の場合と全く同じ式を与えることが示されており、また最近提案されているその他の種々のNon-Ewald法とも共通の性質があることがわかっています8)。周期境界条件を課さない手法には、基盤となる共通の物理があるものと思われます。
 ところで、計算科学では、既に定まったアルゴリズムを基にどれだけ高速に計算するかが競われることも多いと思います。もちろんそれは必要なことですが、アルゴリズムあるいはモデルそのものから考え直すことによって、全く新たな世界が開ける可能性もあると思います。現在、私たちは、上記したZD法を、周期境界系としてではなく3次元トーラス系としての蛋白質やDNAの水溶液などのヘテロな系に対して応用する研究に取り組んでおり、それなりの良い結果が得られ始めています。
 新しいアルゴリズムやモデルを用いる計算の研究では、「天動説」に凝り固まったレビューアとの戦いが必ずおきて、論文を出版する際には苦労します。しかし、研究が成功した場合の波及効果も大きく、実際、周期境界系を使わずに計算ができることは、多くの生体超分子の高速シミュレーションを、より容易にかつ計算資源もより少なく具現化できることにつながります。
 最後に、ここで紹介した研究は、福田育夫博士(理研)、神谷成敏博士(阪大蛋白研)、米澤康滋博士(近畿大)との共同研究によるものです。皆様に感謝いたします。

【参考文献】

1) W. L. Jorgensen et al., J. Chem. Phys . 79, 926 (1983); 2) O. Gereben, L.Pusztai, Chem. Phys. Lett. 507, 80 (2011); 3) Y. Yonetani, J. Chem. Phys. 124,204501 (2006); 4) M. Neumann, Mol. Phys. 50, 841 (1983); 5) H. Frohlich,“Theory of Dielectrics” Clarendon Press (1958); 6) I. Fukuda et al., J. Chem.Phys. 134, 164107 (2011); 7) I. Fukuda et al., J. Chem. Phys. 137, 054314 (2012);8) I. Fukuda, H. Nakamura, Biophys. Rev. 4, 161 (2012); 9) D. Wolf et al., J. Chem.Phys. 110, 8254 (1999)

BioSupercomputing Newsletter Vol.7

SPECIAL INTERVIEW
「京」開発担当者に聞く今後のスパコン戦略とエクサスケールに向けた取り組み
富士通株式会社 テクニカルコンピューティングソリューション事業本部
エグゼクティブ・アーキテクト 奥田 基
「京」への展開で、実利用に向けた最適化とさらなる規模の拡張を進める大規模仮想ライブラリ
東京大学大学院工学系研究科 化学システム工学専攻 教授 船津 公人
研究報告
水の誘電率計算から得られる古くて新しい問題
大阪大学蛋白質研究所 中村 春木(分子スケールWG)
大規模並列計算用流体・構造連成解析プログラムの開発
理化学研究所 情報基盤センター 杉山 和靖(臓器全身スケールWG)
スーパーコンピュータを用いた大規模遺伝子ネットワーク推定ソフトウェア SiGN
東京大学大学院情報理工学系研究科 玉田 嘉紀(データ解析融合WG)
ISLiM研究開発ソフトウェアのソース・コード公開に向けた活動
理化学研究所 次世代計算科学研究開発プログラム 田村 栄悦
SPECIAL INTERVIEW
「京」を用いた大規模シミュレーションによって細胞内分子ダイナミクスの理解と予測を実現する
理化学研究所 基幹研究所 杉田理論分子科学研究室 主任研究員 杉田 有治(課題1 代表)
日本の優れたコンピュータ技術を活かして革新的な分子動力学創薬に挑戦
東京大学 先端科学技術研究センター 特任教授 藤谷 秀章(課題2 代表)
報告
大学の新入生に行った計算生命科学の講義
理化学研究所 HPCI計算生命科学推進プログラム 鎌田 知佐
体制
計算科学技術推進体制
イベント情報