分子シミュレーション装置
【課題】本発明は、アモルファスの電気伝導度を推定する分子シミュレーション装置を提供することを課題とする。
【解決手段】アモルファスの電気伝導度を推定する分子シミュレーション装置であって、アモルファスに存在する異なる波動関数のキャリアの乗り移り易さ（異なる波動関数のエネルギレベルの差及び異なる波動関数の重なり領域の体積に基づく評価）、空間上の広さ（異なる波動関数の各体積の和から重なり領域の体積を引いた体積に基づく評価）及びキャリアが存在する確率をそれぞれ評価し、当該３つの評価に基づいて電気伝導度を推定することを特徴とする。
【解決手段】アモルファスの電気伝導度を推定する分子シミュレーション装置であって、アモルファスに存在する異なる波動関数のキャリアの乗り移り易さ（異なる波動関数のエネルギレベルの差及び異なる波動関数の重なり領域の体積に基づく評価）、空間上の広さ（異なる波動関数の各体積の和から重なり領域の体積を引いた体積に基づく評価）及びキャリアが存在する確率をそれぞれ評価し、当該３つの評価に基づいて電気伝導度を推定することを特徴とする。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、アモルファスの電気伝導度を推定する分子シミュレーション装置に関する。
【背景技術】
【０００２】
従来、半導体などの物質の電気特性を推定する様々な手法が提案されている。特許文献１には、量子力学的な効果を近似的に考慮した半導体の電気特性を推定する手法が開示されている。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開平１１−３９９８号公報
【特許文献２】特開２００７−１３４５７４号公報
【特許文献３】特開平５−１５１１９１号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
上記の提案されている手法は、特定の結晶状態を有する物質に適用した手法である。しかし、アモルファスは、特定の結晶構造を有していないので、このような従来の手法を適用して電気特性を推定することができない。
【０００５】
そこで、本発明は、アモルファスの電気伝導度を推定する分子シミュレーション装置を提供することを課題とする。
【課題を解決するための手段】
【０００６】
本発明に係る分子シミュレーション装置は、アモルファスの電気伝導度を推定する分子シミュレーション装置であって、アモルファスに存在する異なる波動関数のキャリアの乗り移り易さ、空間上の広さ及びキャリアが存在する確率をそれぞれ評価し、当該３つの評価に基づいて電気伝導度を推定することを特徴とする。
【０００７】
アモルファスは、周期的な結晶構造を有さないので、波動関数が局所的に存在している。そのため、アモルファスでキャリアが移動して電気伝導が起こるには、アモルファスに存在する異なる波動関数間をキャリアが乗り移って電気伝導が行われると考えられる。さらに、その異なる波動関数においてそれぞれ異なる空間を広く占めていると、キャリアが他の波動関数に乗り移れる可能性が高くなる。また、波動関数にキャリアが多く存在していると、キャリアが他の波動関数に乗り移る可能性が高くなる。そこで、この分子シミュレーション装置では、アモルファスに存在する異なる波動関数についてのキャリアの乗り移り易さ、空間上の広さ及びキャリアが存在する確率をそれぞれ評価する。そして、分子シミュレーション装置では、その３つの評価に基づいてアモルファスの電気伝導度を推定する。このように、分子シミュレーション装置では、異なる波動関数についてのキャリアの乗り移り易さ、空間上の広さ及びキャリアが存在する確率をそれぞれ評価することにより、アモルファスの電気伝導度を精度良く推定することができる。
【０００８】
本発明の上記分子シミュレーション装置では、異なる波動関数のキャリアの乗り移り易さを、異なる波動関数のエネルギレベルの差及び異なる波動関数の重なり領域の体積に基づいて評価すると好適である。
【０００９】
異なる波動関数間のエネルギレベルの差が小さいと、キャリアがその波動関数間を乗り移り易い。また、異なる波動関数が重なっている領域が大ききと、キャリアはその波動関数間を乗り移り易い。そこで、この分子シミュレーション装置では、異なる波動関数ついてのエネルギレベルの差及び重なり領域の体積を評価することにより、キャリアの乗り移り易さを高精度に評価することができる。
【００１０】
本発明の上記分子シミュレーション装置では、異なる波動関数の空間上の広さを、異なる波動関数の各体積の和から重なり領域の体積を引いた体積に基づいて評価すると好適である。
【００１１】
異なる波動関数がそれぞれ異なる空間（重なり領域以外の空間）を広く占めていると、キャリアがその異なる波動関数間を乗り移った後に他の波動関数に乗り移れる可能性が高くなる。そこで、この分子シミュレーション装置では、異なる波動関数の各体積の和から重なり領域の体積を引いた体積で評価することにより、異なる波動関数における空間上の広さを高精度に評価することができる。
【発明の効果】
【００１２】
本発明は、異なる波動関数についてのキャリアの乗り移り易さ、空間上の広さ及びキャリアが存在する確率をそれぞれ評価することにより、アモルファスの電気伝導度を精度良く推定することができる。
【図面の簡単な説明】
【００１３】
【図１】本実施の形態に係る分子シミュレーション装置の構成図である。
【図２】アモルファス体の構造の一例を示す図である。
【図３】実空間上の異なる２つの波動関数の固有エネルギの一例を示す図である。
【図４】実空間上の異なる２つの波動関数の重なり状態の一例を示す図である。
【図５】アモルファス体の空間分割の一例を示す図である。
【図６】波動関数のψとφの関係の一例を示すグラフである。
【図７】電気伝導度推定処理（基本版）の流れを示すフローチャートである。
【図８】電気伝導度推定処理（改良版）の流れを示すフローチャートである。
【図９】体積算出処理の流れを示すフローチャートである。
【図１０】体積算出処理の流れを示すフローチャートである。
【発明を実施するための形態】
【００１４】
以下、図面を参照して、本発明に係る分子シミュレーション装置の実施の形態を説明する。
【００１５】
本実施の形態では、本発明に係る分子シミュレーション装置を、パーソナルコンピュータ（以下、パソコンと記載）上に構成される分子シミュレーション装置に適用する。本実施の形態に係る分子シミュレーション装置では、アモルファスの電気伝導度（代用評価指標）を推定する。
【００１６】
図１〜図６を参照して、本実施の形態に係る分子シミュレーション装置１について説明する。図１は、本実施の形態に係る分子シミュレーション装置の構成図である。図２は、アモルファス体の構造の一例を示す図である。図３は、実空間上の異なる２つの波動関数の固有エネルギの一例を示す図である。図４は、実空間上の異なる２つの波動関数の重なり状態の一例を示す図である。図５は、アモルファス体の空間分割の一例を示す図である。図６は、波動関数のψとφの関係の一例を示すグラフである。
【００１７】
分子シミュレーション装置１は、記憶装置（ハードディスクなど）に格納されているアモルファス体の電気伝導度推定処理用のプログラムをパソコン本体１０で実行することによって構成され、その推定結果をディスプレイ１１に表示したりあるいは記憶装置に格納する。この際、ユーザは、このアモルファス体の電気伝導度推定に必要な情報をキーボード１２などを利用して入力し、推定結果をディスプレイ１１などで確認する。
【００１８】
分子シミュレーション装置１におけるアモルファス体の電気伝導度推定処理について具体的に説明する前に、アモルファス体の電気伝導度の推定方法及びその推定において必要な波動関数の体積の算出方法について説明する。
【００１９】
アモルファスに限らず、物質中の電子は、物質中に存在する多数の波動関数のいずれかに位置しており、その波動関数の広がりの範囲を自由に移動している。
【００２０】
シリコンなどの半導体は、周期的な結晶構造を有するため、波動関数も周期構造を持っている。このような半導体中の電子は、同じ波動関数に位置したままで自由に移動できるが、実際には他の電子や熱振動、不純物、格子欠陥との衝突によって移動が阻害される。それが、電気伝導度を有限の値としている主な要因である。したがって、周期的な結晶構造の半導体の電気伝導度をシミュレーションによって予測するには、他の電子や熱振動、不純物、格子欠陥との衝突の頻度を計算すればよい。
【００２１】
しかし、アモルファス体（アモルファスの半導体）は、周期的な結晶構造を有さないため、波動関数が局所的に存在している。図２には、アモルファス体ＡＭの構造の一例を示しており、多数の原子ＡＴ，・・・が不規則に存在しており、多数の波動関数ＷＦ，・・・が局所的に分布している。このアモルファス体の電子が移動して電気伝導が起こるには、電子が異なる波動関数間を乗り移ることが必要と考えられる。したがって、アモルファス体の場合、異なる波動関数における電子（キャリア）の乗り移り易さが電気伝導度を決定していると考えられる。
【００２２】
キャリア（電子、ホール）が固有エネルギの異なる波動関数間を乗り移るには、そのエネルギ差を熱振動によるフォノンのやりとりなどで補償しなければならない。図３の例で示すように、固有エネルギＥｍの波動関数と固有エネルギＥｎの波動関数とのエネルギ差が熱エネルギ程度（２６ｍｅＶ程度）のエネルギ内ならやりとりは起こり易いが、それ以上のエネルギ差の場合には変化させることは難しい。したがって、異なる波動関数間のエネルギ差が小さいほど、キャリアが乗り移り易い（評価指標１）。
【００２３】
キャリアが異なる波動関数間を乗り移る場合、その波動関数同士の空間上の重なりが大きいほど乗り移り易い（評価指標２）。図４の例で示すように、波動関数ＡＦｍと波動関数ＡＦｎとは空間上で重なり領域ＯＬを有しており、この重なり領域ＯＬを利用して波動関数ＡＦｍと波動関数ＡＦｎとの間をキャリアが移動する。
【００２４】
また、異なる波動関数間を乗り移り易さだけでなく、それによって伝導が起こり易いことが必要である。それには、異なる波動関数間が乗り移り易かったとしても、その異なる
波動関数が同じ空間（重なり領域が大部分）を占めていては意味がなく、キャリアが波動関数間を乗り移った後にそれ以外の波動関数に移動するためには、異なる波動関数がそれぞれ異なる空間（重なり領域以外の空間）を広く占めていることも重要となる（評価指標３）。
【００２５】
また、キャリアが異なる波動関数間を乗り移る場合、個々の波動関数にキャリアが存在する確率が高いことも重要となる（評価指標４）。波動関数に存在するキャリアが少ないと、キャリアが他の波動関数に乗り移る可能性も低くなる。
【００２６】
この４つの評価指標が、アモルファス体における電気伝導に寄与すると考えられる。したがって、この４つの評価指標に基づいた評価をそれぞれ行い、その４つの評価結果を加味して電気伝導度を推定する。
【００２７】
実際の電気伝導では、キャリアが多数の波動関数をランダムに乗り移って伝導すると考えられる。その組み合わせの数は膨大な数になるので、全ての組み合わせについて評価することは実用的でない。そこで、少数の波動関数の組み合わせで評価し、その全ての組み合わせの評価値の和によって電気伝導のし易さを表わすことにより、実際の電気伝導度を高精度に推定（評価）できる。
【００２８】
ここでは、異なる２つの波動関数の組み合わせで評価する場合について説明する。上記した４つの評価指標による評価を加味して電気伝導度（電気伝導率）の代用評価値Ｆを算出する式（１）を示す。式（１）において、ｎ，ｍは、２つの異なる波動関数を示している。
【数１】
【００２９】
Ｖｎ，Ｖｍが各波動関数の体積であり、Ｖｎｍが２つの波動関数の重なり領域の体積であり、各体積Ｖｎ，Ｖｍ，Ｖｎｍは下記の方法で算出される。Ｅｎ，Ｅｍは、各波動関数の固有エネルギであり、既存の分子シミュレーションによって算出される。Ｅｆは、フェルミエネルギであり、既存の分子シミュレーションによって算出される。Ｅｔｈは、定数であり、熱エネルギ（２６ｍｅＶ）程度である。ＫＢは、ボルツマン定数である。Ｔは、系の絶対温度であり、ユーザによる入力などによって予め決まった値である。
【００３０】
式（１）のΣは、ｎ，ｍで表される異なる２つの波動関数の組み合わせについて、ｎ＝ｍでない全ての組み合わせについての評価値を加算することを表わしている。
【００３１】
式（１）のｅｘｐ｛｜Ｅｎ−Ｅｍ｜／Ｅｔｈ｝は、２つの波動関数の固有エネルギの差を定数Ｅｔｈで除算することにより評価指標１による評価を行うファクタである。式（１）の×Ｖｎｍは、２つの波動関数の重なり領域の体積を乗算することにより評価指標２による評価を行うファクタである。式（１）の（Ｖｎ＋Ｖｍ−Ｖｎｍ）は、２つの波動関数の体積の和から重なり領域の体積を減算することにより評価指標３による評価を行うファクタである。式（１）のｅｘｐ｛（Ｅｆ−（Ｅｎ＋Ｅｍ）／２）／（ＫＢ×Ｔ）｝は、量子力学論に基づいて２つの波動関数の平均エネルギにおける存在確率を乗算することにより評価指標４による評価を行うファクタである。
【００３２】
上記の基本的なアモルファス体の電気伝導度の推定方法の考え方では、アモルファス体に存在する全ての波動関数におけるキャリア（特に、電子）が電気伝導に寄与するとみなしている。しかし、実際の物質では、エネルギの低い波動関数は、全て電子が存在し、フェルミの排他律という法則に従ってそれ以上電子を取り込めない状態になっている。そのため、電子は、身動きが取れないので、伝導に寄与できない。半導体において、実際に伝導に寄与しうるのは、フェルミエネルギよりも高い固有エネルギを持つ波動関数に存在する電子か、逆に、フェルミエネルギよりも低い固有エネルギを持つ波動関数において電子が欠けた状態にあるホールと呼ばれる仮想的な正電荷の粒子である。
【００３３】
そこで、上記の基本的な方法を改良し、２つの波動関数が共にフェルミエネルギよりも高い固有エネルギを持つ場合の電子による伝導と２つの波動関数が共にフェルミエネルギよりも低い固有エネルギを持つ場合のホールによる伝導に分けて評価を行う。具体的には、状態ｎの波動関数の固有エネルギＥｎと状態ｍの波動関数の固有エネルギＥｍが共にフェルミエネルギＥｆよりも高い場合、電子がその２つの波動関数間を乗り移る（電気伝導に寄与する）ことができるので、電気伝導度に電子による伝導を加味する。状態ｎの波動関数の固有エネルギＥｎと状態ｍの波動関数の固有エネルギＥｍが共にフェルミエネルギＥｆよりも低い場合、ホールがその２つの波動関数間を乗り移る（電気伝導に寄与する）ことができるので、電気伝導度にホールによる電気伝導を加味する。状態ｎの波動関数の固有エネルギＥｎと状態ｍの波動関数の固有エネルギＥｍのうちの一方がフェルミエネルギＥｆよりも高く、他方がフェルミエネルギＥｆよりも低い場合（禁制帯を挟んで２つの波動関数が存在する場合）、キャリアがその２つの波動関数間を乗り移ることができない（電気伝導に寄与しない）ので、電気伝導度には加味しない。
【００３４】
各体積Ｖｎ，Ｖｍ，Ｖｎｍ及びＶｎ＋Ｖｍ−Ｖｎｍの算出方法については説明する。波動関数は、連続的な値を持って空間に分布している。そのため、波動関数の広さ（体積）などの概念は、厳密な定義がなく、特別に定める必要がある。そこで、以下に、２つの算出方法について説明する。
【００３５】
まず、１つの目の算出方法について説明する。この方法では、波動関数は連続的な値を持つが、その波動関数の絶対値｜ψ｜が所定の規定値｜ψｒ｜よりも大きければ「そこに波動関数が存在する」とみなし、規定値｜ψｒ｜よりも小さければ「そこに波動関数が存在しない」とみなす。具体的な方法としては、図５に示すように、アモルファス体の空間を同じ体積を有する多数の小さい断片（直方体）に分割する。既存の分子シミュレーションによって、状態ｍの波動関数の絶対値｜ψｍ｜と状態ｎの波動関数の絶対値｜ψｎ｜が求めることができ、さらに、任意の断片ｉについて状態ｍの波動関数の絶対値｜ψｉｍ｜と状態ｎの波動関数の絶対値｜ψｉｎ｜が予め求まることができる。そして、断片ｉの状態ｍの波動関数の絶対値｜ψｉｍ｜が規定値｜ψｒ｜よりも大きいか否かを判定し、｜ψｉｍ｜が｜ψｒ｜よりも大きいと判定した場合には状態ｍの波動関数の体積Ｖｍに断片ｉの体積Ｖｉを加算し、この処理を空間の全ての断片について行い、状態ｍの波動関数の体積Ｖｍを求める。状態ｎの波動関数の体積Ｖｎについても同様に求める。また、断片ｉの状態ｍの波動関数の絶対値｜ψｉｍ｜が規定値｜ψｒ｜よりも大きいかつ状態ｎの波動関数の絶対値｜ψｉｎ｜が規定値｜ψｒ｜よりも大きいか否かを判定し、｜ψｉｍ｜と｜ψｉｎ｜が共に規定値｜ψｒ｜より大きいと判定した場合には状態ｍの波動関数と状態ｎの波動関数との重なり領域の体積Ｖｎｍに断片ｉの体積Ｖｉを加算し、この処理を空間の全ての断片について行い、重なり領域の体積Ｖｎｍを求める。最後に、Ｖｎ＋Ｖｍ−Ｖｎｍを算出する。なお、規定値｜ψｒ｜は、ユーザによって指定される。
【００３６】
次に、２つの目の算出方法について説明する。この方法では、波動関数をそのまま連続的な値として取り扱うために、波動関数が式（２）で示す規格化条件を満たしているものとする。この式（２）において、Ｖは、積分範囲としてアモルファス体の空間の全体積について積分を行うということを意味する。この式（２）を満たさない場合、ψに適当な係数を乗算した値を改めてψと置き直すことにより、式（２）を満たすように調整する。
【数２】
【００３７】
そして、各体積Ｖｎ，Ｖｍ，Ｖｎｍを個別に求めるのではなく、式（３）によってＶｎ＋Ｖｍ−Ｖｎｍをまとめて求める。
【数３】
【００３８】
式（３）において、φは波動関数の値である。その２つのφの差｜φｎ−φｍ｜は、状態ｎの波動関数か状態ｍの波動関数のどちらかが存在するということを満たす指標として用いることができ、その値を全空間（全体積Ｖ）で積分することにより、（Ｖｎ＋Ｖｍ−Ｖｎｍ）の概念を示す指標を求めることができる。
【００３９】
Ｖｎ＋Ｖｍ−Ｖｎｍの意味するところは、状態ｎの波動関数と状態ｍの波動関数の占める空間において重なり領域以外で占める領域の体積である。基本的には、波動関数が有るか無いかの概念であるので、波動関数が非常に大きな値を持っていても「有る」という判定以上のものはない。したがって、図６に示すように、波動関数の値φを頭打ち（上限を１とする）に加工した関数を用いる。図６に示すようなグラフの関係が得られる適切な関数としては、例えば、式（４）がある。
【数４】
【００４０】
さらに、式（１）で電気伝導度を算出するためには、Ｖｎｍが、必要となるので、別個に求める。Ｖｎｍは、状態ｎの波動関数と状態ｍの波動関数の両方に存在する領域（重なり領域）の体積の概念である。したがって、状態ｎの波動関数の値φｎと状態ｍの波動関数の値φｍの積φｎ×φｍの平方根の積分により、その概念に近い指標を求めることができる。その具体的な式としては、式（５）で示す。
【数５】
【００４１】
それでは、分子シミュレーション装置１におけるアモルファス体の電気伝導度の推定処理について説明する。アモルファス体の電気伝導度推定処理用のプログラムをパソコン本体１０で実行することにより、電気伝導度の代用評価値Ｆを求めるためのメインルーチンが実行され、そのメインルーチンによってＶｎｍ及びＶｎ＋Ｖｍ−Ｖｎｍを算出するためのサブルーチンが呼び出されて実行される。ここでは、メインルーチンについては、上記で説明した電気伝導度の推定方法に基づく基本版とその基本版から改良した改良版について説明する。また、サブルーチンについては、上記で説明した２つの算出方法に基づくものをそれぞれ説明する。
【００４２】
まず、電気伝導度の代用評価値Ｆを求めるためのメインルーチンの基本版を図７のフローチャートに沿って説明する。図７は、電気伝導度推定処理（基本版）の流れを示すフローチャートである。ここでは、アモルファス体に存在する波動関数の個数をＮ個とする。ｎ，ｍは、Ｎ個の波動関数中の２つの状態ｎの波動関数と状態ｍの波動関数の組み合わせを更新してゆくためのカウンタである。
【００４３】
系の絶対温度Ｔは、ユーザによる入力などによって予め決められ、パソコンの記憶装置に格納されている。また、Ｎ個の波動関数の固有エネルギＥは、既存の分子シミュレーションによって予め算出され、パソコンの記憶装置に格納されている。また、フェルミエネルギＥｆは、既存の分子シミュレーションによって予め算出され、パソコンの記憶装置に格納されている。また、波動関数の個数Ｎは、既存の分子シミュレーションによって予め算出され、パソコンの記憶装置に格納されている。なお、Ｎ個をユーザ入力などによって規定し、Ｎ個までの波動関数でシミュレーションを行ってもよい。
【００４４】
メインルーチンが起動すると、ｍに１を代入して初期化するとともに、代用評価値Ｆに０を代入して初期化する（Ｓ１０）。
【００４５】
初期化後、ループ毎に、状態ｍの波動関数を次の波動関数に更新するために、ｍをインクリメントする（Ｓ１１）。また、ｎに１を代入して初期化する（Ｓ１２）。そして、今回の組み合わせの状態ｎの波動関数と状態ｍの波動関数についての体積Ｖｎｍ及びＶｎ＋Ｖｍ−Ｖｎｍを算出するためのサブルーチンを呼び出して実行し、Ｖｎｍ及びＶｎ＋Ｖｍ−Ｖｎｍの算出する（Ｓ１３）。さらに、Ｓ１３で算出されたＶｎｍ，Ｖｎ＋Ｖｍ−Ｖｎｍ及び記憶装置から取り出した固有エネルギＥｎ，Ｅｍを用いて、式（６）によって、今回の組み合わせの状態ｎの波動関数と状態ｍの波動関数でのΔＦを算出する（Ｓ１４）。そして、代用評価値Ｆの前回値に今回算出したΔＦを加算し、その加算値を代用評価値Ｆに代入する（Ｓ１５）。そして、状態ｎの波動関数を次の波動関数に更新するために、ｎをインクリメントする（Ｓ１６）。さらに、ｎがｍよりも小さいか否かを判定する（Ｓ１７）。Ｓ１７にてｎがｍより小さいと判定した場合、Ｓ１３の処理に戻り、状態ｎの波動関数を変えた次の組み合わせの処理を行う。一方、Ｓ１７にてｎがｍに等しくなったと判定した場合、ｍがＮより小さいか否かを判定する（Ｓ１８）。Ｓ１８にてｍがＮより小さいと判定した場合、Ｓ１１の処理に戻り、状態ｍの波動関数を次の波動関数に変えた組み合わせの処理を行う。
【数６】
【００４６】
Ｓ１８にてｍがＮと等しくなったと判定した場合、代用評価値Ｆの算出が終了したので、その代用評価値Ｆをアモルファス体の電気伝導度としてディスプレイ１１に表示するかあるいはパソコンの記憶装置に格納する。
【００４７】
次に、電気伝導度の代用評価値Ｆを求めるためのメインルーチンの改良版を図８のフローチャートに沿って説明する。図８は、電気伝導度推定処理（改良版）の流れを示すフローチャートである。
【００４８】
この改良版では、上記の基本版と比較すると、基本版におけるΔＦの算出（Ｓ１４）だけが異なっており、他の処理は同じである。そこで、このΔＦの算出についてのみ説明する。
【００４９】
Ｓ１３でＶｎｍ及びＶｎ＋Ｖｍ−Ｖｎｍを算出すると、状態ｎの波動関数の固有エネルギＥｎがフェルミエネルギＥｆよりも高いかつ状態ｍの波動関数の固有エネルギＥｍがフェルミエネルギＥｆよりも高いか否かを判定する（Ｓ１４ａ）。Ｓ１４ａにてＥｎがＥｆよりも高いかつＥｍがＥｆよりも高いと判定した場合（電子が２つの波動関数間を乗り移る場合）、Ｓ１３で算出されたＶｎｍ，Ｖｎ＋Ｖｍ−Ｖｎｍ及び記憶装置から取り出した固有エネルギＥｎ，Ｅｍを用いて、式（６）によって、今回の組み合わせの状態ｎの波動関数と状態ｍの波動関数でのΔＦを算出する（Ｓ１４ｂ）。
【００５０】
一方、Ｓ１４ａにてＥｎとＥｍの少なくとも一方がＥｆよりも低いと判定した場合、状態ｎの波動関数の固有エネルギＥｎがフェルミエネルギＥｆよりも低いかつ状態ｍの波動関数の固有エネルギＥｍがフェルミエネルギＥｆよりも低いか否かを判定する（Ｓ１４ｃ）。Ｓ１４ｃにてＥｎがＥｆよりも低いかつＥｍがＥｆよりも低いと判定した場合（ホールが２つの波動関数間を乗り移る場合）、Ｓ１３で算出されたＶｎｍ，Ｖｎ＋Ｖｍ−Ｖｎｍ及び記憶装置から取り出した固有エネルギＥｎ，Ｅｍを用いて、式（７）によって、今回の組み合わせの状態ｎの波動関数と状態ｍの波動関数でのΔＦを算出する（Ｓ１４ｄ）。なお、式（７）では、ホールの存在確率とエネルギの関係は電子と逆になるので、（Ｅｍ＋Ｅｎ）／２からフェルミエネルギＥｆを減算している（式（６）と逆）。
【数７】
【００５１】
Ｓ１４ａ、Ｓ１４ｃにてＥｎとＥｍのうちの一方がＥｆよりも高く、他方がＥｆよりも低いと判定した場合（キャリアが２つの波動関数間を乗り移ることができない場合）、ΔＦに０を代入する（Ｓ１４ｅ）。
【００５２】
そして、代用評価値Ｆの前回値に今回算出したΔＦを加算し、その加算値を代用評価値Ｆに代入する（Ｓ１５）。ただし、Ｓ１４ｅでΔＦに０を代入した場合、代用評価値Ｆは前回値のままである。
【００５３】
次に、上記の１つの目の算出方法によるＶｎｍ及びＶｎ＋Ｖｍ−Ｖｎｍを算出するためのサブルーチンを図９のフローチャートに沿って説明する。図９は、体積算出処理の流れを示すフローチャートである。ここでは、アモルファス体の空間を分割する断片の個数をＭ個とする。ｉは、Ｍ個の断片を更新してゆくためのカウンタである。
【００５４】
規定値｜ψｒ｜及び断片の体積Ｖｉは、ユーザによる入力などによって予め決められ、パソコンの記憶装置に格納されている。また、各断片における状態ｎの波動関数の絶対値｜ψｉｎ｜は、既存の分子シミュレーションによって予め算出され、パソコンの記憶装置に格納されている。また、各断片における状態ｍの波動関数の絶対値｜ψｉｍ｜は、既存の分子シミュレーションによって予め算出され、パソコンの記憶装置に格納されている。
【００５５】
サブルーチンが起動すると、状態ｎの波動関数の体積Ｖｎ、状態ｍの波動関数の体積Ｖｍ、重なり領域の体積Ｖｎｍに０をそれぞれ代入して初期化する（Ｓ２０）。また、ｉに０を代入して初期化する（Ｓ２１）。
【００５６】
初期化後、ループ毎に、次の断片に更新するために、ｉをインクリメントする（Ｓ２２）。断片ｉにおける状態ｎの波動関数の絶対値｜ψｉｎ｜が規定値｜ψｒ｜より大きいか否かを判定する（Ｓ２３）。Ｓ２３にて｜ψｉｎ｜が｜ψｒ｜より大きいと判定した場合（状態ｎの波動関数が断片ｉに存在する場合）、状態ｎの波動関数の体積Ｖｎの前回値に断片の体積Ｖｉを加算し、その加算値を状態ｎの波動関数の体積Ｖｎに代入する（Ｓ２４）。さらに、断片ｉにおける状態ｍの波動関数の絶対値｜ψｉｍ｜が規定値｜ψｒ｜より大きいか否かを判定する（Ｓ２５）。Ｓ２５にて｜ψｉｍ｜が｜ψｒ｜より大きいと判定した場合（状態ｎの波動関数と状態ｍの波動関数が共に断片ｉに存在する場合）、重なり領域の体積Ｖｎｍの前回値に断片の体積Ｖｉを加算し、その加算値を重なり領域の体積Ｖｎｍに代入するとともに（Ｓ２６）、状態ｍの波動関数の体積Ｖｍの前回値に断片の体積Ｖｉを加算し、その加算値を状態ｍの波動関数の体積Ｖｍに代入する（Ｓ２７）。Ｓ２３にて｜ψｉｎ｜が｜ψｒ｜以下と判定した場合、断片ｉにおける状態ｍの波動関数の絶対値｜ψｉｍ｜が規定値｜ψｒ｜より大きいか否かを判定する（Ｓ２８）。Ｓ２８にて｜ψｉｍ｜が｜ψｒ｜より大きいと判定した場合（状態ｍの波動関数が断片ｉに存在する場合）、状態ｍの波動関数の体積Ｖｍの前回値に断片の体積Ｖｉを加算し、その加算値を状態ｍの波動関数の体積Ｖｍに代入する（Ｓ２９）。ちなみに、Ｓ２３にて｜ψｉｎ｜が｜ψｒ｜以下と判定しかつＳ２８にて｜ψｉｍ｜が｜ψｒ｜以下と判定した場合（状態ｎの波動関数と状態ｍの波動関数が共に断片ｉに存在しない場合）、Ｖｎ、Ｖｍ、Ｖｎｍ全て更新されない。次に、ｉがＭより小さいか否かを判定する（Ｓ３０）。Ｓ３０にてｉがＭより小さいと判定した場合、Ｓ２２の処理に戻り、次の断片についての処理を行う。
【００５７】
一方、Ｓ３０にてｉがＭに等しくなったと判定した場合、状態ｎの波動関数の体積Ｖｎ、状態ｎの波動関数の体積Ｖｍ、重なり領域の体積Ｖｎｍを算出したので、Ｖｎ，Ｖｍ，Ｖｎｍを用いてＶｎ＋Ｖｍ−Ｖｎｍを算出する（Ｓ３１）。この算出されたＶｎｍ及びＶｎ＋Ｖｍ−Ｖｎｍは、メインルーチンで使用される。
【００５８】
次に、上記の２つの目の算出方法によるＶｎｍ及びＶｎ＋Ｖｍ−Ｖｎｍを算出するためのサブルーチンを図１０のフローチャートに沿って説明する。図１０は、体積算出処理の流れを示すフローチャートである。ここでは、式（３）で示すアモルファス体の全体積Ｖの範囲での積分をＭ個の区間に分割して積分を行う。ｉは、Ｍまでカウントするためのカウンタである。Ｖｎ＋Ｖｍ−Ｖｎｍの算出過程の変数をＶａとする。
【００５９】
各区間における状態ｎの波動関数の絶対値｜ψｉｎ｜は、既存の分子シミュレーションによって予め算出され、パソコンの記憶装置に格納されている。また、各区間における状態ｍの波動関数の絶対値｜ψｉｍ｜は、既存の分子シミュレーションによって予め算出され、パソコンの記憶装置に格納されている。
【００６０】
サブルーチンが起動すると、Ｖａ、重なり領域の体積Ｖｎｍに０をそれぞれ代入して初期化する（Ｓ４０）。また、ｉに０を代入して初期化する（Ｓ４１）。
【００６１】
初期化後、ループ毎に、ｉをインクリメントする（Ｓ４２）。区間ｉにおける状態ｎの波動関数の絶対値｜ψｉｎ｜を用いて、式（４）によって区間ｉにおける状態ｎの波動関数を算出し、その算出値をφｎに代入する（Ｓ４３）。また、区間ｉにおける状態ｍの波動関数の絶対値｜ψｉｍ｜を用いて、式（４）によって区間ｉにおける状態ｍの波動関数を算出し、その算出値をφｍに代入する（Ｓ４４）。さらに、式（３）に基づいて、その算出したφｎとφｍを用いて｜φｎ−φｍ｜×Ｖｉを算出し、その算出値をＶａの前回値に加算し、その加算値をＶａに代入する（Ｓ４５）。また、式（５）に基づいて、その算出したφｎとφｍを用いて（√（φｎ×φｍ））×Ｖｉを算出し、その算出値をＶｎｍの前回値に加算し、その加算値をＶｎｍに代入する（Ｓ４６）。次に、ｉがＭより小さいか否かを判定する（Ｓ４７）。Ｓ４７にてｉがＭより小さいと判定した場合、Ｓ４２の処理に戻り、次の区間についての処理を行う。
【００６２】
一方、Ｓ４７にてｉがＭに等しくなったと判定した場合、重なり領域の体積Ｖｎｍ及びＶａを算出したので、ＶａをＶｎ＋ＶＭ−Ｖｎｍに代入する（Ｓ４８）。この算出されたＶｎｍ及びＶｎ＋ＶＭ−Ｖｎｍは、メインルーチンで使用される。
【００６３】
この分子シミュレーション装置１によれば、２つの異なる波動関数におけるキャリアの乗り移り易さ、空間上の広さ及びキャリアが存在する確率の各評価指標を加味した演算式を構築し、シミュレーションすることにより、アモルファス体の電気伝導度を精度良く推定することができる。
【００６４】
特に、分子シミュレーション装置１では、２つの異なる波動関数のエネルギレベルの差及び重なり領域の体積の各評価指標をその演算式に加味することにより、２つの異なる波動関数におけるキャリアの乗り移り易さを高精度に評価することができる。また、分子シミュレーション装置１では、２つの異なる波動関数の各体積の和から重なり領域の体積を引いた体積をその演算式に加味することにより、２つの異なる波動関数における空間上の広さを高精度に評価することができる。
【００６５】
さらに、分子シミュレーション装置１では、電気伝導度（代用評価値）を算出するときに、２つの異なる波動関数が共にフェルミエネルギより高い固有エネルギを持つ場合には電子による電気伝導を加味し、２つの異なる波動関数が共にフェルミエネルギより低い固有エネルギを持つ場合にはホールによる電気伝導を加味し、それ以外の場合には加味しないことにより、より高精度にアモルファス体の電気伝導度を推定することができる。
【００６６】
また、分子シミュレーション装置１では、アモルファス体を多数の断片に分割し、各断片における波動関数の絶対値が規定値より大きいか否か（その断片に波動関数に存在するか否か）によって波動関数に関する各体積を算出することにより、波動関数に関する各体積を簡単に算出することができる。また、分子シミュレーション装置１では、波動関数を頭打ち加工した関数とし、波動関数を連続的な値として波動関数に関する各体積を算出することにより、波動関数に関する各体積を高精度に算出することができる。
【００６７】
以上、本発明に係る実施の形態について説明したが、本発明は上記実施の形態に限定されることなく様々な形態で実施される。
【００６８】
例えば、本実施の形態ではアモルファス体の電気伝導度推定処理用のプログラムを実行することによってパソコン上に構成される分子シミュレーション装置に適用したが、分子シミュレーション装置の専用機として構成してもよいし、あるいは、ＣＤ−ＲＯＭなどの記憶媒体に格納されたプログラム、インタネットなどのネットワークを介して利用可能なプログラムなどに適用してもよい。
【００６９】
また、本実施の形態では式（１）において異なる波動関数として異なる２つの波動関数の組み合わせで評価を行う一例を示したが、異なる３つの波動関数の組み合わせ、異なる４つの波動関数の組み合わせなどで評価を行ってもよい。
【００７０】
また、本実施の形態では式（１）において評価指標１〜４による各評価を反映したファクタの一例を示したが、評価指標１〜４による各評価を反映したファクタについては他のものでもよい。
【００７１】
また、本実施の形態では波動関数に関する体積を求める手法の二例を示したが、他の手法によって波動関数の体積を求めてもよい。
【符号の説明】
【００７２】
１…分子シミュレーション装置、１０…パーソナルコンピュータ本体、１１…ディスプレイ、１２…キーボード
【技術分野】
【０００１】
本発明は、アモルファスの電気伝導度を推定する分子シミュレーション装置に関する。
【背景技術】
【０００２】
従来、半導体などの物質の電気特性を推定する様々な手法が提案されている。特許文献１には、量子力学的な効果を近似的に考慮した半導体の電気特性を推定する手法が開示されている。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開平１１−３９９８号公報
【特許文献２】特開２００７−１３４５７４号公報
【特許文献３】特開平５−１５１１９１号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
上記の提案されている手法は、特定の結晶状態を有する物質に適用した手法である。しかし、アモルファスは、特定の結晶構造を有していないので、このような従来の手法を適用して電気特性を推定することができない。
【０００５】
そこで、本発明は、アモルファスの電気伝導度を推定する分子シミュレーション装置を提供することを課題とする。
【課題を解決するための手段】
【０００６】
本発明に係る分子シミュレーション装置は、アモルファスの電気伝導度を推定する分子シミュレーション装置であって、アモルファスに存在する異なる波動関数のキャリアの乗り移り易さ、空間上の広さ及びキャリアが存在する確率をそれぞれ評価し、当該３つの評価に基づいて電気伝導度を推定することを特徴とする。
【０００７】
アモルファスは、周期的な結晶構造を有さないので、波動関数が局所的に存在している。そのため、アモルファスでキャリアが移動して電気伝導が起こるには、アモルファスに存在する異なる波動関数間をキャリアが乗り移って電気伝導が行われると考えられる。さらに、その異なる波動関数においてそれぞれ異なる空間を広く占めていると、キャリアが他の波動関数に乗り移れる可能性が高くなる。また、波動関数にキャリアが多く存在していると、キャリアが他の波動関数に乗り移る可能性が高くなる。そこで、この分子シミュレーション装置では、アモルファスに存在する異なる波動関数についてのキャリアの乗り移り易さ、空間上の広さ及びキャリアが存在する確率をそれぞれ評価する。そして、分子シミュレーション装置では、その３つの評価に基づいてアモルファスの電気伝導度を推定する。このように、分子シミュレーション装置では、異なる波動関数についてのキャリアの乗り移り易さ、空間上の広さ及びキャリアが存在する確率をそれぞれ評価することにより、アモルファスの電気伝導度を精度良く推定することができる。
【０００８】
本発明の上記分子シミュレーション装置では、異なる波動関数のキャリアの乗り移り易さを、異なる波動関数のエネルギレベルの差及び異なる波動関数の重なり領域の体積に基づいて評価すると好適である。
【０００９】
異なる波動関数間のエネルギレベルの差が小さいと、キャリアがその波動関数間を乗り移り易い。また、異なる波動関数が重なっている領域が大ききと、キャリアはその波動関数間を乗り移り易い。そこで、この分子シミュレーション装置では、異なる波動関数ついてのエネルギレベルの差及び重なり領域の体積を評価することにより、キャリアの乗り移り易さを高精度に評価することができる。
【００１０】
本発明の上記分子シミュレーション装置では、異なる波動関数の空間上の広さを、異なる波動関数の各体積の和から重なり領域の体積を引いた体積に基づいて評価すると好適である。
【００１１】
異なる波動関数がそれぞれ異なる空間（重なり領域以外の空間）を広く占めていると、キャリアがその異なる波動関数間を乗り移った後に他の波動関数に乗り移れる可能性が高くなる。そこで、この分子シミュレーション装置では、異なる波動関数の各体積の和から重なり領域の体積を引いた体積で評価することにより、異なる波動関数における空間上の広さを高精度に評価することができる。
【発明の効果】
【００１２】
本発明は、異なる波動関数についてのキャリアの乗り移り易さ、空間上の広さ及びキャリアが存在する確率をそれぞれ評価することにより、アモルファスの電気伝導度を精度良く推定することができる。
【図面の簡単な説明】
【００１３】
【図１】本実施の形態に係る分子シミュレーション装置の構成図である。
【図２】アモルファス体の構造の一例を示す図である。
【図３】実空間上の異なる２つの波動関数の固有エネルギの一例を示す図である。
【図４】実空間上の異なる２つの波動関数の重なり状態の一例を示す図である。
【図５】アモルファス体の空間分割の一例を示す図である。
【図６】波動関数のψとφの関係の一例を示すグラフである。
【図７】電気伝導度推定処理（基本版）の流れを示すフローチャートである。
【図８】電気伝導度推定処理（改良版）の流れを示すフローチャートである。
【図９】体積算出処理の流れを示すフローチャートである。
【図１０】体積算出処理の流れを示すフローチャートである。
【発明を実施するための形態】
【００１４】
以下、図面を参照して、本発明に係る分子シミュレーション装置の実施の形態を説明する。
【００１５】
本実施の形態では、本発明に係る分子シミュレーション装置を、パーソナルコンピュータ（以下、パソコンと記載）上に構成される分子シミュレーション装置に適用する。本実施の形態に係る分子シミュレーション装置では、アモルファスの電気伝導度（代用評価指標）を推定する。
【００１６】
図１〜図６を参照して、本実施の形態に係る分子シミュレーション装置１について説明する。図１は、本実施の形態に係る分子シミュレーション装置の構成図である。図２は、アモルファス体の構造の一例を示す図である。図３は、実空間上の異なる２つの波動関数の固有エネルギの一例を示す図である。図４は、実空間上の異なる２つの波動関数の重なり状態の一例を示す図である。図５は、アモルファス体の空間分割の一例を示す図である。図６は、波動関数のψとφの関係の一例を示すグラフである。
【００１７】
分子シミュレーション装置１は、記憶装置（ハードディスクなど）に格納されているアモルファス体の電気伝導度推定処理用のプログラムをパソコン本体１０で実行することによって構成され、その推定結果をディスプレイ１１に表示したりあるいは記憶装置に格納する。この際、ユーザは、このアモルファス体の電気伝導度推定に必要な情報をキーボード１２などを利用して入力し、推定結果をディスプレイ１１などで確認する。
【００１８】
分子シミュレーション装置１におけるアモルファス体の電気伝導度推定処理について具体的に説明する前に、アモルファス体の電気伝導度の推定方法及びその推定において必要な波動関数の体積の算出方法について説明する。
【００１９】
アモルファスに限らず、物質中の電子は、物質中に存在する多数の波動関数のいずれかに位置しており、その波動関数の広がりの範囲を自由に移動している。
【００２０】
シリコンなどの半導体は、周期的な結晶構造を有するため、波動関数も周期構造を持っている。このような半導体中の電子は、同じ波動関数に位置したままで自由に移動できるが、実際には他の電子や熱振動、不純物、格子欠陥との衝突によって移動が阻害される。それが、電気伝導度を有限の値としている主な要因である。したがって、周期的な結晶構造の半導体の電気伝導度をシミュレーションによって予測するには、他の電子や熱振動、不純物、格子欠陥との衝突の頻度を計算すればよい。
【００２１】
しかし、アモルファス体（アモルファスの半導体）は、周期的な結晶構造を有さないため、波動関数が局所的に存在している。図２には、アモルファス体ＡＭの構造の一例を示しており、多数の原子ＡＴ，・・・が不規則に存在しており、多数の波動関数ＷＦ，・・・が局所的に分布している。このアモルファス体の電子が移動して電気伝導が起こるには、電子が異なる波動関数間を乗り移ることが必要と考えられる。したがって、アモルファス体の場合、異なる波動関数における電子（キャリア）の乗り移り易さが電気伝導度を決定していると考えられる。
【００２２】
キャリア（電子、ホール）が固有エネルギの異なる波動関数間を乗り移るには、そのエネルギ差を熱振動によるフォノンのやりとりなどで補償しなければならない。図３の例で示すように、固有エネルギＥｍの波動関数と固有エネルギＥｎの波動関数とのエネルギ差が熱エネルギ程度（２６ｍｅＶ程度）のエネルギ内ならやりとりは起こり易いが、それ以上のエネルギ差の場合には変化させることは難しい。したがって、異なる波動関数間のエネルギ差が小さいほど、キャリアが乗り移り易い（評価指標１）。
【００２３】
キャリアが異なる波動関数間を乗り移る場合、その波動関数同士の空間上の重なりが大きいほど乗り移り易い（評価指標２）。図４の例で示すように、波動関数ＡＦｍと波動関数ＡＦｎとは空間上で重なり領域ＯＬを有しており、この重なり領域ＯＬを利用して波動関数ＡＦｍと波動関数ＡＦｎとの間をキャリアが移動する。
【００２４】
また、異なる波動関数間を乗り移り易さだけでなく、それによって伝導が起こり易いことが必要である。それには、異なる波動関数間が乗り移り易かったとしても、その異なる
波動関数が同じ空間（重なり領域が大部分）を占めていては意味がなく、キャリアが波動関数間を乗り移った後にそれ以外の波動関数に移動するためには、異なる波動関数がそれぞれ異なる空間（重なり領域以外の空間）を広く占めていることも重要となる（評価指標３）。
【００２５】
また、キャリアが異なる波動関数間を乗り移る場合、個々の波動関数にキャリアが存在する確率が高いことも重要となる（評価指標４）。波動関数に存在するキャリアが少ないと、キャリアが他の波動関数に乗り移る可能性も低くなる。
【００２６】
この４つの評価指標が、アモルファス体における電気伝導に寄与すると考えられる。したがって、この４つの評価指標に基づいた評価をそれぞれ行い、その４つの評価結果を加味して電気伝導度を推定する。
【００２７】
実際の電気伝導では、キャリアが多数の波動関数をランダムに乗り移って伝導すると考えられる。その組み合わせの数は膨大な数になるので、全ての組み合わせについて評価することは実用的でない。そこで、少数の波動関数の組み合わせで評価し、その全ての組み合わせの評価値の和によって電気伝導のし易さを表わすことにより、実際の電気伝導度を高精度に推定（評価）できる。
【００２８】
ここでは、異なる２つの波動関数の組み合わせで評価する場合について説明する。上記した４つの評価指標による評価を加味して電気伝導度（電気伝導率）の代用評価値Ｆを算出する式（１）を示す。式（１）において、ｎ，ｍは、２つの異なる波動関数を示している。
【数１】
【００２９】
Ｖｎ，Ｖｍが各波動関数の体積であり、Ｖｎｍが２つの波動関数の重なり領域の体積であり、各体積Ｖｎ，Ｖｍ，Ｖｎｍは下記の方法で算出される。Ｅｎ，Ｅｍは、各波動関数の固有エネルギであり、既存の分子シミュレーションによって算出される。Ｅｆは、フェルミエネルギであり、既存の分子シミュレーションによって算出される。Ｅｔｈは、定数であり、熱エネルギ（２６ｍｅＶ）程度である。ＫＢは、ボルツマン定数である。Ｔは、系の絶対温度であり、ユーザによる入力などによって予め決まった値である。
【００３０】
式（１）のΣは、ｎ，ｍで表される異なる２つの波動関数の組み合わせについて、ｎ＝ｍでない全ての組み合わせについての評価値を加算することを表わしている。
【００３１】
式（１）のｅｘｐ｛｜Ｅｎ−Ｅｍ｜／Ｅｔｈ｝は、２つの波動関数の固有エネルギの差を定数Ｅｔｈで除算することにより評価指標１による評価を行うファクタである。式（１）の×Ｖｎｍは、２つの波動関数の重なり領域の体積を乗算することにより評価指標２による評価を行うファクタである。式（１）の（Ｖｎ＋Ｖｍ−Ｖｎｍ）は、２つの波動関数の体積の和から重なり領域の体積を減算することにより評価指標３による評価を行うファクタである。式（１）のｅｘｐ｛（Ｅｆ−（Ｅｎ＋Ｅｍ）／２）／（ＫＢ×Ｔ）｝は、量子力学論に基づいて２つの波動関数の平均エネルギにおける存在確率を乗算することにより評価指標４による評価を行うファクタである。
【００３２】
上記の基本的なアモルファス体の電気伝導度の推定方法の考え方では、アモルファス体に存在する全ての波動関数におけるキャリア（特に、電子）が電気伝導に寄与するとみなしている。しかし、実際の物質では、エネルギの低い波動関数は、全て電子が存在し、フェルミの排他律という法則に従ってそれ以上電子を取り込めない状態になっている。そのため、電子は、身動きが取れないので、伝導に寄与できない。半導体において、実際に伝導に寄与しうるのは、フェルミエネルギよりも高い固有エネルギを持つ波動関数に存在する電子か、逆に、フェルミエネルギよりも低い固有エネルギを持つ波動関数において電子が欠けた状態にあるホールと呼ばれる仮想的な正電荷の粒子である。
【００３３】
そこで、上記の基本的な方法を改良し、２つの波動関数が共にフェルミエネルギよりも高い固有エネルギを持つ場合の電子による伝導と２つの波動関数が共にフェルミエネルギよりも低い固有エネルギを持つ場合のホールによる伝導に分けて評価を行う。具体的には、状態ｎの波動関数の固有エネルギＥｎと状態ｍの波動関数の固有エネルギＥｍが共にフェルミエネルギＥｆよりも高い場合、電子がその２つの波動関数間を乗り移る（電気伝導に寄与する）ことができるので、電気伝導度に電子による伝導を加味する。状態ｎの波動関数の固有エネルギＥｎと状態ｍの波動関数の固有エネルギＥｍが共にフェルミエネルギＥｆよりも低い場合、ホールがその２つの波動関数間を乗り移る（電気伝導に寄与する）ことができるので、電気伝導度にホールによる電気伝導を加味する。状態ｎの波動関数の固有エネルギＥｎと状態ｍの波動関数の固有エネルギＥｍのうちの一方がフェルミエネルギＥｆよりも高く、他方がフェルミエネルギＥｆよりも低い場合（禁制帯を挟んで２つの波動関数が存在する場合）、キャリアがその２つの波動関数間を乗り移ることができない（電気伝導に寄与しない）ので、電気伝導度には加味しない。
【００３４】
各体積Ｖｎ，Ｖｍ，Ｖｎｍ及びＶｎ＋Ｖｍ−Ｖｎｍの算出方法については説明する。波動関数は、連続的な値を持って空間に分布している。そのため、波動関数の広さ（体積）などの概念は、厳密な定義がなく、特別に定める必要がある。そこで、以下に、２つの算出方法について説明する。
【００３５】
まず、１つの目の算出方法について説明する。この方法では、波動関数は連続的な値を持つが、その波動関数の絶対値｜ψ｜が所定の規定値｜ψｒ｜よりも大きければ「そこに波動関数が存在する」とみなし、規定値｜ψｒ｜よりも小さければ「そこに波動関数が存在しない」とみなす。具体的な方法としては、図５に示すように、アモルファス体の空間を同じ体積を有する多数の小さい断片（直方体）に分割する。既存の分子シミュレーションによって、状態ｍの波動関数の絶対値｜ψｍ｜と状態ｎの波動関数の絶対値｜ψｎ｜が求めることができ、さらに、任意の断片ｉについて状態ｍの波動関数の絶対値｜ψｉｍ｜と状態ｎの波動関数の絶対値｜ψｉｎ｜が予め求まることができる。そして、断片ｉの状態ｍの波動関数の絶対値｜ψｉｍ｜が規定値｜ψｒ｜よりも大きいか否かを判定し、｜ψｉｍ｜が｜ψｒ｜よりも大きいと判定した場合には状態ｍの波動関数の体積Ｖｍに断片ｉの体積Ｖｉを加算し、この処理を空間の全ての断片について行い、状態ｍの波動関数の体積Ｖｍを求める。状態ｎの波動関数の体積Ｖｎについても同様に求める。また、断片ｉの状態ｍの波動関数の絶対値｜ψｉｍ｜が規定値｜ψｒ｜よりも大きいかつ状態ｎの波動関数の絶対値｜ψｉｎ｜が規定値｜ψｒ｜よりも大きいか否かを判定し、｜ψｉｍ｜と｜ψｉｎ｜が共に規定値｜ψｒ｜より大きいと判定した場合には状態ｍの波動関数と状態ｎの波動関数との重なり領域の体積Ｖｎｍに断片ｉの体積Ｖｉを加算し、この処理を空間の全ての断片について行い、重なり領域の体積Ｖｎｍを求める。最後に、Ｖｎ＋Ｖｍ−Ｖｎｍを算出する。なお、規定値｜ψｒ｜は、ユーザによって指定される。
【００３６】
次に、２つの目の算出方法について説明する。この方法では、波動関数をそのまま連続的な値として取り扱うために、波動関数が式（２）で示す規格化条件を満たしているものとする。この式（２）において、Ｖは、積分範囲としてアモルファス体の空間の全体積について積分を行うということを意味する。この式（２）を満たさない場合、ψに適当な係数を乗算した値を改めてψと置き直すことにより、式（２）を満たすように調整する。
【数２】
【００３７】
そして、各体積Ｖｎ，Ｖｍ，Ｖｎｍを個別に求めるのではなく、式（３）によってＶｎ＋Ｖｍ−Ｖｎｍをまとめて求める。
【数３】
【００３８】
式（３）において、φは波動関数の値である。その２つのφの差｜φｎ−φｍ｜は、状態ｎの波動関数か状態ｍの波動関数のどちらかが存在するということを満たす指標として用いることができ、その値を全空間（全体積Ｖ）で積分することにより、（Ｖｎ＋Ｖｍ−Ｖｎｍ）の概念を示す指標を求めることができる。
【００３９】
Ｖｎ＋Ｖｍ−Ｖｎｍの意味するところは、状態ｎの波動関数と状態ｍの波動関数の占める空間において重なり領域以外で占める領域の体積である。基本的には、波動関数が有るか無いかの概念であるので、波動関数が非常に大きな値を持っていても「有る」という判定以上のものはない。したがって、図６に示すように、波動関数の値φを頭打ち（上限を１とする）に加工した関数を用いる。図６に示すようなグラフの関係が得られる適切な関数としては、例えば、式（４）がある。
【数４】
【００４０】
さらに、式（１）で電気伝導度を算出するためには、Ｖｎｍが、必要となるので、別個に求める。Ｖｎｍは、状態ｎの波動関数と状態ｍの波動関数の両方に存在する領域（重なり領域）の体積の概念である。したがって、状態ｎの波動関数の値φｎと状態ｍの波動関数の値φｍの積φｎ×φｍの平方根の積分により、その概念に近い指標を求めることができる。その具体的な式としては、式（５）で示す。
【数５】
【００４１】
それでは、分子シミュレーション装置１におけるアモルファス体の電気伝導度の推定処理について説明する。アモルファス体の電気伝導度推定処理用のプログラムをパソコン本体１０で実行することにより、電気伝導度の代用評価値Ｆを求めるためのメインルーチンが実行され、そのメインルーチンによってＶｎｍ及びＶｎ＋Ｖｍ−Ｖｎｍを算出するためのサブルーチンが呼び出されて実行される。ここでは、メインルーチンについては、上記で説明した電気伝導度の推定方法に基づく基本版とその基本版から改良した改良版について説明する。また、サブルーチンについては、上記で説明した２つの算出方法に基づくものをそれぞれ説明する。
【００４２】
まず、電気伝導度の代用評価値Ｆを求めるためのメインルーチンの基本版を図７のフローチャートに沿って説明する。図７は、電気伝導度推定処理（基本版）の流れを示すフローチャートである。ここでは、アモルファス体に存在する波動関数の個数をＮ個とする。ｎ，ｍは、Ｎ個の波動関数中の２つの状態ｎの波動関数と状態ｍの波動関数の組み合わせを更新してゆくためのカウンタである。
【００４３】
系の絶対温度Ｔは、ユーザによる入力などによって予め決められ、パソコンの記憶装置に格納されている。また、Ｎ個の波動関数の固有エネルギＥは、既存の分子シミュレーションによって予め算出され、パソコンの記憶装置に格納されている。また、フェルミエネルギＥｆは、既存の分子シミュレーションによって予め算出され、パソコンの記憶装置に格納されている。また、波動関数の個数Ｎは、既存の分子シミュレーションによって予め算出され、パソコンの記憶装置に格納されている。なお、Ｎ個をユーザ入力などによって規定し、Ｎ個までの波動関数でシミュレーションを行ってもよい。
【００４４】
メインルーチンが起動すると、ｍに１を代入して初期化するとともに、代用評価値Ｆに０を代入して初期化する（Ｓ１０）。
【００４５】
初期化後、ループ毎に、状態ｍの波動関数を次の波動関数に更新するために、ｍをインクリメントする（Ｓ１１）。また、ｎに１を代入して初期化する（Ｓ１２）。そして、今回の組み合わせの状態ｎの波動関数と状態ｍの波動関数についての体積Ｖｎｍ及びＶｎ＋Ｖｍ−Ｖｎｍを算出するためのサブルーチンを呼び出して実行し、Ｖｎｍ及びＶｎ＋Ｖｍ−Ｖｎｍの算出する（Ｓ１３）。さらに、Ｓ１３で算出されたＶｎｍ，Ｖｎ＋Ｖｍ−Ｖｎｍ及び記憶装置から取り出した固有エネルギＥｎ，Ｅｍを用いて、式（６）によって、今回の組み合わせの状態ｎの波動関数と状態ｍの波動関数でのΔＦを算出する（Ｓ１４）。そして、代用評価値Ｆの前回値に今回算出したΔＦを加算し、その加算値を代用評価値Ｆに代入する（Ｓ１５）。そして、状態ｎの波動関数を次の波動関数に更新するために、ｎをインクリメントする（Ｓ１６）。さらに、ｎがｍよりも小さいか否かを判定する（Ｓ１７）。Ｓ１７にてｎがｍより小さいと判定した場合、Ｓ１３の処理に戻り、状態ｎの波動関数を変えた次の組み合わせの処理を行う。一方、Ｓ１７にてｎがｍに等しくなったと判定した場合、ｍがＮより小さいか否かを判定する（Ｓ１８）。Ｓ１８にてｍがＮより小さいと判定した場合、Ｓ１１の処理に戻り、状態ｍの波動関数を次の波動関数に変えた組み合わせの処理を行う。
【数６】
【００４６】
Ｓ１８にてｍがＮと等しくなったと判定した場合、代用評価値Ｆの算出が終了したので、その代用評価値Ｆをアモルファス体の電気伝導度としてディスプレイ１１に表示するかあるいはパソコンの記憶装置に格納する。
【００４７】
次に、電気伝導度の代用評価値Ｆを求めるためのメインルーチンの改良版を図８のフローチャートに沿って説明する。図８は、電気伝導度推定処理（改良版）の流れを示すフローチャートである。
【００４８】
この改良版では、上記の基本版と比較すると、基本版におけるΔＦの算出（Ｓ１４）だけが異なっており、他の処理は同じである。そこで、このΔＦの算出についてのみ説明する。
【００４９】
Ｓ１３でＶｎｍ及びＶｎ＋Ｖｍ−Ｖｎｍを算出すると、状態ｎの波動関数の固有エネルギＥｎがフェルミエネルギＥｆよりも高いかつ状態ｍの波動関数の固有エネルギＥｍがフェルミエネルギＥｆよりも高いか否かを判定する（Ｓ１４ａ）。Ｓ１４ａにてＥｎがＥｆよりも高いかつＥｍがＥｆよりも高いと判定した場合（電子が２つの波動関数間を乗り移る場合）、Ｓ１３で算出されたＶｎｍ，Ｖｎ＋Ｖｍ−Ｖｎｍ及び記憶装置から取り出した固有エネルギＥｎ，Ｅｍを用いて、式（６）によって、今回の組み合わせの状態ｎの波動関数と状態ｍの波動関数でのΔＦを算出する（Ｓ１４ｂ）。
【００５０】
一方、Ｓ１４ａにてＥｎとＥｍの少なくとも一方がＥｆよりも低いと判定した場合、状態ｎの波動関数の固有エネルギＥｎがフェルミエネルギＥｆよりも低いかつ状態ｍの波動関数の固有エネルギＥｍがフェルミエネルギＥｆよりも低いか否かを判定する（Ｓ１４ｃ）。Ｓ１４ｃにてＥｎがＥｆよりも低いかつＥｍがＥｆよりも低いと判定した場合（ホールが２つの波動関数間を乗り移る場合）、Ｓ１３で算出されたＶｎｍ，Ｖｎ＋Ｖｍ−Ｖｎｍ及び記憶装置から取り出した固有エネルギＥｎ，Ｅｍを用いて、式（７）によって、今回の組み合わせの状態ｎの波動関数と状態ｍの波動関数でのΔＦを算出する（Ｓ１４ｄ）。なお、式（７）では、ホールの存在確率とエネルギの関係は電子と逆になるので、（Ｅｍ＋Ｅｎ）／２からフェルミエネルギＥｆを減算している（式（６）と逆）。
【数７】
【００５１】
Ｓ１４ａ、Ｓ１４ｃにてＥｎとＥｍのうちの一方がＥｆよりも高く、他方がＥｆよりも低いと判定した場合（キャリアが２つの波動関数間を乗り移ることができない場合）、ΔＦに０を代入する（Ｓ１４ｅ）。
【００５２】
そして、代用評価値Ｆの前回値に今回算出したΔＦを加算し、その加算値を代用評価値Ｆに代入する（Ｓ１５）。ただし、Ｓ１４ｅでΔＦに０を代入した場合、代用評価値Ｆは前回値のままである。
【００５３】
次に、上記の１つの目の算出方法によるＶｎｍ及びＶｎ＋Ｖｍ−Ｖｎｍを算出するためのサブルーチンを図９のフローチャートに沿って説明する。図９は、体積算出処理の流れを示すフローチャートである。ここでは、アモルファス体の空間を分割する断片の個数をＭ個とする。ｉは、Ｍ個の断片を更新してゆくためのカウンタである。
【００５４】
規定値｜ψｒ｜及び断片の体積Ｖｉは、ユーザによる入力などによって予め決められ、パソコンの記憶装置に格納されている。また、各断片における状態ｎの波動関数の絶対値｜ψｉｎ｜は、既存の分子シミュレーションによって予め算出され、パソコンの記憶装置に格納されている。また、各断片における状態ｍの波動関数の絶対値｜ψｉｍ｜は、既存の分子シミュレーションによって予め算出され、パソコンの記憶装置に格納されている。
【００５５】
サブルーチンが起動すると、状態ｎの波動関数の体積Ｖｎ、状態ｍの波動関数の体積Ｖｍ、重なり領域の体積Ｖｎｍに０をそれぞれ代入して初期化する（Ｓ２０）。また、ｉに０を代入して初期化する（Ｓ２１）。
【００５６】
初期化後、ループ毎に、次の断片に更新するために、ｉをインクリメントする（Ｓ２２）。断片ｉにおける状態ｎの波動関数の絶対値｜ψｉｎ｜が規定値｜ψｒ｜より大きいか否かを判定する（Ｓ２３）。Ｓ２３にて｜ψｉｎ｜が｜ψｒ｜より大きいと判定した場合（状態ｎの波動関数が断片ｉに存在する場合）、状態ｎの波動関数の体積Ｖｎの前回値に断片の体積Ｖｉを加算し、その加算値を状態ｎの波動関数の体積Ｖｎに代入する（Ｓ２４）。さらに、断片ｉにおける状態ｍの波動関数の絶対値｜ψｉｍ｜が規定値｜ψｒ｜より大きいか否かを判定する（Ｓ２５）。Ｓ２５にて｜ψｉｍ｜が｜ψｒ｜より大きいと判定した場合（状態ｎの波動関数と状態ｍの波動関数が共に断片ｉに存在する場合）、重なり領域の体積Ｖｎｍの前回値に断片の体積Ｖｉを加算し、その加算値を重なり領域の体積Ｖｎｍに代入するとともに（Ｓ２６）、状態ｍの波動関数の体積Ｖｍの前回値に断片の体積Ｖｉを加算し、その加算値を状態ｍの波動関数の体積Ｖｍに代入する（Ｓ２７）。Ｓ２３にて｜ψｉｎ｜が｜ψｒ｜以下と判定した場合、断片ｉにおける状態ｍの波動関数の絶対値｜ψｉｍ｜が規定値｜ψｒ｜より大きいか否かを判定する（Ｓ２８）。Ｓ２８にて｜ψｉｍ｜が｜ψｒ｜より大きいと判定した場合（状態ｍの波動関数が断片ｉに存在する場合）、状態ｍの波動関数の体積Ｖｍの前回値に断片の体積Ｖｉを加算し、その加算値を状態ｍの波動関数の体積Ｖｍに代入する（Ｓ２９）。ちなみに、Ｓ２３にて｜ψｉｎ｜が｜ψｒ｜以下と判定しかつＳ２８にて｜ψｉｍ｜が｜ψｒ｜以下と判定した場合（状態ｎの波動関数と状態ｍの波動関数が共に断片ｉに存在しない場合）、Ｖｎ、Ｖｍ、Ｖｎｍ全て更新されない。次に、ｉがＭより小さいか否かを判定する（Ｓ３０）。Ｓ３０にてｉがＭより小さいと判定した場合、Ｓ２２の処理に戻り、次の断片についての処理を行う。
【００５７】
一方、Ｓ３０にてｉがＭに等しくなったと判定した場合、状態ｎの波動関数の体積Ｖｎ、状態ｎの波動関数の体積Ｖｍ、重なり領域の体積Ｖｎｍを算出したので、Ｖｎ，Ｖｍ，Ｖｎｍを用いてＶｎ＋Ｖｍ−Ｖｎｍを算出する（Ｓ３１）。この算出されたＶｎｍ及びＶｎ＋Ｖｍ−Ｖｎｍは、メインルーチンで使用される。
【００５８】
次に、上記の２つの目の算出方法によるＶｎｍ及びＶｎ＋Ｖｍ−Ｖｎｍを算出するためのサブルーチンを図１０のフローチャートに沿って説明する。図１０は、体積算出処理の流れを示すフローチャートである。ここでは、式（３）で示すアモルファス体の全体積Ｖの範囲での積分をＭ個の区間に分割して積分を行う。ｉは、Ｍまでカウントするためのカウンタである。Ｖｎ＋Ｖｍ−Ｖｎｍの算出過程の変数をＶａとする。
【００５９】
各区間における状態ｎの波動関数の絶対値｜ψｉｎ｜は、既存の分子シミュレーションによって予め算出され、パソコンの記憶装置に格納されている。また、各区間における状態ｍの波動関数の絶対値｜ψｉｍ｜は、既存の分子シミュレーションによって予め算出され、パソコンの記憶装置に格納されている。
【００６０】
サブルーチンが起動すると、Ｖａ、重なり領域の体積Ｖｎｍに０をそれぞれ代入して初期化する（Ｓ４０）。また、ｉに０を代入して初期化する（Ｓ４１）。
【００６１】
初期化後、ループ毎に、ｉをインクリメントする（Ｓ４２）。区間ｉにおける状態ｎの波動関数の絶対値｜ψｉｎ｜を用いて、式（４）によって区間ｉにおける状態ｎの波動関数を算出し、その算出値をφｎに代入する（Ｓ４３）。また、区間ｉにおける状態ｍの波動関数の絶対値｜ψｉｍ｜を用いて、式（４）によって区間ｉにおける状態ｍの波動関数を算出し、その算出値をφｍに代入する（Ｓ４４）。さらに、式（３）に基づいて、その算出したφｎとφｍを用いて｜φｎ−φｍ｜×Ｖｉを算出し、その算出値をＶａの前回値に加算し、その加算値をＶａに代入する（Ｓ４５）。また、式（５）に基づいて、その算出したφｎとφｍを用いて（√（φｎ×φｍ））×Ｖｉを算出し、その算出値をＶｎｍの前回値に加算し、その加算値をＶｎｍに代入する（Ｓ４６）。次に、ｉがＭより小さいか否かを判定する（Ｓ４７）。Ｓ４７にてｉがＭより小さいと判定した場合、Ｓ４２の処理に戻り、次の区間についての処理を行う。
【００６２】
一方、Ｓ４７にてｉがＭに等しくなったと判定した場合、重なり領域の体積Ｖｎｍ及びＶａを算出したので、ＶａをＶｎ＋ＶＭ−Ｖｎｍに代入する（Ｓ４８）。この算出されたＶｎｍ及びＶｎ＋ＶＭ−Ｖｎｍは、メインルーチンで使用される。
【００６３】
この分子シミュレーション装置１によれば、２つの異なる波動関数におけるキャリアの乗り移り易さ、空間上の広さ及びキャリアが存在する確率の各評価指標を加味した演算式を構築し、シミュレーションすることにより、アモルファス体の電気伝導度を精度良く推定することができる。
【００６４】
特に、分子シミュレーション装置１では、２つの異なる波動関数のエネルギレベルの差及び重なり領域の体積の各評価指標をその演算式に加味することにより、２つの異なる波動関数におけるキャリアの乗り移り易さを高精度に評価することができる。また、分子シミュレーション装置１では、２つの異なる波動関数の各体積の和から重なり領域の体積を引いた体積をその演算式に加味することにより、２つの異なる波動関数における空間上の広さを高精度に評価することができる。
【００６５】
さらに、分子シミュレーション装置１では、電気伝導度（代用評価値）を算出するときに、２つの異なる波動関数が共にフェルミエネルギより高い固有エネルギを持つ場合には電子による電気伝導を加味し、２つの異なる波動関数が共にフェルミエネルギより低い固有エネルギを持つ場合にはホールによる電気伝導を加味し、それ以外の場合には加味しないことにより、より高精度にアモルファス体の電気伝導度を推定することができる。
【００６６】
また、分子シミュレーション装置１では、アモルファス体を多数の断片に分割し、各断片における波動関数の絶対値が規定値より大きいか否か（その断片に波動関数に存在するか否か）によって波動関数に関する各体積を算出することにより、波動関数に関する各体積を簡単に算出することができる。また、分子シミュレーション装置１では、波動関数を頭打ち加工した関数とし、波動関数を連続的な値として波動関数に関する各体積を算出することにより、波動関数に関する各体積を高精度に算出することができる。
【００６７】
以上、本発明に係る実施の形態について説明したが、本発明は上記実施の形態に限定されることなく様々な形態で実施される。
【００６８】
例えば、本実施の形態ではアモルファス体の電気伝導度推定処理用のプログラムを実行することによってパソコン上に構成される分子シミュレーション装置に適用したが、分子シミュレーション装置の専用機として構成してもよいし、あるいは、ＣＤ−ＲＯＭなどの記憶媒体に格納されたプログラム、インタネットなどのネットワークを介して利用可能なプログラムなどに適用してもよい。
【００６９】
また、本実施の形態では式（１）において異なる波動関数として異なる２つの波動関数の組み合わせで評価を行う一例を示したが、異なる３つの波動関数の組み合わせ、異なる４つの波動関数の組み合わせなどで評価を行ってもよい。
【００７０】
また、本実施の形態では式（１）において評価指標１〜４による各評価を反映したファクタの一例を示したが、評価指標１〜４による各評価を反映したファクタについては他のものでもよい。
【００７１】
また、本実施の形態では波動関数に関する体積を求める手法の二例を示したが、他の手法によって波動関数の体積を求めてもよい。
【符号の説明】
【００７２】
１…分子シミュレーション装置、１０…パーソナルコンピュータ本体、１１…ディスプレイ、１２…キーボード
【特許請求の範囲】
【請求項１】
アモルファスの電気伝導度を推定する分子シミュレーション装置であって、
アモルファスに存在する異なる波動関数のキャリアの乗り移り易さ、空間上の広さ及びキャリアが存在する確率をそれぞれ評価し、当該３つの評価に基づいて電気伝導度を推定することを特徴とする分子シミュレーション装置。
【請求項２】
異なる波動関数のキャリアの乗り移り易さを、異なる波動関数のエネルギレベルの差及び異なる波動関数の重なり領域の体積に基づいて評価することを特徴とする請求項１に記載する分子シミュレーション装置。
【請求項３】
異なる波動関数の空間上の広さを、異なる波動関数の各体積の和から重なり領域の体積を引いた体積に基づいて評価することを特徴とする請求項１又は請求項２に記載する分子シミュレーション装置。
【請求項１】
アモルファスの電気伝導度を推定する分子シミュレーション装置であって、
アモルファスに存在する異なる波動関数のキャリアの乗り移り易さ、空間上の広さ及びキャリアが存在する確率をそれぞれ評価し、当該３つの評価に基づいて電気伝導度を推定することを特徴とする分子シミュレーション装置。
【請求項２】
異なる波動関数のキャリアの乗り移り易さを、異なる波動関数のエネルギレベルの差及び異なる波動関数の重なり領域の体積に基づいて評価することを特徴とする請求項１に記載する分子シミュレーション装置。
【請求項３】
異なる波動関数の空間上の広さを、異なる波動関数の各体積の和から重なり領域の体積を引いた体積に基づいて評価することを特徴とする請求項１又は請求項２に記載する分子シミュレーション装置。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【公開番号】特開２０１０−１７７２６０（Ｐ２０１０−１７７２６０Ａ）
【公開日】平成２２年８月１２日（２０１０．８．１２）
【国際特許分類】
【出願番号】特願２００９−１５４４２（Ｐ２００９−１５４４２）
【出願日】平成２１年１月２７日（２００９．１．２７）
【出願人】（０００００３２０７）トヨタ自動車株式会社 (59,920)
【公開日】平成２２年８月１２日（２０１０．８．１２）
【国際特許分類】
【出願日】平成２１年１月２７日（２００９．１．２７）
【出願人】（０００００３２０７）トヨタ自動車株式会社 (59,920)
[ Back to top ]