(57)【要約】      （修正有）【課題】第一原理計算に基づく計算結果を少ない個数のＴＢパラメータを用いたＴＢ法で再現するとともに、用いたＴＢパラメータのtransferabilityを確保する。【解決手段】強結合（ＴＢ）法により固体のバンド構造および全エネルギーを計算するＴＢ計算手段２と、第１の構造を有する固体について、ＴＢ計算手段によるバンド構造が第一原理計算によるバンド構造を再現するように、ＴＢ法に含まれる第１のパラメータを決定する第１パラメータ決定手段３と、第２の構造を有する固体について、決定した第１のパラメータを用いてＴＢ計算手段が計算したバンド構造が第一原理計算によるバンド構造を再現するように、ＴＢ法に含まれる第２のパラメータを決定する第２パラメータ決定手段４と、を備えている。【選択図】図１
2016171207.tif 000002
【特許請求の範囲】
【請求項１】
強結合（ＴＢ：Tight-Binding）法により固体のバンド構造および全エネルギーを計算するＴＢ計算手段と、
第１の構造を有する固体について、前記ＴＢ計算手段によるバンド構造が第一原理計算によるバンド構造を再現するように、前記ＴＢ法に含まれる第１のパラメータを決定する第１パラメータ決定手段と、
第２の構造を有する前記固体について、前記決定した第１のパラメータを用いて前記ＴＢ計算手段が計算したバンド構造が前記第一原理計算によるバンド構造を再現するように、前記ＴＢ法に含まれる第２のパラメータを決定する第２パラメータ決定手段と、を備え、
前記第１のパラメータは、原子軌道の広がりおよび軌道エネルギーに関するパラメータであり、
前記第２のパラメータは、前記固体のハミルトニアンの対角項の値が前記固体の構造に応じて変化することに関するパラメータである、
ことを特徴とする、強結合パラメータ決定装置。
【請求項２】
前記第１の構造と前記第２の構造との間では、結晶構造、および／または、格子定数が異なる、
請求項１に記載の強結合パラメータ決定装置。
【請求項３】
前記第２のパラメータは、前記固体のハミルトニアンの対角項を実行的な原子密度の級数で表したときの係数、および／または、前記固体の構造に応じて変化する対角項のエネルギーシフト量を表す係数を含む、
請求項１または２に記載の強結合パラメータ決定装置。
【請求項４】
前記第２の構造を有する前記固体について、前記決定した第１のパラメータを用いて前記ＴＢ計算手段が計算した全エネルギーが前記第一原理計算による全エネルギーを再現するように、前記ＴＢ法に含まれる第３のパラメータを決定する第３パラメータ決定手段を備え、
前記第３のパラメータは、前記固体のオンサイトエネルギーが前記固体の構造に応じて変化することに関するパラメータである、
請求項１ないし３のいずれか１項に記載の強結合パラメータ決定装置。
【請求項５】
前記第１パラメータ決定手段は、前記第１のパラメータに含まれる軌道エネルギーの初期値として、前記第一原理計算により求められた前記原子軌道の部分状態密度（ｐＤＯＳ：partial Density Of State）の重心を用いる、
請求項１ないし４のいずれか１項に記載の強結合パラメータ決定装置。
【請求項６】
コンピュータが、第１の構造を有する固体について、強結合（ＴＢ：Tight-Binding）法によりバンド構造を計算するステップと、
前記第１の構造を有する固体について、前記ＴＢ法によるバンド構造が第一原理計算によるバンド構造を再現するように、前記ＴＢ法に含まれる第１のパラメータを決定するステップと、
第２の構造を有する前記固体について、前記決定した第１のパラメータを用いて前記ＴＢ法によるバンド構造を計算するステップと、
前記第２の構造を有する前記固体について、前記決定した第１のパラメータを用いた前記ＴＢ法によるバンド構造が前記第一原理計算によるバンド構造を再現するように、前記ＴＢ法に含まれる第２のパラメータを決定するステップと、を含み、
前記第１のパラメータは、原子軌道の広がりおよび軌道エネルギーに関するパラメータであり、
前記第２のパラメータは、前記固体のハミルトニアンの対角項の値が前記固体の構造に応じて変化することに関するパラメータである、
ことを特徴とする、強結合パラメータ決定方法。
【請求項７】
前記第１の構造と前記第２の構造との間では、結晶構造、および／または、格子定数が異なる、
請求項６に記載の強結合パラメータ決定方法。
【請求項８】
前記第２のパラメータは、前記固体のハミルトニアンの対角項を実効的な原子密度の級数で表したときの係数、および／または、前記固体の構造に応じて変化する対角項のエネルギーシフト量を表す係数を含む、
請求項６または７に記載の強結合パラメータ決定方法。
【請求項９】
前記第２の構造を有する前記固体について、前記決定した第１のパラメータを用いた前記ＴＢ法による全エネルギーが前記第一原理計算による全エネルギーを再現するように、前記ＴＢ法に含まれる第３のパラメータを決定するステップを含み、
前記第３のパラメータは、前記固体のオンサイトエネルギーが前記固体の構造に応じて変化することに関するパラメータである、
請求項６ないし８のいずれか１項に記載の強結合パラメータ決定方法。
【請求項１０】
前記第１のパラメータに含まれる軌道エネルギーの初期値として、前記第一原理計算により求められた前記原子軌道の部分状態密度（ｐＤＯＳ：partial Density Of State）の重心を用いる、
請求項６ないし９のいずれか１項に記載の強結合パラメータ決定方法。
【請求項１１】
第１の構造を有する固体について、強結合（ＴＢ：Tight-Binding）法によりバンド構造を計算する処理と、
前記第１の構造を有する固体について、前記ＴＢ法によるバンド構造が第一原理計算によるバンド構造を再現するように、前記ＴＢ法に含まれる第１のパラメータを決定する処理と、
第２の構造を有する前記固体について、前記決定した第１のパラメータを用いた前記ＴＢ法によりバンド構造を計算する処理と、
前記第２の構造を有する前記固体について、前記決定した第１のパラメータを用いた前記ＴＢ法によるバンド構造が前記第一原理計算によるバンド構造を再現するように、前記ＴＢ法に含まれる第２のパラメータを決定する処理と、をコンピュータに実行させ、
前記第１のパラメータは、原子軌道の広がりおよび軌道エネルギーに関するパラメータであり、
前記第２のパラメータは、前記固体のハミルトニアンの対角項の値が前記固体の構造に応じて変化することに関するパラメータである、
ことを特徴とする、プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、強結合（ＴＢ：Tight Binding）法により電子状態を記述するためのパラメータを算出する強結合パラメータ決定装置、強結合パラメータ決定方法、および、プログラムに関する。
【背景技術】
【０００２】
固体の電子状態（例えば、バンド構造）を解析する手法として、密度汎関数法に基づく第一原理電子状態計算（以下、「第一原理計算」という。）や、強結合（ＴＢ：Tight Binding）法等が知られている。
【０００３】
関連技術として、非特許文献１には、doubleζ型のスレータ型原子軌道（ＳＴＯ：Slater Type Orbital）を用いてＴＢパラメータを決定する技術が記載されている。
【０００４】
また、非特許文献２には、ハミルトニアンの対角項および非対角項、ならびに、オーバーラップ行列をＴＢパラメータで表し、ＴＢ法で計算したバンド構造および全エネルギーによって、第一原理計算によるバンド構造および全エネルギーを再現するように、ＴＢパラメータを決定する技術が記載されている。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】J. Cerda and F. Soria, "Accurate and Transferable Extended Hueckel-Type Tight-Binding Parameters," Phys. Rev. B61, 7965, (2000)
【非特許文献２】M. J. Mehl and D. A. Papaconstantopoulos, "Applications of a Tight-Binding Total-Energy Method for Transition and Noble Metals: Elastic Constants, Vacancies, and Surfaces of Monatomic Metals," Phys. Rev. B54, 4519 (1996)
【発明の概要】
【発明が解決しようとする課題】
【０００６】
以下の分析は、本発明者によってなされたものである。
【０００７】
非特許文献１に記載された技術によると、同一原子から成る固体であっても構造（例えば、結晶構造、格子定数等）を変化させた場合には、第一原理計算によるバンド構造をＴＢ法に基づいて十分に再現することができない。すなわち、非特許文献１に記載された技術によると、transferabilityを有するＴＢパラメータを求めることができないという問題がある。
【０００８】
また、非特許文献２に記載された技術では、ハミルトニアンの対角項および非対角項、ならびに、オーバーラップ行列がそれぞれ別個のＴＢパラメータで表されるため、多数のＴＢパラメータ（例えば、単一原子結晶についてs, p, d軌道を考慮した場合、６９個のパラメータ）が必要となる。かかる手法によると、導入されたパラメータの物理的意味が不明確となり、ＴＢパラメータのフィッテング作業も困難となる上、得られたパラメータのtransferabilityも保障されないという問題がある。
【０００９】
そこで、第一原理計算に基づく計算結果を少ない個数のＴＢパラメータを用いたＴＢ法で再現するとともに、ＴＢパラメータのtransferabilityを確保することが課題となる。本発明の目的は、かかる課題解決に寄与する強結合パラメータ決定装置、強結合パラメータ決定方法、および、プログラムを提供することにある。
【課題を解決するための手段】
【００１０】
本発明の第１の態様に係る強結合パラメータ決定装置は、強結合（ＴＢ：Tight-Binding）法により固体のバンド構造および全エネルギーを計算するＴＢ計算手段を備えている。また、前記強結合パラメータ決定装置は、第１の構造を有する固体について、前記ＴＢ計算手段によるバンド構造が第一原理計算によるバンド構造を再現するように、前記ＴＢ法に含まれる第１のパラメータを決定する第１パラメータ決定手段を備えている。さらに、前記強結合パラメータ決定装置は、第２の構造を有する前記固体について、前記決定した第１のパラメータを用いて前記ＴＢ計算手段が計算したバンド構造が前記第一原理計算によるバンド構造を再現するように、前記ＴＢ法に含まれる第２のパラメータを決定する第２パラメータ決定手段を備えている。ここで、前記第１のパラメータは、原子軌道の広がりおよび軌道エネルギーに関するパラメータである。また、前記第２のパラメータは、前記固体のハミルトニアンの対角項の値が前記固体の構造に応じて変化することに関するパラメータである。
【００１１】
本発明の第２の態様によると、コンピュータによる強結合パラメータ決定方法が提供される。強結合パラメータ決定方法は、第１の構造を有する固体について、強結合（ＴＢ：Tight-Binding）法によりバンド構造を計算するステップを含む。また、前記強結合パラメータ決定方法は、前記第１の構造を有する固体について、前記ＴＢ法によるバンド構造が第一原理計算によるバンド構造を再現するように、前記ＴＢ法に含まれる第１のパラメータを決定するステップを含む。さらに、前記強結合パラメータ決定方法は、第２の構造を有する前記固体について、前記決定した第１のパラメータを用いて前記ＴＢ法によるバンド構造を計算するステップを含む。また、前記強結合パラメータ決定方法は、前記第２の構造を有する前記固体について、前記決定した第１のパラメータを用いた前記ＴＢ法によるバンド構造が前記第一原理計算によるバンド構造を再現するように、前記ＴＢ法に含まれる第２のパラメータを決定するステップを含む。ここで、前記第１のパラメータは、原子軌道の広がりおよび軌道エネルギーに関するパラメータである。また、前記第２のパラメータは、前記固体のハミルトニアンの対角項の値が前記固体の構造に応じて変化することに関するパラメータである。
【００１２】
本発明の第３の態様によると、第１の構造を有する固体について、強結合（ＴＢ：Tight-Binding）法によりバンド構造を計算する処理をコンピュータに実行させるプログラムが提供される。また、前記プログラムは、前記第１の構造を有する固体について、前記ＴＢ法によるバンド構造が第一原理計算によるバンド構造を再現するように、前記ＴＢ法に含まれる第１のパラメータを決定する処理を前記コンピュータに実行させる。さらに、前記プログラムは、第２の構造を有する前記固体について、前記決定した第１のパラメータを用いた前記ＴＢ法によりバンド構造を計算する処理を前記コンピュータに実行させる。また、前記プログラムは、前記第２の構造を有する前記固体について、前記決定した第１のパラメータを用いた前記ＴＢ法によるバンド構造が前記第一原理計算によるバンド構造を再現するように、前記ＴＢ法に含まれる第２のパラメータを決定する処理を前記コンピュータに実行させる。ここで、前記第１のパラメータは、原子軌道の広がりおよび軌道エネルギーに関するパラメータである。また、前記第２のパラメータは、前記固体のハミルトニアンの対角項の値が前記固体の構造に応じて変化することに関するパラメータである。
【００１３】
なお、プログラムは、非一時的なコンピュータ可読記録媒体（non-transitory computer-readable storage medium）に記録されたプログラム製品として提供することもできる。
【発明の効果】
【００１４】
本発明に係る強結合パラメータ決定装置、強結合パラメータ決定方法、および、プログラムによると、第一原理計算に基づく計算結果を少ない個数のＴＢパラメータを用いたＴＢ法で再現するとともに、ＴＢパラメータのtransferabilityを確保することが可能となる。
【図面の簡単な説明】
【００１５】
【図１】本発明の一実施形態に係る強結合パラメータ決定装置の構成を一例として示すブロック図である。
【図２】ＴＥ−ＴＢ法の概要（パラメータとＬＤＡハミルトニアンとの対応）を示す図である。
【図３】第１の実施形態に係る強結合パラメータ決定装置により決定されるＴＢパラメータを例示する表である。
【図４】第１の実施形態に係る強結合パラメータ決定装置の構成を一例として示すブロック図である。
【図５】第１の実施形態に係る強結合パラメータ決定装置の動作を一例として示すブロック図である。
【図６】ステップＡ１で決定したパラメータによるバンド構造の再現性を示す計算例１（シリコンのダイヤモンド構造、格子定数５．４３Å）である。
【図７】ステップＡ２で決定したパラメータによるバンド構造の再現性を示す計算例２（シリコンのダイヤモンド構造、格子定数４．６１５５Å）である。
【図８】ステップＡ２で決定したパラメータによるバンド構造の再現性を示す計算例２（シリコンのダイヤモンド構造、格子定数４．９Å）である。
【図９】ステップＡ２で決定したパラメータによるバンド構造の再現性を示す計算例２（シリコンのダイヤモンド構造、格子定数５．１５８５Å）である。
【図１０】ステップＡ２で決定したパラメータによるバンド構造の再現性を示す計算例２（シリコンのダイヤモンド構造、格子定数５．４３Å）である。
【図１１】ステップＡ２で決定したパラメータによるバンド構造の再現性を示す計算例２（シリコンのダイヤモンド構造、格子定数５．７０１５Å）である。
【図１２】ステップＡ２で決定したパラメータによるバンド構造の再現性を示す計算例３（シリコンのＦＣＣ構造、格子定数３．７Å）である。
【図１３】ステップＡ２で決定したパラメータによるバンド構造の再現性を示す計算例３（シリコンのＦＣＣ構造、格子定数３．８Å）である。
【図１４】ステップＡ２で決定したパラメータによるバンド構造の再現性を示す計算例３（シリコンのＦＣＣ構造、格子定数３．９Å）である。
【図１５】ステップＡ２で決定したパラメータによるバンド構造の再現性を示す計算例３（シリコンのＳＣ構造、格子定数２．４Å）である。
【図１６】ステップＡ２で決定したパラメータによるバンド構造の再現性を示す計算例３（シリコンのＳＣ構造、格子定数２．５Å）である。
【図１７】ステップＡ２で決定したパラメータによるバンド構造の再現性を示す計算例３（シリコンのＳＣ構造、格子定数２．６Å）である。
【図１８】ステップＡ２で決定したパラメータによるバンドエネルギーの計算例５である。
【図１９】ステップＡ１ないしＡ４で決定したパラメータによる全エネルギーの計算例である。
【発明を実施するための形態】
【００１６】
はじめに、一実施形態の概要について説明する。なお、この概要に付記する図面参照符号は、専ら理解を助けるための例示であり、本発明を図示の態様に限定することを意図するものではない。
【００１７】
図１は、一実施形態に係る強結合パラメータ決定装置の構成を例示するブロック図である。図１を参照すると、強結合パラメータ決定装置１は、ＴＢ（Tight-Binding）計算手段２、第１パラメータ決定手段３、および、第２パラメータ決定手段４を備えている。
【００１８】
ＴＢ計算手段２は、強結合（ＴＢ：Tight-Binding）法により固体のバンド構造および全エネルギーを計算する。また、第１パラメータ決定手段３は、第１の構造（基準とする構造、例えば、シリコンであれば、安定構造に相当する格子定数のダイヤモンド構造等）を有する固体について、ＴＢ計算手段２によるバンド構造が第一原理計算によるバンド構造を再現するように、ＴＢ法に含まれる第１のパラメータを決定する。さらに、第２パラメータ決定手段４は、第２の構造（第１の構造とは、結晶構造、および／または、格子定数が異なるもの）を有する前記固体について、決定した第１のパラメータを用いてＴＢ計算手段２が計算したバンド構造が第一原理計算によるバンド構造を再現するように、ＴＢ法に含まれる第２のパラメータを決定する。ここで、第１のパラメータは、原子軌道の広がりおよび軌道エネルギーに関するパラメータである（例えば、ＳＴＯ波動関数に含まれる波動関数の広がりを表すパラメータ(ζa1,l, ζa2,l, θal)および軌道エネルギーεal等）。また、第２のパラメータは、前記固体のハミルトニアンの対角項の値が前記固体の構造（結晶構造、格子定数等）に応じて変化することに関するパラメータを考慮する。（例えば、固体のハミルトニアンの対角項を実効的な原子密度の級数で表したときの係数aal，bal, cal, dal等，及び結晶構造により異なる対角項のエネルギーシフト量を表した係数γal等）
【００１９】
かかる強結合パラメータ決定装置１によると、原子軌道関数に含まれるパラメータおよび軌道エネルギーを用いてハミルトニアンおよびオーバーラップ行列を求めることができ、さらに、ハミルトニアンの対角項には構造変化の影響を取り込むことが可能となる。また、強結合パラメータ決定装置１によると、固体の構造（例えば、結晶構造、格子定数等）を変化させた場合においても、第一原理計算によるバンド構造をＴＢ法に基づいて十分に再現することが可能となる。すなわち、強結合パラメータ決定装置１によると、第一原理計算に基づく計算結果を少ない個数のＴＢパラメータ（例えば、単一原子結晶について３８個）を用いたＴＢ法で再現することができると同時に、ＴＢパラメータのtransferabilityを確保することが可能となる。
【００２０】
＜実施形態１＞
次に、第１の実施形態に係る強結合パラメータ決定装置について、図面を参照して説明する。本発明の強結合パラメータ決定装置は、本発明者等の研究によって導き出されたＴＥ−ＴＢ（Total Energy associated Tight Binding）法に基づく。そこで、はじめに、ＴＥ−ＴＢ法の理論的な枠組みについて説明する。
【００２１】
［ＴＥ−ＴＢ法］
強結合（ＴＢ：Tight-Binding）法は、パラメータ（以下、「ＴＢパラメータ」という。）を用いることにより高速な量子力学計算を可能にする方法である。ＴＢ法を利用するには、物理的な素性が明確なＴＢパラメータを導入した上で、計算精度を保証するようにＴＢパラメータを決定する必要がある。
【００２２】
本発明者等は、局所密度汎関数法（ＬＤＡ：Local Density Approximation)のハミルトニアンをベースにしてＴＢパラメータを自動算出するアルゴリズムを開発した。すなわち、第一原理計算で求められたバンド構造および全エネルギーを参照し、これらを再現するようにＴＢパラメータを決定する。今回開発したアルゴリズムでは、局所的な周辺環境の変化に対応してＴＢパラメータを決定している。この算出方法は、バンド構造のみならず全エネルギーも同時に再現するものであり、かかる手法をＴＥ−ＴＢ法と称する。ＴＥ−ＴＢ法で算出したパラメータを用いることで、結晶中の空孔、格子間原子（interstitial atom）、不純物等の動的挙動、結晶転位、非晶質のＭＤ（Molecular Dynamics）等の量子力学的挙動に加えて、未知物質の構造・機能が精度良く計算されるものと期待される。
【００２３】
図２は、ＴＥ−ＴＢ法の概要を示す。図２に示すように、局所密度汎関数（ＬＤＡ：Local Density Approximation）法によると、全エネルギーELDATotalは、バンドエネルギーEband、格子クーロンエネルギーΦlattice、および、オンサイト（On-site）エネルギーΦon-siteの和として表される。以下、図２に示す全エネルギーELDATotalの表式を導出するとともに、全エネルギーELDATotalを構成する各項と、ＴＢパラメータとの関係について詳述する。
【００２４】
１．強結合（ＴＢ）法と局所密度汎関数（ＬＤＡ）法のハミルトニアン
（１．１）スレータ型軌道（ＳＴＯ：Slater-Type Orbital）によるＴＢ法ハミルトニアン
（１．１．１）ＳＴＯ波動関数
以下の式（１）のように、分子軌道(ＭＯ：Molecular Orbital)の波動関数を原子軌道φ(a)nlmの線形結合で与えるものとする。ここで、(n, l, m)は、それぞれ、主量子数n、方位量子数l、および、磁気量子数mを表す。
【００２５】
2016171207.tif 000003
（１）
【００２６】
式（１）における係数Cαa;nlmは、一般化固有値方程式Hψα = εαSψαを解くことにより決定される。原子軌道基底{φ(a)nlm}がＴＢ法の枠組みで規定されている場合には、ハミルトニアン行列Hとオーバーラップ行列Sを求めることができる。
【００２７】
ここでは、原子aの軌道として、一例として、スレータ型軌道（ＳＴＯ：Slater-Type Orbital）φSTO(a)nlm(r - Ra)を採用する。球面調和関数Ymlと動径方向関数R(a)nlから成るＳＴＯ波動関数は、φSTO(a)nlm(r - Ra) = R(a)nl(r) Yml(θ, φ)と表すことができる。
【００２８】
スクリーニングパラメータζを用いると、ＳＴＯ波動関数の動径方向関数を式（２）のように表すことができる。
【００２９】
2016171207.tif 000004
（２）
【００３０】
なお、式（２）の代わりに、２つのスクリーニングパラメータζ1, ζ2を用いてＳＴＯの動径方向関数を式（３）のように表すこともできる。
【００３１】
2016171207.tif 000005
（３）
【００３２】
以下では、ＳＴＯ波動関数の動径方向関数として、式（３）を用いる場合について説明する。式（３）で表されるdouble-ζタイプのＳＴＯ波動関数は、結晶のバンド構造を定量的に再現するのに適している。ただし、本発明における原子軌道関数は、式（３）の関数形に限定されるものではない。
【００３３】
ここで、係数C1, C2は、２種類の波動関数の混合の度合いと規格化で決まる。式（３）を用いたときの規格化条件は、式（４）によって表される。
【００３４】
2016171207.tif 000006
（４）
【００３５】
したがって、係数C1, C2は独立ではない。そこで、θal = C2/C1とパラメタライズし、係数C1, C2の代わりにパラメータθalを導入する。
【００３６】
（１．１．２）ＴＢ法ハミルトニアン
強結合（ＴＢ）法のハミルトニアンを拡張ヒュッケル近似に従って評価する。すなわち、ハミルトニアンの対角成分および非対角成分を、それぞれ、式（５）および式（６）のように表す。
【００３７】
2016171207.tif 000007
（５）
【００３８】
2016171207.tif 000008
（６）
【００３９】
式（５）、（６）において、l, l'は、それぞれ原子軌道の方位量子数、すなわち、s軌道（l = 0）、p軌道（l = 1）、d軌道（l = 2）を表す。また、式（６）におけるオーバーラップ行列Sal,bl'については、上述のＳＴＯ波動関数を用いて計算する。なお、式（６）のKwhはWolfsberg-Helmholzの定数を表し、原子間距離Rabの関数である。
【００４０】
（１．２）ＬＤＡによる全エネルギー
局所密度近似汎関数法（ＬＤＡ：Local Density Approximation）による全エネルギーELDATotalは、式（７）のように表される。
【００４１】
2016171207.tif 000009
（７）
【００４２】
式（７）の右辺第１項のεαはn番目のバンドαについての固有エネルギーを表し、和は占有バンドについて取る。EXCとμXC(r)は、ＬＤＡの表式における交換相関エネルギーおよび交換相関ポテンシャルである。n(r)は価電子帯（Valence Band）の電子密度、Zaは実効的原子番号（原子番号-コア電子の数）、和Σ{a≠b}は、原子同士のすべての組み合わせについて和を取ることを表す。
【００４３】
（１．２．１）バンドエネルギー
式（７）の第１項はバンドエネルギーEbandと呼ばれ、次の式（８）ように表される。
【００４４】
2016171207.tif 000010
（８）
【００４５】
（１．２．２）交換相互作用エネルギー
式（７）の右辺第３項および第４項は、全エネルギーのうちの交換相互作用の寄与を表す部分である。これらは式（８）のバンドエネルギーには含まれない。ここでは、これらをまとめて、式（９）のように交換相互作用エネルギーΦXCを定義する。
【００４６】
2016171207.tif 000011
（９）
【００４７】
交換相互作用エネルギーについては、原子個々に局所的な扱いをする。
【００４８】
（１．２．３）格子クーロンエネルギー
静電エネルギー項には、核子‐核子間相互作用、核子‐電子間相互作用、電子‐電子間相互作用が含まれる。このうち、核子‐電子間および電子‐電子間相互作用は、バンドエネルギーに含まれている。ただし、電子‐電子間相互作用はダブルカウントされているので、式（７）の第２項においてその分を差し引いている。式（７）のうちのバンドエネルギーに含まれない静電エネルギー項を足し合わせたものをクーロンエネルギーΦCoul.とすると、クーロンエネルギーΦCoul.は式（１０）となる。
【００４９】
2016171207.tif 000012
（１０）
【００５０】
ここで、価電子帯（バレンス）の電子密度を各原子の価電子帯（バレンス）電荷密度の和で表すと、式（１１）となる。
【００５１】
2016171207.tif 000013
（１１）
【００５２】
式（１１）を式（１０）に代入すると、クーロンエネルギーΦCoul.は式（１２）のように表される。
【００５３】
2016171207.tif 000014
2016171207.tif 000015
（１２）
【００５４】
ここで、式（１２）において、式（１３）で定義格子クーロンエネルギーΦlatticeを導入した。
【００５５】
2016171207.tif 000016
（１３）
【００５６】
ここで、格子クーロンエネルギーΦlatticeは、クーロン相互作用のうちの構造に起因するもののみで構成され、ある原子の周りの電子の電荷分布na(r)が他の原子と重なりがなく、系全体が電気的に中性である場合にはゼロとなる。一方、ある原子から他の原子への電荷移動がある場合には、格子クーロンエネルギーΦlatticeの寄与を正確に評価する必要がある。
【００５７】
（１．２．４）オンサイト（On-Site）エネルギー
上記定式化において、式（９）の交換相互作用エネルギーΦXC、および、式（１２）で表されるクーロンエネルギーΦCoul.の第２項は、原子固有の局所的なエネルギーである。これらをまとめて、式（１４）で表されるオンサイト（On-site）エネルギーΦon-siteとする。
【００５８】
2016171207.tif 000017
（１４）
【００５９】
オンサイトエネルギーΦon-siteの値は、原子の密度や電子密度na(r)等の局所的な環境のみで決まる。
【００６０】
（１．２．５）ＬＤＡでの全エネルギー（まとめ）
最終的に、ＬＤＡによる全エネルギー式（７）は、式（８）のバンドエネルギーEbandと、式（１３）の格子クーロンエネルギーΦlatticeと、式（１４）のオンサイトエネルギーΦon-siteの和として、式（１５）のように表すことができる。
【００６１】
2016171207.tif 000018
（１５）
【００６２】
２．ＴＥ−ＴＢ（Total Energy associated Tight Binding）法によるＴＢパラメータ決定手順
次に、ＴＥ−ＴＢ法によるパラメータの決定手順について説明する。局所密度汎関数法に基づく全エネルギーは式（１５）の右辺に含まれる３つのエネルギー項の和として表される。以下では、式（１５）の右辺の各項について、対応するＴＢパラメータを決定する手順を説明する。
【００６３】
（２．１）基準とする結晶構造に対するバンド構造を再現（図２のステップＡ１）
ステップＡ１は、非特許文献１に記載された手法と同様である。基準とする構造とは、例えば、シリコン（Silicon）であれば、その安定な結晶構造であるダイヤモンド構造（格子定数5.43Å）である。基準とする構造に対して、ＴＢ法プログラムと第一原理計算でそれぞれバンド構造を計算し、その結果を比較して両者の差が小さくなるように、ＳＴＯパラメータ(ζa1,l, ζa2,l, θal)と、基準とする結晶構造での軌道エネルギーεal(ρ0al)を決定する。ここで、aは結晶を構成する原子を示す。各パラメータは１原子当たり、軌道毎（l = s, p, d）の３つの値をとるため、波動関数について９つ(ζa1,l, ζa2,l, θal)、軌道エネルギーεalについて３つの合計１２個のパラメータを決定する。
【００６４】
（２．２）異なる格子定数、結晶構造のバンド構造を再現（図２のステップＡ２）
固体のダイナミクスを考える場合にはハミルトニアンの対角項が周辺の環境変化に影響を受け増減することを考慮しておく必要がある。このために、同一組成の結晶に関して、異なる格子定数、異なる結晶構造についてもバンド構造再現する必要がある。この点については、非特許文献２で用いられている考え方を採用した。すなわち、ハミルトニアンの対角項を局所的な周辺環境に依存する式（１６）ないし式（１８）で表し、５つのパラメータa, b, c, d, λ、さらにこれに結晶構造のトポロジーによる影響を考慮するパラメーラγを加えエネルギーが最適化されるように決定する。
【００６５】
2016171207.tif 000019
（１６）
【００６６】
ここで、ρalは実効的な原子の密度であり、原子aの周りで、次式（１７）によって与えられるものとする。
【００６７】
2016171207.tif 000020
（１７）
【００６８】
2016171207.tif 000021
（１８）
【００６９】
式（１７）の和はすべての近隣の原子bについて取る。lは、原子の軌道の方位量子数を表し、λは近隣の原子の寄与を調整するフィッティングパラメータである。FC(R)は、カットオフ関数で、一例として、原子単位（atomic unit）でL = 0.5, R0 = 14.0に設定することができる。なお、式（１６）のパラメータaは安定な結晶構造に対する計算で定まるため、実際に決定すべきパラメータはb, c, d,λ,γの５つである。これらのパラメータの決定方法については後述する。
【００７０】
ＳＴＯ波動関数と軌道エネルギーのＴＢパラメータ（１原子種当り１２個）については、ステップＡ１で決定した値をそのまま用いる。
【００７１】
次に、基準となる構造に加えて、異なる結晶構造を考慮する。さらに、各結晶構造ごとに、格子定数が異なるものを複数考慮する。格子定数の取り方は任意であるが、一例として、安定状態に対して±10％程度の範囲で数点を取ることができる。例えば、シリコン（Silicon）の場合であれば、ダイヤモンド構造で格子定数が異なるものを数点選び、さらに、ＦＣＣ（Face Centered Cubic）構造、および、ＳＣ（Simple Cubic）構造について、安定構造を含むようにしつつ、異なる格子定数を有する構造を数点選択する。
【００７２】
これらの結晶構造に対して、ＴＢ法プログラムと第一原理計算法のバンド構造を比較しその差の合計が小さくなるように、パラメータを決定する。その際、ステップＡ１で決定したＳＴＯ波動関数に含まれるパラメータ(ζa1,l, ζa2,l, θal)と軌道エネルギーεalをそのまま計算に使用する。また、ρal = ρ0alのときにεal(ρal)がεal(ρ0al)に一致するように、aalを決定する。したがって、式（１６）、（１７）の６つのパラメータλal, aal, bal, cal, dal, γalのうちの独立なパラメータはλal, bal, cal, dal, γalの５つになる。ここで得られたパラメータに基づいて、用意した構造毎にバンドエネルギーETE-TBbandを計算する。
【００７３】
ただし，γは実効的な原子密度だけでは表現できない第2の構造の第１の構造（基準とする構造）に対する局所電子密度のトポロジカルな相違を考慮するためのパラメータであり、したがって第一の構造に対しては考慮する必要はない。
【００７４】
なお、ステップＡ１、Ａ２において、バンド構造の最適化は完了し，後述のステップＡ４では、全エネルギーの最適化のみを行う。
【００７５】
（２．３）格子クーロンエネルギーの計算（図２のステップＡ３）
格子クーロンエネルギーΦlatticeは、慎重に決定する必要がある。なぜなら、格子クーロンエネルギーΦlatticeは一般にゼロではないものの、その値が非常に小さいからである。具体的には、ＭＯ（Molecular Orbital）を決定した後、局所電荷密度na(r)をＳＴＯ波動関数で計算し、これを用いて格子クーロンエネルギーΦlatticeを式（１９）（式（１３））に基づいて計算する。
【００７６】
2016171207.tif 000022
（１９）
【００７７】
すなわち、上記ステップＡ１で決定したパラメータを用いて電荷密度を計算し、計算した電荷密度を用いることにより、式（１９）の格子クーロンエネルギーΦlatticeを計算することができる。
【００７８】
２．４オンサイト（On-site）エネルギーの計算（図２のステップＡ４）
式（１４）で与えられるオンサイト（On-site）エネルギーΦon-siteは、自己エネルギーと交換相互作用エネルギーの和であり、ＴＢ法の枠組みでは計算することができないエネルギー項である。そこで、第一原理計算法とＴＢ法で全エネルギーを一致させるために、式（１６）と同様の形式で表現されるパラメータをもう一組導入する。
【００７９】
2016171207.tif 000023
（２０）
【００８０】
2016171207.tif 000024
（２１）
【００８１】
ここで、
2016171207.tif 000025
2016171207.tif 000026
2016171207.tif 000027
【００８２】
Φaon-siteは局所電子密度na(r)の汎関数である。式(２１)は、これを原子配置で決まる局所的な実効原子密度で表される項(Paの展開部分)と電子密度na(r)の結晶構造依存性示す項との一般的な級数展開で表したものである。
【００８３】
パラメータA, B, C, Dは、局所電子電荷密度na(r)の汎関数と考えることができ、原子毎に決定すべきものである。実効的な密度Paは原子と電子の密度の関数であり、ここでは、式（２２）で表されるものと仮定する。
【００８４】
2016171207.tif 000028
（２２）
【００８５】
電子密度の結晶構造に依存した分布については、ウィグナーサイツ半径の1/rsの展開式で考慮している。式（２１）は1/rs4の３次の項まで表しているが、より高次の項まで取ることもできる。
【００８６】
上記パラメータを、異なる結晶構造と格子定数についてＬＤＡで計算されたELDATotalを再現するように決定する。
【００８７】
以上の４つのステップＡ１〜Ａ４で決定されるＴＢパラメータを図３に例示する。なお、図３の表において、軌道レベルεal(ρ0al)の表式に含まれる原子密度ρ0alは、基準となる構造における原子密度ρalである。
【００８８】
以上の手順を踏んで得られたパラメータによると、異なる構造に対するバンド構造と全エネルギーが、ＴＥ−ＴＢと第一原理計算結果との間でよく一致することが発明者等によって確認されている。したがって、パラメータを求めるために採用した構造以外の構造を有する固体についても、得られたパラメータを用いてバンド構造を計算することができるものと期待される。実際、シリコン（Si）の場合には、ダイヤモンド構造、ＦＣＣ構造、および、ＳＣ構造を用いて決定したパラメータに基づいてＴＥ−ＴＢ法で計算したβ-tin構造のバンド構造と全エネルギーは、格子定数の大きさを変えた第一原理計算で得られるバンド構造と全エネルギーと良く一致することが確かめられている。さらに、ＴＥ−ＴＢ法と第一原理計算との間で、全エネルギーも良く一致することが確認されている。すなわち、ＴＥ−ＴＢ法によると、transferabilityを有するＴＢパラメータを求めることができる。
【００８９】
３．軌道エネルギーの初期値
第一原理計算のバンド構造を再現するためのＴＢパラメータとして、各元素の軌道エネルギーεal（l = s, p, d）がある。軌道エネルギーパラメータの初期値をどうように選ぶかは、バンド構造の再現性に大きな影響を及ぼす可能性がある。ここでは、軌道エネルギーεalの初期値として、第一原理計算で求めた系の局所状態密度ｐＤＯＳ（partial Density Of State）の重心に相当するエネルギーを採用するものとする。
【００９０】
具体的には、第一原理計算プログラムを用いて各元素aの軌道l毎に局所状態密度pDOSa,l(E)を求めておく。局所状態密度を用いると、重心に相当するエネルギーεa,lは次（２３）で定義される。
【００９１】
2016171207.tif 000029
（２３）
【００９２】
式（２３）において、積分範囲[Elowera,l, Euppera,l]は、各元素および軌道について物理的に重要な領域をとる。
【００９３】
このように軌道エネルギーεalの初期値を選ぶことにより、高精度のパラメータフィッティングが実現される。特に、化合物系のバンド構造を再現する際には、このように初期値を適切に選択することが重要となる。
【００９４】
４．パラメータを最適化する評価関数
第一原理計算結果とＴＥ−ＴＢ法での計算結果のずれを評価関数と呼ぶ所定の関数で表し、この関数の値が小さくなるように最適化し、パラメータを決定する。ここでは、評価関数を、一例として、式（２４）ないし（２６）とした。ただし、本発明の評価関数は、かかる式に限定されない。
【００９５】
2016171207.tif 000030
（２４）
【００９６】
2016171207.tif 000031
（２５）
【００９７】
2016171207.tif 000032
（２６）
【００９８】
式（２５）において、ETE-TBband,i, EFPband,iは、それぞれＴＥ−ＴＢ法および第一原理計算で求めたバンドエネルギーを表す。また、式（２６）において、ETE-TBtotal, EFPtotalは、それぞれＴＥ−ＴＢ法および第一原理計算で求めた全エネルギーを表す。さらに、式（２５）、（２６）において、Ndは比較する構造の数を表す。また、式（２５）において、Nkは比較するバンド図のk点の個数を表し、Nbは比較するバンドの本数を表す。
【００９９】
［構成］
次に、本実施形態の強結合パラメータ決定装置の構成について、図面を参照して説明する。図４は、本実施形態に係る強結合パラメータ決定装置１の構成を例示するブロック図である。図４を参照すると、強結合パラメータ決定装置１は、ＴＢ計算手段２、第１パラメータ決定手段３、第２パラメータ決定手段４、第３パラメータ決定手段５、および、第一原理計算ＤＢ（Database）６を備えている。
【０１００】
ＴＢ計算手段２は、上述のＴＥ−ＴＢ（Total Energy associated Tight Binding）法により固体のバンド構造および全エネルギーを計算する。
【０１０１】
第一原理計算ＤＢ６は、ＴＥ−ＴＢ法によりＴＢパラメータを決定する対象とされる固体について、局所密度汎関数（ＬＤＡ：Local Density Approximation）法に基づいて計算したバンド構造、全エネルギー、部分状態密度等の計算結果を保持するデータベースである。
【０１０２】
第１パラメータ決定手段３は、第１の構造を有する固体（例えば、シリコンであれば、格子定数5.43Åのダイヤモンド構造等）について、ＴＢ計算手段２によるバンド構造が第一原理計算によるバンド構造を再現するように、ＴＥ−ＴＢ法に含まれる第１のパラメータを決定する（図２のステップＡ１）。具体的には、第１パラメータ決定手段３は、ＴＢ計算手段２によって求めたバンド構造が、第一原理計算ＤＢ６に格納されたバンド構造を再現するように、例えば、式（２４）の右辺第１項から成る評価関数に基づいて、第１のパラメータを最適化する。なお、第１パラメータ決定手段３が用いる最適化手法には、公知の手法を用いることができる。
【０１０３】
ここで、第１のパラメータは、図３の表に示す原子軌道パラメータであり、原子軌道の広がりを表すＳＴＯ波動関数パラメータζa1,l, ζa2,l, θal、および、軌道エネルギーパラメータεalを含む。また、第１パラメータ決定手段３は、第１のパラメータに含まれる軌道エネルギーパラメータεalの初期値として、第一原理計算により求められた原子軌道の部分状態密度（ｐＤＯＳ：partial Density Of State）の重心（式（２３）参照）を用いる。
【０１０４】
第２パラメータ決定手段４は、第２の構造を有する前記固体について、第１パラメータ決定手段３により決定された第１のパラメータを用いてＴＢ計算手段２が計算したバンド構造が、第一原理計算によるバンド構造を再現するように、ＴＥ−ＴＢ法に含まれる第２のパラメータを決定する（図２のステップＡ２）。具体的には、第２パラメータ決定手段３は、第１パラメータ決定手段３により決定された第１のパラメータを用いてＴＢ計算手段２によって求めたバンド構造が、第一原理計算ＤＢ６に格納されたバンド構造を再現するように、例えば、式（２４）の右辺第１項から成る評価関数に基づいて、第２のパラメータを最適化する。
【０１０５】
ここで、第１の構造と第２の構造との間では、結晶構造、および／または、格子定数が異なる。具体的には、第２の構造として、例えば、シリコンであれば、第１の構造とは格子定数が異なるダイヤモンド構造、または、ＦＣＣ構造もしくはＳＣ構造等を用いることができる。また、第２のパラメータは、図３の表に示すバンドエネルギーパラメータであり、固体のハミルトニアンの対角項の値が前記固体の構造に応じて変化することに関するパラメータである。例えば、第２のパラメータは、固体のハミルトニアンの対角項を実行的な原子密度の級数で表したときの係数を含んでいてもよい。より具体的には、第２のパラメータは、式（１６）ないし（１８）に含まれるパラメータλal, bal, cal, dal, γalとしてもよい。
【０１０６】
第３パラメータ決定手段５は、第２の構造を有する前記固体について、第１パラメータ決定手段３により決定された第１のパラメータを用いてＴＢ計算手段２が計算した全エネルギーが、第一原理計算による全エネルギーを再現するように、ＴＢ法に含まれる第３のパラメータを決定する（図２のステップＡ４）。具体的には、第３パラメータ決定手段３は、第１パラメータ決定手段３により決定された第１のパラメータを用いてＴＢ計算手段２によって求めた全エネルギーが、第一原理計算ＤＢ６に格納された全エネルギーを再現するように、例えば、式（２４）の右辺第２項から成る評価関数に基づいて、第３のパラメータを最適化する。
【０１０７】
ここで、第３のパラメータは、図３の表に示すオンサイトエネルギーパラメータであり、固体のオンサイトエネルギーが固体の構造に応じて変化することに関するパラメータである。具体的には、第３のパラメータは、式（２１）および（２２）に含まれるパラメータA, B, C, D,およびΛ，Γ，Θとしてもよい。
【０１０８】
［動作］
次に、本実施形態の強結合パラメータ決定装置１（図４）の動作について、図面を参照して説明する。図５は、強結合パラメータ決定装置１の動作を例示するフロー図である。
【０１０９】
図５を参照すると、ＴＢ計算手段２は、ＴＥ−ＴＢ法により第１の構造を有する固体のバンド構造を計算する（ステップＳ１）。
【０１１０】
第１パラメータ決定手段３は、第１の構造を有する固体について、ＴＢ計算手段２によるバンド構造が第一原理計算によるバンド構造を再現できているか否かを判定する（ステップＳ２）。例えば、第１パラメータ決定手段３は、式（２４）の評価関数の右辺第１項が所定の閾値以下となった場合、ＴＥ−ＴＢ法によるバンド構造が第一原理計算によるバンド構造を再現したとみなしてもよい。
【０１１１】
第１パラメータ決定手段３はＴＢ計算手段２によるバンド構造が第一原理計算によるバンド構造を再現できていない場合（ステップＳ２のＮｏ）、第１のパラメータを調整する（ステップＳ３）。一方、第１パラメータ決定手段３はＴＢ計算手段２によるバンド構造が第一原理計算によるバンド構造を再現できている場合（ステップＳ２のＹｅｓ）、第１のパラメータをステップＳ１で使用した値に決定する（ステップＳ４）。
【０１１２】
次に、ＴＢ計算手段２は、ステップＳ４で第１パラメータ決定手段３により決定された第１のパラメータを用いて、ＴＥ−ＴＢ法により、第２の構造を有する前記固体のバンド構造を計算する（ステップＳ５）。
【０１１３】
第２パラメータ決定手段４は、第２の構造を有する前記固体について、第１パラメータ決定手段３により決定された第１のパラメータを用いてＴＢ計算手段２が計算したバンド構造が、第一原理計算によるバンド構造を再現できているか否かを判定する（ステップＳ６）。例えば、第２パラメータ決定手段４は、式（２４）の評価関数の右辺第１項が所定の閾値以下となった場合、ＴＥ−ＴＢ法によるバンド構造が第一原理計算によるバンド構造を再現したとみなしてもよい。
【０１１４】
第２パラメータ決定手段４は、第２の構造を有する固体について、ＴＢ計算手段２によるバンド構造が第一原理計算によるバンド構造を再現できていない場合（ステップＳ６のＮｏ）、第２のパラメータを調整する（ステップＳ７）。一方、第２パラメータ決定手段４はＴＢ計算手段２によるバンド構造が第一原理計算によるバンド構造を再現できている場合（ステップＳ６のＹｅｓ）、第２のパラメータをステップＳ５で使用した値に決定する（ステップＳ８）。
【０１１５】
次に、ＴＢ計算手段２は、ステップＳ４で第１パラメータ決定手段３により決定された第１のパラメータを用いて、ＴＥ−ＴＢ法により、第２の構造を有する前記固体のバンド構造を計算する（ステップＳ９）。
【０１１６】
第３パラメータ決定手段５は、第２の構造を有する前記固体について、第１パラメータ決定手段３により決定された第１のパラメータを用いてＴＢ計算手段２が計算した全エネルギーが、第一原理計算による全エネルギーを再現できているか否かを判定する（ステップＳ１０）。例えば、第３パラメータ決定手段５は、式（２４）の評価関数の右辺第２項が所定の閾値以下となった場合、ＴＥ−ＴＢ法による全エネルギーが第一原理計算による全エネルギーを再現したとみなしてもよい。
【０１１７】
第３パラメータ決定手段５は、第２の構造を有する固体について、ＴＢ計算手段２による全エネルギーが第一原理計算による全エネルギーを再現できていない場合（ステップＳ１０のＮｏ）、第３のパラメータを調整する（ステップＳ１１）。一方、第３パラメータ決定手段５はＴＢ計算手段２による全エネルギーが第一原理計算による全エネルギーを再現できている場合（ステップＳ１０のＹｅｓ）、第３のパラメータをステップＳ９で使用した値に決定する（ステップＳ１２）。
【０１１８】
なお、第２パラメータ決定手段４および第３パラメータ決定手段５は、上記のように別個にパラメータ最適化を行う代わりに、同時にパラメータ最適化を行ってもよい。例えば、式（２４）の評価関数が所定の閾値以下となった場合、ＴＥ−ＴＢ法によるバンド構造および全エネルギーが第一原理計算によるバンド構造および全エネルギーを再現したとみなすことにより、第２パラメータ決定手段４による第２のパラメータの最適化（ステップＳ５〜Ｓ８）と、第３パラメータ決定手段５による第３のパラメータの最適化（ステップＳ９〜Ｓ１２）を、同時に実行してもよい。
【０１１９】
［計算例］
次に、シリコン結晶を例として、上記実施形態に係る強結合パラメータ決定装置で決定したパラメータが様々な格子定数、結晶構造の結晶に対して、高い精度で第一原理計算を再現することを計算例に基づいて示す。
【０１２０】
（計算例１）
図６は、ステップＡ１で決定したパラメータによるバンド構造の再現性を示す計算例である。ここでは、第１の構造として、シリコンのダイヤモンド結晶の安定構造（格子定数５．４３Å）を考え、これついてステップＡ１で決定したパラメータが第一原理計算で得られたバンド構造を再現している様子を示したものである。ここで、ＴＥ−ＴＢ法の計算はプログラムＥＬＳＥＳ（Extra-Large Scale Electronic Structure calculation program）で行い、第一原理計算はプログラムＡＢＩＮＩＴで行った（以下の計算例２〜６も同様である。）。
【０１２１】
（計算例２）
図７〜図１１は、ステップＡ２で決定したパラメータによるバンド構造の再現性を示す計算例である。図７〜図１１において、第２の構造としてシリコンのダイヤモンド構造を採用し、格子定数をそれぞれ、４．６１５５Å、４．９Å、５．１５８５Å、５．４３Å、５．７０１５Åとした。
【０１２２】
（計算例３）
図１２〜図１４は、ステップＡ２で決定したパラメータによるバンド構造の再現性を示す計算例である。図１２〜図１４において、第２の構造としてシリコンのＦＣＣ構造を採用し、格子定数をそれぞれ、３．７Å、３．８Å、３．９Åとした。
【０１２３】
（計算例４）
図１５〜図１７は、ステップＡ２で決定したパラメータによるバンド構造の再現性を示す計算例である。図１５〜図１７において、第２の構造としてシリコンのＳＣ構造を採用し、格子定数をそれぞれ、２．４Å、２．５Å、２．６Åとした。
【０１２４】
（計算例５）
図１８は、ステップＡ２で決定したパラメータによるバンドエネルギーの再現性を示す計算例である。ここでは、第２の構造としてシリコンのダイヤモンド構造、ＦＣＣ構造、および、ＳＣ構造を採用し、各構造について格子定数を数点選び、ＴＥ−ＴＢ法（ＥＬＳＥＳ）と第一原理計算（ＡＢＩＮＩＴ）によるバンドエネルギーの計算結果を比較して示す。
【０１２５】
（計算例６）
図１９は、ステップＡ１ないしＡ４で決定したパラメータによるシリコンの全エネルギーの計算例である。図１９の左側はＳＣ構造、中央はＦＣＣ構造、右側はタイヤモンド構造の全エネルギーの計算結果を、ＴＥ−ＴＢ法（ＥＬＳＥＳ）と第一原理計算(ＡＢＩＮＩＴ)について比較して示す。
【０１２６】
本実施形態の強結合パラメータ決定装置１によると、原子の局所的な環境変化（格子定数、結晶構造多形）させた場合に対しても第一原理計算で得られるバンド構造、全エネルギーを同時に再現することが可能となる。また、強結合パラメータ決定装置１によると、ＴＢパラメータを自動的に算出可能となり、結晶構造、原子間距離の変化に対応してユニバーサルに電子状態と全エネルギーを求めることができる。
【０１２７】
非特許文献２では、式（１６）の形をハミルトニアンの対角項に用い、また、ハミルトニアンの非対角項とオーバーラップ行列を、各軌道の結合ごとに別の形式のパラメタリゼーションで表している。これらのパラメタリゼーションでは単一原子結晶についてs, p, d軌道を考慮した場合に、６９個ものパラメータが必要であり、かつ、transferabilityが保障されない。また、非特許文献２では、全エネルギーを考慮しておらず、非特許文献２に記載された技術を動的なシミュレーションに適用することは不可能である。
【０１２８】
一方、本実施形態の強結合パラメータ決定装置１によると、ＳＴＯ波動関数と軌道エネルギーで、ハミルトニアン／オーバーラップ行列が計算でき、さらに対角項には構造変化の影響を取り込んでいる。すなわち、本実施形態の手法では、パラメータの物理的素性が明確であり、同一原子を含む異なる物質での計算にも用いることが可能である（すなわち、ＴＢパラメータはtransferabilityを有する）等、非特許文献２に記載された技術と比較して実用性の高いパラメタリゼーションとなっている。また、ＴＥ−ＴＢ法に必要なパラメータ数は、例えば単一原子結晶について３８個と少ない数で最適化を実現している。
【０１２９】
また、本実施形態の強結合パラメータ決定装置１では、全エネルギーがＬＤＡの表式に基づいてパラメタリゼーションしているため、ダイナミックな構造変化を伴う量子論的分子動力学計算において各瞬間の全エネルギーを、レファレンスとしたＬＤＡと同様の精度で求めることが可能である。
【０１３０】
さらに、強結合パラメータ決定装置１によって算出したＴＢパラメータを用いた計算は、第一原理計算と比較して圧倒的に高速な計算方法であり、第一原理計算と比較して超大規模かつ長時間の計算を実行することができる。すなわち、かかるＴＢパラメータを用いた計算によると、これまで困難であった固体中の原子拡散等の遅い動的性質に対する量子力学的計算が可能となり、新材料の設計開発に新しい手段を提供することが可能となる。
【０１３１】
なお、本発明において、下記の形態が可能である。
［形態１］
上記第１の視点に係る強結合パラメータ決定方法のとおりである。
［形態２］
前記強結合パラメータ決定方法において、前記第１の構造と前記第２の構造との間では、結晶構造、および／または、格子定数が異なっていてもよい。
［形態３］
前記強結合パラメータ決定方法において、前記第２のパラメータは、前記固体のハミルトニアンの対角項を実効的な原子密度の級数で表したときの係数、および／または、前記固体の構造に応じて変化する対角項のエネルギーシフト量を表す係数を含んでいてもよい。
［形態４］
前記強結合パラメータ決定方法は、前記第２の構造を有する前記固体について、前記決定した第１のパラメータを用いた前記ＴＢ法による全エネルギーが前記第一原理計算による全エネルギーを再現するように、前記ＴＢ法に含まれる第３のパラメータを決定するステップを含んでいてもよい。また、前記第３のパラメータは、前記固体のオンサイトエネルギーが前記固体の構造に応じて変化することに関するパラメータであってもよい。
［形態５］
前記強結合パラメータ決定方法において、前記第１のパラメータに含まれる軌道エネルギーの初期値として、前記第一原理計算により求められた前記原子軌道の部分状態密度（ｐＤＯＳ：partial Density Of State）の重心を用いてもよい。
［形態６］
上記第２の態様に係る強結合パラメータ決定装置のとおりである。
［形態７］
前記強結合パラメータ決定装置において、前記第１の構造と前記第２の構造との間では、結晶構造、および／または、格子定数が異なっていてもよい。
［形態８］
前記強結合パラメータ決定装置において、前記第２のパラメータは、前記固体のハミルトニアンの対角項を実行的な原子密度の級数で表したときの係数、および／または、前記固体の構造に応じて変化する対角項のエネルギーシフト量を表す係数を含んでいてもよい。
［形態９］
前記強結合パラメータ決定装置は、前記第２の構造を有する前記固体について、前記決定した第１のパラメータを用いて前記ＴＢ計算手段が計算した全エネルギーが前記第一原理計算による全エネルギーを再現するように、前記ＴＢ法に含まれる第３のパラメータを決定する第３パラメータ決定手段を備えていてもよい。また、前記第３のパラメータは、前記固体のオンサイトエネルギーが前記固体の構造に応じて変化することに関するパラメータであってもよい。
［形態１０］
前記強結合パラメータ決定装置において、前記第１パラメータ決定手段は、前記第１のパラメータに含まれる軌道エネルギーの初期値として、前記第一原理計算により求められた前記原子軌道の部分状態密度（ｐＤＯＳ：partial Density Of State）の重心を用いてもよい。
［形態１１］
上記第３の態様に係るプログラムのとおりである。
［形態１２］
前記プログラムにおいて、前記第１の構造と前記第２の構造との間では、結晶構造、および／または、格子定数が異なっていてもよい。
［形態１３］
前記プログラムにおいて、前記第２のパラメータは、前記固体のハミルトニアンの対角項を実効的な原子密度の級数で表したときの係数を含んでいてもよい。
［形態１４］
前記プログラムは、前記第２の構造を有する前記固体について、前記決定した第１のパラメータを用いた前記ＴＢ法による全エネルギーが前記第一原理計算による全エネルギーを再現するように、前記ＴＢ法に含まれる第３のパラメータを決定する処理を、前記コンピュータに実行させてもよい。また、前記第３のパラメータは、前記固体のオンサイトエネルギーが前記固体の構造に応じて変化することに関するパラメータであってもよい。
［形態１５］
前記プログラムにおいて、前記第１のパラメータに含まれる軌道エネルギーの初期値として、前記第一原理計算により求められた前記原子軌道に対する部分状態密度（ｐＤＯＳ：partial Density Of State）の重心を用いてもよい。
【０１３２】
なお、本発明の全開示（請求の範囲を含む）の枠内において、さらにその基本的技術思想に基づいて、実施形態の変更・調整が可能である。また、本発明の全開示の枠内において種々の開示要素（各請求項の各要素、各実施形態の各要素、各図面の各要素等を含む）の多様な組み合わせ、ないし、選択が可能である。すなわち、本発明は、請求の範囲を含む全開示、技術的思想にしたがって当業者であればなし得るであろう各種変形、修正を含むことは勿論である。特に、本書に記載した数値範囲については、当該範囲内に含まれる任意の数値ないし小範囲が、別段の記載のない場合でも具体的に記載されているものと解釈されるべきである。
【符号の説明】
【０１３３】
１    強結合パラメータ決定装置
２    ＴＢ計算手段
３    第１パラメータ決定手段
４    第２パラメータ決定手段
５    第３パラメータ決定手段
６    第一原理計算ＤＢ
【図１】
2016171207.tif 000033
【図２】
2016171207.tif 000034
【図３】
2016171207.tif 000035
【図４】
2016171207.tif 000036
【図５】
2016171207.tif 000037
【図６】
2016171207.tif 000038
【図７】
2016171207.tif 000039
【図８】
2016171207.tif 000040
【図９】
2016171207.tif 000041
【図１０】
2016171207.tif 000042
【図１１】
2016171207.tif 000043
【図１２】
2016171207.tif 000044
【図１３】
2016171207.tif 000045
【図１４】
2016171207.tif 000046
【図１５】
2016171207.tif 000047
【図１６】
2016171207.tif 000048
【図１７】
2016171207.tif 000049
【図１８】
2016171207.tif 000050
【図１９】
2016171207.tif 000051