JP6626665B2 - 輸送係数を算出する方法、装置、及びプログラム - Google Patents

Description

本発明は、分子動力学シミュレーションを用いて粘性係数、拡散係数及び熱伝導係数等の輸送係数を算出する方法、装置、及びプログラムに関する。

解析対象となる物質の拡散係数や粘性係数、熱伝導係数などの輸送係数は、揺動散逸定理によれば、分子動力学計算から得られる物理量の自己相関関数を時間積分することにより得られることが知られている。例えば特許文献１では、輸送係数として熱伝導係数を求める例を挙げている。特許文献１では、解析対象となる物質を原子レベルで再現したモデルデータを用いて分子動力学計算を行い、平衡状態の分子の位置及びその挙動を計算する。その計算結果に基づき熱流束の時間変化を表す時系列データを算出し、熱流束の自己相関関数から熱伝導係数を算出する。

特許文献２は、輸送係数を算出するにあたり、現実には存在しない寄与が分子動力学計算に混入してしまうことに着目して、これを除去する手段を提案している。除去後の自己相関関数をラプラス変換してフィティング関数とし、指数関数で表現される当該フィティング関数を時間積分することで輸送係数を計算することが記載されている（段落００３６参照）。

特開２０１０−１３９５００号公報 特開２０１４−１０６５５４号公報

輸送係数は、物理量の時間相関関数の積分値で与えられる。例えば、せん断応力Ｐｘｙとすれば、粘性係数ηは下記式（１）で表現される。ｋ_Ｂはボルツマン定数、Ｖは体積、Ｔは温度である。＜…＞は自己相関関数を示している。

式（１）に示すように、輸送係数の算出には、自己相関関数の値を時点０から無限大∞まで積分する必要がある。しかし、現実には無限時間の計算は不可能であるので、有限時間で計算を打ち切る必要がある。図４に示すように、自己相関関数は、長時間領域にて揺らぎが生じ、この揺らぎが輸送係数の値を算出する精度に影響を与える。図７に示すように、輸送係数の計算を進めると、図中にて四角で示す部分で概ね一定の値に収束し、その後、揺らぎの影響により悪化（decay）する。そこで、自己相関関数の長時間領域での揺らぎの影響を低減するために、図７に示すように、輸送係数（この例は粘性係数）の値の計算を続け、値が概ね一定値に収束し、その後、悪化することを確認できるまで計算を続け、計算を終了する。輸送係数として採用する値は、或る程度一定値に収束した部分の平均値を採用することが考えられる。

しかしながら、上記方法では、自己相関関数の長時間領域での揺らぎの影響を低減するために、輸送係数の値が或る程度収束し、その後悪化することを確認できるまで積分計算を継続しなくてはならず、十分に長い計算時間が必要で、計算コストが大きくなる。

また、特許文献２のように、自己相関関数の値をラプラス変換して、指数関数で表現されるフィティング関数を用いてフィティングした場合には、指数関数の形状に起因して短時間領域の値を指数関数で精度良く表現できず、その結果、輸送係数の算出精度に影響を与えることが判明した。

本発明は、このような課題に着目してなされたものであって、その目的は、計算コストを低減するとともに、算出精度を向上させた輸送係数を算出する方法、装置、及びプログラムを提供することである。

本発明は、上記目的を達成するために、次のような手段を講じている。

すなわち、本発明の輸送係数を算出する方法は、
解析対象となる物質を表す原子モデルデータを用いた分子動力学計算に基づき、予め定めた解析温度及び解析圧力下における平衡状態にある分子の挙動を算出し、輸送係数に対応する物理量を時点ｔ_１から時点ｔ_ｋまで表すｋ個の時系列データを算出するステップと、
前記時点ｔ_１から時点ｔ_ｋをｍ個（ｍ＜ｋ）のグループに区分し、グループ毎に算出した物理量の平均値に基づきｍ個の自己相関関数の値を算出するステップと、
前記ｍ個の自己相関関数の値を、ｅｘｐ｛−（ｔ／τ）^β｝で表されるＫＷＷ関数（Kohlausch-Williams-Watts）を含む近似式で近似し、τ及びβを含むパラメータを決定するステップと、
決定された前記パラメータを用いて前記近似式を時間積分して輸送係数を算出するステップと、
を含む。

本発明の輸送係数を算出する装置は、
解析対象となる物質を表す原子モデルデータを用いた分子動力学計算に基づき、予め定めた解析温度及び解析圧力下における平衡状態にある分子の挙動を算出し、輸送係数に対応する物理量を時点ｔ_１から時点ｔ_ｋまで表すｋ個の時系列データを算出する物理量算出部と、
前記時点ｔ_１から時点ｔ_ｋをｍ個（ｍ＜ｋ）のグループに区分し、グループ毎に算出した物理量の平均値に基づきｍ個の自己相関関数の値を算出する自己相関関数算出部と、
前記ｍ個の自己相関関数の値を、ｅｘｐ｛−（ｔ／τ）^β｝で表されるＫＷＷ関数（Kohlausch-Williams-Watts）を含む近似式で近似し、τ及びβを含むパラメータを決定する近似部と、
決定された前記パラメータを用いて前記近似式を時間積分して輸送係数を算出する輸送係数算出部と、
を備える。

このように、時点ｔ_１から時点ｔ_ｋをｍ個（ｍ＜ｋ）のグループに区分し、グループ毎に算出した物理量の平均値に基づき自己相関関数の値を算出するので、平均化により長時間領域の揺らぎを抑制でき、揺らぎの影響を確認するために必要以上の時間を計算する必要がないので計算コストを低減することが可能となる。その反面、物理量の平均値を取ることで、時系列データの数が減少するので、線形スケールでの自己相関関数が分からず、ログスケールでの自己相関関数となり、指数関数による近似が必要となる。そこで、近似式としてＫＷＷ関数を用いるので、単なる指数関数では表現できない短時間領域を適切に表現でき、計算精度を向上させることができる。したがって、計算コストを低減するとともに、算出精度を向上させることが可能となる。

本発明の輸送係数を算出する装置を模式的に示すブロック図。 本発明の輸送係数を算出する方法を示すフローチャート。 Multiple-tau相関法で算出した緩和弾性率の時間変化を示す図。 Multiple-tau相関法を用いずにｋ個の時系列データから算出した自己相関関数の値を示す図。 指数関数を重ね合わせた近似式に関する説明図。 ＫＷＷ関数を重ね合わせた近似式に関する説明図。 従来の平均化により輸送係数を求める方法に関する説明図。

以下、本発明の一実施形態を、図面を参照して説明する。

［輸送係数を算出する装置］
本実施形態の装置は、予め定めた解析温度及び解析圧力下における解析対象となる物質の輸送係数を、分子動力学シミュレーションを用いて算出する装置である。輸送係数としては、粘性係数、拡散係数及び熱導電率などが挙げられる。本実施形態では、粘性係数を算出する例を用いて説明する。

図１に示すように、装置１は、初期設定部１０と、分子動力学算出部１１と、物理量算出部１２と、自己相関関数算出部１３と、近似部１４と、輸送係数算出部１５と、を有する。これら各部１０〜１５は、ＣＰＵ、メモリ、各種インターフェイス等を備えたパソコン等の情報処理装置において予め記憶されている図示しない処理ルーチンをＣＰＵが実行することによりソフトウェア及びハードウェアが協働して実現される。

図１に示す初期設定部１０は、キーボードやマウス等の既知の操作部を介してユーザからの操作を受け付け、解析対象となる物質を表す原子モデルや解析温度、解析圧力など分子動力学計算に必要な設定等の各種設定を実行し、これら設定値をメモリに記憶する。

図１に示す分子動力学算出部１１は、初期設定部１０により設定された原子モデルデータ、解析温度、解析圧力などの各種パラメータを用い、原子モデルの分子動力学計算に基づき平衡状態における各分子モデルの位置に関する時系列データを算出する。これに伴い、図１に示す物理量算出部１２は、分子動力学算出部１１の算出結果に基づき輸送係数に対応する物理量を時点ｔ_１から時点ｔ_ｋまで表すｋ個の時系列データを算出する。ｋ個は、原子モデルの挙動を再現する時点の個数であり、解析時間の長さでもある。解析時間の長さは初期設定部１０により予め設定される。ｋが大きいほど、計算時間が長くなる。図１に示すように、算出されたｋ個の時系列データはメモリに記憶される。輸送係数が粘性係数の場合はせん断応力テンソルとなり、輸送係数が拡散係数の場合は速度となり、輸送係数が熱伝導係数の場合は熱流束となる。

本字実施形態では、揺動散逸定理に従い、輸送係数は、物理量の時間相関関数の積分値で与えられる。せん断応力Ｐｘｙとすれば、粘性係数ηは下記式（１）で表現され、速度ｖとすれば、拡散係数Ｄは下記式（２）で表現され、熱流束密度Ｊとすれば、熱伝導係数λは下記式（３）で表現される。ここで、ｋ_Ｂはボルツマン定数、Ｖは体積、Ｔは温度である。＜…＞は自己相関関数を示している。

図１に示す自己相関関数算出部１３は、図１に示すように、前記時点ｔ_１から時点ｔ_ｋをｍ個（ｍ＜ｋ）のグループに分け、グループ毎に物理量の平均値を算出し、当該平均値に基づきｍ個の自己相関関数の値Ｇ（ｔ_ｉ）を算出する。ただし、ｉ＝１〜ｍである。本実施形態では、Multiple-tau相関法を用いてｍ個の自己相関関数の値（例えば応力緩和時系列データ）を算出している。勿論、物理量を何らかの手法により平均化する手法であれば、Multiple-tau相関法に限定されない。例えば、最大エントロピー法、ブロックアベレージ法などが挙げられる。なお、図１のメモリに、物理量の時系列データ（ｋ個）から平均値を求め、ｍ個の自己相関関数の値を算出する様子を模式的に示しているが、これは理解を容易にするための説明図であり、厳密ではない。

図１に示す近似部１４は、自己相関関数算出部１３が算出した自己相関関数の値Ｇ（ｔ_ｉ）を、ｅｘｐ｛−（ｔ／τ）^β｝で表されるＫＷＷ関数（Kohlausch-Williams-Watts）を含む近似式で近似し、τ及びβを含むパラメータを決定する。

本実施形態では、次のように決定している。
近似式は、Ｎ個のＫＷＷ関数の重ね合わせ展開式である式（４）で表現される。

まず、近似部１４は、Ｎ＝１として、最小二乗法によってパラメータａ_ｉ、τ_ｉ、β_ｉ（ｉ＝１〜Ｎ）を決定する。具体的には、応力緩和時系列データ（ｍ個の自己相関関数の値）に対する二乗誤差の総和Δ^２が最小となるように、パラメータａ_ｉ、τ_ｉ、β_ｉ（ｉ＝１〜Ｎ）を決定する。前記総和Δ^２は式（５）で表される。

次に、近似部１４は、決定したパラメータａ_ｉ、τ_ｉ、β_ｉ及び近似式（４）で算出される値と、ｍ個の自己相関関数の値との誤差が所定の許容値以下となるまで、Ｎを１ずつ大きくして最小二乗法によるパラメータａ_ｉ、τ_ｉ、β_ｉの決定を繰り返す。
具体的には、展開式の長さＮは、Ｎ＝１から始めて、式（６）が満たされるまで１ずつ増やしている。

式（６）は、一致度Ｒ^２が９９．９％以上であればよいことを示している。誤差としては、０．１％以下である。所定の許容値に対応する値０．９９９は、予め初期設定部１０を介して設定される。なお、本実施形態では、Ｎ＝３で式（６）の条件を満たした。

図１に示す輸送係数算出部１５は、決定されたパラメータを用いて近似式を時間積分し、輸送係数を算出する。粘性係数であれば、式（７）で表される。Γはガンマ関数を示す。パラメータａ_ｉ、τ_ｉ、β_ｉ、Ｎは、近似部１４により決定されている。

本実施形態の効果を示すために、液体スチレンの粘性係数を算出した。

１）図３は、本実施形態の方法（Multiple-tau相関法）で緩和弾性率の時間変化を算出した図である。図４は、Multiple-tau相関法を用いずにｋ個の時系列データから自己相関関数の値を算出した図である。図４には、長時間領域において生じる揺らぎが存在するが、図３では揺らぎが低減されていることが分かる。なお、緩和弾性率は自己相関関数の値にＶ／ｋ_ＢＴを乗じて算出されるので、図４の値にＶ／ｋ_ＢＴを乗じれば、図３と図４は比較可能になる。

２）図５は、本実施形態の方法（Multiple-tau相関法）で算出した緩和弾性率の時間変化（応力緩和）を、指数関数を重ね合わせた近似式で近似した値を示す図である。指数関数は、ｅｘｐ｛−（ｔ／τ）｝である。指数関数の重ね合わせの数を３つとした。図６は、本実施形態の方法（Multiple-tau相関法）で算出した緩和弾性率の時間変化（応力緩和）を、ＫＷＷ関数を重ね合わせた近似式で近似した図である。ＫＷＷ関数の重ね合わせの数を３つとした。図５において、短時間領域では、近似式と応力緩和の値の誤差が大きく、近似式が適切とはいえないことが分かる。図６は、図５に比べて、近似式と応力緩和の値の誤差が小さく、近似式が適切といえる。

３）図７は、粘性係数を、式（１）の時間積分により算出した図である。図中に示すように、粘性係数の値が概ね一定値に収束し、その後、悪化することを確認できるまで計算を続け、計算を終了する。計算時間は、時点５．Ｅ＋０４まで計算した。粘性係数は、値が収束したと考えられる箇所（図中の四角）部分の平均値を採用した。

４）液体スチレンの粘性係数の算出値及び計算時間は、次の通りである。
実験値：６９６μＰａＳ
図７に示す平均値を用いた方法：７０７μＰａＳ；計算時間５．Ｅ＋０４
図５に示す指数関数を用いた近似式：７２９μＰａＳ；計算時間１．Ｅ＋０４
図６に示すＫＷＷ関数を用いた近似式：６９２μＰａＳ；計算時間１．Ｅ＋０４
上記の通り、図７の従来方法は、実験値に対して或る程度の精度で輸送係数を算出できるが、計算コストが高いというデメリットが存在する。
図５に示す方法は、計算コストが少なく好ましいが、短時間領域の物理量を近似式でうまく表現できていないために、輸送係数の算出精度が悪い。
図６に示す方法は、計算コストが少なく、更に、輸送係数の算出精度が他の方法に比して高い。これは、Multiple-tau相関法による平均化の効果で長時間領域の揺らぎが低減し、更に、短時間領域の値を近似式で適切に表現できているためと考えられる。
したがって、本実施形態の方法であれば、計算コストの低減と、算出精度の向上とを両立することができることが理解できる。本実施形態では、計算誤差を５％から１％に低減し、計算時間を約５分の１に低減している。

［輸送係数を算出する方法］
上記装置１を用いた輸送係数を算出する方法について、図２を用いて説明する。

まず、ステップＳＴ１において、図１に示す初期設定部１０は、解析対象となる物質を表す原子モデルや解析温度、解析圧力など分子動力学計算に必要な設定等の各種設定を実行し、これら設定値をメモリに記憶する。

次のステップＳＴ２において、図１に示す分子動力学算出部１１は、予め設定された各種パラメータを用い、原子モデルの分子動力学計算に基づき平衡状態にある各分子モデルの位置に関する時系列データを算出する。これに伴い、図１に示す物理量算出部１２は、ステップＳＴ３において、分子動力学計算の算出結果に基づき輸送係数に対応する物理量（本実施形態ではせん断応力テンソル）を時点ｔ_１から時点ｔ_ｋまで表すｋ個の時系列データを算出する。

次のステップＳＴ４において、図１に示す自己相関関数算出部１３は、時点ｔ_１から時点ｔ_ｋをｍ個（ｍ＜ｋ）のグループに区分し、グループ毎に算出した物理量の平均値に基づきｍ個の自己相関関数の値を算出する。本実施形態では、Multiple-tau相関法を用いているが、これに限定されない。例えば、最大エントロピー法、又はブロックアベレージ法が利用可能である。

ステップＳＴ５〜８において、図１に示す近似部１４は、ｍ個の自己相関関数の値Ｇ（ｔ_ｉ）を、ｅｘｐ｛−（ｔ／τ）^β｝で表されるＫＷＷ関数（Kohlausch-Williams-Watts）を含む近似式で近似し、τ及びβを含むパラメータを決定する。

具体的には、ステップＳＴ５において、近似部１４は、Ｎ＝１とする。次のステップＳＴ６において、近似部１４は、最小二乗法を用いて、ｍ個の自己相関関数の値Ｇ（ｔ_ｉ）を式（４）で近似し、パラメータａ_ｉ、τ_ｉ、β_ｉ（ｉ＝１〜Ｎ）を決定する。次のステップＳＴ７において、近似部１４は、決定したパラメータａ_ｉ、τ_ｉ、β_ｉ及び近似式（４）で算出される値と、ｍ個の自己相関関数の値との誤差が所定の許容値以下であるかを判定する。ステップＳＴ７において、誤差が所定の許容値以下であると判断されるまで（ＳＴ７：ＮＯ）、Ｎを１ずつ大きくし（ＳＴ８）、ステップＳＴ６を繰り返す。勿論、このようにＮを徐々に増大させて検証する方法を採用せずに、複数のＫＷＷ関数を重ね合わた近似式を定義しておいてもよい。

ステップＳＴ７において、誤差が所定の許容値以下であると判断されると（ＳＴ７：ＹＥＳ）、次のステップＳＴ９において、輸送係数算出部１５は、決定されたパラメータａ_ｉ、τ_ｉ、β_ｉを用いて近似式（４）を時間積分して輸送係数を算出する。なお、近似式（４）を時間積分した式は（７）である。

以上のように、本実施形態の輸送係数を算出する方法は、
解析対象となる物質を表す原子モデルデータを用いた分子動力学計算に基づき、予め定めた解析温度及び解析圧力下における平衡状態にある分子の挙動を算出し、輸送係数に対応する物理量を時点ｔ_１から時点ｔ_ｋまで表すｋ個の時系列データを算出するステップＳＴ３と、
時点ｔ_１から時点ｔ_ｋをｍ個（ｍ＜ｋ）のグループに区分し、グループ毎に算出した物理量の平均値に基づきｍ個の自己相関関数の値Ｇ（ｔ_ｉ）［ｉ＝１〜ｍ］を算出するステップＳＴ４と、
ｍ個の自己相関関数の値Ｇ（ｔ_ｉ）を、ｅｘｐ｛−（ｔ／τ）^β｝で表されるＫＷＷ関数（Kohlausch-Williams-Watts）を含む近似式で近似し、τ及びβを含むパラメータを決定するステップ（ＳＴ５〜８）と、
決定されたパラメータを用いて近似式を時間積分して輸送係数を算出するステップＳＴ９と、
を含む。

本実施形態の輸送係数を算出する装置は、
解析対象となる物質を表す原子モデルデータを用いた分子動力学計算に基づき、予め定めた解析温度及び解析圧力下における平衡状態にある分子の挙動を算出し、輸送係数に対応する物理量を時点ｔ_１から時点ｔ_ｋまで表すｋ個の時系列データを算出する物理量算出部１２と、
前記時点ｔ_１から時点ｔ_ｋをｍ個（ｍ＜ｋ）のグループに区分し、グループ毎に算出した物理量の平均値に基づきｍ個の自己相関関数の値Ｇ（ｔ_ｉ）［ｉ＝１〜ｍ］を算出する自己相関関数算出部１３と、
前記ｍ個の自己相関関数の値Ｇ（ｔ_ｉ）を、ｅｘｐ｛−（ｔ／τ）^β｝で表されるＫＷＷ関数（Kohlausch-Williams-Watts）を含む近似式で近似し、τ及びβを含むパラメータを決定する近似部１４と、
決定されたパラメータを用いて近似式を時間積分して輸送係数を算出する輸送係数算出部１５と、
を備える。

このように、時点ｔ_１から時点ｔ_ｋをｍ個（ｍ＜ｋ）のグループに区分し、グループ毎に算出した物理量の平均値に基づき自己相関関数の値Ｇ（ｔ_ｉ）［ｉ＝１〜ｍ］を算出するので、平均化により長時間領域の揺らぎを抑制でき、揺らぎの影響を確認するために必要以上の時間を計算する必要がないので計算コストを低減することが可能となる。その反面、物理量の平均値を取ることで、時系列データの数が減少するので、線形スケールでの自己相関関数が分からず、ログスケールでの自己相関関数となり、指数関数による近似が必要となる。そこで、近似式としてＫＷＷ関数を用いるので、単なる指数関数では表現できない短時間領域を適切に表現でき、計算精度を向上させることができる。したがって、計算コストを低減するとともに、算出精度を向上させることが可能となる。

本実施形態では、自己相関関数の値Ｇ（ｔ_ｉ）は、Multiple-tau相関法、最大エントロピー法、又はブロックアベレージ法を用いて算出する。

本実施形態では、近似式は、式（４）で表され、
Ｎ＝１として、最小二乗法によってパラメータａ_ｉ、τ_ｉ、β_ｉ（ｉ＝１〜Ｎ）を決定し、
決定したパラメータａ_ｉ、τ_ｉ、β_ｉ及び前記近似式で算出される値Ｇ’（ｔ）と、前記ｍ個の自己相関関数の値との誤差が所定の許容値以下となるまで、Ｎを１ずつ大きくして、前記最小二乗法によるパラメータａ_ｉ、τ_ｉ、β_ｉの決定を繰り返す。

この方法によれば、近似式の展開数を増大させて、誤差が所定の許容値以下にしているので、近似の精度が目標の許容値を満たし、その結果、輸送係数の精度を更に向上させることが可能となる。

本実施形態のプログラムは、上記方法を構成する各ステップをコンピュータに実行させるプログラムである。このプログラムを実行することによっても、上記方法の奏する作用効果を得ることが可能となる。言い換えると、上記算出方法を使用しているともいえる。

以上、本発明の実施形態について図面に基づいて説明したが、具体的な構成は、これらの実施形態に限定されるものでないと考えられるべきである。本発明の範囲は、上記した実施形態の説明だけではなく特許請求の範囲によって示され、さらに特許請求の範囲と均等の意味および範囲内でのすべての変更が含まれる。

例えば、図１に示す各部１０〜１５は、所定プログラムをコンピュータのＣＰＵで実行することで実現しているが、各部を専用メモリや専用回路で構成してもよい。

上記の各実施形態で採用している構造を他の任意の実施形態に採用することは可能である。各部の具体的な構成は、上述した実施形態のみに限定されるものではなく、本発明の趣旨を逸脱しない範囲で種々変形が可能である。

１１…分子動力学算出部
１２…物理量算出部
１３…自己相関関数算出部
１４…近似部
１５…輸送係数算出部

Claims (7)

1. コンピュータが実行する方法であって、
   解析対象となる物質を表す原子モデルデータを用いた分子動力学計算に基づき、予め定めた解析温度及び解析圧力下における平衡状態にある分子の挙動を算出し、輸送係数に対応する物理量を時点ｔ_１から時点ｔ_ｋまで表すｋ個の時系列データを算出するステップと、
   前記時点ｔ_１から時点ｔ_ｋをｍ個（ｍ＜ｋ）のグループに区分し、グループ毎に算出した物理量の平均値に基づきｍ個の自己相関関数の値を算出するステップと、
   前記ｍ個の自己相関関数の値を、ｅｘｐ｛−（ｔ／τ）^β｝で表されるＫＷＷ関数（Kohlausch-Williams-Watts）を含む近似式で近似し、τ及びβを含むパラメータを決定するステップと、
   決定された前記パラメータを用いて前記近似式を時間積分して輸送係数を算出するステップと、
   を含む、輸送係数を算出する方法。
2. 前記自己相関関数の値は、Multiple-tau相関法、最大エントロピー法、又はブロックアベレージ法を用いて算出する、請求項１に記載の方法。
3. 前記近似式は、式（４）で表され、
   
   Ｎ＝１として、最小二乗法によってパラメータａ_ｉ、τ_ｉ、β_ｉ（ｉ＝１〜Ｎ）を決定し、
   決定したパラメータａ_ｉ、τ_ｉ、β_ｉ及び前記近似式で算出される値Ｇ’（ｔ）と、前記ｍ個の自己相関関数の値との誤差が所定の許容値以下となるまで、Ｎを１ずつ大きくして、前記最小二乗法によるパラメータａ_ｉ、τ_ｉ、β_ｉの決定を繰り返す、請求項１又は２に記載の方法。
4. 解析対象となる物質を表す原子モデルデータを用いた分子動力学計算に基づき、予め定めた解析温度及び解析圧力下における平衡状態にある分子の挙動を算出し、輸送係数に対応する物理量を時点ｔ_１から時点ｔ_ｋまで表すｋ個の時系列データを算出する物理量算出部と、
   前記時点ｔ_１から時点ｔ_ｋをｍ個（ｍ＜ｋ）のグループに区分し、グループ毎に算出した物理量の平均値に基づきｍ個の自己相関関数の値を算出する自己相関関数算出部と、
   前記ｍ個の自己相関関数の値を、ｅｘｐ｛−（ｔ／τ）^β｝で表されるＫＷＷ関数（Kohlausch-Williams-Watts）を含む近似式で近似し、τ及びβを含むパラメータを決定する近似部と、
   決定された前記パラメータを用いて前記近似式を時間積分して輸送係数を算出する輸送係数算出部と、
   を備える、輸送係数を算出する装置。
5. 前記自己相関関数の値は、Multiple-tau相関法、最大エントロピー法、又はブロックアベレージ法を用いて算出する、請求項４に記載の装置。
6. 前記近似式は、式（４）で表され、
   
   前記近似部は、Ｎ＝１として、最小二乗法によってパラメータａ_ｉ、τ_ｉ、β_ｉ（ｉ＝１〜Ｎ）を決定し、
   決定したパラメータａ_ｉ、τ_ｉ、β_ｉ及び前記近似式で算出される値Ｇ’（ｔ）と、前記ｍ個の自己相関関数の値との誤差が所定の許容値以下となるまで、Ｎを１ずつ大きくして、前記最小二乗法によるパラメータａ_ｉ、τ_ｉ、β_ｉの決定を繰り返す、請求項４又は５に記載の装置。
7. 請求項１〜３のいずれかに記載の方法をコンピュータに実行させるプログラム。