JP2014106554A - 輸送係数の算出方法、輸送係数の算出装置、及び輸送係数の算出プログラム - Google Patents

Abstract

【課題】量子論を考慮した適切な輸送係数を算出する輸送係数の算出方法を提供する。
【解決手段】解析対象となる物質を表す原子モデルを用いた分子動力学計算に基づき、予め定めた解析温度及び解析圧力下における平衡状態にある分子の輸送係数に対応する物理量の時系列データを算出するステップＳＴ２〜３と、解析温度において確率１／２で励起される振動の周波数の逆数ｔ_maxを次の式（４）により算出するステップＳＴ４と、上記逆数ｔ_maxをパラメータとする窓関数と物理量の時系列データとの畳み込み演算を行うステップＳＴ５と、畳み込み演算後の物理量の時系列データから自己相関関数を算出するステップＳＴ６と、自己相関関数に基づき輸送係数を算出するステップＳＴ７〜８と、を有する。
ｔ_max＝ｈ／（ｋ_Ｂ Ｔ ln３） …（４）
ただし、ｈはプランク定数、ｋ_Ｂはボルツマン定数、Ｔは解析温度。
【選択図】図１

Description

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

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

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

しかしながら、分子動力学では物質を構成する分子の運動に対する量子効果を無視した近似を用いるため、現実には存在しない寄与が物理量に混入してしまい、算出される輸送係数に大きな誤差が生じてしまう。具体的には、或る温度Ｔ[K]において確率１／２で励起される振動運動の周波数ｖは、ｖ＝（ｋ_ＢＴ）／ｈ（ln ３）で表される。ｈはプランク定数、ｋ_Ｂはボルツマン定数である。この周波数ｖよりも大きな周波数で振動する運動は、量子論に従う限りほとんど励起されない。一方、分子動力学において粒子の運動は量子論ではなく古典力学に基づき計算されるため、変位（エネルギー）が与えられれば、フェムト秒（１．０×１０^−１５秒）程度の微小時間スケールであっても、現実には存在しない仮想的な振動運動が含まれることになる。物理量は分子の運動と共に時間変化するので、分子動力学計算の結果を用いると、本来存在しない高振動数領域にパワースペクトルが生じる。輸送係数は流束量の自己相関関数の全時間領域にわたる積分値であることから、上記非現実的な分子振動に由来するパワーが加算されると、その値は過大評価され、求められる輸送係数の誤差も大きくなることが判明した。

特許文献１には、分子動力学から得られた物理量の自己相関関数から輸送係数を計算するにあたり、分子振動に起因する振動成分を低周波透過フィルタによって目分量で除去し、カットオフ周波数ゼロの極限値を求めるべき輸送係数の値としているようである。

しかしながら、特許文献１では、上記量子論を考慮すべきことの記載がなく、特許文献１の図５及び図６に示すように、カットオフ周波数の下限を目分量で決定し、プロット結果を近似した直線を用いてカットオフ周波数ゼロの極限値を外挿により求めているので、直線で近似できない場合は算出できず、算出できたとしても輸送係数の過大評価につながってしまうおそれがある。

本発明は、このような課題に着目してなされたもので、その目的は、量子論を考慮した適切な輸送係数を算出する輸送係数の算出方法、算出装置及び算出プログラムを提供することである。

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

すなわち、本発明の輸送係数の算出方法は、解析対象となる物質の輸送係数をコンピュータが算出するための輸送係数の算出方法であって、
解析対象となる物質を表す原子モデルを用いた分子動力学計算に基づき、予め定めた解析温度及び解析圧力下における平衡状態にある分子の位置に関する時系列データを算出するステップと、
前記分子動力学計算の結果に基づき輸送係数に対応する物理量の時系列データを算出するステップと、
前記解析温度において確率１／２で励起される振動の周波数の逆数ｔ_maxを次の式（４）により算出するステップと、
ｔ_max＝ｈ／（ｋ_Ｂ Ｔ ln３） …（４）
ただし、ｈはプランク定数、ｋ_Ｂはボルツマン定数、Ｔは解析温度である。
前記逆数ｔ_maxをカットオフパラメータとする窓関数と前記物理量の時系列データとの畳み込み演算を行うステップと、
前記畳み込み演算後の物理量の時系列データから自己相関関数を算出するステップと、
前記自己相関関数に基づき輸送係数を算出するステップと、
を有することを特徴とする。

このように、解析温度において確率１／２で励起される振動の周波数の逆数をボルツマン分布に基づく式（４）を用いて算出し、当該値をカットオフパラメータとする窓関数を用いた畳み込み演算（積分）を行うので、量子論に基づき除去すべき振動成分を決定し、当該振動成分を除去することができ、量子論を考慮してより適切な輸送係数を算出することが可能となる。

輸送係数を精度良く算出するためには、自己相関関数をラプラス変換でフィッティングし、フィティング関数を時間積分することで輸送係数を算出することが好ましい。

本発明は、装置としても特定可能である。
すなわち、本発明の送係数の算出装置は、解析対象となる物質を表す原子モデルを用いた分子動力学計算に基づき、予め定めた解析温度及び解析圧力下における平衡状態にある分子の位置に関する時系列データを算出する分子動力学算出部と、前記分子動力学算出部の算出結果に基づき輸送係数に対応する物理量の時系列データを算出する物理量算出部と、
前記解析温度において確率１／２で励起される振動の周波数の逆数ｔ_maxを次の式（４）により算出するカットオフパラメータ決定部と、
ｔ_max＝ｈ／（ｋ_Ｂ Ｔ ln３） …（４）
ただし、ｈはプランク定数、ｋ_Ｂはボルツマン定数、Ｔは解析温度である。
前記カットオフパラメータ決定部が算出した逆数ｔ_maxをカットオフパラメータとする窓関数と前記物理量算出部が算出した物理量の時系列データとの畳み込み演算を行う畳み込み演算部と、前記畳み込み演算後の物理量の時系列データから自己相関関数を算出する自己相関関数算出部と、前記自己相関関数に基づき輸送係数を算出する輸送係数算出部と、
を備えることを特徴とする。
この装置を使用することによっても、上記算出方法が奏する作用効果を得ることができる。

輸送係数を精度良く算出するためには、前記輸送係数算出部は、前記自己相関関数をラプラス変換でフィッティングし、フィティング関数を時間積分することで輸送係数を算出するように構成されていることが好ましい。

本発明は、上記算出方法を構成するステップを、プログラムの観点から特定することも可能である。
すなわち、本発明の算出プログラムは、解析対象となる物質を表す原子モデルを用いた分子動力学計算に基づき、予め定めた解析温度及び解析圧力下における平衡状態にある分子の位置に関する時系列データを算出するステップと、前記分子動力学計算の結果に基づき輸送係数に対応する物理量の時系列データを算出するステップと、前記解析温度において確率１／２で励起される振動の周波数の逆数ｔ_maxを次の式（４）により算出するステップと、
ｔ_max＝ｈ／（ｋ_Ｂ Ｔ ln３） …（４）
ただし、ｈはプランク定数、ｋ_Ｂはボルツマン定数、Ｔは解析温度。
前記逆数ｔ_maxをカットオフパラメータとする窓関数と前記物理量の時系列データとの畳み込み演算を行うステップと、前記畳み込み演算後の物理量の時系列データから自己相関関数を算出するステップと、前記自己相関関数に基づき輸送係数を算出するステップと、
をコンピュータに実行させることを特徴とする。
このプログラムを実行することによっても、上記算出方法が奏する作用効果を得ることができる。

本発明に係る輸送係数の算出方法を示すフローチャート。 本発明に係る輸送係数の算出装置を模式的に示すブロック図。 輸送係数の算出結果を示す説明図。

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

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

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

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

図２に示す分子動力学算出部１１は、初期設定部１０により設定された原子モデル、解析温度、解析圧力などの各種パラメータを用い、原子モデルの分子動力学計算に基づき平衡状態における各分子モデルの位置に関する時系列データを算出する。これに伴い、図２に示す物理量算出部１２は、分子動力学算出部１１の算出結果に基づき輸送係数に対応する物理量の時系列データを算出する。輸送係数が粘性係数の場合はせん断応力テンソルとなり、輸送係数が拡散係数の場合は速度となり、輸送係数が熱伝導係数の場合は熱流束となる。

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

図２に示すカットオフパラメータ決定部１３は、初期設定部１０により設定された解析温度において確率１／２で励起される振動の周波数の逆数ｔ_maxを量子論に基づき算出する。具体的には、ボルツマン分布に基づき下記式（４）を用いて算出する。
ｔ_max＝ｈ／（ｋ_Ｂ Ｔ ln３） …（４）
ただし、ｈはプランク定数、ｋ_Ｂはボルツマン定数、Ｔは解析温度である。
この値ｔ_maxを算出するのは、解析温度において確率１／２で励起される振動の周波数よりも大きな周波数の振動は、量子論に従う限りほとんど励起されないため、かかる周波数以上の振動成分を分子動力学に基づく算出結果から除去するためである。

図２に示す畳み込み演算部１４は、カットオフパラメータ決定部１３で算出した値ｔ_maxをカットオフパラメータとする窓関数と、物理量算出部１２で算出した物理量（せん断応力テンソル）の時系列データとを畳み込み演算（積分）し、物理量の時系列データからボルツマン分布に基づき振動成分を除去する。具体的には、下記式（５）のような畳み込み変換により行う。Ｆは物理量、Ｗが窓関数である。本実施形態で用いる窓関数は矩形窓である。

図２に示す自己相関関数算出部１５は、畳み込み演算部１４による畳み込み演算後の物理量の時系列データから自己相関関数を算出する。

図２に示す輸送係数算出部１６は、自己相関関数算出部１５が取得した自己相関関数をラプラス変換でフィッティングし、フィッティング関数を時間積分することで輸送係数（本実施形態では粘性係数）を算出する。

［輸送係数の算出方法］
上記算出装置１を用いた粗視化モデルの非結合ポテンシャルの算出方法について説明する。

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

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

次のステップＳＴ４（図１参照）において、図２に示すカットオフパラメータ決定部１３は、カットオフパラメータであるｔ_maxを、上記式（４）を用いて算出する。

次のステップＳＴ５（図１参照）において、図２に示す畳み込み演算部１４は、ｔ_maxをカットオフパラメータとする窓関数と物理量（せん断応力テンソル）の時系列データとを畳み込み演算（積分）する。

次のステップＳＴ６（図１参照）において、図２に示す自己相関関数算出部１５は、演算後の物理量の時系列データから自己相関関数を算出する。次のステップＳＴ７（図１参照）において、図２に示す輸送係数算出部１６は、自己相関関数をラプラス変換する。次のステップＳＴ８において、図２に示す輸送係数算出部１６は、ラプラス変換でフィッティングした関数を時間積分して輸送係数（本実施形態では粘性係数）を算出する。

上記装置及び方法による算出結果の一例を図３に示す。図３の例では、温度３００Ｋ及び１気圧にて液体イソプレンの粘性係数をカットオフ周波数を変えて算出した例である。式（４）を用いて算出すれば、図３にて矢印で示すように、ｔ_maxは250fs[フェムト秒；１０^−１５秒]で、カットオフ周波数は、４THz[テラヘルツ；１０^１２ヘルツ]である。このときの計算値は、2453μＰ[マイクロポアズ]となり、実験値2620μＰに近い値となる。一方、図３を見れば明らかであるが、カットオフ周波数を４THzよりも大きく設定すると、量子論では励起されない非現実的な振動成分が加算され、過大評価されて、輸送係数が大きくなり、誤差が大きくなる。一方、カットオフ周波数を４THzよりも小さく設定すると、上記非現実的な振動成分以外の要因により輸送係数に誤差が生じると考えられる。これは、窓関数の畳み込み演算により、長時間平均がゼロにならない余計な成分が残り、これが長時間積分によって支配的となるのが要因と考えられる。カットオフ周波数が４THzであれば、適切な区間で積分でき、長時間平均がゼロにならない余計な成分が残らず、適切に演算できると考えられる。

以上のように、本実施形態の輸送係数の算出方法は、解析対象の物質の輸送係数をコンピュータが算出するための方法であって、解析対象となる物質を表す原子モデルを用いた分子動力学計算に基づき、予め定めた解析温度及び解析圧力下における平衡状態にある分子の位置に関する時系列データを算出するステップＳＴ２と、前記分子動力学計算の結果に基づき輸送係数に対応する物理量の時系列データを算出するステップＳＴ３と、前記解析温度において確率１／２で励起される振動の周波数の逆数ｔ_maxを次の式（４）により算出するステップＳＴ４と、前記逆数ｔ_maxをカットオフパラメータとする窓関数と前記物理量の時系列データとの畳み込み演算を行うステップＳＴ５と、前記畳み込み演算後の物理量の時系列データから自己相関関数を算出するステップＳＴ６と、前記自己相関関数に基づき輸送係数を算出するステップＳＴ７、ＳＴ８と、を有する。
ｔ_max＝ｈ／（ｋ_Ｂ Ｔ ln３）…（４）
ただし、ｈはプランク定数、ｋ_Ｂはボルツマン定数、Ｔは解析温度である。

本実施形態の輸送係数の算出装置は、解析対象となる物質を表す原子モデルを用いた分子動力学計算に基づき、予め定めた解析温度及び解析圧力下における平衡状態にある分子の位置に関する時系列データを算出する分子動力学算出部１１と、分子動力学算出部１１の算出結果に基づき輸送係数に対応する物理量の時系列データを算出する物理量算出部１２と、解析温度において確率１／２で励起される振動の周波数の逆数ｔ_maxを式（４）により算出するカットオフパラメータ決定部１３と、カットオフパラメータ決定部１３が算出した逆数ｔ_maxをカットオフパラメータとする窓関数と物理量算出部１２が算出した物理量の時系列データとの畳み込み演算を行う畳み込み演算部１４と、前記畳み込み演算後の物理量の時系列データから自己相関関数を算出する自己相関関数算出部１５と、前記自己相関関数に基づき輸送係数を算出する輸送係数算出部１６と、を備える

このように、解析温度において確率１／２で励起される振動の周波数の逆数をボルツマン分布に基づく式（４）を用いて算出し、当該値をカットオフパラメータとする窓関数を用いた畳み込み演算（積分）を行うので、量子論に基づき除去すべき振動成分を決定し、当該振動成分を除去することができ、量子論を考慮してより適切な輸送係数を算出することが可能となる。

本実施形態では、輸送係数算出部１６は、自己相関関数をラプラス変換でフィッティングし、フィティング関数を時間積分することで輸送係数を算出するように構成されている。
解析対象の物質が流体の場合の自己相関関数は、減衰が著しいという特性を有するため、時間軸を均等に区画する通常の積分を行えば、自己相関関数の特性により誤差が生じてしまう問題がある。そこで、本実施形態では、自己相関関数をラプラス変換してフィッティング関数とし、指数関数の重ね合わせで表現される当該フィティング関数を時間積分して輸送係数を算出するようにすれば、ｌｏｇスケールで表現し、且つ個々の指数関数が厳密に時間積分可能であるので、その展開係数の和で計算することにより、精度良く自己相関関数の時間積分を算出でき、輸送係数を精度良く算出することが可能となる。

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

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

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

また、畳み込み演算で用いる窓関数は、上記のような矩形窓でなくてもよく、種々の窓関数を利用可能である。例えば、ガウス関数やハニング関数などを利用可能である。
上記の各実施形態で採用している構造を他の任意の実施形態に採用することは可能である。各部の具体的な構成は、上述した実施形態のみに限定されるものではなく、本発明の趣旨を逸脱しない範囲で種々変形が可能である。

１１…分子動力学算出部
１２…物理量算出部
１３…カットオフパラメータ決定部
１４…畳み込み演算部
１５…自己相関関数算出部
１６…輸送係数算出部

Claims (5)

1. 解析対象となる物質の輸送係数をコンピュータが算出するための輸送係数の算出方法であって、
   解析対象となる物質を表す原子モデルを用いた分子動力学計算に基づき、予め定めた解析温度及び解析圧力下における平衡状態にある分子の位置に関する時系列データを算出するステップと、
   前記分子動力学計算の結果に基づき輸送係数に対応する物理量の時系列データを算出するステップと、
   前記解析温度において確率１／２で励起される振動の周波数の逆数ｔ_maxを次の式（４）により算出するステップと、
   ｔ_max＝ｈ／（ｋ_Ｂ Ｔ ln３） …（４）
   ただし、ｈはプランク定数、ｋ_Ｂはボルツマン定数、Ｔは解析温度。
   前記逆数ｔ_maxをカットオフパラメータとする窓関数と前記物理量の時系列データとの畳み込み演算を行うステップと、
   前記畳み込み演算後の物理量の時系列データから自己相関関数を算出するステップと、
   前記自己相関関数に基づき輸送係数を算出するステップと、
   を有することを特徴とする輸送係数の算出方法。
2. 前記自己相関関数をラプラス変換でフィッティングし、フィティング関数を時間積分することで輸送係数を算出する請求項１に記載の輸送係数の算出方法。
3. 解析対象となる物質を表す原子モデルを用いた分子動力学計算に基づき、予め定めた解析温度及び解析圧力下における平衡状態にある分子の位置に関する時系列データを算出する分子動力学算出部と、
   前記分子動力学算出部の算出結果に基づき輸送係数に対応する物理量の時系列データを算出する物理量算出部と、
   前記解析温度において確率１／２で励起される振動の周波数の逆数ｔ_maxを次の式（４）により算出するカットオフパラメータ決定部と、
   ｔ_max＝ｈ／（ｋ_Ｂ Ｔ ln３） …（４）
   ただし、ｈはプランク定数、ｋ_Ｂはボルツマン定数、Ｔは解析温度。
   前記カットオフパラメータ決定部が算出した逆数ｔ_maxをカットオフパラメータとする窓関数と前記物理量算出部が算出した物理量の時系列データとの畳み込み演算を行う畳み込み演算部と、
   前記畳み込み演算後の物理量の時系列データから自己相関関数を算出する自己相関関数算出部と、
   前記自己相関関数に基づき輸送係数を算出する輸送係数算出部と、
   を備えることを特徴とする輸送係数の算出装置。
4. 前記輸送係数算出部は、前記自己相関関数をラプラス変換でフィッティングし、フィティング関数を時間積分することで輸送係数を算出するように構成されている請求項３に記載の輸送係数の算出装置。
5. 解析対象となる物質を表す原子モデルを用いた分子動力学計算に基づき、予め定めた解析温度及び解析圧力下における平衡状態にある分子の位置に関する時系列データを算出するステップと、
   前記分子動力学計算の結果に基づき輸送係数に対応する物理量の時系列データを算出するステップと、
   前記解析温度において確率１／２で励起される振動の周波数の逆数ｔ_maxを次の式（４）により算出するステップと、
   ｔ_max＝ｈ／（ｋ_Ｂ Ｔ ln３） …（４）
   ただし、ｈはプランク定数、ｋ_Ｂはボルツマン定数、Ｔは解析温度。
   前記逆数ｔ_maxをカットオフパラメータとする窓関数と前記物理量の時系列データとの畳み込み演算を行うステップと、
   前記畳み込み演算後の物理量の時系列データから自己相関関数を算出するステップと、
   前記自己相関関数に基づき輸送係数を算出するステップと、
   をコンピュータに実行させることを特徴とする輸送係数の算出プログラム。