説明

電子状態計算方法、電子状態計算プログラム及び電子状態計算装置

【課題】 TDDFT法による物質の電子状態の計算を効率的に行う。
【解決手段】 電子状態計算方法は、TDDFT法によって物質の電子状態を計算する方法であって、物質の電子状態の初期状態を示す情報等を入力する入力ステップ(S01)と、時間ステップnにおいて、時間ステップn−1から時間ステップn−I−1までのハートレー項V(n−1)〜V(n−I−1)を用いて、ハートレー項V(n)の初期値Vinit(n)を算出する初期値算出ステップ(S03)と、初期値Vinit(n)からハートレー項V(n)を算出するハートレー項算出ステップ(S05)と、算出対象の1つ前の時間ステップでの物質の電子状態を示す情報とハートレー項V(n−1)とから、物質の電子状態を算出する電子状態算出ステップ(S02)と、を含む。

【発明の詳細な説明】
【技術分野】
【0001】
本発明は、TDDFT法によって物質の電子状態を計算する電子状態計算方法、電子状態計算プログラム及び電子状態計算装置に関する。
【背景技術】
【0002】
電子状態によって定まる時間依存の物質の物性量を算出する方法として、その物性量の要因となっている電子状態を記述する基本方程式を用い、外部刺激に対する応答の時間発展を求める方法が広く用いられている。この種の算出方法では、シミュレーションを行なう空間内に外部刺激を与え、その電子状態の揺らぎを記述する方法として、例えば、時間依存密度汎関数法(以下、TDDFT(Time-Dependent Density Functional Theory)法と呼ぶ)がある(例えば、非特許文献1参照)。
【0003】
この計算手法の大きな特徴として、励起状態の情報を提供できる特徴をもっていることが知られている。このため、物質の物性量の中でも励起情報が必要な発光・吸収に関するスペクトルあるいは光学定数の算出には、重要な手法として使われている。
【0004】
この種の方法では電子状態は一般に波動関数で表され、Kohn-Sham方程式を用いセルフコンシステントな状態を求めておき、これにある時間において微小な摂動として外部刺激を与え、時間依存のKohn-Sham方程式を用い、その応答をその解である波動関数の変化として、時間を追って追跡する。一般に、これらの時間変化を記述する方法として、空間的、時間的に差分化して計算する方法が採られている。
【非特許文献1】K. Yabana and G.F. Bertsch,“Time-dependent local-density approximation in real time”, Phys. Rev B54,4484(1996)
【発明の開示】
【発明が解決しようとする課題】
【0005】
この時間発展の各ステップで波動関数を求めるために、その時間ステップより以前のステップでの電子状態を初期値として用い安定な状態となったところで、次の時間ステップへ進めることにより、時間発展計算が行なわれる。この安定な状態とするために、特にハートレー項などのハミルトニアンの各項を、数値的に繰り返し計算することが行われている。計算対象となる物質によっては、この繰り返し計算の回数が多くなるため、計算時間が非常に長くなるという問題があった。繰り返し計算の回数は、物質の分子形状に大きく影響され、特に高分子では鎖長が長くなればなるほど、繰り返し回数が増えていた。
【0006】
本発明は、以上の問題点を解決するためになされたものであり、TDDFT法による物質の電子状態の計算を効率的に行うことができる電子状態計算方法、電子状態計算プログラム及び電子状態計算装置を提供することを目的とする。
【課題を解決するための手段】
【0007】
上記の目的を達成するために、本発明に係る電子状態計算方法は、コンピュータにより、TDDFT法によって物質の電子状態を計算する電子状態計算方法であって、物質の電子状態の初期状態を示す情報を含むTDDFT法によって物質の電子状態を計算するためのパラメータを入力して、TDDFT法における時間ステップj毎の、物質の電子状態のハミルトニアンH(j)におけるハートレー項V(j)を含む当該物質の電子状態を示す情報を記憶する記憶手段に記憶させる入力ステップと、時間ステップnにおいて、時間ステップn−1から時間ステップn−I−1までのハートレー項V(n−1)〜V(n−I−1)を記憶手段から読み出して、予め設定されるI個の係数α〜αを用いて、時間ステップnにおけるハートレー項V(n)の初期値Vinit(n)を予め記憶した式
【数1】

により演算をして算出する初期値算出ステップと、初期値算出ステップにおいて算出されたハートレー項V(n)の初期値Vinit(n)から、予め記憶したハートレー項算出ルールに基づき差分法により当該ハートレー項V(n)を算出すると共に、算出したハートレー項V(n)を記憶手段に記憶させるハートレー項算出ステップと、ハートレー項算出ステップにおいて算出された物質のハートレー項V(n)及び記憶手段により記憶された情報から、予め記憶した電子状態算出ルールに基づき物質の電子状態を算出すると共に、算出した電子状態を示す情報を記憶手段に記憶させる電子状態算出ステップと、を含むことを特徴とする(但し、Iは予め定められる1以上の整数であり、n,jは時間ステップを示すインデックスである)。
【0008】
本発明に係る電子状態計算方法では、時間ステップnにおけるハートレー項V(n)の初期値Vinit(n)が、当該時間ステップの1つ前の時間ステップn−1から時間ステップn−I−1までのハートレー項V(n−1)〜V(n−I−1)に基づいて算出される。このように算出された初期値Vinit(n)は算出するハートレー項V(n)の値に近くなるため、ハートレー項V(n)の算出を従来と比べてより少ない繰り返し計算で行うことが可能になる。これにより、本発明に係る電子状態計算方法によれば、TDDFT法による物質の電子状態の計算を効率的に行うことができる。
【0009】
電子状態計算方法は、電子状態算出ステップにおいて算出された物質の電子状態を示す情報から、予め記憶した物性量算出ルールに基づき当該物質の物性量を算出する物性量算出ステップと、物性量算出ステップにおいて算出された物質の物性量を示す情報を出力する出力ステップと、を更に備えることが望ましい。この構成によれば、効率的に物質の物性量を算出することができる。
【0010】
入力ステップにおいて、物質の電子状態を計算するためのパラメータとして、物質の電子状態の初期状態を示す情報に与えられる外部刺激を示す情報を入力することが望ましい。この構成によれば、時間変化に応じた物質の電子状態あるいは物性量を適切に算出することができる。
【0011】
初期値算出において、時間ステップnにおけるハートレー項V(n)の初期値Vinit(n)を予め記憶した式
【数2】

により演算をして算出する、ことが望ましい。この構成によれば、簡易かつ確実に物質の電子状態の計算を効率的に行うことができる。
【0012】
電子状態計算方法は、記憶手段からハートレー項V(j)を含む物質の電子状態を示す情報を読み出して、当該ハートレー項V(j)を含む当該物質の電子状態を示す情報から予め記憶した係数値算出ルールに基づき、係数αを算出して当該記憶手段に記憶させる係数算出ステップを更に含み、初期値算出ステップにおいて、係数αを記憶手段から読み出して、初期値Vinit(n)の算出に用いるものとして設定する、ことが望ましい。この構成によれば、ハートレー項V(n)の初期値Vinit(n)を、より算出するハートレー項V(n)の値に近づけることができ、ハートレー項V(n)の算出を従来と比べて更に少ない繰り返し計算で行うことが可能になる。
【0013】
ところで、本発明は、上記のように電子状態計算方法の発明として記述できる他に、以下のように電子状態計算プログラム及び電子状態計算装置の発明としても記述することができる。これはカテゴリが異なるだけで、実質的に同一の発明であり、同様の作用及び効果を奏する。
【0014】
即ち、本発明に係る電子状態計算プログラムは、コンピュータに、TDDFT法によって物質の電子状態を計算させる電子状態計算プログラムであって、TDDFT法における時間ステップj毎の、物質の電子状態のハミルトニアンH(j)におけるハートレー項V(j)を含む当該物質の電子状態を示す情報を記憶する記憶手段と、物質の電子状態の初期状態を示す情報を含むTDDFT法によって物質の電子状態を計算するためのパラメータを入力して、記憶手段に記憶させる入力手段と、時間ステップnにおいて、時間ステップn−1から時間ステップn−I−1までのハートレー項V(n−1)〜V(n−I−1)を記憶手段から読み出して、予め設定されるI個の係数α〜αを用いて、時間ステップnにおけるハートレー項V(n)の初期値Vinit(n)を予め記憶した式
【数3】

により演算をして算出する初期値算出手段と、初期値算出手段により算出されたハートレー項V(n)の初期値Vinit(n)から、予め記憶したハートレー項算出ルールに基づき差分法により当該ハートレー項V(n)を算出すると共に、算出したハートレー項V(n)を記憶手段に記憶させるハートレー項算出手段と、ハートレー項算出手段により算出された物質のハートレー項V(n)及び記憶手段により記憶された情報から、予め記憶した電子状態算出ルールに基づき物質の電子状態を算出すると共に、算出した電子状態を示す情報を記憶手段に記憶させる電子状態算出手段と、をコンピュータに機能させることを特徴とする。
【0015】
本発明に係る電子状態計算装置は、TDDFT法によって物質の電子状態を計算する電子状態計算装置であって、TDDFT法における時間ステップj毎の、物質の電子状態のハミルトニアンH(j)におけるハートレー項V(j)を含む当該物質の電子状態を示す情報を記憶する記憶手段と、物質の電子状態の初期状態を示す情報を含むTDDFT法によって物質の電子状態を計算するためのパラメータを入力して、記憶手段に記憶させる入力手段と、時間ステップnにおいて、時間ステップn−1から時間ステップn−I−1までのハートレー項V(n−1)〜V(n−I−1)を記憶手段から読み出して、予め設定されるI個の係数α〜αを用いて、時間ステップnにおけるハートレー項V(n)の初期値Vinit(n)を予め記憶した式
【数4】

により演算をして算出する初期値算出手段と、初期値算出手段により算出されたハートレー項V(n)の初期値Vinit(n)から、予め記憶したハートレー項算出ルールに基づき差分法により当該ハートレー項V(n)を算出すると共に、算出したハートレー項V(n)を記憶手段に記憶させるハートレー項算出手段と、ハートレー項算出手段により算出された物質のハートレー項V(n)及び記憶手段により記憶された情報から、予め記憶した電子状態算出ルールに基づき物質の電子状態を算出すると共に、算出した電子状態を示す情報を記憶手段に記憶させる電子状態算出手段と、を備えることを特徴とする。
【発明の効果】
【0016】
本発明によれば、ハートレー項V(n)の算出を従来と比べてより少ない繰り返し計算で行うことが可能になるため、TDDFT法による物質の電子状態の計算を効率的に行うことができる。
【発明を実施するための最良の形態】
【0017】
以下、図面と共に本発明に係る電子状態計算方法、電子状態計算プログラム及び電子状態計算装置の好適な実施形態について詳細に説明する。なお、図面の説明においては同一要素には同一符号を付し、重複する説明を省略する。
【0018】
本実施形態に係る電子状態計算方法は、TDDFT法(TDDFTシミュレーション)によって物質の電子状態を計算するものである。即ち、電子状態計算方法は、時間ステップを設定して、当該時間ステップ毎の物質の電子状態(電子状態の時間発展)を求めるものである。計算された物質の電子状態は、後述するように当該物質の物理的あるいは化学的物性量の算出に用いられる。
【0019】
本実施形態において電子状態を示す波動関数ψ(n)は、以下の式のように算出対象の時間ステップn(n番目の時間ステップ)の前の時間ステップn−1の波動関数ψ(n−1)に対する時間発展演算子T(n−1)により算出される(なお、上記では、時間発展演算子を便宜上Tと記載しているが、実際には以下の式に記載されるようにハット記号(^)が付けられたTである。以降についても同様)。
【数5】

なお、nは整数であり時間ステップを示すインデックスである。ここで時間発展演算子T(n)は、以下の式により求められる。
【数6】

ここで、iは虚数を、H(n)は時間ステップnのハミルトニアンを、Δtは時間ステップの間隔をそれぞれ示している。上記の時間発展は、Δtに関して4次以上の展開を行なっても良いし、Suzuki-Trotter形式で演算を行なっても良い(参考文献:Masuo Suzuki, “General NonsymmetricHigher-Order Decomposition of Exponential Operators and Symplectic Integrator”、J.Phys.Soc. Jpn. 61, p3015(1992))。
【0020】
時間ステップnにおけるハミルトニアンH(n)は、以下の式のように時間ステップnにおける電子状態(例えば波動関数ψ(n))に依存している。
【数7】

ここで、ψ(n)は、時間ステップnの物質におけるi番目の電子状態を、{ψ(n)}は、その扱う電子状態一式をそれぞれ示している。例えば、擬ポテンシャルを使った密度汎関数法では、ハミルトニアンH(n)は、
【数8】

と表される。第1項は電子の運動エネルギー項であり、第2項はイオンと電子の相互作用ポテンシャルである。また、第3項がハートレー項であり、第4項は交換相関相互作用を表す。この場合は、時間ステップnのハミルトニアンH(n)は時間ステップnの電子密度(波動関数から算出可能)に依存している。
【0021】
時間発展演算の中で一番変化の大きいハートレー項V(n)は、実空間解析では次の式
【数9】

を空間差分により差分法を用いて連立方程式に変換し、共役勾配法等の繰り返し法を用いて求める。従来の方法では、この際、初期値Vinit(n)として、以下の式
【数10】

に示すように、1つ前の時間ステップn−1のハートレー項の値が用いられて、上述の繰り返し計算を行いV(n)が求められていた。
【0022】
引き続いて、本実施形態に係る電子状態計算方法が実行される電子状態計算装置10を説明する。図1に電子状態計算装置10を示す。電子状態計算装置10は、具体的には、ワークステーションやPC(Personal Computer)等のコンピュータである。電子状態計算装置10は、例えばCPU(Central ProcessingUnit)やメモリ等のハードウェアにより構成されており、これらの構成要素が動作することにより後述する電子状態計算装置10としての機能が発揮される。なお、本実施形態に係る電子状態計算方法をコンピュータに実行させるプログラムが電子状態計算装置10において実行されることにより本方法が行われてもよい。
【0023】
図1に示すように電子状態計算装置10は、記憶部11と、入力部12と、初期値算出部13と、ハートレー項算出部14と、電子状態算出部15と、物性量算出部16と、出力部17と、係数算出部18と、計算制御部19とを備えて構成される。また、電子状態計算装置10は、外部装置20と接続されており、外部装置20から情報が入力される。
【0024】
記憶部11は、物質の電子状態を計算するための情報を記憶する記憶手段である。記憶部11は、例えば、電子状態計算装置10が備えるメモリにより実現される。具体的には、物質を示す情報である原子配置やTDDFT法の計算に用いるためのパラメータ(例えば、時間発展のステップ数、ステップ間隔、汎関数の選択、計算空間の大きさ(差分格子間隔と格子点数))、及び物質の電子状態を示す情報等である。
【0025】
物質の電子状態を示す情報としては、例えば、波動関数及び電子密度等である。記憶部11に記憶される物質の電子状態を示す情報として、本実施形態に係る電子状態計算方法による算出の初期状態と、本実施形態に係る電子状態計算方法により算出される時間ステップ毎の電子状態を示す情報とを含む。時間ステップj毎の物質の電子状態を示す情報には、当該物質の電子状態のハミルトニアンH(j)におけるハートレー項V(j)を含む。特に、記憶部11は、TDDFT法における時間ステップj毎の、物質の電子状態のハミルトニアンH(j)におけるハートレー項V(j)を含む当該物質の電子状態を示す情報を記憶する。ここで、jは1以上の整数であり、時間ステップを示すインデックスである。なお、記憶部11は、電子状態を示す情報を格納する際には、何れの時間ステップにおける電子状態であるかが識別できるように、時間ステップを示す情報に対応付けておく。
【0026】
入力部12は、物質の電子状態の初期状態を示す情報を含むTDDFT法によって物質の電子状態を計算するためのパラメータを入力して、記憶部11に記憶させる入力手段である。入力部12に入力される物質の電子状態の初期状態としては、物質の電子状態を記述する基本方程式を満たすセルフコンシステントなもの(波動関数)とするのが望ましい。
【0027】
また、入力部12は、物質の電子状態を計算するためのパラメータとして、物質の電子状態の初期状態を示す情報に与えられる外部刺激を示す情報を入力することとしてもよい。外部刺激を示す情報が入力されると、外部刺激が与えられた物質の電子状態が計算される。ここで外部刺激は、計算対象の物質を構成する原子あるいは分子の位置を時間変化させるものである。入力される外部刺激は、摂動の範囲の刺激であり、物質を構成する原子あるいは分子の位置の時間変化が電子の運動に比べて十分に小さいと見なせるものとするのがよい。セルフコンシステントな初期状態に対して、上記の外部刺激を与えると、物質の時間変化に応じた電子状態、及び物性量が適切に計算できるからである。
【0028】
外部刺激は、セルフコンシステントな電子状態を示す波動関数ψに対して、時間ステップj=0にて、
【数11】

とすることで与えられる。ここで、ψ´は外部刺激が与えられた波動関数を、Eは十分小さい電場を、rは刺激を与える方向をそれぞれ示している。なお、参考文献K. Yabana and G.F. Bertsch, “Time-dependent local-densityapproximation in real time”, Phys. Rev B54, 4484(1996)では、解析系の大きさの逆数を指標として、これよりも十分小さな波数とすることが必要であると記述されている。
【数12】

これは全ての電子に均一な速度場を与えることに等しく、これによって生じる電場を扱うといっても良いので上記標識を用いている。
【0029】
より具体的には、原子単位系を用いて、rはx方向、y方向、z方向の何れかの座標である。ここで、原子単位系はauで表され、h/2π=m=e=1(但し、ここではhはプランク定数を表す)で求められる単位系である。ボーア半径を1とするので、1au=0.529177Å、エネルギーは1ハートレーを1とするため1au=27.12eVとなる。表式が簡便となるため量子力学計算では良く用いる単位系である。
【0030】
上記の情報の入力は、例えば、ユーザによる操作によって外部装置20から行われる。また、所定の演算を行って上記の情報(例えば、セルフコンシステントな電子状態を表す情報等)が生成されてもよい。入力部12は、入力した各情報を記憶部11に記憶させる。また、外部装置20から入力されたパラメータを予め電子状態計算装置10の記憶部11に記憶させておき、電子状態計算方法の実行時に記憶された情報が用いられてもよい。入力部12は、各情報を入力するとその旨を計算制御部19に通知する。
【0031】
初期値算出部13は、時間ステップnにおいて、ハートレー項V(n)の初期値Vinit(n)を算出する初期値算出手段である。また、初期値算出部13は、算出した初期値Vinit(n)を記憶部11に格納する。具体的には、初期値算出部13は、時間ステップn−1から時間ステップn−I−1までのハートレー項V(n−1)〜V(n−I−1)を記憶部11から読み出して、予め設定されるI個の係数α〜αを用いて、初期値Vinit(n)を予め記憶した式
【数13】

により演算をして算出する。ここで、Iは予め定められる1以上の整数である。即ち、初期値算出部13は、初期値Vinit(n)の算出に、1つ前の時間ステップn−1でのハートレー項V(n−1)の値だけでなく、2つ以上前の時間ステップでのハートレー項の値を用いる。ハートレー項を用いる時間ステップの数は、I+1である。
【0032】
特に、初期値算出部13は、初期値Vinit(n)の算出に、2つ前までの時間ステップでのハートレー項を用いることとするのが好ましい。即ち、初期値算出部13は、初期値Vinit(n)を予め記憶した式
【数14】

により演算をして算出する(上記の一般化した式(11)におけるI=1)。なお、上記の式(12)における係数αは、式(11)における係数αに対応する。これは、計算の容易さ、記憶容量、及び確実に本発明を実施する等の観点によるものである。
【0033】
上記のように、初期値Vinit(n)の算出に、2つ前までの過去の時間ステップでのハートレー項の値を用いることとすると、初期値Vinit(n)は、算出するハートレー項V(n)の値に近くなる。これは、図2に示すように、2つ前までの過去の時間ステップでのハートレー項の値を用いることによって、ハートレー項V(n)の変動割合(図2においてはα[V(n−1)−V(n−2)]が相当する)が考慮された初期値Vinit(n)となるためである。
【0034】
また、各[V(n−j)−V(n−j−1)]の係数αは、予め設定されておいてもよいし、後述するように係数算出部18により算出されてもよい。係数αは記憶部11に記憶されており、初期値算出部13は、係数αを記憶部11から読み出して、初期値Vinit(n)の算出に用いるものとして設定する。係数αは、計算に用いるハートレー項V(n)の数に応じて、適切に設定されることが望ましい。上記のように2つ前までの過去の時間ステップでのハートレー項の値を用いる場合には、例えば、係数αを1とするのが望ましい。また、計算機実験等によって、経験的に適切な値を設定するのが望ましい。
【0035】
ハートレー項算出部14は、初期値算出部13により算出されたハートレー項V(n)の初期値Vinit(n)から、予め記憶したハートレー項算出ルールに基づき差分法により当該ハートレー項V(n)を算出するハートレー項算出手段である。また、ハートレー項算出部14は、算出したハートレー項V(n)を記憶部11に記憶させる。ハートレー項算出部14は、上記の算出のための演算に応じたパラメータ等を記憶部11から読み出して、演算に用いる。
具体的には、ハートレー項算出部14は、記憶部11から算出するハートレー項V(n)の初期値Vinit(n)を読み出して、上述したように、ハートレー項V(n)を、次の式
【数15】

を空間差分により差分法を用いて連立方程式に変換し、共役勾配法等の繰り返し法を用いて求める。ハートレー項算出部14は、算出したハートレー項V(n)の値を記憶部11に記憶させる。
【0036】
電子状態算出部15は、算出対象の1つ前の時間ステップでの物質の電子状態を示す情報を記憶部11から読み出して、当該物質の電子状態を示す情報とハートレー項算出部14によって算出されたハートレー項V(n−1)とから、予め記憶した電子状態算出ルールに基づき物質の電子状態を算出する電子状態算出手段である。また、電子状態算出部15は、算出した電子状態を示す情報を記憶部11に記憶させる。また、電子状態算出部15は、上記の算出のための演算に応じたパラメータ等を記憶部11から読み出して、演算に用いる。
【0037】
具体的には、電子状態算出部15は、算出対象の1つ前の時間ステップn−1において、電子状態を示す波動関数ψ(n−1)からハミルトニアンH(n−1)を算出する。続いて、電子状態算出部15は、ハミルトニアンH(n−1)から時間発展演算子T(n−1)を算出する。続いて、電子状態算出部15は、時間ステップnにおける波動関数ψ(n)を、以下の式により算出する。
【数16】

【0038】
また、電子状態算出部15は、電子状態として波動関数ψ(n)だけでなく、例えば、以下の式により電子密度ρ(n)を算出することとしてもよい。
【数17】

【0039】
ここで、wは、重み係数を示す値であり、予め記憶部11に記憶されており、電子密度ρ(n)の算出の際に記憶部11から読み出されて算出に用いられる。
【0040】
物性量算出部16は、電子状態算出部15により算出された物質の電子状態を示す情報から、予め記憶した物性量算出ルールに基づき当該物質の物性量を算出する物性量算出手段である。物性量算出部16は、算出した物性量を示す情報を記憶部11に記憶させる。
【0041】
具体的には例えば、物性量算出部16は、物質の物性量として、発光・吸収スペクトルを算出するための双極子モーメントμ(n)を、物性量算出部16に算出された電子密度ρ(n)を記憶部11から読み出して、以下の式
【数18】

により算出する。ここで、rは、空間上の位置を示すベクトルである。物性量算出部16は、双極子モーメントμ(n)を各時間ステップ毎に算出する。
【0042】
また、物性量算出部16は、物質の物性量として、上記以外にも、ベクトルポテンシャルA(n)を用いて以下の式により表面電荷σ(n)を算出することとしてもよい。
【数19】

【0043】
出力部17は、物性量算出部16によって算出された物質の物性量を示す情報を、記憶部11から読み出して、出力する出力手段である。出力部17は、例えば、電子状態計算装置10が備えるディスプレイ等の表示装置に物性量の情報を表示することによって出力を行う。また、電子状態計算装置10から他の装置に出力することとしてもよい。また、出力部17は、電子状態算出部15によって算出された物質の電子状態を示す情報を、記憶部11から読み出して、出力することとしてもよい。
【0044】
係数算出部18は、初期値Vinit(n)の算出に用いるハートレー項を含む物質の電子状態を示す情報を記憶部11から読み出して、当該ハートレー項を含む物質の電子状態を示す情報から予め記憶した係数値算出ルールに基づき、I個の係数αを算出して記憶部11に記憶させる係数算出手段である。ここで算出されるαは適切なα、即ち、算出するハートレー項V(n)により近い初期値Vinit(n)を算出させるものである。係数算出部18は、算出した係数αを記憶部11に記憶させる。
【0045】
具体的には、係数算出部18は、以下のような最適化解を求める処理により、係数αを求める。以下では、初期値Vinit(n)の算出に、2つ前までの時間ステップでのハートレー項を用いる場合の例を示す。まず、空間をm個に差分化し表現すると、ハートレー項V(n)はm次元のベクトルで表される。差分化に伴う定数をm×mの正定値対称行列Aとし、電子密度ρ(n)に関連した既知量をm次元のベクトルb(n)とすると、式(7)は、
【数20】

という連立一次方程式の形に表すことができる。このとき、共役勾配法と同様に重みつき誤差(Vinit(n)−V(n),A(Vinit(n)−V(n)))を最小とする最適な係数αは、関数
【数21】

の最小化を
【数22】

により計算することで、求めることができる。
【0046】
具体的には、表式を簡単にするため時間ステップjにおけるハートレー項V(j)は、x、また、算出対象の初期値Vinit(n)をxと書き表すと、2つ前までの時間ステップを用いる場合x=xn−1+α(xn−1−xn−2)(ここで、α=αである)であるので、この関数f(x)は、
【数23】

と書き表せる。この関数f(x)を最小化する係数αは式(18)から、
【数24】

と求めることができる。ここで、ハートレー項算出部14により時間ステップnより前の時間ステップにおいて算出されたハートレー項V(j)は、連立一次方程式(16)の近似解であってもよいとし、以下の関係式
【数25】

が成り立つものとした。
【0047】
また、3つ前までの時間ステップを用いる場合その係数をα(=α)及びβ(=α)とすると、x=xn−1+α(xn−1−xn−2)+β(xn−2−xn−3)であるので、この関数f(x)は、
【数26】

となる。係数α及びβについて、この関数f(x)の停留点は、
【数27】

及び
【数28】

から、
【数29】

と求めることができる。また、上記算出された係数α及びβにおいて、関数f(x)が最小値を取る条件は、
【数30】

を満たす場合である。
【0048】
更に、最適な係数αを求める別の方法として、残差ベクトルb(n)−AVinit(n)の大きさを最小化しても良い。即ち、関数
【数31】

を用い、
【数32】

の計算から、残差を最小とする係数αを求めることも可能である。
【0049】
計算制御部19は、本実施形態に係る電子状態計算方法のTDDFT法による処理を全体的に統括する手段である。具体的には、上記の各処理を行う構成要素12〜19からの通知を受けて、各処理を行う構成要素12〜19に対して次に行うべき処理を指示する。計算制御部19は、予め設定されたステップ数の計算が完了したか否かを判断して、計算を終了させる等を行う。ここで、計算制御部19は、設定されたステップ数の情報を記憶部11から読み出して、上記の判断に用いる。
【0050】
上記の各処理を行う構成要素12〜19は、例えば、電子状態計算装置10が備えるCPUにより実現される。以上が電子状態計算装置10の構成である。
【0051】
引き続いて、図3のフローチャートを用いて、電子状態計算装置10により実行される本実施形態に係る電子状態計算方法を説明する。
【0052】
電子状態計算装置10では、まず、入力部12によって、物質の電子状態の算出に必要な、上述したパラメータ等の情報が入力される(S01、入力ステップ)。入力された各情報は、記憶部11に記憶される。入力される情報には、物質の電子状態を示す波動関数の初期値ψ(0)及び外部刺激が含まれる。ここで、TDDFT法における時間ステップnは0(初期値)とされる。
【0053】
続いて、時間ステップnは、それまでの時間ステップnに1が加えられたものとされる。続いて、電子状態算出部15によって、時間ステップnにおける波動関数ψ(n)が算出される(S02、電子状態算出ステップ)。具体的には、前の時間ステップn−1において算出されたハミルトニアンH(n−1)から時間発展演算子T(n−1)が算出される。続いて、前の時間ステップn−1において算出された波動関数ψ(n−1)と時間発展演算子T(n−1)とが上述したように演算されることにより、時間ステップnにおける波動関数ψ(n)が算出される。なお、時間ステップnが1の場合は、時間ステップn−1(即ち時間ステップが0の場合)における各値は、算出された値でなく、記憶部11に初期値として記憶されているものが用いられる。算出された波動関数ψ(n)は、記憶部11に記憶される。
【0054】
続いて、初期値算出部13によって、時間ステップnにおける、ハートレー項V(n)の初期値Vinit(n)が算出される(S03、初期値算出ステップ)。初期値Vinit(n)は、上述したように、時間ステップnの1つ前の時間ステップn−1からI+1前の時間ステップn−I−1までのハートレー項V(n−1)〜V(n−I−1)から算出される。なお、時間ステップnがI以下である場合は、各ハートレー項V(n−1)〜V(n−I−1)が取得できないため、上記の方法による初期値Vinit(n)の算出は行わない。その場合、従来と同様に、前の時間ステップn−1のハートレー項V(n−1)を初期値Vinit(n)とすることとしてもよい。算出された初期値Vinit(n)は、記憶部11に記憶される。
【0055】
ここで、初期値Vinit(n)の算出の前に、係数算出部18によって、初期値Vinit(n)の算出に用いられる係数αの算出が行われることとしてもよい(S04、係数算出ステップ)。係数αの算出が行われた場合、係数αは記憶部11に記憶されて、初期値算出部13に初期値Vinit(n)の算出の際に読み出されて用いられる。なお、係数αの算出が行われない場合、予め係数αが記憶部11に記憶されている。
【0056】
続いて、ハートレー項算出部14によって、初期値算出部13により算出されたハートレー項V(n)の初期値Vinit(n)から、予め記憶したハートレー項算出ルールに基づき差分法により当該ハートレー項V(n)が算出される(S05、ハートレー項算出ステップ)。算出されたハートレー項V(n)は、記憶部11に記憶される。
【0057】
続いて、電子状態算出部15によって、算出された波動関数ψ(n)からハミルトニアンH(n)が算出される(S06、電子状態算出ステップ)。算出されたハミルトニアンH(n)は、記憶部11に記憶される。
【0058】
続いて、電子状態算出部15によって、算出された波動関数ψ(n)から電子密度ρ(n)が算出される(S07、電子状態算出ステップ)。算出された電子密度ρ(n)は、記憶部11に記憶される。
【0059】
続いて、物性量算出部16によって、算出された電子密度ρ(n)から、物質の物性量である双極子モーメントμ(n)が、予め記憶した物性量算出ルールに基づき算出される(S08、物性量算出ステップ)。なお、算出される物質の物性量は、必ずしも双極子モーメントである必要は無く、物質の電子状態から算出されるものであればよい。算出された双極子モーメントμ(n)は、記憶部11に記憶される。
【0060】
続いて、出力部17によって、物性量算出部16によって算出された物質の物性量を示す情報が出力される(S09、出力ステップ)。また、併せて、算出された波動関数ψ(n)等の物質の電子状態を示す情報が出力されてもよい。なお、この出力は、各時間ステップnにおいて行われる必要は必ずしも無く、本処理の最後にまとめて行われてもよい。
【0061】
続いて、計算制御部19によって、現在の時間ステップnが、処理を終了させるために予め設定された時間ステップ数であるか否かが判断される(S10)。現在の時間ステップnが、設定された時間ステップ数であった場合は、ここで処理は終了される。一方、現在の時間ステップnが、設定された時間ステップ数でなかった場合は、上述したように、それまでの時間ステップnに1が加えられたものして、即ち、時間ステップを更新して、S02〜S10の処理を繰り返し行う。以上が、本実施形態に係る電子状態計算方法である。
【0062】
上述したように、本実施形態では、時間ステップnにおけるハートレー項V(n)の初期値Vinit(n)が、当該時間ステップの前の複数の時間ステップにおけるハートレー項に基づいて算出される。このように算出された初期値Vinit(n)は算出するハートレー項V(n)の値に近くなるため、ハートレー項V(n)の算出を従来と比べてより少ない繰り返し計算で行うことが可能になる。これにより、本発明に係る電子状態計算方法によれば、TDDFT法による物質の電子状態の計算を効率的に行うことができる。
【0063】
また、本実施形態のように、物質の物性値を算出することとすれば、効率的に物質の物性量を算出することができる。
【0064】
また、本実施形態のように、物質のセルフコンシステントな電子状態に対して外部刺激を与えることにより、当該電子状態の外部刺激に対する応答を計算することとすれば、時間変化に応じた物質の電子状態あるいは物性量を適切に算出することができる。
【0065】
また、本実施形態のように、ハートレー項V(n)を算出する時間ステップnの2つ前までのハートレー項を用いて、初期値Vinit(n)を算出することとすれば、簡易かつ確実に物質の電子状態の計算を効率的に行うことができる。
【0066】
また、本実施形態のように初期値Vinit(n)を算出するための式に含まれる係数αを算出することとすれば、ハートレー項V(n)の初期値Vinit(n)を、より算出するハートレー項V(n)の値に近づけることができ、ハートレー項V(n)の算出を従来と比べて更に少ない繰り返し計算で行うことが可能になる。しかしながら、この場合、係数αを算出する処理が増えることとなるので、係数αの算出によるハートレー項V(n)の算出のための繰り返し計算を減少させる効果が顕著な場合に用いることが望ましい。
【0067】
引き続いて、上述した一連の物質の電子状態の計算を行う処理をコンピュータに実行させるための電子状態計算プログラムを説明する。図4に示すように、電子状態計算プログラム31は、コンピュータが備える記録媒体30に形成されたプログラム格納領域30a内に格納されている。
【0068】
電子状態計算プログラム31は、記憶モジュール31aと、入力モジュール31bと、初期値算出モジュール31cと、ハートレー項算出モジュール31dと、電子状態算出モジュール31eと、物性量算出モジュール31fと、出力モジュール31gと、係数算出モジュール31hと、計算制御モジュール31iとを備えて構成される。記憶モジュール31a、入力モジュール31b、初期値算出モジュール31c、ハートレー項算出モジュール31d、電子状態算出モジュール31e、物性量算出モジュール31f、出力モジュール31g、係数算出モジュール31h及び計算制御モジュール31iを実行させることにより実現される機能は、上述した電子状態計算装置10の記憶部11、入力部12、初期値算出部13、ハートレー項算出部14、電子状態算出部15、物性量算出部16、出力部17、係数算出部18及び計算制御部19の機能とそれぞれ同様である。
【0069】
なお、電子状態計算プログラム31は、その一部若しくは全部が、通信回線等の伝送媒体を介して伝送され、他の機器により受信されて記録(インストールを含む)される構成としてもよい。
【実施例】
【0070】
次に本発明を実施例によって更に詳しく解説するが、本発明は実施例に限定されるものではない。
【0071】
[実施例1]
フルオレン1分子を本発明における電子状態の算出対象とした。本実施例では、密度汎関数法を用いた実空間実時間TDDFT法において、初期値Vinitの算出に、式(12)を用いた。また、初期値Vinitとして式(8)を用いて算出した場合を比較例とする。各場合について、1,000ステップまでの各時間ステップにおけるハートレー項Vの算出に要する繰り返し回数、及び10,000ステップまでの平均繰り返し回数を計測した。
【0072】
本実施例で用いる入力等として、密度汎関数法としては交換相関相互作用にP.J.Perdew and Alex Zunger, ”Self-interaction correction to density-functionalapproximation for many-electron systems”, Phys. Rev. B23, p5048(1981)(以後、Perdew-Zungerとして引用)により与えられるパラメータを用い、また擬ポテンシャルはN.Troullior and Jose Luis Martins, ”Efficient pseudopotentials for plane-wavecalculations”, Phys. Rev. B43, p1993(1991)(以後、Troullior-Martinsとして引用)により示された手法により予め生成しておいたものを用いた。原子配置は最適化された位置に固定し、空間の差分格子間隔は0.3Å、格子点数は、65×53×47である。これによりセルフコンシステントになるまでKohn-Sham方程式に従って電子状態を予め生成しておいた。これに外部刺激として0.01V/Åのδ関数型の電場をt=0において印加した。
【0073】
また、波動関数の時間発展計算は、時間間隔2πΔt/h(但し、ここではhはプランク定数を表す)を0.002eV−1として、時間発展演算子をΔtについて4次まで展開することにより行われた。また、上述したように電子状態として電子密度ρ(n)の算出を行った。また、それに基づき、物質の物性量として、物理量である光の吸収発光に関る双極子モーメントμ(n)を算出した。また、これらの計算は時間ステップが10、000ステップとなることを判定条件として計算を終了させることとした。
【0074】
ハミルトニアンH(n)の更新は式(6)により行なわれるが、ここでのハートレー項V(n)を求めるにあたり式(7)を実空間での差分法により連立方程式に置き換え、共役勾配法を用いた。このときの初期値として、本発明における予測手法を用い、最適な初期値Vinit(n)を提供することとした。初期値Vinit(n)の予測には式(12)を用い、α=1とした。
【0075】
この結果、10,000ステップまでのハートレー項Vの算出に要する平均繰り返し回数は、比較例では17回であるのに対し、本発明では5回であった。本発明は、従来の方法に比べて繰り返し回数が非常に小さくて済み、計算時間の短縮に大きく寄与していることが認識される。
【0076】
また初期1,000ステップまでの各時間ステップにおいてハートレー項Vの算出に要する繰り返し回数の時間発展の様子を図5に示す。図5において横軸はTDDFT法による時間ステップを示す値であり、縦軸は各時間ステップにおいてハートレー項Vの算出に要する繰り返し回数である。本発明による方法が、初期の時間ステップではより効果があることが認識される。
【0077】
[実施例2]
実施例1と同様に、本発明による計算と、比較例による計算とを行い、10,000ステップまでの平均繰り返し回数を計測した。本実施例では、フルオレン1量体、2量体、3量体、4量体、6量体及び8量体に対して、それぞれ計算を行った。
【0078】
本実施例で用いる入力及び処理内容等は、実施例1とほぼ同様である。交換相関相互作用にPerdew-Zungerにより与えられるパラメータを用い、また擬ポテンシャルはTroullior-Martinsにより示された手法により予め生成しておいたものを用いた。原子配置はいずれも最適化された位置に固定し、空間の差分格子間隔は0.3Åで、格子点数はフルオレン1量体の場合は65×53×47、2量体の場合は91×55×51、3量体の場合は119×55×53、4量体の場合は147×55×55、6量体の場合は203×55×55、8量体の場合は257×55×55である。これによりセルフコンシステントになるまでKohn-Sham方程式に従って電子状態を予め生成しておいた。
【0079】
また、実施例1と同じ外部刺激を加え、時間間隔2πΔt/h(但し、ここではhはプランク定数を表す)を0.002eV−1とし、物理量として双極子モーメントμ(n)を算出した。
【0080】
特に本発明で重要なハミルトニアンH(n)の更新は、実施例1と同様に式(6)により行なわれるが、ここでのハートレー項Vを求めるにあたり式(7)を実空間での差分法により連立方程式に置き換え、共役勾配法を用いた。このときの初期値として、本発明における予測手法を用い、本発明による初期値Vinit(n)を提供することとした。初期値Vinit(n)の予測には式(12)を用い、α=1とした。
【0081】
この結果、10,000ステップまでのハートレー項Vの算出に要する平均繰り返し回数と量数との関係について、従来法と本発明とによる予測手法を用いることにより得られた結果を、図6に示す。
【0082】
従来法では分子形状(例えば高分子鎖長)が長くなるにつれて、計算にかかるコストが非常に大きくなるが、本発明を用いると共役勾配法の平均繰り返し回数は鎖長に依らずほぼ一定の値、この場合5回程度となっていることがわかる。このように、本発明は、分子形状による影響を最小限にとどめる効果があることもわかる。
【図面の簡単な説明】
【0083】
【図1】本発明の実施形態に係る電子状態計算装置の構成図である。
【図2】本発明の実施形態において算出されるハートレー項の初期値を概念的に示した図である。
【図3】本発明の実施形態に係る電子状態計算方法を示すフローチャートである。
【図4】本発明の実施形態に係る電子状態計算プログラムを示す図である。
【図5】本発明と比較例とによる電子状態の算出の際の、各時間ステップにおけるハートレー項Vの算出に要する共役勾配法の繰り返し回数を示したグラフである。
【図6】本発明と比較例とによる電子状態の算出の際の、フルオレンの量体数に応じた、ハートレー項Vの算出に要する共役勾配法の平均繰り返し回数を示したグラフである。
【符号の説明】
【0084】
10…電子状態計算装置、11…記憶部、12…入力部、13…初期値算出部、14…ハートレー項算出部、15…電子状態算出部、16…物性量算出部、17…出力部、18…係数算出部、19…計算制御部、20…外部装置、30…記録媒体、30a…プログラム格納領域、31…電子状態計算プログラム、31a…記憶モジュール、31b…入力モジュール、31c…初期値算出モジュール、31d…ハートレー項算出モジュール、31e…電子状態算出モジュール、31f…物性量算出モジュール、31g…出力モジュール、31h…係数算出モジュール、31i…計算制御モジュール。

【特許請求の範囲】
【請求項1】
コンピュータにより、TDDFT法によって物質の電子状態を計算する電子状態計算方法であって、
前記物質の電子状態の初期状態を示す情報を含む前記TDDFT法によって前記物質の電子状態を計算するためのパラメータを入力して、前記TDDFT法における時間ステップj毎の、前記物質の電子状態のハミルトニアンH(j)におけるハートレー項V(j)を含む当該物質の電子状態を示す情報を記憶する記憶手段に記憶させる入力ステップと、
時間ステップnにおいて、時間ステップn−1から時間ステップn−I−1までの前記ハートレー項V(n−1)〜V(n−I−1)を前記記憶手段から読み出して、予め設定されるI個の係数α〜αを用いて、前記時間ステップnにおけるハートレー項V(n)の初期値Vinit(n)を予め記憶した式
【数1】

により演算をして算出する初期値算出ステップと、
前記初期値算出ステップにおいて算出された前記ハートレー項V(n)の初期値Vinit(n)から、予め記憶したハートレー項算出ルールに基づき差分法により当該ハートレー項V(n)を算出すると共に、算出したハートレー項V(n)を前記記憶手段に記憶させるハートレー項算出ステップと、
前記ハートレー項算出ステップにおいて算出された前記物質の前記ハートレー項V(n)及び前記記憶手段により記憶された情報から、予め記憶した電子状態算出ルールに基づき前記物質の電子状態を算出すると共に、算出した電子状態を示す情報を前記記憶手段に記憶させる電子状態算出ステップと、
を含む電子状態計算方法(但し、Iは予め定められる1以上の整数であり、n,jは時間ステップを示すインデックスである)。
【請求項2】
前記電子状態算出ステップにおいて算出された前記物質の電子状態を示す情報から、予め記憶した物性量算出ルールに基づき当該物質の物性量を算出する物性量算出ステップと、
前記物性量算出ステップにおいて算出された前記物質の物性量を示す情報を出力する出力ステップと、
を更に備える請求項1に記載の電子状態計算方法。
【請求項3】
前記入力ステップにおいて、前記物質の電子状態を計算するためのパラメータとして、前記物質の電子状態の初期状態を示す情報に与えられる外部刺激を示す情報を入力することを特徴とする請求項1又は2に記載の電子状態計算方法。
【請求項4】
前記初期値算出において、前記時間ステップnにおけるハートレー項V(n)の初期値Vinit(n)を予め記憶した式
【数2】

により演算をして算出する、
ことを特徴とする請求項1〜3の何れか一項に記載の電子状態計算方法。
【請求項5】
前記記憶手段から前記ハートレー項V(j)を含む前記物質の電子状態を示す情報を読み出して、当該ハートレー項V(j)を含む当該物質の電子状態を示す情報から予め記憶した係数値算出ルールに基づき、係数αを算出して当該記憶手段に記憶させる係数算出ステップを更に含み、
前記初期値算出ステップにおいて、係数αを前記記憶手段から読み出して、前記初期値Vinit(n)の算出に用いるものとして設定する、
ことを特徴とする請求項1〜4の何れか一項に記載の電子状態計算方法。
【請求項6】
コンピュータに、TDDFT法によって物質の電子状態を計算させる電子状態計算プログラムであって、
前記TDDFT法における時間ステップj毎の、前記物質の電子状態のハミルトニアンH(j)におけるハートレー項V(j)を含む当該物質の電子状態を示す情報を記憶する記憶手段と、
前記物質の電子状態の初期状態を示す情報を含む前記TDDFT法によって前記物質の電子状態を計算するためのパラメータを入力して、前記記憶手段に記憶させる入力手段と、
時間ステップnにおいて、時間ステップn−1から時間ステップn−I−1までの前記ハートレー項V(n−1)〜V(n−I−1)を前記記憶手段から読み出して、予め設定されるI個の係数α〜αを用いて、前記時間ステップnにおけるハートレー項V(n)の初期値Vinit(n)を予め記憶した式
【数3】

により演算をして算出する初期値算出手段と、
前記初期値算出手段により算出された前記ハートレー項V(n)の初期値Vinit(n)から、予め記憶したハートレー項算出ルールに基づき差分法により当該ハートレー項V(n)を算出すると共に、算出したハートレー項V(n)を前記記憶手段に記憶させるハートレー項算出手段と、
前記ハートレー項算出手段により算出された前記物質の前記ハートレー項V(n)及び前記記憶手段により記憶された情報から、予め記憶した電子状態算出ルールに基づき前記物質の電子状態を算出すると共に、算出した電子状態を示す情報を前記記憶手段に記憶させる電子状態算出手段と、
を前記コンピュータに機能させるための電子状態計算プログラム(但し、Iは予め定められる1以上の整数であり、n,jは時間ステップを示すインデックスである)。
【請求項7】
TDDFT法によって物質の電子状態を計算する電子状態計算装置であって、
前記TDDFT法における時間ステップj毎の、前記物質の電子状態のハミルトニアンH(j)におけるハートレー項V(j)を含む当該物質の電子状態を示す情報を記憶する記憶手段と、
前記物質の電子状態の初期状態を示す情報を含む前記TDDFT法によって前記物質の電子状態を計算するためのパラメータを入力して、前記記憶手段に記憶させる入力手段と、
時間ステップnにおいて、時間ステップn−1から時間ステップn−I−1までの前記ハートレー項V(n−1)〜V(n−I−1)を前記記憶手段から読み出して、予め設定されるI個の係数α〜αを用いて、前記時間ステップnにおけるハートレー項V(n)の初期値Vinit(n)を予め記憶した式
【数4】

により演算をして算出する初期値算出手段と、
前記初期値算出手段により算出された前記ハートレー項V(n)の初期値Vinit(n)から、予め記憶したハートレー項算出ルールに基づき差分法により当該ハートレー項V(n)を算出すると共に、算出したハートレー項V(n)を前記記憶手段に記憶させるハートレー項算出手段と、
前記ハートレー項算出手段により算出された前記物質の前記ハートレー項V(n)及び前記記憶手段により記憶された情報から、予め記憶した電子状態算出ルールに基づき前記物質の電子状態を算出すると共に、算出した電子状態を示す情報を前記記憶手段に記憶させる電子状態算出手段と、
を備える電子状態計算装置(但し、Iは予め定められる1以上の整数であり、n,jは時間ステップを示すインデックスである)。

【図1】
image rotate

【図2】
image rotate

【図3】
image rotate

【図4】
image rotate

【図5】
image rotate

【図6】
image rotate


【公開番号】特開2009−26140(P2009−26140A)
【公開日】平成21年2月5日(2009.2.5)
【国際特許分類】
【出願番号】特願2007−189829(P2007−189829)
【出願日】平成19年7月20日(2007.7.20)
【国等の委託研究の成果に係る記載事項】(出願人による申告)文部科学省、平成18年度「先端大型研究施設戦略活用プログラム[地球シミュレータ戦略活用プログラム]」)の成果に係る出願
【出願人】(000002093)住友化学株式会社 (8,981)
【出願人】(504194878)独立行政法人海洋研究開発機構 (110)