(57)【要約】【課題】臨界点近傍の温度・圧力領域においても、精度良い物性値評価を実現する。【解決手段】複数の粒子からなる計算体系を模擬した分子動力学計算を用いて、評価温度および評価圧力における物性値を評価する際に、評価温度および評価圧力に基づいて解析条件温度および解析条件圧力を求める解析条件補正工程（工程Ｓ４）と、計算体系の温度を解析条件温度とし、計算体系の圧力を解析条件圧力とする分子動力学計算を行う分子動力学計算工程（工程Ｓ５）と、分子動力学計算で得られた結果に基づいて評価温度および評価圧力における評価対象物性値を算出する物性値評価工程（工程Ｓ９）と、を実行する。たとえば、解析条件温度は、評価温度に実際の臨界温度を乗じたものを分子動力学計算で求めた臨界温度で除して求め、解析条件圧力は、評価圧力に実際の臨界圧力を乗じたものを分子動力学計算で求めた臨界圧力で除して求める。【選択図】図１
2007248273.tif 000002
【特許請求の範囲】
【請求項１】
複数の粒子からなる計算体系を模擬した分子動力学計算を用いて、評価温度および評価圧力における物性値を評価する物性値評価方法において、
前記評価温度および前記評価圧力に基づいて解析条件温度および解析条件圧力を求める解析条件補正工程と、
前記計算体系の温度を前記解析条件温度とし、前記計算体系の圧力を前記解析条件圧力とする分子動力学計算を行う分子動力学計算工程と、
前記分子動力学計算で得られた結果に基づいて前記評価温度および前記評価圧力における評価対象物性値を算出する物性値評価工程と、
を有することを特徴とする物性値評価方法。
【請求項２】
前記粒子の実際の臨界温度および臨界圧力である臨界温度真値および臨界圧力真値を入力する工程と、
分子動力学計算によって前記粒子の臨界温度および臨界圧力である臨界温度計算値および臨界圧力計算値を計算する工程と、
を有し、
前記解析条件補正工程は、前記評価温度、前記臨界温度真値並びに前記臨界温度計算値、および、前記評価圧力、前記臨界圧力真値並びに前記臨界圧力計算値に基づいて前記解析条件温度および前記解析条件圧力を求めるものである、
ことを特徴とする請求項１記載の物性値評価方法。
【請求項３】
前記解析条件補正工程は、前記評価温度に前記臨界温度真値を乗じたものを前記臨界温度計算値で除して前記解析条件温度を求め、前記評価圧力に前記臨界圧力真値を乗じたものを前記臨界圧力計算値で除して前記解析条件圧力を求めるものである、
ことを特徴とする請求項２記載の物性値評価方法。
【請求項４】
前記評価対象物性値は、拡散係数を含むことを特徴とする請求項１ないし請求項３いずれか記載の物性値評価方法。
【請求項５】
前記評価対象物性値は、粘度を含むことを特徴とする請求項１ないし請求項４いずれか記載の物性値評価方法。
【請求項６】
前記評価対象物性値は、熱伝導度を含むことを特徴とする請求項１ないし請求項５いずれか記載の物性値評価方法。
【請求項７】
複数の粒子からなる計算体系を模擬した分子動力学計算を用いて、評価温度および評価圧力における物性値を評価する物性値評価装置において、
前記評価温度および前記評価圧力に基づいて解析条件温度および解析条件圧力を求める解析条件補正手段と、
前記計算体系の温度を前記解析条件温度とし、前記計算体系の圧力を前記解析条件圧力とする分子動力学計算を行う分子動力学計算手段と、
前記分子動力学計算で得られた結果に基づいて前記評価温度および前記評価圧力における評価対象物性値を算出する物性値評価手段と、
を有することを特徴とする物性値評価装置。
【請求項８】
複数の粒子からなる計算体系を模擬した分子動力学計算を用いて、評価温度および評価圧力における物性値を評価させる物性値評価プログラムにおいて、コンピュータに、
前記評価温度および前記評価圧力に基づいて解析条件温度および解析条件圧力を求めさせる解析条件補正手順と、
前記計算体系の温度を前記解析条件温度とし、前記計算体系の圧力を前記解析条件圧力とする分子動力学計算を行わせる分子動力学計算手順と、
前記分子動力学計算で得られた結果に基づいて前記評価温度および前記評価圧力における評価対象物性値を算出させる物性値評価手順と、
を実行させることを特徴とする物性値評価プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、分子動力学法を用いて物性値を評価する物性値評価方法、物性値評価装置および物性値評価プログラムに関する。
【背景技術】
【０００２】
計算機の性能向上に伴い、系の構成要素である原子や分子レベル、あるいは超微粒子レベルの観点に立ち、微視的および巨視的な物理量等の評価を行なう、分子シミュレーションが盛んに行なわれるようになってきた。連続体としての支配方程式を差分法や有限要素法などで離散化して、その方程式の解を求める数値解析法とは異なり、分子シミュレーションでは、系の構成要素である分子の運動を直接追跡する。よって、実験室で作り出すことが困難な極限状態における諸問題を解決できる可能性も期待されている。
【０００３】
分子シミュレーションには大別してモンテカルロ法と分子動力学法がある。モンテカルロ法は、分子の配置をある確率法則の下に乱数を用いて作成していく確率論的な方法であり、分子動力学法は、分子の運動方程式を連立して解く決定論的方法である。モンテカルロ法、分子動力学法ともに統計力学が基本となっており、それぞれの粒子の位置や速度、エネルギーを統計的に処理することによって、巨視的な物性値を算出することができる。
【０００４】
分子動力学法では、原子や分子を球形や楕円体状の粒子として表現する。それぞれの粒子間にはポテンシャル関数に基づく相互作用が設定されており、それぞれの粒子はそのポテンシャルに従った力を受けて計算体系の中を運動する。
【０００５】
通常の分子動力学法では、粒子間のポテンシャルとして
（ａ）ＬＪポテンシャル
（ｂ）クーロンポテンシャル
（ｃ）結合ポテンシャル
（ｄ）結合角ポテンシャル
などが考慮される。
【０００６】
（ａ）ＬＪポテンシャルφＬＪは、
φＬＪ（ｒ）＝４ε｛（σ／ｒ）ｍ−（σ／ｒ）ｎ｝                      （１）
の形で知られるポテンシャルである。ｒは対象とする２粒子間の距離である。εとσは実数であり、対象とする粒子種毎に定められる。ｍとｎは正の整数で、ｍ＝１２、ｎ＝６とされることが多い。
【０００７】
（ｂ）クーロンポテンシャルφＱは、粒子の持つ電荷の間に働くポテンシャルで、
φＱ（ｒ）＝ｑ１ｑ２／ｒ                                              （２）
の形で与えられる。ｒは対象とする２粒子間の距離であり、ｑ１とｑ２は対象とする２粒子のそれぞれが持つ電荷である。
【０００８】
（ｃ）結合ポテンシャルφＢは、分子内の化学結合を模擬するもので、対象とする２つの原子あるいは粒子間の距離をｒ、安定な原子あるいは粒子間の距離をｒ０としたときに、たとえば、
φＢ（ｒ）＝ｋＢ（ｒ−ｒ０）２                                        （３）
の形で与えられる。ここでｋＢは定数である。
【０００９】
（ｄ）結合角ポテンシャルφＡは、２つの原子あるいは粒子のなす角度に対するポテンシャルで、安定な角度をθ０、現在の角度をθとしたときに、たとえば、
φＡ（θ）＝ｋＡ（θ−θ０）２                                        （４）
の形で与えられる。ｋＡは定数である。
【００１０】
分子動力学計算では設定されたタイムステップ毎に、すべての粒子対に対して上記ポテンシャル値を計算し、そこから算出される各粒子に働く力を用いて運動方程式を数値的に解いて、すべての粒子の位置、速度、加速度などの値を時系列データとして出力する。しかしながら、分子等の粒子の挙動を、すべての温度・圧力領域に対して正確に再現できるポテンシャルモデルは見つかっておらず、沸点など興味のある対象を再現するようなポテンシャルモデルを使用するのが通常である。
【００１１】
物性値は、出力された各粒子の位置、速度、加速度などについての時系列データに対して統計処理を施して算出する。一般的な、拡散係数、粘度、熱伝導度の計算方法を以下に示す（たとえば特許文献１および非特許文献１参照）。
【００１２】
粒子の自己拡散係数Ｄは、次に示す式（５）あるいは式（６）によって計算される。
【数１】
2007248273.tif 000003
【数２】
2007248273.tif 000004
【００１３】
ここで、ｒ（ｔ）とｖ（ｔ）はそれぞれ時刻ｔにおける粒子の位置と速度、ｔ＝０は時系列データの最初の時刻である。粒子が複数個の場合には、各粒子の自己拡散係数を平均することにより、計算精度は向上する。
【００１４】
粘度ηは、次に示す式（７）あるいは式（８）によって計算される。
【数３】
2007248273.tif 000005
【数４】
2007248273.tif 000006
【００１５】
ここで、ｋＢはボルツマン定数、Ｔは温度、Ｖは分子動力学計算で対象としている系の体積、ｚｊ（ｔ）は時刻ｔにおける粒子のｚ方向の位置であり、Ｐｊｘは運動量であって、質量ｍを用いて、
【数５】
2007248273.tif 000007
で表されるものである。
【００１６】
また、Ｊｚｘ（ｔ）は、次式で表される。
【数６】
2007248273.tif 000008
【００１７】
ここで、Ｎは粒子数、Ｚｉｊ＝Ｚｉ−Ｚｊ、ｆｉｊｘ（ｔ）は粒子ｊから粒子ｉに作用する力のｘ成分である。
【００１８】
分子動力学法を用いた熱伝導度評価では、まず、ある物質あるいは化合物を想定し、その物質の分子モデルを表現するパラメータを入力データとして設定する。それら入力データに対して分子動力学計算を、平衡状態のある一定時間について行なう。分子動力学計算から得られた各粒子の位置と速度から、ｊ番目の粒子の全エネルギー（＝運動エネルギーとポテンシャルエネルギーの和）Ｅｊを次式によって算出する。
【数７】
2007248273.tif 000009
【００１９】
ここで、ｍｊ、ｖｊはそれぞれｊ番目の粒子の質量と速度であり、Ｎは全粒子数である。またφは、式（１）ないし式（４）で与えられる各ポテンシャルの和である。
【００２０】
これらの諸量を用いて熱伝導度λを算出する評価式としては、たとえば次のようなものがある。
【数８】
2007248273.tif 000010
【００２１】
ここで、ｔは分子動力学でシミュレーションした時間、ｋＢはボルツマン定数、Ｔは温度、Ｎは計算に利用した原子あるいは分子の個数、Ｅ＊ｊ、Ｅ＊ｋはそれぞれｊ番目、ｋ番目の原子あるいは分子の持つエネルギーの時間平均値からのずれ、ｚｊ（ｔ）はｊ番目の粒子の時刻ｔにおける位置、ｚｋ（０）はｋ番目の粒子の初期位置である。
【００２２】
熱伝導度λを算出する評価式の他の例としては、次式がある。
【数９】
2007248273.tif 000011
【００２３】
ここで、ｔは分子動力学計算を実施した時間、ｋＢはボルツマン定数、Ｔは温度、Ｎは計算に利用した原子あるいは分子の個数、Ｅ＊ｊ（ｔ）、Ｅ＊ｊ（０）はそれぞれｊ番目の原子あるいは分子の時刻ｔおよび初期時間(時刻０)におけるエネルギーの時間平均値からのずれ、ｚｊ（ｔ）、ｚｊ（０）はｊ番目の粒子の時刻ｔおよび初期時刻における位置である。
【００２４】
また、次式によって熱伝導度λを求めることもできる。
【数１０】
2007248273.tif 000012
【００２５】
ここで、ｔは分子動力学計算を実施した時間、Ｖは分子動力学計算で用いた計算体系の体積、ｋＢはボルツマン定数、Ｔは温度であり、＜＞は時間平均を表す。Ａ（ｔ）は、次式で定義される。
【数１１】
2007248273.tif 000013
【数１２】
2007248273.tif 000014
【００２６】
ここでｖｉ（ｔ）はｉ番目の原子あるいは分子の時刻ｔにおける速度、ｘｉｊ（ｔ）は時刻ｔにおけるｉ番目とｊ番目の原子あるいは分子の相対位置、ｆｉｊ（ｔ）は時刻ｔにおけるｉ番目とｊ番目の原子あるいは分子間に働く力である。Ａ２（ｔ）についての時間平均＜Ａ２（ｔ）＞を評価する際には、ｔ０についての平均を計算する。
【００２７】
図４は、従来の分子動力学計算に基づいて評価した水の自己拡散係数を、実験値のフィッティング式によって求めたものと併せて示したグラフである。
【００２８】
図４において、四角（□：７００Ｋ）および三角（△：８００Ｋ）で示した自己拡散係数は、従来の分子動力学計算に基づいて評価した自己拡散係数である。また、実線および一点鎖線で示した自己拡散係数は、非特許文献１に示された次式で表される自己拡散係数の実験値のフィッティング式によって求めたものである。
【数１３】
2007248273.tif 000015
【００２９】
ここで、Ｔは温度、ρは想定する温度および圧力に対応する水の密度である。
【００３０】
図４に示されるように、従来の分子動力学計算に基づいて評価した自己拡散係数は、臨界点（６４７Ｋ、２２ＭＰａ）の近傍で、実験値のフィッティング式で求めた自己拡散係数と大きく乖離していることがわかる。
【特許文献１】特開２００５−２３３７５２号公報
【非特許文献１】神山新一、外１名、「分子動力学シミュレーション」、朝倉書店、１９９７年
【非特許文献２】W. J. Lamb、外２名、“Self-diffusion in compressed supercritical water”、J. Chem. Phys.、Vol. 74、１９８１年、ｐ.６８７５−６８８０
【発明の開示】
【発明が解決しようとする課題】
【００３１】
分子動力学計算では、すべての温度・圧力領域に対して正確に再現できるポテンシャルモデルが現時点では見つかっていないため、特定の事象を再現できるポテンシャルモデルを用いることが一般的である。このため、分子動力学計算の結果を用いて評価された物性値は、臨界点近傍の温度・圧力領域では、評価誤差が大きくなる場合がある。
【００３２】
上述の課題を解決するため、本発明は、臨界点近傍の温度・圧力領域においても、精度良い物性値評価を実現することを目的とする。
【課題を解決するための手段】
【００３３】
上記目的を解決するため、本発明は、複数の粒子からなる計算体系を模擬した分子動力学計算を用いて、評価温度および評価圧力における物性値を評価する物性値評価方法において、前記評価温度および前記評価圧力に基づいて解析条件温度および解析条件圧力を求める解析条件補正工程と、前記計算体系の温度を前記解析条件温度とし、前記計算体系の圧力を前記解析条件圧力とする分子動力学計算を行う分子動力学計算工程と、前記分子動力学計算で得られた結果に基づいて前記評価温度および前記評価圧力における評価対象物性値を算出する物性値評価工程と、を有することを特徴とする。
【００３４】
また、本発明は、複数の粒子からなる計算体系を模擬した分子動力学計算を用いて、評価温度および評価圧力における物性値を評価する物性値評価装置において、前記評価温度および前記評価圧力に基づいて解析条件温度および解析条件圧力を求める解析条件補正手段と、前記計算体系の温度を前記解析条件温度とし、前記計算体系の圧力を前記解析条件圧力とする分子動力学計算を行う分子動力学計算手段と、前記分子動力学計算で得られた結果に基づいて前記評価温度および前記評価圧力における評価対象物性値を算出する物性値評価手段と、を有することを特徴とする。
【００３５】
また、本発明は、複数の粒子からなる計算体系を模擬した分子動力学計算を用いて、評価温度および評価圧力における物性値を評価させる物性値評価プログラムにおいて、コンピュータに、前記評価温度および前記評価圧力に基づいて解析条件温度および解析条件圧力を求めさせる解析条件補正手順と、前記計算体系の温度を前記解析条件温度とし、前記計算体系の圧力を前記解析条件圧力とする分子動力学計算を行わせる分子動力学計算手順と、前記分子動力学計算で得られた結果に基づいて前記評価温度および前記評価圧力における評価対象物性値を算出させる物性値評価手順と、を実行させることを特徴とする。
【発明の効果】
【００３６】
本発明によれば、臨界点近傍の温度・圧力領域においても、精度良い物性値評価を実現できる。
【発明を実施するための最良の形態】
【００３７】
本発明に係る物性値評価装置の実施の形態を、図面を参照して説明する。
【００３８】
図２は、本発明に係る一実施の形態の物性値評価装置のブロック図である。
【００３９】
物性値評価装置は、本体１、入力手段２および出力手段３を有している。本体１は、解析条件入力部４、解析条件補正部５、分子動力学計算部６、物性値評価部７、および、記憶手段８を有している。
【００４０】
本体１は、たとえば１台のコンピュータであるが、複数のコンピュータを組み合わせたものであってもよい。なお、本体１に含まれる解析条件入力部４、解析条件補正部５、分子動力学計算部６、および、物性値評価部７は、それぞれ独立して存在している必要はない。入力手段２は、たとえばキーボードやマウスである。出力手段３は、たとえば液晶ディスプレーやプリンタである。
【００４１】
図１は、本発明に係る一実施の形態の物性値評価の流れ図である。
【００４２】
この物性値評価は、解析条件入力、解析条件補正、分子動力学計算および物性値評価の４つの段階を有している。
【００４３】
解析条件入力の段階では、入力手段２を介して解析条件入力部４に、原子・分子データが入力される（工程Ｓ１）。また、同様に、評価対象の温度および圧力などの解析条件が入力される（工程Ｓ２）。
【００４４】
解析条件補正の段階では、まず、入力手段２を介して解析条件補正部５に臨界温度および臨界圧力などの臨界点データが入力される（工程Ｓ３）。解析条件補正部５は、この臨界点データを使用して、温度および圧力などの解析条件を補正する（工程Ｓ３）。
【００４５】
この補正は、分子動力学計算の入力となる解析条件温度をＴ’＝（Ｔ’ｃ／Ｔｃ）Ｔとして求め、解析条件圧力をＰ’＝（Ｐ’ｃ／Ｐｃ）Ｐとして求めるものである。ここで、ＴｃとＰｃは、物性値評価の対象とする粒子の系の実際の臨界温度と臨界圧力である。また、Ｔ’ｃとＰ’ｃは、分子動力学を用いて表現している系の臨界温度と臨界圧力である。ＴおよびＰは、物性値を評価したい温度および圧力の条件である。
【００４６】
分子動力学計算の段階において、分子動力学計算部６は、原子・分子データと解析条件温度および解析条件圧力を用いて、分子動力学法により原子・分子の挙動を計算する（工程Ｓ５）。このとき、粒子間のポテンシャルとして、上述の式（１）ないし式（４）で表される、ＬＪポテンシャル、クーロンポテンシャル、結合ポテンシャル、結合角ポテンシャルなどが用いられる。
【００４７】
分子動力学計算では、設定されたタイムステップ毎に、すべての粒子対に対して上記ポテンシャル値を計算し、そこから算出される各粒子に働く力を用いて運動方程式を数値的に解いて、すべての粒子の位置、速度、加速度などの値を時系列データとして出力する。分子動力学計算の結果は、記憶手段８に出力されて保存される。
【００４８】
物性値評価段階において、物性値評価部７は、分子動力学計算で得られた計算データから物性値を計算する。物性値評価部７は、まず、記憶手段７から粒子の位置・速度等のデータを読み込み（工程Ｓ７）、全粒子について運動エネルギー、位置エネルギー、相互作用などを計算する（工程Ｓ８）。
【００４９】
その後、物性値評価部７は、計算された各粒子の位置・速度・加速度についての時系列データに対して統計処理を施して、拡散係数、粘度および熱伝導度などの物性値を評価する（工程Ｓ９）。拡散係数、粘度および熱伝導度などの物性値は、上述の通り式（５）ないし式（１６）を用いて評価できる。
【００５０】
評価された物性値は、記憶手段８に保存されるとともに、出力手段３に出力される（工程Ｓ１０）。
【００５１】
臨界点付近で、状態方程式からのずれを表す圧縮率因子を、臨界温度と臨界圧力で規格化した温度と圧力の関数として求めると、気体の種類によらず、一致する傾向がある。臨界温度と臨界圧力で規格化した温度と圧力の関数として圧縮率因子を求めると、実在の物質と分子動力学のモデルで表現された物質の間においても、一致すると考えられる。これは、両者の状態方程式が一致することを意味する。
【００５２】
状態方程式は、温度と圧力を与えたときの密度を求める式とみなすことができる。また、拡散係数、粘度、熱伝導度などの物性値は、密度の関数として表すことができる。したがって、臨界温度と臨界圧力を用いて補正することにより、分子動力学による評価結果と実際の物性値の間によい一致が期待できる。
【００５３】
図３は、本実施の形態の物性評価装置によって評価した拡散係数の圧力依存性を示すグラフである。
【００５４】
図３の横軸は圧力を臨界圧力で除した換算圧力であり、縦軸は自己拡散係数である。図４の横軸は圧力であり、縦軸は自己拡散係数である。また、図３において、実線、破線および一点鎖線で示した自己拡散係数は、式（１７）で表される自己拡散係数の実験値のフィッティング式によって求めたものである。
【００５５】
また、図３において、実線、破線および一点鎖線で示した自己拡散係数は、換算温度がそれぞれ１．０００、１．０５３および１．２０３のときの値である。ここで、換算温度とは、温度を臨界温度で除した値である。すなわち、実線、破線および一点鎖線においては、温度を臨界温度の実測値である６４７．０９６Ｋで除したものであり、換算温度１．０００、１．０５３および１．２０３の実際の温度は、それぞれ、６４７．０９６Ｋ、６８２Ｋおよび７７８Ｋとなる。
【００５６】
また、図３において、丸（○）、四角（□）および三角（△）で示した自己拡散係数は、本実施の形態の物性評価装置によって評価したものであって、換算温度がそれぞれ１．０００、１．０５３および１．２０３に相当する。本実施の形態の物性評価装置によって評価した自己拡散係数においては、換算温度は、本実施の形態の物性評価装置によって評価した臨界温度である６６５Ｋで除したものであり、換算温度１．０００、１．０５３および１．２０３に対応する解析条件の温度は、それぞれ、６６５Ｋ、７００Ｋおよび８００Ｋとなる。
【００５７】
図４から、解析条件温度および解析条件圧力の補正を行わない分子動力学計算による結果に基づいて評価した自己拡散係数は、臨界点（６４７Ｋ、２２ＭＰａ）の近傍で、実験値のフィッティング式で求めた自己拡散係数と大きく乖離していることがわかる。一方、図３に示されるように、本実施の形態の物性値評価装置によって求めた自己拡散係数は、換算圧力が１となる臨界点近傍であっても、実験値のフィッティング式で求めた自己拡散係数とよく一致していることがわかる。
【００５８】
このように、臨界点の実測値および分子動力学計算による評価値によって、解析条件温度および解析条件圧力を補正することによって、臨界点近傍であっても、高い精度で自己拡散係数を評価することができる。
【００５９】
粘度、熱伝導度などの密度に依存する物性値も、解析条件温度と解析条件圧力を補正することにより、同様の効果が期待できる。
【００６０】
なお、以上の説明は単なる例示であり、本発明は上述の実施の形態に限定されず、様々な形態で実施することができる。
【図面の簡単な説明】
【００６１】
【図１】本発明に係る一実施の形態の物性評価方法の流れ図である。
【図２】本発明に係る一実施の形態の物性評価装置のブロック図である。
【図３】本発明に係る一実施の形態の物性評価装置によって評価した拡散係数の圧力依存性を示すグラフである。
【図４】従来の分子動力学計算に基づいて評価した水の自己拡散係数を、実験値のフィッティング式によって求めたものと併せて示したグラフである。
【符号の説明】
【００６２】
１…本体、２…入力手段、３…出力手段、４…解析条件入力部、５…解析条件補正部、６…分子動力学計算部、７…物性値評価部、８…記憶手段
【図１】
2007248273.tif 000016
【図２】
2007248273.tif 000017
【図３】
2007248273.tif 000018
【図４】
2007248273.tif 000019