JP2005050321A - プログラム及びそれを用いた計算方法 - Google Patents

Abstract

【課題】 高分子中のプロトンの拡散係数を計算するに際し、高精度かつ効率のよい計算手法を提供する。
【解決手段】 本発明のプログラムは、原子の種類、座標、結合情報、ポテンシャルパラメータタイプ、電荷、質量に関する情報をメモリに格納する手順と、原子対の原子の番号を表に登録する第１の登録手順と、プロトンと水分子の番号を表に登録する第２の登録手順と、プロトンと水の結合パターンを表に登録する第３の登録手順と、原子間に働く古典的な相互作用を計算する第１の計算手順と、半量子的な相互作用を計算する第２の計算手順と、前記第１の計算手順の結果および第２の計算手順の結果から全ての原子の座標を更新する第３の計算手順と、プロトンと水分子中の水素を区別しプロトンの座標を決定する第４の計算手順と、拡散係数を計算する第５の計算手順を少なくとも実行させることを特徴とするものである。
【選択図】なし

Description

本発明は、電子計算機を利用したＭＤシミュレーションを用いての高速かつ高精度に高分子中のプロトンの拡散係数を計算するプログラム及びそれを用いた計算方法に関する。

近年の電子計算機の速度の飛躍的向上と方法論の発展により、分子動力学法（Ｍｏｌｅｃｕｌａｒ Ｄｙｎａｍｉｃｓ ｓｉｍｕｌａｔｉｏｎ：以下ＭＤと略す）においては、液体、高分子、タンパク質等のモデルについて、実験では困難な、原子レベルでの詳細な構造や運動に関する信頼性のある知見を得ることに成功を納めつつある。このＭＤ法とは、ニュートンの運動方程式を数値積分することにより、各原子の位置を時間に対して追跡するというシミュレーション手法である。

高分子中のプロトンの拡散係数をＭＤ法を用いて予測する方法としては、古典力学に基づく古典ＭＤ計算を行う方法と、量子力学に基づく量子ＭＤ計算（第一原理ＭＤ計算）を行う方法の２種類の方法を挙げることができる。

古典ＭＤ計算は原子に働く力をパラメータとして入力する手法であり、計算速度および汎用性に優れた手法であるが、この方法を用いてプロトンの拡散係数を計算すると、許容し難い計算誤差が生じるという問題点があった。古典ＭＤ計算を用いて高分子中のプロトンの拡散係数を予測した例としては、Ｅｎｎａｒｉらによる研究（非特許文献１〜３）が挙げられるが、発明者らが追試したところではプロトンの拡散係数について実験値を全く再現できないことがわかった。このことから、プロトンの拡散係数を精度よく計算するためには、量子効果（プロトンホッピング機構での輸送）を考慮する必要があるといえる。

一方、量子ＭＤ計算はシュレディンガー方程式を分子軌道法や密度汎関数法等を用いて解くことにより、経験的パラメータを使用することなく、高精度の物性予測を可能とする手法である。この方法は非常に小さな系を詳細に解析するには適しているが、高分子／液体／イオン混合系のような複雑系を対象とした場合、莫大な計算時間を必要とし、全く実用的ではなかった。

近年、量子効果を取り入れたＭＤ計算を効率化、高速化するためのアルゴリズム（半量子ＭＤ計算）がいくつかの研究グループによって研究されている。半量子ＭＤ計算とはまず最初に原子座標（反応座標）に対するポテンシャル曲面を分子軌道法や密度汎関数法によって求めておき、得られたポテンシャル曲面に沿って多体効果を考慮して原子を運動させる手法である。あらかじめ、ポテンシャル曲面を求めておくため、量子ＭＤ法による全電子計算と比較して計算速度の点で有利である。水中のプロトンの拡散性を半量子ＭＤ計算で扱った例が既にいくつかの研究グループによって報告されている。

例えば、Stillingerらのグループは水をＯ^2-とＨ⁺に分けて考慮し、経験ポテンシャルを当てはめたＰＭ６モデルを提案しているが（非特許文献４〜６）、この手法は水中のプロトン移動をシミュレートすることに特化した手法であるため、第３成分を含む系の計算をする場合にはその都度、ポテンシャルを再構築する必要があり、汎用性に劣るという問題点があった。

一方、Borgisらのグループ、VothらのグループはＨ₅Ｏ₂ ⁺クラスターについてモデルポテンシャルを構築し （ＥＶＢモデル：empirical valence bond）、プロトン移動のメカニズムを調べている（非特許文献７〜１０）。この手法は、計算アルゴリズムに非効率的な部分が多く、さらに、ポテンシャル関数が複雑であるため、古典ＭＤ法と比較した場合、計算速度および汎用性が劣っており、この方法を高分子中のプロトン拡散係数の計算に適用することは困難であった。

最近になって、Tuckermanらによって開発された手法（非特許文献１１）はBorgis 、Vothらの手法に類似の方法であるが、ポテンシャル関数が単純であるため、Borgis 、Vothらの手法と比較して拡張性に優れた手法である。また、TuckermanらはＨ₅Ｏ₂ ⁺クラスターについて、高精度のＭＰ２／６−３１Ｇ^**レベルの分子軌道計算を行ってポテンシャル曲面を決定し、このポテンシャル曲面を再現するように、ＥＶＢモデルのパラメータを決定しているため、この手法は計算精度にも優れているといえる。しかしながら、この手法はＨ₅Ｏ₂ ⁺用に開発されているため、Ｈ⁺（Ｈ₂Ｏ）_nクラスター／高分子の混合系のダイナミクスの解析には適していないという問題点があった。
J. Ennari, M. Elomaa, F. Sundholm, Polymer, 40, 5035 (1999). J. Ennari, M. Elomaa, I. Neelov, F. Sundholm, Polymer, 41, 985 (2000). J. Ennari, Lars-Olof Pietil, V. Virkkunen, F. Sundholm, Polymer, 43, 5427 (2002). F. H. Stillinger and C. W. David, J. Chem. Phys., 69, 1473 (1978). F. H. Stillinger, J. Chem. Phys., 71, 1647 (1979). T. A. Weber and F. H. Stillinger, J. Phys. Chem., 86, 1314 (1982). U. W. Schmitt and G. A. Voth, J. Phys. Chem. B, 102, 5547 (1998). U. W. Schmitt and G. A. Voth, J. Chem. Phys., 111, 9361 (1999). R. Vuilleumier, D. Borgis, Chem. Phys. Lett., 284, 71 (1998). R. Vuilleumier, D. Borgis, J. Phys. Chem. B, 102, 4261 (1998). D. E. Sagnella, M. E. Tuckerman, J. Chem. Phys., 108, 2073 (1998).

前記のように、古典ＭＤ法は計算精度に問題があり、量子ＭＤ法は計算速度に問題がある。また、半量子ＭＤ法は汎用性に乏しく、Ｈ⁺（Ｈ₂Ｏ）_nクラスター／高分子の系を取扱うことが困難であった。さらに、半量子ＭＤ法は量子ＭＤ法と比較すると計算速度に優れているが、古典ＭＤ法と比較すると大幅に劣っているため、さらなる計算速度向上を図る必要があった。

本発明は、従来技術の背景に鑑み、高分子中のプロトンの拡散係数を計算する際に量子効果を取り入れて高精度計算を可能とし、かつ効率のよいプログラム及びそれを用いた計算方法を提供せんとするものである。

本発明は、上記の課題を解決するために、次のような手段を採用するものである。すなわち、本発明のプログラムは、原子の種類、座標、結合情報、ポテンシャルパラメータタイプ、電荷、質量に関する情報をメモリに格納する手順と、一定距離内にある原子対の原子の番号を表に登録する第１の登録手順と、前記原子対表を用いて一定距離内にあるプロトンと水分子の番号を表に登録する第２の登録手順と、前記プロトン・水分子表からプロトンと水の結合パターンを表に登録する第３の登録手順と、前記原子対表から原子間に働く古典的な相互作用を計算する第１の計算手順と、前記プロトンと水の結合パターン表から半量子的な相互作用を計算する第２の計算手順と、前記第１の計算手順の結果および第２の計算手順の結果から全ての原子の座標を更新する第３の計算手順と、更新された原子の座標からプロトンと水分子中の水素を区別しプロトンの座標を決定する第４の計算手順と、プロトンの座標の軌跡から拡散係数を計算する第５の計算手順を少なくとも実行させることを特徴とするものである。また、本発明の計算方法は、かかるプログラムを格納した電子計算機を用いて高分子中のプロトンの拡散係数を求める高分子中のプロトンの拡散係数を計算することを特徴とするものである。

本発明の新規なプロトンの拡散係数計算方法により、実験を行うことなく、高精度で高速に高分子中のプロトンの拡散係数を計算することができ、高分子電解質膜等の研究開発においてスピードアップ、ブレークスルーが期待できる。

高分子中のプロトンの拡散係数を計算するにあたって、系を構成する原子に作用する力を全て半量子的に扱うと多大な計算コストを必要とするが、本発明者らは、プロトンとプロトン近傍にある水分子との相互作用は半量子的に取り扱い、それ以外の相互作用は古典的に取扱うことにより、計算速度を向上させることができることに気付き、さらに、計算アルゴリズムの効率化、計算に必要となるポテンシャルパラメータの開発を行って本発明に到達したものである。

図１に本発明のフローを示した。図１において、１は初期構造入力処理ステップ、２は原子対表作成処理ステップ、３はプロトン・水分子表作成処理ステップ、４は結合パターン表作成処理ステップ、５は古典的相互作用計算処理ステップ、６は半量子的相互作用計算処理ステップ、７は原子座標更新処理ステップ、８はプロトン判別処理ステップ、９は拡散係数計算処理ステップを表わす。これらの処理ステップはそれぞれ図１に示したように第１〜３の登録手順および第１〜５の計算手順を実行させるプログラムである。

初期構造入力処理ステップは原子の種類、座標、結合情報、ポテンシャルパラメータタイプ、電荷、質量に関する情報をメモリに格納する手順を実行させるものである。この読み込みには当該情報が記述されたファイルを読み込むことなどで行うことができる。

原子対表作成処理ステップは、図２の点線で示したような一定距離以内にある原子対の原子の番号をメモリ上の表に登録する第１の登録手順を実行させるプログラムである。このプログラムによって実行させる処理を以下に詳述する。原子ｉと原子ｊの距離Ｒ_ijを計算し、Ｒ_ijが一定距離Ｒ_listよりも小さければ、原子の番号ｉとｊをセットで原子対表に登録する。この操作を全ての原子対の組み合わせ（ｉ＝１〜Ｎ−１、ｊ＝ｉ＋１〜Ｎ）について行う。ただし、Ｎは系の全原子数を表わす。ここで、原子対登録の判断基準となる距離Ｒ_listは計算精度を考慮すると１１オングストローム以上であることが好ましく、計算速度を考慮すると１６オングストローム以下であることが好ましい。

プロトン・水分子表作成処理ステップは、前記原子対表を用いて、図２の実線で示したような一定距離内にあるプロトンと水分子の番号をメモリ上の表に登録する第２の登録手順を実行させるものである。このプログラムによって実行させる処理を以下に詳述する。プロトンｉと水分子ｊの酸素原子の距離Ｒ_ijを計算し、Ｒ_ijが一定距離Ｒ^W _listよりも小さければ、プロトンの番号ｉと水分子の番号ｊをセットで水分子表に登録する。この操作を前記原子対表に登録されたプロトンと水の全ての組み合わせについて行う。ここで、登録の判断基準となる距離Ｒ^W _listは前記Ｒ_listよりも小さい必要があるが、この距離Ｒ^W _listは計算精度を考慮すると６オングストローム以上であることが好ましく、計算速度を考慮すると１２オングストローム以下であることが好ましい。

結合パターン表作成処理ステップは、前記プロトン・水分子表を用いて、図３で例を示したようにプロトンと水の結合パターンをメモリ上の表に登録する第３の登録手順を実行させるものである。登録手順３においてはプロトンの拡散係数を計算するためには少なくとも３つの酸素原子間をホッピングするプロトンを考慮する必要があるが、対称性を考慮するとプロトン近傍にある６個以上の水分子について結合パターンを登録することが好ましい。一方、計算速度を考慮すると半量子的に取扱う水分子数が２０個以下について結合パターンを登録することが好ましい。また、系に複数のプロトンが存在する場合には、計算速度を考慮すると１個のプロトンのみを半量子的に取り扱うとして結合パターン表に登録し、それ以外のプロトンは古典的に取り扱うとして結合パターン表には登録しないことが好ましい。

古典的相互作用計算処理ステップは、前記原子対表を用いて、古典的な相互作用を計算する第１の計算手順を実行させるものである。古典的相互作用の計算には数式（８）で示したポテンシャル関数を用いる。

数式（８）

ここで、数式（８）の右辺第１項は結合長のポテンシャル関数を表わし、第２項は結合角のポテンシャル関数を表わし、第３項は二面角のポテンシャル関数を表わし、第４項はｉｎｖｅｒｓｉｏｎのポテンシャル関数を表わし、第５項はｖｄＷポテンシャル関数を表わし、第６項は静電ポテンシャル関数を表わす。なお、数式（８）の第１項から第４項を総称して分子内ポテンシャル関数と呼び、第５項、第６項を分子間ポテンシャル関数と呼ぶ。数式（８）中のパラメータは高分子については、ＡＭＢＥＲ［W. D. Cornell, P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz Jr, D. M. Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwell and P. A. Kollman, J. Am. Chem. Soc., 117, 5179 (1995)]、ＤＲＥＩＤＩＮＧ[S.L. Mayo, B.D. Olafson, W.A. Goddard III, J. Phys. Chem., 94, 8897 (1990)］等の学術的に信頼性の高い公知のパラメータを用いることが好ましいが、計算精度について考慮した場合、分子内ポテンシャル関数のパラメータについては分子軌道法によって決定し、分子間ポテンシャル関数のパラメータについては密度と蒸発熱の実験値を再現するように決定することがさらに好ましい。特に、第３項の二面角ポテンシャルパラメータは高分子の構造を決定する重要な因子であるため、ＲＨＦ／６−３１Ｇ^**レベルの分子軌道計算［R. Ditchfield, W. J. Hehre, J. A. Pople, J. Chem. Phys., 54, 724 (1971)］によって決定したパラメータを用いて計算を行うことがさらに好ましい。また、第６項の電荷パラメータについては、ＱＥｑ法［A. K. Rappe and W. A. Goddard III, J. Phys. Chem., 95, 3358-3363 (1991)］等の学術的に信頼性の高い手法を用いて決定することが好ましく、計算精度を考慮すると、Ｂ３ＬＹＰ／６−３１Ｇ^**レベルの分子軌道法によって決定したパラメータを用いて計算を行うことがさらに好ましい。Ｈ₅Ｏ₂ ⁺のｖｄＷポテンシャルパラメータについてはＴｕｃｋｅｒｍａｎらによって報告されているパラメータを用いることが好ましい。Ｈ₅Ｏ₂ ⁺の電荷パラメータについては、プロトン移動に伴って値が変化するパラメータを用いることが好ましいが、例えば、ＷｅｉとＳａｌａｈｕｂ［D. Wei and D. R. Salahub, J. Chem. Phys., 101, 7633 (1994)］によって開発された分極の効果を取り入れた方法やＶｕｉｌｌｅｕｍｉｅｒとＢｏｒｇｉｓ［R. Vuilleumier, D. Borgis, J. Mol. Struct., 552, 117 (2000)］の重み付き平均電荷等を用いることで好適に電荷の値を変化させることができる。また、数式（８）の第６項で示した静電相互作用はＥｗａｌｄの方法［M. P. Tosi, In F. Seitz and D. Turnbull, editors, "Solid St. Phys. Advances in Research and Applications", Vol. 16 1-120. Academic Press New York and London (1964)］により無限遠までの計算を行うことが好ましい。これは、静電気力が遠達力であり、静電ポテンシャルが非常に重要で全ポテンシャルをほとんど支配しているため、より正確な計算が必要とされるからである。なお、Ｅｗａｌｄ法を用いる場合には、ガウス型曲線の広がり幅を決定するパラメータαおよび逆格子ベクトルの最大値｜ｎ｜² _maxを入力する必要があるが、千〜数千原子程度の系を計算対象とした場合には、計算精度を考慮するとαは０．２０オングストローム^-1以上、｜ｎ｜² _maxは２０以上とすることが好ましく、計算速度を考慮するとαは０．３０オングストローム^-1以下、｜ｎ｜² _maxは５０以下とすることが好ましい。ｖｄＷ相互作用および静電相互作用の実空間のカットオフ半径は計算精度を考慮するとｒ_cは１０オングストローム以上とすることが好ましく、計算速度を考慮するとｒ_cは１５オングストローム以下とすることが好ましい。

半量子的相互作用計算処理ステップは、前記結合パターン表を用いて、原子に働く半量子的相互作用を計算する第２の計算手順を実行させるものである。半量子的相互作用の計算には数式（１）〜（７）で示したポテンシャル関数を用いる。

数式（１）

数式（２）

数式（３）

数式（４）

数式（５）

数式（６）

数式（７）

数式（１）において、係数ａ_iはＨの固有ベクトル、Ｖ_ijは行列Ｈのｉｊ成分であるが、対角成分Ｖ_ijを数式（２）〜（６）式で定義し、非対角成分Ｖ_ijを数式（７）で定義した。ここで、数式（３）の第１項はＨ₃Ｏ⁺分子内の酸素−水素の結合長ポテンシャル関数を表わし、第２項はＨ₂Ｏ分子内の酸素−水素の結合長ポテンシャル関数を表わし、数式（４）の第１項はＨ₃Ｏ⁺分子内の水素−酸素−水素の結合角ポテンシャル関数を表わし、第２項はＨ₂Ｏ分子内の水素−酸素−水素の結合角ポテンシャル関数を表わし、数式（５）の第１項はＨ₃Ｏ⁺中の酸素原子とＨ₂Ｏ中の酸素原子間のｖｄＷポテンシャル関数を表わし、第２項はＨ₂Ｏ中の酸素原子間のｖｄＷポテンシャル関数を表わし、数式（６）の第１項はＨ₃Ｏ⁺中の酸素および水素原子とＨ₂Ｏ中の酸素および水素原子の静電ポテンシャル関数を表わし、第２項はＨ₂Ｏ中の酸素および水素原子とＨ₂Ｏ中の酸素および水素原子の静電ポテンシャル関数を表わし、数式（７）はプロトンおよび水分子の相互作用を表わす。なお、数式（３）〜（７）は半量子的に取扱うプロトンが１個である場合に用いる関数であるが、複数のプロトンを半量子的に取扱う場合には、数式（３）〜（７）について半量子的に取扱うプロトンの数でさらに和をとる必要がある。ここで、数式（３）の第１項、数式（４）の第１項、数式（５）の第１項のポテンシャル関数中のパラメータはＴｕｃｋｅｒｍａｎらによって開発されたパラメータを用い、数式（５）の第２項のパラメータおよび数式（６）のパラメータ（ｑ）は水分子集合体に対して最適化されたＴＩＰ３Ｐパラメータ［W. L. Jorgensen, R. W. Impey, J. Chandrasekhar, J. D. Madura, and M. L. Klein, J. Chem. Phys., 79, 926 (1983)]またはＳＰＣ／Ｅ[H. J. C. Berendsen, J. R. Grigera, T. P. Staatsma, J. Phys. Chem., 91, 6269 (1987)］等を用いることが好ましい。また、数式（３）の第２項のパラメータＫ_bおよび数式（４）の第２項のパラメータＫθは表１に示される範囲の値で規定して計算を行うことが好ましい。この範囲のパラメータを用いることにより、水のＩＲスペクトルを好適に再現することができる。なお、数式（３）の第２項のパラメータｒ₀および数式（４）の第２項のパラメータθ₀は前記のＴＩＰ３Ｐパラメータ、ＳＰＣ／Ｅパラメータ等を用いることが好ましい。また、数式（６）の第１項および数式（７）の第２項のパラメータは表１に示される範囲の値を用いて計算を実行することが好ましい。この範囲のパラメータを用いることで、Ｔｕｃｋｅｒｍａｎらが分子軌道法を用いて計算したＨ₅Ｏ₂ ⁺のポテンシャル曲面を好適に再現することができる。なお、数式（７）の第１項のパラメータはＴｕｃｋｅｒｍａｎらによって開発されたパラメータを用いることが好ましい。

原子座標更新処理ステップは前記計算手順１および計算手順２の結果を用いて、微小時間後の原子の座標を数値積分によって計算する第３の計算手順を実行させるものである。数値積分には異なる時定数の力ごとに大きさの異なる時間刻みを設定することのできる多重時間刻み数値積分法［ＲＥＳＰＡ法：G. J. Martyna, Mol. Phys., 87, 1117(1996)］）を用いて計算を行うことが好ましい。この方法を用いることによって計算効率を大幅に改善することができるが、時間刻みを表２に示される範囲の値で規定して計算を行うことがさらに好ましい。この範囲の時間刻みを用いることによって、運動性が速く計算コストの比較的小さい分子内力の計算頻度と運動性が遅く計算コストの大きい分子間力の計算頻度をそれぞれ最適にすることができるため、計算精度を維持したまま計算時間を短縮することができる。

プロトン判別処理ステップは原子の座標データを用いて、水分子中の水素原子とプロトンの座標を決定する第４の計算手順を実行させるものである。ある時刻ｔにおけるプロトンと水素を区別する方法としては以下のような方法が好適である。まず、１シミュレーションステップ前の時刻ｔ−Δｔにおいて、原子間の距離を計算し、図４に示したような便宜上の番号を原子に付け、原子を区別しておく。図４に示した例では、プロトンを１番とし、プロトンから最も近い水分子の酸素原子を２番、水素原子を３、４番とした。また、プロトンから２番目に近い水分子の酸素原子を５番とした。さらに、３番の水素原子から最も近い水分子の酸素原子を６番とし、４番の水素原子から最も近い水分子の酸素原子を７番とした。次に、時刻ｔにおいて、ｒ１（１番−５番の距離）、ｒ２（３番−６番の距離）、ｒ３（４番−７番の距離）を比較して最も距離の小さいものをプロトンとみなし、プロトンの座標を決定する。

拡散係数計算処理ステップは前記プロトンの座標の軌跡から、プロトンの拡散係数を計算する第５の計算手順を実行させるものである。拡散係数の計算は数式（９）に従って平均二乗変位の時間に対する傾きから、拡散係数を計算する。ここで、平均二乗変位の計算は計算精度を考慮すると、プロトン座標のデータ量として３００ｐｓ以上のデータを用いて計算を行うことが好ましい。なお、拡散係数の計算は、計算精度を考慮すると、１ｐｓ以上の平均二乗変位のデータを用いることが好ましく、モデルのサイズ依存性を回避するためには２０ｐｓ以下の平均二乗変位のデータを用いることが好ましい。

数式（９）

（ただし、式中のＤは拡散係数：単位ｃｍ²／ｓ、ｒは原子の座標：単位ｃｍ、ｔは時刻：単位ｓを表わす）
本発明の新規なプロトンの拡散係数計算方法により、実験を行うことなく、高精度で高速に高分子中のプロトンの拡散係数を計算することができる。さらに、計算したプロトンの拡散係数を数式（１０）に代入することでプロトン伝導度の計算が可能である。

数式（１０）

（ただし、式中のκは電気伝導度：単位Ｓ／ｃｍ、ｃはプロトン密度：単位ｍｍｏｌ／ｃｍ³、ｚはプロトンの電荷：単位ｅ、Ｄは拡散係数：単位ｃｍ²／ｓ、Ｔは温度：単位Ｋ、Ｆはファラデー定数：単位Ｃ／ｍｏｌ、Ｒは気体定数：単位Ｊ／ｍｏｌ・Ｋを表わす）
以上説明したように、本発明によってプロトンの拡散係数やプロトン伝導度をシミュレーションで予測することが可能となったことによって、高分子電解質膜等の研究開発においてスピードアップ、ブレークスルーが期待できる。

本発明の具体的実施態様を以下に実施例を持って述べるが、本発明はこれに限定されるものではない。

実施例１
高分子、水、プロトンを混合したモデルを構築し、高分子中のプロトンの拡散係数を本発明の半量子的相互作用を取り入れた分子動力学法によって計算した。高分子モデルには現在、燃料電池用電解質ポリマーとして主流となっているＤｕＰｏｎｔ社のＮａｆｉｏｎを用いた。なお、分子動力学計算は東レ株式会社で独自に開発したＦＯＲＴＲＡＮプログラムを用いて行い、分子軌道計算は汎用プログラムＧＡＵＳＳＩＡＮ９８を用いて行った。また、計算にはＰｅｎｔｉｕｍ（登録商標）４／２．４ＧＨｚを用いた。

以下に計算モデルおよび計算方法について詳述する。プロトンの拡散係数は以下のステップＡ〜Ｅの一連の処理によって計算した。

ステップＡ：計算モデル作成
Ａｃｃｅｌｒｙｓ社製の分子設計システムＣｅｒｉｕｓ２上で３Ｄ−Ｓｋｅｔｃｈｅｒを用いて構造式（１）に示すポリマーを作成した。

構造式（１）

なお、構造式（１）において、ｎ＝５とした（分子量＝４７５３、スルホン酸基の割合＝１．０５：単位ｍｍｏｌ／ｇ）。次いで、Ａｍｏｒｐｈｏｕｓ Ｂｕｉｌｄｅｒを用いてモンテカルロ計算を行い、水、プロトンを３次元周期境界条件下にランダムに配置して系を構成した。なお、表３に示した組成で複数のモデルを作成し、Ｎａｆｉｏｎ・水・プロトンで構成される分子集合体構造ファイルを出力した。

表３

次に、Ｎａｆｉｏｎ・水・プロトンで構成される分子集合体構造ファイルを前記初期座標入力処理ステップにて読み込み、原子の種類、座標、結合情報、ポテンシャルパラメータタイプ、電荷、質量をメモリに格納した。

ステップＢ．各表の作成
前記登録手順１〜３に従ってそれぞれ原子対表、プロトン・水分子表、結合パターン表をメモリに登録した。ここで、登録手順１においては、原子対表を作成するための原子間距離は１１オングストロームとし、登録手順２においてはプロトン・水分子表を作成するための原子間距離は７オングストロームとした。また、登録手順３においては、１個のプロトンのみを対象として、該プロトンとその近傍にある６個の水分子について結合パターンをメモリに格納した。

ステップＣ：相互作用計算
計算手順１に従って古典的相互作用を計算した。ここで、ｖｄＷ相互作用および実空間の静電相互作用の計算はカットオフ半径ｒ_c＝１０オングストロームとし、逆空間の静電相互作用はα＝０．２１オングストローム^-1、｜ｎ｜² _max＝３０としてＥｗａｌｄ法を用いて計算を行った。古典的相互作用計算に用いたポテンシャルパラメータは、Ｎａｆｉｏｎ分子については結合長および結合角の力の定数（Ｋ^bond，Ｋ^angle）はＤＲＥＩＤＩＮＧパラメータを用い、平衡位置（Ｒ⁰，θ⁰）はＢ３ＬＹＰ／６−３１Ｇ^**レベルの分子軌道計算によって決定したパラメータを用い、主鎖の二面角ポテンシャルパラメータ（Ｋ^torsion，ｎ，φ⁰）はＲＨＦ／６−３１Ｇ^**レベルの分子軌道計算によって決定したパラメータを用い、側鎖にはＤＲＥＩＤＩＮＧパラメータを用い、ｉｎｖｅｒｓｉｏｎポテンシャルパラメータ（Ｋ^inversion，ω⁰）にはＤＲＥＩＤＩＮＧパラメータを用いた。また、Ｎａｆｉｏｎ分子のｖｄＷパラメータ（Ｄ⁰，ｒ⁰）にはＣＦ₂部分に文献［P. Jedlovszky, M. Mezei. , J. Chem. Phys., 110, 2991 (1999)］のパラメータを用い、ＳＯ₃ ^-部分には文献［W. R. Cannon, B. M. Pettitt, J. A. McCammon, J. Phys. Chem., 98, 6225 (1994)］のパラメータ、エーテル酸素にはＡＭＢＥＲパラメータを用い、電荷（ｑ）はＢ３ＬＹＰ／６−３１Ｇ^**レベルの分子軌道計算によって決定した値を用いた。水分子については分子間ポテンシャルパラメータにはＴＩＰ３Ｐパラメータを用い、分子内ポテンシャルパラメータは力の定数（Ｋ^bond，Ｋ^angle）には表４に記載の（Ｋ_b，Ｋθ）の値を用い、平衡位置（Ｒ⁰，θ⁰）にはＴＩＰ３Ｐパラメータを用いた。プロトンについては系に存在する１０個のプロトンのうち、１個を半量子的に取り扱い、残りの９個を古典的に取扱った。なお、古典的に取扱うプロトンにはＥｎｎａｒｉらのパラメータを用いた。Ｈ₅Ｏ₂ ⁺の電荷についてはＷｅｉとＳａｌａｈｕｂによって開発された方法を用いて計算した。

表４

（ただし、ｑ_O ⁺、ｑ_H ⁺はそれぞれＨ₃Ｏ⁺中の酸素と水素の電荷を表わす）
次に、計算手順２に従って半量子的相互作用を計算した。半量子的相互作用計算に用いたポテンシャルパラメータについて以下に詳述する。数式（３）の第１項のパラメータ（Ｋ_b ⁺，ａ，ｒ₀ ⁺）および数式（４）の第１項のパラメータ（Ｋθ⁺，θ₀ ⁺）はＴｕｃｋｅｒｍａｎらによって開発されたパラメータを用い、数式（３）の第２項のパラメータ（ｒ₀）、数式（４）の第２項のパラメータ（θ₀）、数式（５）の第２項のパラメータ（Ｄ₀，Ｒ₀）、数式（６）のパラメータ（ｑ）はＴＩＰ３Ｐパラメータを用いた。数式（３）の第２項のパラメータ（Ｋ_b）、数式（４）の第２項のパラメータ（Ｋθ）、数式（６）の第１項のパラメータ（ｑ⁺）、数式（７）の第２項のパラメータ（Ｖ_end）は表４に示したパラメータを用いた。また、数式（５）の第１項のパラメータ（Ｄ₀ ⁺，Ｒ₀ ⁺）および数式（７）の第１項のパラメータ（Ｖ⁰，ｂ，ｒ_eq）はＴｕｃｋｅｒｍａｎらによって開発されたパラメータを用いた。

ステップＤ：原子座標更新
次に、計算手順３に従って微小時間後の原子の座標を計算した。ここで、温度をＮｏｓｅ−Ｈｏｏｖｅｒ法［Ｓ． Ｎｏｓｅ ａｎｄ Ｍ． Ｌ． Ｋｌｅｉｎ， Ｍｏｌ． Ｐｈｙｓ．， ５０， １０５５ （１９８３）／Ｓ． Ｎｏｓｅ， Ｍｏｌ． Ｐｈｙｓ．， ５２， ２５５ （１９８４）／Ｓ． Ｎｏｓｅ， Ｊ． Ｃｈｅｍ． Ｐｈｙｓ．， ８１， ５１１ （１９８４）］を用いて２５℃に制御し、圧力を斜交系セルを用いるＰａｒｒｉｎｅｌｌｏ−Ｒａｈｍａｎ法［Ｍ． Ｐａｒｒｉｎｅｌｌｏ ａｎｄ Ａ． Ｒａｈｍａｎ， Ｊ． Ａｐｐｌ． Ｐｈｙｓ．， ５２， ７１８２ （１９８１）］を用いて１ａｔｍに制御した。数値積分には多重時間刻み数値積分法を用い、時間刻みは表５に示した値を用いた。次に、計算手順４に従って更新された原子の座標データを用いて、水分子中の水素原子とプロトンを判別し、プロトンの座標をファイルに出力した。

表５

ステップＥ：拡散係数計算
ステップＢ〜Ｄを繰り返し実行し、半量子的に取扱ったプロトンの座標の軌跡（５００ｐｓ）を得た。ここで、計算時間を短縮するために、ステップＡにおいて原子対表、プロトン・水分子表を作成する頻度は１００ステップに１回とした。次に、計算手順５に従って前記半量子的に取扱ったプロトンの座標の軌跡から、プロトンの拡散係数を計算した。ここで、平均二乗変位を計算するためのプロトンの座標の軌跡は１００ｐｓ〜５００ｐｓのデータを用い、拡散係数の計算に用いる平均二乗変位のデータは、２ｐｓ〜１０ｐｓのデータを用いた。また、同様の方法で水の拡散係数を計算した。

比較例１
実施例１と同一のモデルについて、古典的相互作用のみによる分子動力学法によってプロトンの拡散係数を計算した。なお、プロトンのポテンシャルパラメータにはＥｎｎａｒｉらのパラメータを用いた。

比較例２
実施例１と同一のモデルについて、多重時間刻み数値積分法を用いることなく、プロトンの拡散係数を計算した。なお、時間刻みは０．０６２５ｆｓに設定した。

比較例３
構造式（１）において、ｎ＝１であるＮａｆｉｏｎモデル分子（分子量＝９８１、スルホン酸基の割合＝１．０５：単位ｍｍｏｌ／ｇ）と水１０分子で構成されるモデルをＡｃｃｅｌｒｙｓ社製の分子設計システムＣｅｒｉｕｓ２を用いて作成し、分子力学計算を行って構造を緩和させた後、汎用プログラムＧＡＵＳＳＩＡＮ９８を用いてＨＦ／３−２１Ｇレベルの分子軌道計算を行い、原子に働く力を計算した。

実施例１の方法による分子動力学計算から算出したプロトンと水の拡散係数を図５（左）にそれぞれ□と○で示し、比較例１の方法によって算出したプロトンの拡散係数を図５（左）に△で示し、これらの計算に対応するＧｏｔｔｅｓｆｅｌｄらのプロトンと水の拡散係数の実験値［X. Ren, T. E. Springer, T. A. Zawodzinski, and S. Gottesfeld, J. Electrochem. Soc., 147, 466 (2000)］を図５（右）にそれぞれ□と○で示した。なお、比較例１の水の拡散係数は実施例１とほぼ等しい値であったため、記載を省略した。実施例１については膜中の水の量を増加させると、水およびプロトン拡散係数が大きくなるが、プロトンの拡散係数の増加率は水と比べると大きいことがわかる。この結果はＧｏｔｔｅｓｆｅｌｄらの実験事実と定性的に一致した。この結果はプロトンがホッピングによって移動するメカニズムが計算で正確に再現されていることを示す。定量性に欠ける原因は系の大きさ（Ｎａｆｉｏｎの分子量、分子数）が小さいことや、水やＮａｆｉｏｎのポテンシャルパラメータにあると考えられるが、電解質膜の構造と低分子の拡散性の相関を解析する際に要求される精度としては十分に高いと判断できる。一方、比較例１ではプロトンの拡散係数は水の拡散係数に比較して小さな値となり、プロトンの拡散係数が水よりも大きいという実験事実を再現することができなかった。これらの結果から、本発明の量子効果を取り入れたプロトンの拡散係数評価方法は従来の方法より優れていることが示された。

実施例１、比較例１〜３の方法によるプロトンの拡散係数の計算に要したＣＰＵ時間を表６に示した。この表から、実施例１は比較例１の古典ＭＤ計算と比べて、約２倍の計算時間が必要であることがわかるが、これは十分に実用可能な計算時間であるといえる。また、比較例２と比較すると約１／５の計算時間であり、多重時間刻み数値積分法で用いた時間刻みが計算の効率化に有効であることが示された。また、比較例３は実施例１と比べて系を構成する原子数が大幅に小さいことにもかかわらず、計算時間は約２００００倍となっていた。このことから、本発明は従来の量子計算法よりも、非常に効率よい計算方法であることが示された。

本発明のプログラムのフローを説明するための図である。 各原子間に働く相互作用を説明するための図である。 プロトンと水の結合パターンを示す図である。 プロトンと周辺の水分子の相互作用を説明するための図である。 プロトン、水の拡散係数の計算値および実験値を示す図である。

Claims (6)

1. 原子の種類、座標、結合情報、ポテンシャルパラメータタイプ、電荷、質量に関する情報をメモリに格納する手順と、一定距離内にある原子対の原子の番号を表に登録する第１の登録手順と、前記原子対表を用いて一定距離内にあるプロトンと水分子の番号を表に登録する第２の登録手順と、前記プロトン・水分子表からプロトンと水の結合パターンを表に登録する第３の登録手順と、前記原子対表から原子間に働く古典的な相互作用を計算する第１の計算手順と、前記プロトンと水の結合パターン表から半量子的な相互作用を計算する第２の計算手順と、前記第１の計算手順の結果および第２の計算手順の結果から全ての原子の座標を更新する第３の計算手順と、更新された原子の座標からプロトンと水分子中の水素を区別しプロトンの座標を決定する第４の計算手順と、プロトンの座標の軌跡から拡散係数を計算する第５の計算手順を少なくとも実行させることを特徴とするプログラム。
2. 第３の登録手順においては、プロトンの近傍にある６分子以上２０分子以下の水分子を対象として、プロトンと水の結合パターンを表に登録することを特徴とする請求項１記載のプログラム。
3. 第３の登録手順においては、系に複数のプロトンが存在する場合に、１個のプロトンのみを対象として、プロトンと水の結合パターンを表に登録することを特徴とする請求項１または２に記載のプログラム。
4. 第２の計算手順においては、数式（１）〜（７）で示したポテンシャル関数を用い、かつ、数式（１）〜（７）において用いるポテンシャルパラメータは、メモリー中に格納された表１に示される範囲の値を用いることを特徴とする請求項１〜３のいずれかに記載のプログラム。
   数式（１）
   数式（２）
   数式（３）
   数式（４）
   数式（５）
   数式（６）
   数式（７）
   表１
   
5. 第３の計算手順は多重時間刻み数値積分法を用い、時間刻みとしてメモリー中に格納された表２に示される範囲の値を用いることを特徴とする請求項１〜４のいずれかに記載のプログラム。
   表２
   
6. 請求項１〜５のいずれか記載のプログラムを格納した電子計算機を用いて、高分子中のプロトンの拡散係数を計算することを特徴とする計算方法。