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

Description

本発明は、複雑な分子構造を有する材料系、とりわけ樹脂材料系の平衡分子動力学法による計算から、高精度で安定して熱伝導率を求める算出方法、算出装置に関し、さらにこの方法をコンピュータに実行させるための算出プログラムに関する。

絶縁性材料の熱伝導率は、格子振動によって熱の移動が行われることにより規定される。そして、この格子振動は、材料を構成する各々の原子の動きで表現でき、材料の熱伝導率を算出するためには、材料を構成する多数の原子の動きを計算する必要がある。多数の原子の動きを計算する手法として、分子動力学法が知られている。

分子動力学法には、大別して非平衡分子動力学法（ＮＥＭＤ、Ｎｏｎ−Ｅｑｕｉｌｉｂｒｉｕｍ Ｍｏｌｅｃｕｌａｒ Ｄｙｎａｍｉｃｓ）、および平衡分子動力学法（ＥＭＤ、Ｅｑｕｉｌｉｂｒｉｕｍ Ｍｏｌｅｃｕｌａｒ Ｄｙｎａｍｉｃｓ）に基づく研究開発が行われている。ＮＥＭＤ法による熱伝導率シミュレーションは、計算モデルに対して温度勾配を与えた条件で熱流束を計算し、フーリエの法則に基づき、温度勾配と熱流束の関係から熱伝導率を算出する。一方、ＥＭＤ法による熱伝導率シミュレーションは、熱平衡系で熱流束を計算し、グリーン−久保公式に基づき、熱流束の自己相関関数から熱伝導率を算出する。

ＮＥＭＤ法は、系に直接温度勾配を与えるために、系の熱移動を直接評価できることがメリットである。一方、ＮＥＭＤ法では、系に温度勾配を課すため、複雑な材料では境界条件の取り扱いが難しい。例えば、ＮＥＭＤ法では、熱浴系と測定系の境界が非物理的な接触面となる。ＮＥＭＤ法では、この非物理的な境界効果を除去する必要があるため、基本的に計算モデルを大きくとる必要がある。それでもなお、実際に計算可能なモデルサイズでは、非現実的な境界条件（超薄膜に大きな温度差を印加した条件）に相当する問題が残される。

また、系の温度勾配を安定化させるために長時間計算を要することも、ＮＥＭＤ法の弱点である。これらの問題があるために、ＮＥＭＤ法による熱伝導率シミュレーションは、比較的簡単な分子構造であるカーボンナノチューブの１本鎖等への適用に限定されている。このように、現状では、ＮＥＭＤ法は、複雑な樹脂材料のマクロな熱伝導率を算出するためには課題が多いといえる。

これに対して、ＥＭＤ法は、周期境界条件下の熱平衡系でシミュレーションできるため、材料のマクロな熱伝導率を評価できることが主な特徴である。また、ＥＭＤ法は、ＮＥＭＤに比して小さなモデルサイズでも熱伝導率が評価できるため、現実的な計算時間内で熱伝導率を評価できることがメリットである。従って、複雑な構造を有する樹脂材料の熱伝導を評価するためには、ＥＭＤ法によるシミュレーション手法が適していると考えられる。

ただし、ＥＭＤ法で実際に熱伝導率を評価する場合には、熱流束の自己相関関数を計算し、その時間挙動を評価解析することによって、熱伝導率を算出する必要がある。ここで、ＥＭＤ法を用いた各時刻ｔに対する熱流束Ｊ（ｔ）の出力データから熱伝導率を算出する式について簡単に説明する。熱伝導率λは、熱流束Ｊ（ｔ）の自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞から、グリーン−久保公式に基づき、下式（１４）の関係で結びつけられる。

上式（１４）において、Ｖは体積、Ｔは温度、ｋＢはボルツマン定数、＜＞はサンプル平均を表している。熱流束Ｊ（ｔ）は、ＥＭＤ法によって計算することができる。熱流束Ｊ（ｔ）は、下式（１５）で表されるベクトル量である。

上式（１５）において、Σｉは、粒子数をＮとしたときの、ｉ＝１からＮまでの総和を表す。Ｅ_ｉ（ｔ）は、時刻ｔにおける粒子ｉのエネルギー、ｒ_ｉｊ（ｔ）＝ｒ_ｉ（ｔ）−ｒ_ｊ（ｔ）は、時刻ｔにおける粒子ｉと粒子ｊの位置ベクトル差、ｖ_ｉ（ｔ）は時刻ｔにおける粒子ｉの速度ベクトル、Ｆ_ｉｊ（ｔ）＝Ｆ_ｉ（ｔ）−Ｆ_ｊ（ｔ）は時刻ｔにおける粒子ｉと粒子ｊの力ベクトル差である。

しかし、ＥＭＤ法において、実際に計算できるステップ数は、有限であり、上式（１４）の積分範囲も有限な区間で行うことしかできない。このため、上式（１４）の算出方法では、積分範囲に熱伝導率λの計算値が依存することが問題となる。特に、熱伝導率λは、自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞の長時間の挙動に支配されるため、時間の大きい側の積分区間の区切りに対して、計算値は敏感になる。

そこで、実際のＥＭＤ計算では、有限時間しか取り扱えないことに起因して、積分区間の選択を試行錯誤して熱伝導率λを算出することを想定し、本問題を作業性の観点から課題と捉えている。これを解決するために、コンピュータ上の表示処理を工夫し、積分範囲を効率的に選択できる算出装置がある（例えば、特許文献１参照）。

また、ＥＭＤ法によって出力される自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞は、全ての振動モード（フォノン）の情報が含まれるため、複雑で激しい振動形状を示すことが知られている。

例えば、ダイヤモンド結晶のＥＭＤ計算を行い、自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞を算出している従来技術がある（例えば、非特許文献１参照）。自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞のｔ依存性は、非常に激しい振動を呈する結果が示されている。自己相関関数の振動は、音響フォノンと光学フォノンの熱振動が自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞に混合して現れているためであり、激しい振動は、ダイヤモンド結晶の光学フォノンに相当する。

熱伝導率λは、音響フォノンに起因しており、光学フォノンに基づく激しい振動成分を除去する必要がある。非特許文献１では、自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞を、音響フォノンと光学フォノンを独立した緩和過程で表現する簡易なモデル化を行い、自己相関関数から音響フォノンを選択的に抽出し、ダイヤモンド結晶の熱伝導率を定量的に求める算出方法を開示している。

特開平５−１５１１９１号公報

Ｊ．Ｃｈｅ、Ｔ．Ｃａｇｉｎ、Ｗ．Ｄｅｎｇ、ａｎｄ Ｗ．Ａ．Ｇｏｄｄａｒｄ ＩＩＩ、Ｊ．Ｃｈｅｍ．Ｐｈｙｓ．、１１３ ６８８８ （２０００）．

しかしながら、従来技術には、以下のような課題がある。
特許文献１で問題として取り上げられている通り、実際のＥＭＤ計算では、有限時間しか取り扱えないため、上式（１４）では、積分範囲の取り方に熱伝導率の計算結果が強く依存する。この特許文献１では、コンピュータ上の表示処理を工夫することで、本問題の解決を作業性の観点から図っている。しかしながら、熱伝導率を定量的に算出する算出方法として、上式（１４）の算出方法に係わる抜本的な問題解決は、図られていない。

また、非特許文献１で示されたダイヤモンド結晶のような比較的単純な分子構造であっても、自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞は、複雑な時間依存性を示す。特に、自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞の複雑な時間依存性は、計算対象の分子構造の複雑さと相関している。従って、複雑な内部構造を有する樹脂材料系においては、自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞は、ダイヤモンド結晶と比べて、さらに複雑な挙動を示す。このため、自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞の簡易なモデル化手法、およびこのようなモデルを複雑化した手法では、熱伝導率を定量的に算出することが困難と考えられる。

従って、実際のＥＭＤ計算では、有限のステップ数しか取り扱えないことに起因する問題と、計算対象の分子構造に起因する自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞の複雑な挙動を評価解析する困難さの問題、の２つの課題がある。

本発明は、前記のような課題を解決するためになされたものであり、複雑な分子構造を有する材料系において、精度よく短時間で熱伝導率を定量的に評価することができる輸送係数の算出方法、輸送係数の算出装置、および輸送係数の算出プログラムを得ることを目的とする。

本発明に係る輸送係数の算出方法は、原子あるいは分子の構造データおよび温度計算条件を入力するステップと、所定の時間ｔに対して力学量Ｊ（ｔ）を算出するステップと、力学量Ｊ（ｔ）を用いて下式（１６）

に基づいて周波数スペクトルΠ（ω、ω_ｃ）を算出するステップと、周波数スペクトルΠ（ω、ω_ｃ）を用いて下式（１７）

に基づいて重み付き周波数成分Ｓ（ω_ｃ）を算出するステップと、重み付き周波数成分Ｓ（ω_ｃ）を用いて下式（１８）

に基づいて輸送係数Ｋを算出するステップとを備えたものである。

なお、輸送係数として熱伝導率を算出する場合、力学量Ｊ（ｔ）は上式（１５）で現れる力学量を用い、上式（１７）の定数Ｃは、下式（１９）を用いる。

本発明に係る輸送係数の算出方法、輸送係数の算出装置、および輸送係数の算出プログラムによれば、平衡分子動力学法から得られる熱流束の自己相関関数に対して、先ず、低周波数透過フィルターを適用し、熱流束の自己相関関数の中から熱伝導特性に寄与しない高周波の格子振動を自己相関関数から除去し、次に、低周波数透過フィルターを適用した自己相関関数に関する周波数積分値を求め、さらに、低周波数透過フィルターのカットオフ周波数のゼロ極限を求めることによって外挿値を求めることにより、複雑な分子構造を有する材料系において、精度よく短時間で熱伝導率を定量的に評価することができる輸送係数の算出方法、輸送係数の算出装置、および輸送係数の算出プログラムを得ることができる。

本発明の実施の形態２における輸送係数の算出方法のフローチャートである。 本発明の実施の形態２における輸送係数の算出装置の構成図である。 本発明の実施の形態３における平衡分子動力学シミュレータの処理ステップＳ１０１によって出力されたポリエチレンのモデル結晶における熱流束Ｊ（ｔ）の自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞の計算結果である。 本発明の実施の形態３において、図３の自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞を、非特許文献１の算出方法に従ってフィティングした結果を示した図である。 本発明の実施の形態３における熱伝導率の算出方法を、ポリエチレンのモデル結晶における結晶状態での自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞へ適用した結果を示した図である。 本発明の実施の形態３における熱伝導率の算出方法を、ポリエチレンのモデル結晶におけるアモルファス状態での自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞へ適用した結果を示した図である。 本発明の実施の形態３において用いた低周波数透過フィルターの周波数特性を示した図である。 本発明の実施の形態４における輸送係数の算出方法のフローチャートである。 本発明の実施の形態４における輸送係数の算出装置の構成図である。 本発明の実施の形態６における輸送係数の算出方法のフローチャートである。 本発明の実施の形態６における輸送係数の算出装置の構成図である。 本発明の実施の形態７における輸送係数の算出方法のフローチャートである。 本発明の実施の形態７における輸送係数の算出装置の構成図である。

以下、本発明の輸送係数の算出方法、輸送係数の算出装置、および輸送係数の算出プログラムの好適な実施の形態につき図面を用いて説明する。

実施の形態１．
本実施の形態１では、輸送係数として熱伝導率λを算出する例を示す。熱伝導率λを算出する手法として、上式（１６）のΠ（ω、ω_ｃ）について、単純にωに関するゼロ極限を求める手法を用いることも考えられる。しかしながら、複雑な内部構造を有する樹脂材料では、低周波透過フィルターｆω_ｃ（ｔ）だけでは、熱伝導率λに寄与しない格子振動を完全には除去できない。このため、Π（ω、ω_ｃ）について単純にωに対するゼロ極限を求める方法は、数値的に安定しないことを、我々は確認している。

また、実際のＥＭＤ計算では有限時間しか取り扱えないため、実際に周波数ωをゼロに十分に近づけることが難しい。この問題は、上式（１４）に基づく算出方法において、実際のＥＭＤ計算では有限の積分区間しか取扱えないことに起因した熱伝導率λの数値的不安定さと本質的には同じ問題である。

これを解決するためには、長時間計算が必要となるが、それでもなお、実際には有限の計算時間しか取り扱えない。このため、ＥＭＤ計算部の長時間計算は、特許文献１と同様に、本質的な解決手段とはならない。このため、Π（ω、ω_ｃ）に対して単純にωのゼロ極限を求める手法だけでは、熱伝導率を定量的に求めることは容易でない。

一方、本発明の場合、上式（１７）において積分されるΠ（ω、ω_ｃ）は、低周波透過フィルターによってωの高周波数成分が減じられている。このため、上式（１７）において有限な積分区間しか取れないことに起因するＳ（ω_ｃ）の不定性は、ほとんどない。つまり、上式（１７）によるＳ（ω_ｃ）は、数値的に安定して算出することができる。

また、本実施の形態１では、低周波透過フィルターでも完全には除去できない熱運動による格子振動も含めた積分値Ｓ（ω_ｃ）を算出し、ω_ｃについて外挿を行う。すなわち、本工程は、熱伝導に寄与しないが容易に自己相関関数から取り除けない格子振動を徐々に除去していくことを意味している。従って、本発明に係わる上式（１７）の算出方法においては、計算モデルの複雑な内部構造に起因する自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞の複雑な挙動が、Ｓ（ω_ｃ）に多少含まれていてもよい。

さらに、本実施の形態１は、Ｓ（ω_ｃ）のカットオフ周波数ω_ｃについて外挿する。このため、実際のＥＭＤ計算では有限時間しか取り扱えない問題に、強く影響されない。この結果、比較的短時間のＥＭＤ計算により求められた自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞から、精度よく熱伝導率λを算出することができる。

以上のように、実施の形態１によれば、実際のＥＭＤ計算において有限時間しか取り扱えないことに起因する問題と、計算対象の複雑な分子構造に起因する自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞の複雑な挙動を評価解析する困難さの問題とを、同時に解決することができる。

実施の形態２．
先の実施の形態１では、数式に基づいて、本発明の輸送係数の算出方法、輸送係数の算出装置、および輸送係数の算出プログラムを説明した。これに対して、本実施の形態２では、図面に基づいて、本発明の輸送係数の算出方法、輸送係数の算出装置、および輸送係数の算出プログラムについて、詳細に説明する。

図１は、本発明の実施の形態２における輸送係数の算出方法のフローチャートであり、輸送係数として熱伝導率λを算出する場合を示している。この図１のフローチャートは、平衡分子動力学法シミュレータによる処理ステップＳ１０１と、熱伝導率計算シミュレータの処理ステップＳ２０１とで構成されている。

さらに、処理ステップＳ１０１は、Ｓ１１〜Ｓ１４に細分化され、処理ステップＳ２０１は、Ｓ２１〜Ｓ２６に細分化されている。処理ステップＳ１０１における平衡分子動力学法シミュレータは、原子・分子の構造データと温度などの計算条件を入力し（ステップＳ１１、Ｓ１２）、平衡分子動力学計算を実施する（ステップＳ１３）。そして、平衡分子動力学計算によって算出される熱流束Ｊ（ｔ）の時系列データを熱伝導算出シミュレータの処理ステップＳ２０１に対して出力する（ステップＳ１４）。この時系列データを入力することで、処理ステップＳ２０１により熱伝導率が算出される。

このように、処理ステップＳ１０１における平衡分子動力学法シミュレータは、原子の初期座標、初期速度、原子間ポテンシャル関数、原子の質量、原子の電荷等を入力し、平衡分子動力学法に基づいて、原子・分子の座標や速度などの時間発展（時系列データ）を出力するものである。なお、平衡分子動力学法の計算手法は、一般的な教科書に記述されているので説明を省略する。

次に、処理ステップＳ２０１における熱伝導算出シミュレータ内の各ステップＳ２１〜Ｓ２６について説明する。ステップＳ２１は、熱流束Ｊ（ｔ）の自己相関関数のサンプル平均値＜Ｊ（ｔ）・Ｊ（０）＞を算出するステップである。また、ステップＳ２２は、低周波数透過フィルターのカットオフ周波数ω_ｃを入力するステップである。そして、ステップＳ２３は、本発明に係わる上式（１６）を用いて低周波数透過フィルターを適用された自己相関関数の周波数スペクトルΠ（ω、ω_ｃ）を算出するステップである。

次に、ステップＳ２４は、本発明に係わる上式（１７）を用いて低周波数透過フィルターのカットオフ周波数ω_ｃに対してＳ（ω_ｃ）を算出するステップである。次に、ステップＳ２５は、本発明に係わる上式（１８）を用いて、複数のカットオフ周波数ω_ｃに対して計算されたＳ（ω_ｃ）から、低周波数透過フィルターのカットオフ周波数ω_ｃに対するゼロ極限値を外挿するステップである。さらに、ステップＳ２６は、求められた熱伝導率λを出力するステップである。

図２は、本発明の実施の形態２における輸送係数の算出装置の構成図であり、図１のフローチャートに対応した処理を行う装置構成を示している。この図２における輸送係数の算出装置は、入力部１０、処理部２０、および出力部３０で構成されている。さらに、処理部２０は、４つの算出部２１〜２４（第１算出部〜第４算出部に相当）で構成されている。

入力部１０は、原子あるいは分子の少なくとも座標データおよび計算条件等の入力部であり、図１のステップＳ１１、Ｓ１２に相当する。処理部２０内の第１算出部２１は、力学量Ｊ（ｔ）を算出するための平衡分子動力学法に基づく熱流束Ｊ（ｔ）を求める算出部であり、図１のステップＳ１３、Ｓ１４に相当する。

第２算出部２２は、熱流束Ｊ（ｔ）の自己相関関数のサンプル平均値＜Ｊ（ｔ）・Ｊ（０）＞の算出結果と、低周波数透過フィルターのカットオフ周波数ω_ｃとを用いて、上式（１６）に基づく周波数スペクトルΠ（ω、ω_ｃ）を求める算出部であり、図１のステップＳ２１〜Ｓ２３に相当する。

第３算出部２３は、周波数スペクトルΠ（ω、ω_ｃ）を用いて上式（１７）に基づく重み付き周波数成分Ｓ（ω_ｃ）を求める算出部であり、図１のステップＳ２４に相当する。第４算出部２４は、重み付き周波数成分Ｓ（ω_ｃ）を用いて上式（１８）に基づいて熱伝導率λを求める算出部であり、図１のステップＳ２５に相当する。さらに、出力部３０は、熱伝導率λの出力部であり、図１のステップＳ２６に相当する。

なお、本発明に係る熱伝導算出装置は、図２に示した詳細な構成に限定されるものではない。例えば、平衡分子動力学法に基づく熱流束Ｊ（ｔ）の第１算出部２１から得られる熱流束の時系列データを、外部記憶媒体に出力し、出力された熱流束Ｊ（ｔ）の時系列データを改めて入力する入力処理部を備えることでも実現できる。

この場合には、入力した時系列データを用いて、処理部２０内の第２算出部２２〜第４算出部２４で演算処理することで、出力部３０から熱伝導率λを出力できる。

以上のように、実施の形態２によれば、実際のＥＭＤ計算において有限時間しか取り扱えないことに起因する問題と、計算対象の複雑な分子構造に起因する自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞の複雑な挙動を評価解析する困難さの問題とを、同時に解決することができる。

実施の形態３．
本実施の形態３では、本発明を汎用樹脂材料であるポリエチレン材料に適用した例を述べる。ポリエチレン材料は、マクロな形態として、結晶状態と非晶質（アモルファス）状態を取ることが知られている。このため、結晶状態とアモルファス状態の両者について、本発明を適用した。なお、本発明の算出方法を適用する材料系としては、高分子化合物、分子集合体、高分子集合体、結晶性化合物、および非結晶性化合物を挙げることができ、分子の形態に係わらず、本発明に係わる算出方法を適用することができる。

そこで、まず始めに、ポリエチレン材料の結晶状態について本発明を適用した例を述べる。ポリエチレン結晶のモデル結晶は、周期境界条件を課したセルサイズ（１６２．３Å×９．６Å×１４．８Å）に原子数３０７２、ポリマー数８のポリエチレン鎖を配置、質量密度１．０ｇ／ｃｃの計算モデルを用いた。

これらの原子・分子データを先の図１におけるステップＳ１１を通じて平衡分子動力学計算のステップＳ１３へ入力した。平衡分子動力学計算のステップＳ１３は、ＮＶＴアンサンブル（粒子数（Ｎ）、体積（Ｖ）、温度（Ｔ）を一定とした熱力学的条件）において、温度３００Ｋ、時間刻み幅は１ｆｓの計算条件で実行した。この計算条件において、系を熱平衡状態にするために、平衡分子動力学計算をステップ数５×１０^４ステップ、熱伝導率を算出するために平衡分子動力学計算をステップ数３×１０^５ステップ実行した。

図３は、本発明の実施の形態３における平衡分子動力学シミュレータの処理ステップＳ１０１によって出力されたポリエチレンのモデル結晶における熱流束Ｊ（ｔ）の自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞の計算結果である。図３から、先に説明したように、自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞は、複雑な挙動を示すことが確認できる。

図３のように、激しい振動を呈する自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞を、上式（１４）に基づいて単純に時間積分したとしても、激しい振動積分自体が精度よく求められない。このことからも、上式（１４）に基づいた熱伝導率の評価解析手法（例えば、特許文献１）では、熱伝導率を定量的に算出できないことが容易に理解できる。

また、図３のような激しい振動を有する自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞を、何らかのモデル式で近似し、自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞をモデル式に対してフィティングする計算方法を用いた場合にも、熱伝導率の定量的算出の困難さが容易に想像できる。

図４は、本発明の実施の形態３において、図３の自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞を、非特許文献１の算出方法に従ってフィティングした結果を示した図である。図４に示すように、フィティングの対象となるデータのばらつきが大きいため、熱伝導率の定量的算出が困難であることが確認できる。実際に、非特許文献１の方法を適用して熱伝導率を算出すると、０．９Ｗ／（ｍ・Ｋ）となる。ポリエチレン結晶の実験値は、測定サンプルの作成方法等に強く依存するが、概ね１０Ｗ／（ｍ・Ｋ）以上であり、非特許文献１の方法による計算値は、実験値とよく一致しない結果となる。

これに対して、図５は、本発明の実施の形態３における熱伝導率の算出方法を、ポリエチレンのモデル結晶における結晶状態での自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞へ適用した結果を示した図である。より具体的には、上式（１６）を用いた先の図１におけるステップＳ２３の算出を行った後、上式（１７）を用いた先の図１におけるステップＳ２４の算出を行った結果であり、横軸に低周波数透過フィルターのカットオフ周波数ωｃ、縦軸に上式（１７）の左辺Ｓ（ω_ｃ）をプロットしたものである。

先の図１におけるステップＳ２５に示したように、上式（１８）から、図５においてカットオフ周波数ω_ｃのゼロ極限を外挿する算出方法を適用することで、熱伝導率λが算出される。このため、実際にステップＳ２５の算出方法を適用して熱伝導率λを算出するためには、図５に示される計算結果のばらつきが小さいことが重要である。図５に示すとおり、Ｓ（ω_ｃ）は、ω_ｃに対してばらつきが少ないデータとなっていることが確認できる。

実際に、図５の計算結果をω_ｃに対して外挿、すなわち、Ｓ（ω_ｃ）のω_ｃゼロ極限を求めた結果、ポリエチレンのモデル結晶の熱伝導率は、１１．８Ｗ／（ｍ・Ｋ）となった。本発明に係わる算出方法を適用した計算値は、概ね１０Ｗ／（ｍ・Ｋ）以上となった実験値と良好に一致している。このように、先の実施の形態１、２で説明した本発明の算出方法により、ポリエチレン結晶の熱伝導率を定量的に算出できることが確認できた。

次に、ポリエチレン材料のアモルファス状態に本発明を適用した例を述べる。ポリエチレンアモルファスのモデル非晶質は、周期境界条件を課したセルサイズ（３６．８Å×３６．８Å×３６．８Å）に原子数１５３６、ポリマー数４のポリエチレン鎖を配置、質量密度０．９５ｇ／ｃｃの計算モデルを用いた。

これらの原子・分子データを先の図１におけるステップＳ１１を通じて平衡分子動力学計算Ｓ１３へ入力した。平衡分子動力学計算のステップＳ１３は、ＮＶＴアンサンブル（粒子数（Ｎ）、体積（Ｖ）、温度（Ｔ）を一定とした熱力学的条件）において、温度３００Ｋ、時間刻み幅は１ｆｓの計算条件で実行した。この計算条件において、系を熱平衡状態にするために、平衡分子動力学計算をステップ数５×１０^４ステップ、熱伝導率を算出するために平衡分子動力学計算をステップ数３×１０^５ステップ実行した。

図６は、本発明の実施の形態３における熱伝導率の算出方法をポリエチレンのモデル結晶におけるアモルファス状態での自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞へ適用した結果を示した図である。より具体的には、上式（１６）を用いた先の図１におけるステップＳ２３の算出を行った後、上式（１７）を用いた先の図１におけるステップＳ２４の算出を行った結果であり、横軸に低周波数透過フィルターのカットオフ周波数ωｃ、縦軸に（１７）式の左辺Ｓ（ω_ｃ）をプロットしたものである。

図６から、Ｓ（ω_ｃ）は、ω_ｃに対してばらつきが少ないデータとなっていることが確認できる。従って、先の図１におけるステップＳ２５に示したように、カットオフ周波数ω_ｃのゼロ極限を外挿する算出方法を適用することにより、精度よく外挿値を求めることが実現可能であることが確認できた。

実際に、図６の計算結果をω_ｃに対してフィティングして、Ｓ（ω_ｃ）のω_ｃゼロ極限を求めた結果、ポリエチレンのアモルファスモデルの熱伝導率は、０．３６Ｗ／（ｍ・Ｋ）となった。ポリエチレンアモルファスの熱伝導率の実験値は、測定サンプルの作成方法等に依存して多少のばらつきがあるが、概ね０．４Ｗ／（ｍ・Ｋ）程度であり、計算値は、実験値と良好に一致している。ポリエチレンアモルファスについても、本発明に係わる算出方法により熱伝導率を定量的に算出できることが確認できた。

以上のように、実施の形態３によれば、本発明の算出方法を用いることで、ポリエチレン材料のマクロな形態、すなわち、結晶状態やアモルファス状態の熱伝導特性の差を明確に再現きることが確認できた。

なお、図７は、本発明の実施の形態３において用いた低周波数透過フィルターの周波数特性を示した図である。しかしながら、本発明に係る算出方法は、低周波数透過フィルターのゲイン特性の詳細に依存するものではなく、図７に示した低周波数透過フィルターに限定されるものではない。

実施の形態４．
上述した実施の形態１〜３は、熱伝導率を算出するために時系列データとして熱流束Ｊ（ｔ）の時系列データを想定している。このため、熱流束Ｊ（ｔ）の時系列データから熱伝導率を算出する処理ステップＳ２０１を、便宜上、熱伝導率計算シミュレータと称していた。しかしながら、本発明に係る算出方法を含む処理ステップＳ２０１は、熱伝導率の算出方法だけに限定されるものではない。

この処理ステップＳ２０１は、例えば、体積粘性率、電気伝導率、ずれ粘性係数などの輸送係数の定量的算出にも、一般的に適用できるものである。この場合、分子動力学法シミュレータから出力される時系列データは、算出する輸送係数と対応した力学量となる。算出する輸送係数と力学量との対応関係は、一般的な教科書に記述されている内容であるが、例えば、電気伝導率の場合、時系列でデータＪ（ｔ）は、下式（２０）となる。

図８は、本発明の実施の形態４における輸送係数の算出方法のフローチャートである。また、図９は、本発明の実施の形態４における輸送係数の算出装置の構成図であり、図８のフローチャートに対応した処理を行う装置構成を示している。

本実施の形態４の図８における輸送係数を算出するフローチャート、および図９における算出装置の構成と、先の図１における熱伝導率を算出するフローチャート、および図２における算出装置の構成との差異は、図８、図９において添字ａが付された部分である。より具体的には、本実施の形態４における図８、図９では、力学量に応じた輸送係数Ｋの定数Ｃを扱っている点が異なっている。

定数Ｃは、平衡分子動力学シミュレータにおける温度や体積などの原子・分子データや、計算条件からの数値と自然定数で構成される。輸送係数の算出装置も同様に、求めたい輸送係数Ｋに対応する力学量を計算する力学量算出処理部が有ればよい。

以上のように、実施の形態４によれば、体積粘性率、電気伝導率、ずれ粘性係数などの輸送係数の定量的算出においても、実際のＥＭＤ計算において有限時間しか取り扱えないことに起因する問題と、計算対象の複雑な分子構造に起因する自己相関関数＜Ｊ（ｔ）・Ｊ（０）＞の複雑な挙動を評価解析する困難さの問題とを、同時に解決することができ、熱伝導率の場合と同様の効果を得ることができる。

実施の形態５．
先の実施の形態１〜４は、系全体の輸送係数、すなわち、系全体の熱伝導率あるいは電機伝導率の算出方法について説明した。これに対して、本実施の形態５では、輸送係数と分子構造との関係を把握するために、原子毎もしくは原子団毎における輸送係数の割合（原子毎の輸送係数）を算出する計算方法を、数式に基づいて説明する。なお、図面を用いた具体的な説明は、実施の形態６以降で行う。

平衡分子動力学計算によって算出される原子毎の力学量Ｊ_ｉ（ｔ）の時系列データを、熱伝導率算出シミュレータＳ３０１に入力することで輸送係数を算出する。例えば、原子毎の力学量Ｊ_ｉ（ｔ）が熱流束の場合、下式（２１）となる。

熱伝導率算出シミュレータＳ３０１では、原子毎の熱流束Ｊ_ｉ（ｔ）から、下式（２２）に基づいて、原子間周波数スペクトルΠ_ｉｊ（ω；ω_ｃ ^ｍ）を算出する。ここで、ω_ｃ ^ｍは、系全体の計算で用いたカットオフ周波数の中のｍ個目（ｍは、１以上の整数）のカットオフ周波数を意味している。

次に、上式（２２）の原子間周波数スペクトルを用いて、下式（２３）に基づいて、原子間重み付き原子間周波数成分Ｓ_ｉｊ（ω_ｃ ^ｍ）を算出する。

上式（２３）におけるＣは、計算したい輸送係数と力学量を関係づける定数であり、熱伝導率の場合は、上式（１９）を用いる。ここまでの工程は、先の実施の形態１における系全体の輸送係数を算出する方法において、対応する物理量を原子間および原子間毎に算出しているのみの違いである。

原子毎の輸送係数を求めるためには、系全体の輸送係数を算出する際の自己相関関数の時間相関抽出に加え、空間相関を各原子に対して射影する必要がある。

一般的には、空間相関を分離するためには、複数のカットオフ周波数に対して、上式（２３）の原子間について空間相関を分離すればよい。すなわち、複数のカットオフ周波数に対して、上式（２３）の原子ｉと原子ｊについて構成された行列を数値的に行列対角化すればよい。しかしながら、樹脂材料のような複雑な構造を有する材料では、輸送係数の計算に必要な原子数が大きいため、行列が大規模となる。このため、行列の数値対角化計算に膨大な計算量が必要となることが予想される。

実際に、本願発明者らは、全てのカットオフ周波数に対する空間相関を各原子に射影するための数値対角化計算が、計算量の観点から困難であることを見出した。このため、本実施の形態５では、膨大な計算量を必要とせず、同時に、定量的に原子毎の輸送係数を求める手段を説明する。

上式（２３）の原子間周波数スペクトルにおいて、系全体の計算で用いたカットオフ周波数の中から（Ｋ個の中から）、カットオフ周波数ω_ｃ ^ｋにおける原子間相関行列Ａ_ｉｊ（ω_ｃ ^ｋ）を、下式（２４）に基づいて算出する（ただし、ｋ＝１〜Ｋで表される整数）。

カットオフ周波数ω_ｃ ^ｋは、系全体の計算で用いた複数のカットオフ周波数の中の１つのカットオフ周波数を選択することが望ましく、好ましくは、系全体の計算で用いた複数のカットオフ周波数の平均値に最も近い１つを選択することが望ましい。カットオフ周波数ω_ｃ ^ｋが系全体の計算で用いた複数のカットオフ周波数と著しく異なる場合には、以下に説明する周波数相関・原子相関射影係数の近似精度が悪化する。

上式（２４）の原子間相関行列の数値対角化計算によって得られる固有ベクトルｕ_１〜ｕ_Ｎから構成される下式（２５）の射影演算子Ｕ^ｋを算出することができる。

行列の数値対角化法は、一般的な教科書に記述されているので説明を省略する。
上式（２５）の射影演算子Ｕ^ｋを用いて、複数のカットオフ周波数に対して周波数相関・原子相関射影係数を、下式（２６）に基づいて算出する。

上式（２６）のβ_ｉ（ω_ｃ ^ｍ）は、周波数相関において空間相関が混ざりあった状態から、空間相関に係わる射影演算子Ｕ^ｋを通じて、ある特定の原子ｉ成分の割合を抽出することができることを意味している。

上式（１８）の系全体の熱伝導率λのカットオフ周波数依存性、すなわち、重み付き周波数成分Ｓ（ω_ｃ ^ｍ）から、上式（２６）の周波数相関・原子相関射影係数β_ｉ（ω_ｃ ^ｍ）を用いて、下式（２７）に基づき、カットオフ周波数ゼロの外挿値を求めることによって、原子毎の熱伝導率λ_ｉを求めることができる。

以上のように、実施の形態５によれば、空間相関の射影演算子の計算を１回だけで処理し、さらに系全体の熱伝導率の計算結果を利用してカットオフ周波数相関ゼロの外挿値を求めている。この結果、定量的に、かつ、効率よく原子毎の熱伝導率を算出することが可能となる。

実施の形態６．
先の実施の形態５では、数式に基づいて、本発明の原子毎の輸送係数の算出方法、原子毎の輸送係数の算出装置、および輸送係数の算出プログラムを説明した。これに対して、本実施の形態６では、図面に基づいて、本発明の原子毎の輸送係数の算出方法、原子毎の輸送係数の算出装置、および原子毎の輸送係数の算出プログラムについて、詳細に説明する。

図１０は、本発明の実施の形態６における輸送係数の算出方法のフローチャートであり、原子毎の輸送係数として原子毎の熱伝導率λ_ｉを算出する場合を示している。この図１０のフローチャートは、平衡分子動力学法シミュレータによる処理ステップＳ１０１と、原子毎の熱伝導率計算シミュレータの処理ステップＳ３０１とで構成されている。

さらに、処理ステップＳ１０１は、Ｓ１１〜Ｓ１５に細分化され、処理ステップＳ３０１は、Ｓ３１〜Ｓ３８に細分化されている。処理ステップＳ１０１における平衡分子動力学法シミュレータは、原子・分子の構造データと温度などの計算条件を入力し（ステップＳ１１、Ｓ１２）、平衡分子動力学計算を実施する（ステップＳ１３）。

そして、平衡分子動力学計算によって算出される熱流束Ｊ（ｔ）の時系列データを熱伝導算出シミュレータの処理ステップＳ２０１に対して出力する（ステップＳ１５）。この時系列データを入力することで、処理ステップＳ３０１により原子毎の熱伝導率が算出される。

このように、処理ステップＳ１０１における平衡分子動力学法シミュレータは、原子の初期座標、初期速度、原子間ポテンシャル関数、原子の質量、原子の電荷等を入力し、平衡分子動力学法に基づいて、原子・分子の座標や速度などの時間発展（時系列データ）を出力するものである。本実施の形態６の図１０におけるＳ１０１は、先の実施の形態２の図１におけるＳ１０１と比較すると、時系列データを出力するステップが、ステップＳ１４からステップＳ１５に代わっている点が異なっている。

次に、処理ステップＳ３０１における原子毎の熱伝導算出シミュレータ内の各ステップＳ３１〜Ｓ３８について説明する。ステップ３１は、原子毎の熱流束Ｊ_ｉ（ｔ）の自己相関関数のサンプル平均値＜Ｊ_ｉ（ｔ）・Ｊ_ｊ（０）＞を算出するステップである。また、ステップＳ２２は、低周波数透過フィルターのカットオフ周波数ω_ｃを入力するステップである。そして、ステップＳ３３は、本発明に係わる上式（２２）を用いて低周波数透過フィルターを適用された自己相関関数の原子間周波数スペクトルΠ_ｉｊ（ω；ω_ｃ ^ｍ）を算出するステップである。

次に、ステップＳ３４は、本発明に係わる上式（２３）を用いて低周波数透過フィルターのカットオフ周波数ω_ｃ ^ｍに対してＳ_ｉｊ（ω_ｃ ^ｍ）を算出するステップである。次に、ステップＳ３５は、本発明に係わる上式（２４）を用いて、Ａ_ｉｊ（ω_ｃ ^ｋ）を算出するステップである。さらに、ステップＳ３５は、本発明に係わるＡ_ｉｊ（ω_ｃ ^ｋ）を数値的に行列対角化することより、射影演算子Ｕ^ｋを算出するステップでもある。次に、ステップＳ３６は、本発明に係わる上式（２６）を用いて、複数のカットオフ周波数に対して、β_ｉ（ω_ｃ ^ｍ）を算出するステップである。

次に、ステップＳ３７は、本発明に係わる上式（２７）を用いて複数のカットオフ周波数に対して計算されたβ_ｉ（ω_ｃ ^ｍ）と、本発明に係わる上式（１８）を用いて複数のカットオフ周波数に対して計算された重み付き周波数成分Ｓ（ω_ｃ ^ｍ）から、低周波数透過フィルターのカットオフ周波数ω_ｃ ^ｍに対するゼロ極限値を外挿するステップである。さらに、ステップＳ３８は、求められた原子毎の熱伝導率λ_ｉを出力するステップである。

図１１は、本発明の実施の形態６における輸送係数の算出装置の構成図であり、図１０のフローチャートに対応した処理を行う装置構成を示している。この図１１における輸送係数の算出装置は、入力部１０、処理部４０、および出力部５０で構成されている。さらに、処理部４０は、７つの算出部４１〜４７（第１算出部〜第７算出部に相当）で構成されている。

入力部１０は、原子あるいは分子の少なくとも座標データおよび計算条件等の入力部であり、図１０のステップＳ１１、Ｓ１２に相当する。処理部４０内の第１算出部４１は、原子毎の力学量Ｊ_ｉ（ｔ）を算出するための平衡分子動力学法に基づく熱流束Ｊ_ｉ（ｔ）を求める算出部であり、図１０のステップＳ１３、Ｓ１５に相当する。本実施の形態６の図１１における入力部１０は、先の実施の形態２における図２の入力部１０と比較すると、時系列データを出力するステップが、ステップＳ１４からステップＳ１５に代わっている点が異なっている。

第２算出部４２は、原子毎の熱流束Ｊ_ｉ（ｔ）の原子毎の自己相関関数のサンプル平均値＜Ｊ_ｉ（ｔ）・Ｊ_ｉ（０）＞の算出結果と、低周波数透過フィルターのカットオフ周波数ω_ｃ ^ｍとを用いて、上式（２２）に基づく周波数スペクトルΠ_ｉｊ（ω；ω_ｃ ^ｍ）を求める算出部であり、図１０のステップＳ３１、Ｓ３３、およびＳ２２に相当する。

第３算出部４３は、周波数スペクトルΠ_ｉｊ（ω；ω_ｃ ^ｍ）を用いて上式（２３）に基づく原子間重み付き周波数成分Ｓ_ｉｊ（ω_ｃ ^ｍ）を求める算出部であり、図１０のステップＳ３４に相当する。第４算出部４４は、重み付き周波数成分Ｓ_ｉｊ（ω_ｃ ^ｍ）を用いて上式（２４）に基づいて原子間相関行列Ａ_ｉｊ（ω_ｃ ^ｋ）を求める算出部であり、図１０のステップＳ３５に相当する。

第５算出部４５は、原子間相関行列Ａ_ｉｊ（ω_ｃ ^ｋ）の数値対角化を行うことによって上式（２５）の射影演算子Ｕ^ｋを求める算出部であり、図１０のステップＳ３５に対応する。第５算出部４５は、上式（２６）に基づく波数相関・原子相関射影係数β_ｉ（ω_ｃ ^ｍ）を求める算出部であり、図１０のステップＳ３６に相当する。第６算出部４６は、上式（２７）に基づく原子毎の熱伝導率λを求める算出部であり、図１０のステップＳ３７に相当する。さらに、出力部５０は、熱伝導率λの出力部であり、図１０のステップＳ３８に相当する。

なお、本発明に係る熱伝導算出装置は、図１１に示した詳細な構成に限定されるものではない。例えば、平衡分子動力学法に基づく熱流束Ｊ_ｉ（ｔ）の第１算出部４１から得られる熱流束の時系列データを、外部記憶媒体に出力し、出力された熱流束Ｊ_ｉ（ｔ）の時系列データを改めて入力する入力処理部を備えることでも実現できる。

この場合には、入力した時系列データを用いて、処理部４０内の第２算出部４２〜第７算出部４７で演算処理することで、出力部５０から熱伝導率λを出力できる。

以上のように、実施の形態６によれば、空間相関の射影演算子の計算を１回だけで処理し、さらに系全体の熱伝導率の計算結果を利用してカットオフ周波数相関ゼロの外挿値を求めている。この結果、定量的に、かつ、効率よく原子毎の熱伝導率を算出することが可能となる。

実施の形態７．
上述した実施の形態６は、原子毎の熱伝導率を算出するために時系列データとして、原子毎の熱流束Ｊ_ｉ（ｔ）の時系列データを想定している。このため、原子毎の熱流束Ｊ_ｉ（ｔ）の時系列データから原子毎の熱伝導率を算出する処理ステップＳ３０１を、便宜上、原子毎の熱伝導シミュレータと称していた。しかしながら、本発明に係る算出方法を含む処理ステップＳ３０１は、原子毎の熱伝導率の算出方法だけに限定されるものではない。

この処理ステップＳ３０１は、例えば、原子毎の体積粘性率、原子毎の電気伝導率、原子毎のずれ粘性係数などの、原子毎の輸送係数の定量的算出にも、一般的に適用できるものである。この場合、分子動力学法シミュレータから出力される時系列データは、算出する輸送係数と対応した力学量となる。算出する輸送係数と力学量との対応関係は、一般的な教科書に記述されている内容であるが、例えば、電気伝導率の場合、原子毎の時系列でデータＪ_ｉ（ｔ）は、下式（２８）となる。

図１２は、本発明の実施の形態６における輸送係数の算出方法のフローチャートである。また、図１３は、本発明の実施の形態６における輸送係数の算出装置の構成図であり、図１２のフローチャートに対応した処理を行う装置構成を示している。

本実施の形態６の図１２における輸送係数を算出するフローチャート、および図１３における算出装置の構成と、先の図１０における熱伝導率を算出するフローチャート、および図１１における算出装置の構成との差異は、図１２、図１３において添字ａが付された部分である。より具体的には、本実施の形態５における図１２、図１３では、力学量に応じた輸送係数Ｋ_ｉの定数Ｃを扱っている点が異なっている。

定数Ｃは、平衡分子動力学シミュレータにおける温度や体積などの原子・分子データや、計算条件からの数値と自然定数で構成される。輸送係数の算出装置も同様に、求めたい輸送係数Ｋに対応する力学量を計算する力学量算出処理部が有ればよい。

以上のように、実施の形態７によれば、空間相関の射影演算子の計算を１回だけで処理し、さらに系全体の熱伝導率の計算結果を利用してカットオフ周波数相関ゼロの外挿値を求めている。この結果、定量的に、かつ、効率よく原子毎の熱伝導率を算出することが可能となる。

なお、本発明は、平衡分子動力学シミュレータの処理ステップＳ１０１の詳細に依存しない。すなわち、本発明は、入力する原子・分子データ（計算材料）に限定されるものではない。さらに、本発明は、平衡分子動力学法シミュレータに入力する計算条件にも限定されるものではない。

さらに、本発明は、上述した実施の形態に限定されるものではなく、本発明の要旨を逸脱しない範囲で設計変更などがあっても本発明に含まれる。例えば、熱伝導率を計算する工程を、コンピュータにより実現するためのプログラム自体、および、このプログラムをコンピュータで読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータに読み込ませることにより、各機能を実現してもよい。

１０ 入力部、２０ 処理部、２１ 第１算出部、２２ 第２算出部、２３ 第３算出部、２４ 第４算出部、３０ 出力部、４０ 処理部、４１ 第１算出部、４２ 第２算出部、４３ 第３算出部、４４ 第４算出部、４５ 第５算出部、４６ 第６算出部、４７ 第７算出部、５０ 処理部。

Claims (10)

1. 原子あるいは分子の構造データおよび温度計算条件を入力するステップと、
   所定の時間ｔに対して力学量Ｊ（ｔ）を算出するステップと、
   前記力学量Ｊ（ｔ）を用いて下式（１）
   

   

   に基づいて周波数スペクトルΠ（ω、ω_ｃ）を算出するステップと、
   前記周波数スペクトルΠ（ω、ω_ｃ）を用いて下式（２）
   

   

   に基づいて重み付き周波数成分Ｓ（ω_ｃ）を算出するステップと、
   前記重み付き周波数成分Ｓ（ω_ｃ）を用いて下式（３）
   

   

   に基づいて輸送係数Ｋを算出するステップと
   を備えたことを特徴とする輸送係数の算出方法。
2. 請求項１に記載の輸送係数の算出方法において、
   前記力学量Ｊ（ｔ）として、熱流束を示す下式（４）
   

   

   を用いることを特徴とする輸送係数の算出方法。
3. 請求項１に記載の輸送係数の算出方法において、
   前記力学量Ｊ（ｔ）として、電流密度を示す下式（５）
   

   

   を用いることを特徴とする輸送係数の算出方法。
4. 原子あるいは分子の構造データおよび温度計算条件の入力部と、
   所定の時間ｔに対する力学量Ｊ（ｔ）を算出する第１算出部と、
   前記力学量Ｊ（ｔ）を用いて下式（１）
   

   

   に基づいて周波数スペクトルΠ（ω、ω_ｃ）を算出する第２算出部と、
   前記周波数スペクトルΠ（ω、ω_ｃ）を用いて下式（２）
   

   

   に基づいて重み付き周波数成分Ｓ（ω_ｃ）を算出する第３算出部と、
   前記重み付き周波数成分Ｓ（ω_ｃ）を用いて下式（３）
   

   

   に基づいて輸送係数Ｋを算出する第４算出部と
   を備えたことを特徴とする輸送係数の算出装置。
5. 輸送係数の算出装置に用いられる輸送係数の算出プログラムであって、
   コンピュータに、
   原子あるいは分子の構造データおよび温度計算条件を入力する入力部と、
   所定の時間ｔに対する力学量Ｊ（ｔ）を算出する第１算出部と、
   前記力学量Ｊ（ｔ）を用いて下式（１）
   

   

   に基づいて周波数スペクトルΠ（ω、ω_ｃ）を算出する第２算出部と、
   前記周波数スペクトルΠ（ω、ω_ｃ）を用いて下式（２）
   

   

   に基づいて重み付き周波数成分Ｓ（ω_ｃ）を算出する第３算出部と、
   前記重み付き周波数成分Ｓ（ω_ｃ）を用いて下式（３）
   

   

   に基づいて輸送係数Ｋを算出する第４算出部と
   を機能させるための輸送係数の算出プログラム。
6. 原子あるいは分子の構造データおよび温度計算条件を入力するステップと、
   所定の時間ｔに対して原子毎の力学量Ｊ_ｉ（ｔ）を算出するステップと、
   前記力学量Ｊ_ｉ（ｔ）を用いて下式（６）
   

   

   に基づいて原子間周波数スペクトルΠ_ｉｊ（ω；ω_ｃ ^ｍ）を算出するステップと、
   前記原子間周波数スペクトルΠ_ｉｊ（ω；ω_ｃ ^ｍ）を用いて下式（７）
   

   

   に基づいて原子間重み付き原子間周波数成分Ｓ_ｉｊ（ω_ｃ ^ｍ）を算出するステップと、
   前記原子間重み付き原子間周波数成分Ｓ_ｉｊ（ω_ｃ ^ｍ）を用いて下式（８）
   

   

   に基づいて重み付き原子間相関行列算出部Ａ_ｉｊ（ω_ｃ ^ｋ）を算出するステップと、
   前記重み付き原子間相関算出部Ａ_ｉｊ （ω_ｃ ^ｋ）を対角化することによって算出される下式（９）
   

   

   の射影演算子Ｕ^ｋを算出する計算ステップと、
   前記射影演算子Ｕ^ｋと前記原子間重み付き原子間周波数成分Ｓ_ｉｊ（ω_ｃ ^ｍ）を用いて下式（１０）
   

   

   に基づいて周波数相関・原子間相関射影係数β_ｉ（ω_ｃ ^ｍ）を算出する計算ステップと、
   前記周波数相関・原子相関射影係数β_ｉ（ω_ｃ ^ｍ）、および複数のカットオフ周波数に対して計算された重み付き周波数成分Ｓ（ω_ｃ ^ｍ）を用いて下式（１１）
   

   

   に基づいて原子毎の輸送係数Ｋ_ｉを算出する計算ステップと
   を備えたことを特徴とする原子毎の輸送係数の算出方法。
7. 請求項６に記載の輸送係数の算出方法において、
   前記原子毎の力学量Ｊ_ｉ（ｔ）として、原子毎の熱流束を示す下式（１２）
   

   

   を用いることを特徴とする原子毎の輸送係数の算出方法。
8. 請求項６に記載の輸送係数の算出方法において、
   前記原子毎の力学量Ｊ_ｉ（ｔ）として、原子毎の電流密度を示す下式（１３）
   

   

   を用いること特徴とする原子毎の輸送係数の算出方法。
9. 原子あるいは分子の少なくとも座標データおよび所定の温度の入力部と、
   所定の時間ｔに対する原子毎の力学量Ｊ_ｉ（ｔ）の算出部と、
   前記力学量Ｊ_ｉ（ｔ）を用いて下式（６）
   

   

   に基づいて原子間周波数スペクトルルΠ_ｉｊ（ω；ω_ｃ ^ｍ）を算出する第２算出部と、
   前記原子間周波数スペクトルルΠ_ｉｊ（ω；ω_ｃ ^ｍ）を用いて下式（７）
   

   

   に基づいて原子間重み付き周波数成分Ｓ_ｉｊ（ω_ｃ ^ｍ）を算出する第３算出部と、
   前記原子間重み付き周波数成分Ｓ_ｉｊ（ω_ｃ ^ｍ）を用いて下式（８）
   

   

   に基づいて重み付き原子間相関行列Ａ_ｉｊ（ω_ｃ ^ｋ）を算出する算４出部と、
   前記重み付き原子間相関行列Ａ_ｉｊ（ω_ｃ ^ｋ）の行列対角化によって得られる固有ベクトルから構成される下式（９）
   

   

   の射影演算子Ｕ^ｋを算出する第５算出部と
   前記射影演算子Ｕ^ｋと前記重み付き周波数成分Ｓ（ω_ｃ ^ｍ）を用いて下式（１０）
   

   

   に基づいて周波数相関・原子間相関射影係数β_ｉ（ω_ｃ ^ｍ）を算出する第６算出部と、
   前記周波数相関・原子間相関射影係数β_ｉ（ω_ｃ ^ｍ）を用いて下式（１１）
   

   

   に基づいて原子毎の輸送係数Ｋ_ｉを算出する第７算出部と
   を備えた原子毎の輸送係数の算出装置。
10. 原子毎の輸送係数の算出装置に用いられる原子毎の輸送係数の算出プログラムであって、
    コンピュータに、
    原子あるいは分子の少なくとも座標データおよび温度計算条件を入力する入力部と、
    所定の時間ｔに対する原子毎の力学量Ｊ_ｉ（ｔ）を算出する第１算出部と、
    前記原子毎の力学量Ｊ_ｉ（ｔ）を用いて下式（６）
    

    

    に基づいて原子間周波数スペクトルΠ_ｉｊ（ω、ω_ｃ ^ｍ）を算出する第２算出部と、
    前記原子間周波数スペクトルΠ_ｉｊ（ω、ω_ｃ）を用いて下式（７）
    

    

    に基づいて原子間重み付き周波数成分Ｓ_ｉｊ（ω_ｃ）を算出する第３算出部と、
    前記原子間重み付き周波数成分Ｓ_ｉｊ（ω_ｃ）を用いて下式（８）
    

    

    に基づいて重み付き原子間相関行列Ａｉｊ（ω_ｃ ^ｋ）を算出する第４算出部と、
    前記重み付き原子間相関行列Ａ_ｉｊ（ω_ｃ ^ｋ）の行列対角化によって得られる固有ベクトルから構成される下式（９）
    

    

    の射影演算子Ｕ^ｋを算出する第５算出部と、
    前記原子間重み付き周波数成分Ｓ_ｉｊ（ω_ｃ ^ｍ）と前記射影演算子Ｕ^ｋを用いて下式（１０）
    

    

    に基づいて周波数相関・原子相関射影係数β_ｉ（ω_ｃ ^ｍ）を算出する第６算出部と、
    複数のカットオフ周波数に対して計算された重み付き周波数成分Ｓ（ω_ｃ ^ｍ）と前記周波数相関・原子間相関射影係数β_ｉ（ω_ｃ ^ｍ）を用いて下式（１１）
    

    

    に基づいて原子毎の輸送係数Ｋ_ｉを算出する第７算出部と
    を機能させるための原子毎の輸送係数の算出プログラム。

Images