3E01
多配置展開法による一次元水素分子の non-Born-Oppenheimer 波動関数
（東大院・理） 井手 善広，加藤 毅，山内 薫

Non-Born-Oppenheimer wave function of 1D hydrogen moleculeby time-dependent multiconfiguration expansion(School of Science, The University of Tokyo) Yoshihiro Ide, Tsuyoshi Kato, Kaoru Yamanouchi【序】

強光子場の炭化水素分子においては，水素原子またはプロトンが分子内で高速に移動す

る水素マイグレーションと呼ばれる現象が起こることが知られている[1,2]．この水素マイグレー
ション過程における分子内のプロトン移動は，電子の運動と強く相関していると考えられる．この異種粒子間の量子力学的相関運動を取り扱うために，電子動力学理論である時間依存多配置波
動関数(Multiconfiguration time-dependent Hartree-Fock, MCTDHF)理論が多成分系の取り扱いへと拡張された[3]．
本研究では，拡張された MCTDHF 理論が Schrödinger 方程式の厳密解法と整合した結果を与えることを，数値的に検証する．具体的には，一次元水素分子の電子—核基底状態を計算対象として以下の点に着目して解析を行った．まず，MCTDHF 理論によって与えられる基底状態の波動関数
  を用いて，エネルギー期待値   を計算した．また，別途求められている断熱電子状態を使って，各断熱電子状態上の核の相対運動の波動関数    を   から抽出した．ここで， は核間距離を表す．これら   および   を，厳密解波動関数 ることのできる および から同様にして求めと比較し，拡張された MCTDHF 理論が一次元水素分子の基底状態における異種粒子間の量子力学的相関運動を記述できることを確認する．このことは，拡張された MCTDHF 理論に基づいた強光子場中の分子ダイナミクスの解析手法が妥当であることを
示す．
【理論】

拡張された MCTDHF 理論において，一次元水素分子の電子—核波動関数   は次のように展開される．
   ,  , ,   ∑  ,    ,  , (1)

ここで，    ,  ,   1, 2 は電子  の位置座標  およびスピン座標  をまとめて表す．

は電子配置を表し，配置  に対する Slater 行列式を  とし，配置の重みを表す関数を  としている． が核間距離  に依存しない点が，通常の断熱電子状態   ,  ;  による展開（Born-Huang 展開）とは異なる．ここで，! は断熱電子状態を指定する量子数である．
配置の重み関数  ,  および   ,  ,  を構成する " 番目の電子軌道 #$   ,  に対する運動方程式は，次の Dirac-Frenkel の時間依存変分原理を式(1)へ適用することによって得られる．
%& ' () * +,【結果と考察】.'  /  0-(2)

式(1), (2)から得られた， ,  および #$   ,  に対する運動方程式を虚時

間領域で時間発展させ，基底状態の波動関数   を求めた．以下，基底状態について解析を行うので，時刻の引数  を省略する．

  を用いて，エネルギー期待値   1 を計算した．ここで，1 は   の記述に用いた

電子軌道の総数である．計算値を図 1 に示す．1 が増加するに従って，  1 が厳密解 へ向かって単調に収束することが分かった．このことは，電子間および電子と核の間の粒子間相関
エネルギーの評価が，1 の増加と共に，徐々に正確化されたことを意味する（1  5 のとき，99.9 %の粒子間相関エネルギーが計算された）．
次に，波動関数の物理的な内容を比較するため，核の相対運動の波動関数を比較する．通常のBorn-Huang 展開の場合，! 番目の断熱電子状態  上の，核の相対運動の波動関数 34  は次のように定義される．
6  (3)

ただし，Born-Huang 展開された電子—核波動関数  34 は次の形で与えられる．
 34  ,  ,   ∑ 34    ,  ; 

式(3)と同様に考えて， 



および 

(4)

に対しても，断-1.41
Energy / a.u.

34   5   ,  ;   34  ,  ,  6

-1.43
-1.45
1 2 3 4 5
Number of Orbitals

熱電子状態上の，核の波動関数を抽出することができる．
   5   ,  ;     ,  ,  6
   5   ,  ;   , ,  6

6  (5)6 (6)式(5), (6)に基づいて計算された，断熱電子基底状態  上

図 1. エネルギー期待値．MCTDHF
理論（赤），厳密解法（黒）．

の，核の相対運動の波動関数を図 2 に示す．1  1 （電

1.0

厳密解 

0.5

子—核相関に対する Hartree 近似）のとき，  の形状はの形状と比べると幅が狭く，より短い平衡核

間距離を与える．しかし，1  5 のとき，  と は，

図 2 のスケールでは互いに重なる．すなわち，1 の増加に

より，






が，

を構成する  上の核の相対運動の波動関数を構成する核の相対運動の波動関数 に漸近したことが分かる．このことは，電子と核の相対運動の相関の記述が 1 の増加と共に改良されていくことを，核の波動関数の収束によって示したことになる．
【結論】

0.0
1.0 2.0 3.0
Internuclear Distance / a0

図 2. 電子基底状態  上の，核の相対運動の波動関数  ．
MCTDHF 理 論 （ 1  1 , 青；
，厳密解法 （灰色）．
1  5, 赤）

拡張された MCTDHF 理論を使って一次元水素分子の電子—核基底状態   を計算した．  を記述するために用いる電子軌道の数 1 を増加させることによって，エネルギー期待値   1 が厳密解 へ収束していくことが示された．このことから，異種粒子間の相関エネルギーを系統的に改良させることが可能であることが分かった．また同時に，1 の増加と共に，
断熱電子基底状態上の核の波動関数   が，厳密解 ち，電子運動と相関する核の運動が へと漸近することが示された．すなわ

によって記述されることが示された．以上のことから，拡張された MCTDHF 理論の枠内で，電子波動関数を記述するための電子軌道の数を増加させることによって，電子—核波動関数を系統的に改良できることが数値的に確認された．

[1] A. Hishikawa, H. Hasegawa, and K. Yamanouchi, J. Electron Spectrosc. Relat. Phenom. 141, 195 (2004). [2]T. Okino, Y. Furukawa, P. Liu, T. Ichikawa, R. Itakura, K. Hoshina, K. Yamanouchi, and H. Nakano, Chem. Phys.
Lett. 423, 220 (2006). [3] T. Kato and K. Yamanouchi, J. Chem. Phys. 131, 164118 (2009).

3E02

電子ダイナミクスにおける断熱・非断熱経路の
電子配置・分子軌道エネルギー解析（東北大院理 1，横浜市大院生命ナノ 2，東大院理 3，大阪府大院理 4） 河野小山田隆行 1，大村周 1，加藤毅 3，小関史朗 4，菅野裕彦 1，

学1

Energy Analysis of Time-dependent Electron Configurations and MolecularOrbitals in Adiabatic and Nonadiabatic Electron Dynamics(Tohoku Univ. 1, Yokohama City Univ. 2, The Univ. of Tokyo3, Osaka Prefecture Univ. 4)
Hirohiko Kono1, Takayuki Oyamada2, Shu Ohmura1, Tsuyoshi Kato3, Shiro Koseki4,
Manabu Kanno1【序】時間に依存したレーザー光などの外場と分子中の電子との非摂動論的相互作用が引き起こす動力学過程の理論的記述は、分子分光のみならず化学反応における電子のアト秒スケールの実時間運動を理解する上でも極めて有用である。電磁場と電子との１体の相互作用をとおして分子によって獲得されたエネルギーは、多電子間のエネルギー交換を経て、振動な
どの自由度に流れていく。多電子の全系にどのようにエネルギーが流れていくかには、電子間の相互作用、つまり、電子相関が深く関わっている。最近、近赤外レーザーによるイオン化において、一つの電子が連続的に多光子を吸収してイオン化するか、多くの電子がエネル
ギーを交換しながらイオン化するかを区別する実験が報告され[1]、多電子間あるいは軌道間のエネルギー交換を定量化する手法の開発が必要とされている。本研究では，電子ダイナミ
クスを支配する要因を明らかにするため，電子配置や自然軌道に対して化学ポテンシャルなどを提案し，強い近赤外光と相互作用するH2やLiHに適用した結果について報告する。
【Multiconfiguration Time-Dependent Hartree-Fock（MCTDHF）法の開発】本研究では、多電子系ダイナミクスを評価するために、時間t に依存するSlater 行列式 Φ I t  で展開された多配置波動関数 Φt  を用いた[2]。

Φ  t    I 1 CI  t  Φ I  t 
M(1)

ここで、C I t  は配置間相互作用係数である（M は配置の数）。Slater 行列式のなかの分子軌

道  k t は、イオン化の連続状態を記述できるように空間の格子点上の振幅（グリッド基底）で表現した。(1)式をDirac-Frenkel の時間依存の変分原理に代入すると、 C I t  と  k t が従う運動方程式が求まる。このMCTDHF 法では、 k t も時間発展させるので、少ない電子配置で多電子系の時間発展を効率的に近似できる。
【電子配置の特性エネルギー】次に，電磁場と電子との１体の相互作用から多電子の全系にどのようにエネルギーが流れていくかを定量化する。(1)式を使うと，外場と相互作用する系の全エネルギー E  t  は次式で与えられる。
E  t   Φ  t  Hˆ  t  Φ  t   I CI  t   Φ I  t  Hˆ  t  Φ I  t  CI  t 
M

2



1



M
I (  I )CI   t  Φ I  t  Hˆ  t  Φ I   t  (2)

第1項は各軌道に対する平均場近似エネルギーの和で，第2項は軌道間の対相関エネルギーの和である。ある瞬間の決まった軌道に対する各電子配置 Φ I t  のエネルギーを EI (t ) E (t )  CI (t ) で定義すると， E  t  は次式のように  EI  t  の和で与えられる。
M
2
E (t )   I CI EI (t )(3)2

この定義を使えば，配置間のエネルギー交換を定量化できる。また，電子ダイナミクスを外

場に対して断熱か非断熱かに分類できる。つまり，すべての EI  t  が同じ値 E  t  をとる場
合，電子状態変化は外場に対して断熱的であり，そうでなければ，非断熱的である。
一方， E  t  は「電子相関エネルギーを含んだ瞬間的な自然軌道  j (t ) の時間依存化学ポテ
ンシャル  j (t ) 」にも分割できる[3,4]。

E t    j wj t   j t 
No

(4)

誘起双極子

ここで， No は軌道の数， w j t  は  j t  の占有数
で，  j (t )  E (t ) w j (t ) である。 No=M の場合，

 j (t ) は EI  t  を使ってユニークに求まる。また，
 j t  に対しては、レーザー電場  t  からの１体相

レーザー電場

互作用で直接得られるエネルギーSj(t)を定義できる（ d j は  j の双極子モーメント）[3]。

  t    d d  t  dt 

S j  t    dt t0j

図1 ピーク強度 3  10 W cm 、波長12(5)E  t   w j (t )  j (t ) であり，軌道間のエネルギー交換

が無ければ， S j  t    j (t ) である。

【MCTDHFによる計算結果：高強度近赤外光による分子のイオン化】ピーク強度 1014 W cm 2 程度の近赤外光と相互作用する H 2 分子のイオン化を調べた（核間距離21522 nmの２サイクルパルスと相互作用する LiH （R=3 Bohr）の2自然軌

道 2  t  に入ったエネルギー S 2  t （点線）と  2  t  （赤線）。両者とも双極子エネルギー d 2  t  (t ) を基準値として加えている。偏光方向は分子軸に平行で，LiからHに向かう方向

はR=1.6 Bohrに固定）。9つの分子軌道を考慮して，81を正と定義している。誘起双極子モー

個のSlater行列式の時間発展を求めた。偏光軸が分子軸

メントの最大値は0.55デバイである。

zと平行な場合，z軸に沿って大振幅の往復運動をする
が，最も空間的に広がる自然軌道は 2σ g 軌道である。  j  t    j  t    j  0  で定義される化学ポテンシャルの変化量を使うと， j (t ) > Sj(t)を満たすエネルギーアクセプタ軌道（例えば，

H2では， 2g軌道）に分類される軌道が分子のイオン化において重要であることがわかった。

 j (t ) > Sj(t)を生み出す他軌道とのエネルギーの交換には電子相関の変化，つまり，(2)式の第２項の変化が大きく寄与している。
３つの軌道1, 2, 3を使ったLiHの例を図１に示す。2の誘起双極子 d 2  t  は，光の1.5サ
イクル以降は電場の変化に対して非断熱的に振る舞っているが，他軌道とのエネルギー交換を反映した化学ポテンシャルの変化  2  t  とはきわめて良い相関を示している（ d 2  t  と2  t  は同じ時刻でゼロになる）。各軌道の化学ポテンシャルが非定常過程に対しても定義でき，その変化が自然軌道を占有した電子の動きの指標になっていることが明らかになった。
[1] I.V. Hertel et al., Phys. Rev. Lett. 102, 023003 (2009).
[2] T. Kato and H. Kono, J. Chem.Phys. 128, 184102 (2008).
[3] T. Kato and H. Kono, Chem. Phys. 366, 46 (2009).
[4] T. Kato, T. Oyamada, H. Kono, and S. Koseki, Prog. Theor. Phys. Suppl. (2012) in press.

3E03
強レーザー場中の多電子ダイナミクス：多配置波動関数理論の開発（東京大学） 佐藤 健、石川 顕一

Multielectron dynamics in intense laser fields: Multiconfigurational wavefunction
theories(University of Tokyo) Takeshi Sato, Kenichi L. Ishikawa
超短パルス高強度光源を用いて、電子の運動を電子固有の時間スケールで直接観測・操作するアト秒技術が急速に発展している【1】。とくに近年、実験の精密化に伴い、有効一電子描像を超える多電子ダイナミクスや電子相関の効果に関心が集まっている。これらは時間依存 Schrödinger
方程式 (TDSE) によって厳密に記述されるが、多電子系への応用は困難である。TDSE と一電子モデルの間のギャップを埋める、近似的な多電子波動関数理論が必要である。
高強度外場による電離過程や高次高調波発生を記述するには、束縛状態と電離状態をシームレスに記述できなければならない。しかし、もっとも簡単な時間依存 Hartree-Fock (TDHF) 法はレーザー電離を正しく記述できない。二電子系を例に考えると、HF 波動関数が各時刻で閉殻構造をとり、電離電子とコア電子の空間的に異なる運動を記述できないことが原因である。私達はこの問題を解決するため、Generalized Valence-Bond (GVB) 波動関数【2】に基づく時間依存理論を提案した【3】。一重項二電子系に対する GVB 波導関数の空間部分は次式で与えられる：

ΨGVB (r1 , r2 ) ∝ ψ1 (r1 ) ψ2 (r2 ) + ψ2 (r1 ) ψ1 (r2 )(1)= A1 ϕ1 (r1 ) ϕ1 (r2 ) + A2 ϕ2 (r1 ) ϕ2 (r2 ).

(2)

ψ1 , ψ2 は非直交である。一方、直交系 ϕ1 , ϕ2 は一体の密度行列を対角化する自然軌道である。二つの表現は可逆な線形変換で結ばれ、変分的に等価である。多電子系への拡張を見据え後者を APSG
波動関数と呼ぶ【4】。APSG では、軌道の数 n を増やして計算精度を向上させられる。

Ψ(r1 , r2 , t) =n∑Ai (t) ϕi (r1 , t)ϕi (r2 , t) =

i=1

n ∑
n
∑

Bi j (t) ϕ′i (r1 , t)ϕ′j (r2 , t).

(3)

i=1 j=1

上式は n 軌道を用いた変分的に最良の解を与える（自然展開）。従って二電子系では、同式最右辺に基づく時間依存多配置 HF (MCTDHF) 法【5】と等価である。式 (1) 自体にもおもしろい応用がある。すなわち一方の軌道を HF 占有軌道に凍結し他方を変分自由度とすれば、時間依存一電子励起配置間相互作用 (TD-CIS) 法【6】と等価な理論を非占有軌道を参照せずに構築できる。

ΨCIS (r1 , r2 ) = C1 (t)ϕ1 (r1 )ϕ1 (r2 ) +∑{}C j (t) ϕ1 (r1 )ϕ j (r2 ) + ϕ j (r1 ) ϕ1 (r2 )j≥2= ϕ1 (r1 )ψ(r2 , t) + ψ(r1 , t)ϕ1 (r2 ),ここで {ϕi } は占有 (i = 1) および非占有 (i ≥ 2) HF 軌道、ψ = C1 ϕ1 /2 +(4)∑j≥2 C j ϕ j である。

これらの表現に基づき、時間依存変分原理を用いて TD-GVB、TD-APSG、TD-CIS 方程式を導出し、厳密な TDSE 計算が可能な一次元ヘリウムモデル【7】に適用した。Fig. 1 は、波長 750

nm、パルス長（半値全幅）4 fs のレーザーパルスによる一重および二重電離確率を示している。
二配置の TD-APSG (TD-GVB) 法によって、特に高強度領域で、妥当な一重電離確率が与えられている。また非逐次二重電離 (NSDI) に特徴的な強度 1015 W/cm2 近傍の “knee structure”【7】を定性的に再現できている。しかし NSDI 領域の二重電離確率を精度良く求めるには巨大な活性空間が必要であり、再散乱過程における動的電子相関の重要性を示唆している。一方、TD-CIS 法
は強度 1.5 × 1015 W/cm2 程度以下の一重電離確率を極めて精度よく再現している。多重電離は原理的に記述できない。
以上の二電子系での定式化を多電子系に拡張する厳密な手法が MCTDHF 法であり【5】、全電子・軌道を “active” とした CASSCF【8】波導関数の運動方程式を与える。

ΨMCTDHF (r1 , r2 , ··, rN , t) = ΨCASSCF (r1 , r2 , ··, rN , t).

(5)電離過程において重要な電子のみを相関させる理論、あるいは重要性を吟味できる理論を創るために、“core” 軌道の概念を導入すべきである。

Ψ(r1 , r2 , ··, rN , t) = Â [Ψcore (r1 , r2 , ··, rn , t) Ψactive (rn+1 , rn+2 , ··, rN , t)] ,(6)

ここで Â は反対称化演算子である。Ψcore は n/2 個の規格直交軌道からなる閉殻構造である。Ψactiveには、例えば、CASSCF や APSG の active 空間を採用できる。後者は式 (1) または (3) の構造を持つ二電子波動関数の反対称化積である。式 (4) の CIS 波動関数を用いれば、一重電離を描ける
簡便な一電子理論を構築できる。当日の発表では、二電子系における定式化と数値検証の結果を報告する。また開発の状況に応じて、多電子系への拡張を提案する。
100

100

10-1

YieldsYields

10-1

10-2

10-3

TDSE
HF
TD-CISAPSG(02)
APSG(04)
APSG(08)
APSG(12)APSG(20)
1015
Intensity (W/cm2)

TDSE
HF
TD-CISAPSG(02)
APSG(04)
APSG(08)
APSG(12)APSG(20)

10-2

10-3

10-4

1015
Intensity (W/cm2)

Fig. 1：Single (left plot) and double (right) ionization yields as a function of intensity, obtained by TDSE, TDHF, TD-CIS, and TD-APSG methods.
【1】F. Krausz and M. Ivanov, Rev. Mod. Phys., 81, 163 (2009). 【2】F. B. Bobrowicz and W. A. Goddard III, in Methods of Electronic Structure Theory, H. F. Schaefer III, Ed., Plenum Press, 79(1977). 【3】佐藤健, 石川顕一, 第 5 回分子科学討論会, 4A07 (2011). 【4】Antisymmetrized Productof Strongly Orthogonal Geminals; Theoret. Chim. Acta, 1, 327 (1963). 【5】T. Kato and H. Kono,
Chem. Phys. Lett., 392, 533 (2004); J. Caillat et al., Phys. Rev. A, 71, 012712 (2005). 【6】N. Rohringer,
A. Gordon, and R. Santra, Phys. Rev. A, 74, 043420 (2006); L. Greenman et al, Phys. Rev. A,
74, 043420 (2010). 【7】N. E. Dahlen and R. van Leeuwen, Phys. Rev. A, 64, 023405 (2001).
【8】Complete Active Space Self-Consistent Field; B. O. Roos, P. R. Taylor, and P. E. Siegbahn,
Chem. Phys. 48, 157 (1980).

Phase space approach to intense-laser-driven
electronic wavepacket propagation(Weizmann Institute) Norio Takemoto, Asaf Shimshovitz, David J. Tannor
We propose a new method to solve the time-dependent Schrödinger equation based ona phase space perspective [1]. The method employs the periodic von Neumann basis withbiorthogonal exchange (pvb basis) recently introduced for the calculation of the energy
eigenstates of time-independent quantum systems [2]. While the individual elements inthis basis set are time-independent, a small subset is chosen in a time-dependent mannerto adapt to the evolution of the wavepacket in phase space. The resulting method is
simple, eﬃcient, accurate, ﬂexible, and stable. Our long-term goal is to apply the method
to simulate the correlated dynamics of multiple electrons in atoms and molecules in intenseand ultrashort laser pulses. As a ﬁrst demonstration, we present calculations of electron
dynamics in a 1D model atom interacting with the combined ﬁeld of intense near-infrared(NIR) and attosecond extreme-ultraviolet (XUV) laser pulses.
In Fig. 1, we show the simulated evolution of the 1D electronic wavepacket in phasespace. The NIR laser pulse (5 × 1013 W/cm2 , 800 nm, 4 fs) and the XUV laser pulse(1 × 1012 W/cm2 , 15 nm, 250 as) were applied to the electron in a 1D soft-core Coulomb
potential. In the pvb basis, the expansion coeﬃcients of a wavepacket are given by the
projection of the wavepacket to phase space Gaussians imposed with periodic boundary
conditions. Squared moduli of these coeﬃcients are shown by the red color scale on theellipses, which have the same center and aspect ratio as the corresponding phase space
Gaussians. A lattice of 64 × 64 Gaussians is required to span the rectangular boundary ofeach panel in Fig. 1, but only those Gaussians indicated by the ellipses have signiﬁcant
overlap with the wavepacket. The active Gaussians are predetermined by running classicaltrajectories for a simpliﬁed model. The representation with the full set of Gaussians isequivalent to the representation using the same number of Fourier grid points. However,the localized nature of the Gaussians allowed us to propagate the wavepacket using only
about 10 % of the pvb basis with negligible loss of accuracy.
Figure 2 shows the photoelectron momentum distribution that emerges from the simulation in Fig. 1 (blue solid line). Comparison with a simulation using the full pvb basis(red dashed line) indicates that the small active set not only reproduces the qualitative
features, such as the direct (N1 and N1’) and rescattered (N2 and N2’) NIR-photoelectron
cut-oﬀs as well as the NIR-streaked single-XUV-photon ionization peaks (X1 and X1’),
but also has quantitative accuracy.

[1] N. Takemoto, A. Shimshovitz, and D.J. Tannor, J. Chem. Phys. 137, 011102 (2012).
[2] A. Shimshovitz and D.J. Tannor, Phys. Rev. Lett. (in press) [e-print arXiv:1201.2299v1].

8 (a)
6
4
2
0
−2
−4
−6
−8
8 (b)
6
4
2
0
−2
−4
−6
−8
8 (c)
6
4
2
0
−2
−4
−6
−8
−600 −400 −200

Momentum p (a.u.)

Momentum p (a.u.)

Momentum p (a.u.)

10−12 10−9 10−6 10−30200

Position q (a.u.)

400

600

Probability Density (a.u.)

FIG. 1. Snapshots of the wavepacket in phase space. The squared moduli of the overlaps of the wavepacket tothe phase space Gaussians are shown by the red color scale. The dark- and light-blue dots as well as dark-blue +marks represent the classical trajectories used to predetermine which Gaussians will have significant overlap.

10-4

Full
Reduced

10-6

X1'
N2'

N1'

N1

X1

N2

10-8

10-10
10-12
10-14

−3

−2

−1

0

1

2

Photoelectron Momentum p (a.u.)3

FIG. 2. Photoelectron momentum distribution. The results with the reduced (blue solid line) and full (red dashedline) pvb bases are compared.

3E05
位相緩和制御のためのレーザーパルス列の理論的設計（慶大院・理工）○菅原道彦

Theoretical designing of laser pulse sequence for controlling dephasing dynamics(Department of Fundamental Science and Technology, Graduate School of Science and Technology, Keio University) M. Sugawara【序】レーザーによる多原子分子の振動ダイナミクスの量子制御を困難にしている主な要因として、振動位相緩和の存在が挙げられる。本研究では、位相緩和過程の代表例である分子内振動エネルギー再分配（Intramolecular Vibrational Relaxation = IVR）の初期過程において、一旦乱れた固有状態間の位相関係の回復を可能とするレーザーパルス列の設計を試みる。この際、2準位系動力学の解釈に使用されるブロッホベクトルモデルを多準位系に拡張し、フォトンエコ
ー等における位相回復手法の描像を設計の指針とする。
【理論】図１の様な Bixon-Jortner (BJ)モデル準位系を考える。
初期状態 i と光学遷移が許容である中間状態 m がレーザ
ー場 u (t ) と結合し、 m は等間隔 ε で分布しているバックグ{ } と強度 v で結合している。また、照射

ラウンド準位 B jするレーザー、 u (t ) = u0 cos[ωt + δ ] （初期位相 δ ）の角振
動数 ω は i ⇔ m 遷移に共鳴している。以下の手順に従って位相緩和制御を実現するパルス列を設計する。

図１《手順１：ブロッホベクトルの定義》

Bixon-Jortner モデル準位系[1]{ } を用いて、系の状態を Ψ (t ) = a i + ∑ c (t ) φ と表す。この時、 b = 1 − a 、

分子固有状態 φ j2

α j = arg[a j ] 、 β j = arg[c j ] を 用 い て 、 j

j

j

j

番 目 の ブ ロ ッ ホ ベ ク ト ル を

B j = ( ab cos[α − β i ], −ab sin[α − βi ], (a − b ) / 2 ) と定義する。
22《手順２：制御レーザーパルスの設計とその最適化》n 個のパルス列照射下における系の時間n
発 展 演 算 子 U ( n ) (t ) を 、 U ( n ) (t ) = U (t ')U (τ (i ,i −1) )U ( ∆t (i ) , δ (i ) ) と 定 め る 。 こ こ で 、
seq

seq

free

∏

i =1

free

opt

U opt (∆t , δ ) と U free (τ ) はそれぞれ、パルス照射中、パルス間の時間発展演算子である。また、∆t (i ) , δ (i ) は i 番目のパルス持続時間及び初期位相、 τ (i +1,i ) は i − 1 番目と i 番目のパルス間の時間 (t ) ≡ U ( n ) (t ) i で与えられる間隔を表す。初期状態 i にパルス列を照射した後の系の状態は Ψseq
ため、任意の目標状態 f との差 δΨ ≡ f − Ψ を関数 I (t ',"τ ( i ,i −1) , ∆t ( i ) ") = δΨ δΨ ≥ 0 で評価し、I を最小化することにより制御パルス列のパラメータ (t ',"τ (i ,i −1) , ∆t (i ) ") の最適化を実行する。
【結果】BJ モデル系のパラメータを v = 0.01 、 ε = 0.01 、 α = 0.0 、光学過程の相互作用強度
を u0 µ ge = 1 とした。上述のブロッホベクトルの定義によると状態 i は z 軸方向の単位ベクトル(0, 0,1) の集合に対応し、 i に π / 2 パルスを照射することによって生成する重ね合わせ状態
f ≡ ( i + m ) / 2 は (0, −1, 0) の集合として表現される（図 2(a)参照)。この状態 f は状態
m に起因する位相緩和を伴い、準位分布は π / 2 パルス照射後時間と共にバックグラウンド状態へ流出する。この過程は、ブロッホベクトル {B j } が xy 平面上に扇状に広がることに対応する（図 2(b)参照）。そこで、この状態に制御パルス列を照射し、 π / 2 パルスを照射直後の重ね合わせ状態 f を再生させることを試みた。１つ目の例である《タイプ１》のパルス列は、フォトンエコーに類似の制御描像に基づいたもので、図 2(b)の状態にパルス持続時間が 2π である矩形パルスを照射することによって実現される。この 2π パルスによって、扇状に広がったブロッホベクトルは x 軸まわりに 360 度回転するが、回転後の扇の広がり方は照射前と比較す
ると y 軸に対して対称に位置している。このため、 2π パルス照射後一定時間（図２(b)と同じ時間）経過したところで、ベクトルは再び収束し緩和過程が起こる前の状態 f を回復する。

m

m

i

(a)

{B }

{B }

{B }

j

j

i

i

m

i

(c )

(b)

m

j

(d )

図２ 《タイプ１》のパルス列照射：初期状態 i に順番に以下の演算子(a)～(d)を演算した場合の準位分布の時
間 変 化 （ 上 段 ） 、 及 び 各 ス テ ー ジ の 終 時 刻 に お け る ブ ロ ッ ホ ベ ク ト ル （ 下 段 ） ： (a) U opt (π / 2, 0) , (b)U free (π / 4) , (c) U opt (2π , 0) , (d) U free (π / 4)

《タイプ２：初期位相、パルス間隔を数値的に最適化したパルス列》
パルス列として、互いに逆位相mm {B }{B }j

である２つのパルスをパルス時

i

j

i

間幅 p、パルス間隔 f で照射し、その後 t ' だけ時間発展させる時

(a )

m

i
(c )

(b)

間 発 展 演 算 子 と し てU free (t ')U opt ( p, 0)U free ( f )U opt ( p, π )を 採 用 し 、 I (t ', p, f ) の 値 が 最小になる様に、パラメータを決定した。図３(a)~(f)に最適化さ{B }れたパルス列照射中の準位分布

jm及びブロッホベクトルの時間発m展を示す。最適化された２つのパルスによって xy 平面上に扇状

(d )

(e )

{B }

m

j

i

{B }
j

i

(f)

に集められたブロッホベクトル（図３(e)参照）は《タイプ１》の場合と同様に自由時間発展U (t ') によって再び収束し目
free

標状態の f にほぼ到達している。最適化の結果、目標状態への到達度は《タイプ１》よりも良好な結果が得られている。

図３ 《タイプ２》のパルス列照射：初期状態 i に順番に以下の演算子(a)～(f)を演算した場合（パルスシークエンス照射下）での準位分布の時間変
化 、 及 び 各 ス テ ー ジ の 終 時 刻 に お け る ブ ロ ッ ホ ベ ク ト ル ： (a)U opt (π / 2, 0) , (b) U free (π / 4) , (c) U opt ( p, π ) , (d) U free ( f ) , (e) U opt ( p, 0) , (f)U free (t ') ：(a)~ (f)は各時間発展の終時刻でのブロッホベクトル【参考文献】[1] M. Bixon and J. Jortner, J. Chem. Phys., 48, 715 (1968).

3E06
非平面キラル芳香分子(P)-2,2'-biphenol における光誘起電子運動の量子スイッチ（国立台湾大學 1，国立交通大学 2，東北大・理 3）峯尾 浩文 1，八巻 昌弘 2，寺西 慶
哲 2，林 倫年 1，林 聖賢 2，藤村 勇一 3

Quantum switching of –electron dynamics in optically induced
nonplanar chiral aromatic molecule (P)-2,2’-biphenol（ Institute of Applied Mechanics, National Taiwan University, Taipei 106
Taiwan1，Department of Applied Chemistry, Institute of Molecular Science andCenter for Interdisciplinary Molecular Science, National Chiao-Tung University,
Hsin-Chu 300 Taiwan2 ，Department of Chemistry, Graduate School, Tohoku
University, Sendai 980-85783 ）Hirofumi Mineo1, Masahiro Yamaki2, Yoshiaki
Teranisi2, Michitoshi Hayashi1, Sheng H. Lin2, and Yuichi Fujimura3[序] 芳香族分子のπ電子制御は、ナノエレクトロニクス・オプトエレクトロニクス領域において次世代の超高速スイッチングとして注目されている。これに関する量子化学研究については、Barth 等[1]はＭｇポルフィリンをＵＶ円偏光短時間パルス励起したときに生成されるπ電子環電流を計算している。この系では、電子回転は光の角運動量によっておこる。菅野等[2]は、ピラジン環をもつキラル分子の擬縮重状態
を直線偏光ＵＶパルスで励起することにより環電流を生成出来ることを示した。この場合、パルスは角運動量を持っていない。環電流は、２つの擬縮重電子状態によってつくられるコヒーレント状態が持つ角運動量によって生成される。いづれの分子も平面環状分子であり、生成されるπ電子の角運動量は分子面に垂直（Ｚ－軸）方向、その環電流ベクトルは一次元である。本研究では、二次元のスイッチングデバイスの可能性をもつ非平面分子(P)-2,2’-biphenol を取り上げ、ＵＶ直線偏光によるπ電子回転のパルス制御法を構築することが目的である。ここでは、その概略を述べる。光誘起π電子動力学の詳しい理論的展開は次の講演で述べる。
[本論] 図 1(a)は (P)-2,2’-biphenol のπ電子動力学を誘起する３つの電子状態（a,b1,b2）から、a と b1、及び b1 と b2 の
図１．
２つのコヒーレント状態が直線偏光パルスによってつくられることを示す。図 1(b)は関係する遷移双曲子モーメントを示す。 eb(1b)2 は b1 と b2 の線形結合が同位相で励起される、一方、eb(1b)2 は逆位相の線形結合がらを合成したベクトル(lX, lZ)も示されている。２つの芳香環のπ電子の位相を保った
回転には４種類がある。CA-type 回転はL(R) 環では clockwise(anticlockwise)回転を示す。孤立系であるのでπ電子が数回回転した後には逆回りの回転が起こる。
２次元カレント制御を目的として、ポンプーダンプ型の量子制御法で得られたπ電子
回転結果を図３に示す。ポンプパルスにより
電子回転コヒーレンスが生成し、逆回転が起こる前にダンプパルスにより電励起状態を消去する。ここでは、CC→AC→AA→CA の段階的制御を取上げた。図３(a)は目的とす
る角運動量 lZ と lX が数十フェムト秒内で得られたことがわかる。図３(b)には４段階の制御がそれぞれの段階でポンプとダンプパ
ルスがどのように働いていたかをみるために、パルス形を図示した。ポンプとダンプパ
ルスがお互い重なりあっている。これは、直線偏光のパルスが重なり円偏光的な性質をもち、電子の回転を加速させることを示している。
[結論] 非平面芳香環分子、(P)-2,2’-biphenol
は二次元π電子カレントを生成する。ポンプ
ーダンプ制御法によりカレントの量子

Z
x

L

X

Y

R

Z
x

L

X

X

Z

Z
lL

lR

R

lR

X

Y

lX

lL

lZ

iii) CA-type rotation(+)
by eab1 excitationi) CC-type rotationby e(+)
b1b2 excitation

Z
x

L

Y

R

Z
x

L

X

R

Y

X

Z lZ
lR
lL

l
lX L

X

X

Z

lR

ii) AA-type rotationby e(-)
b1b2 excitationiv) AC-type rotationby e(-)
excitation
ab1

Rotations in the same
direction

Rotations in the opposite
direction

図２．
(a)
0.1
lZ 0.0
-0.1
300

0.1
0.0
l X -0.1

200
100
0

e (ab1) e (ab1)(b)e (ab1) e (ab1)

1.0
1.0

E(t)/GVm-1生成されることを示す。生成条件は既に報告している[2]。
図２は、直線偏光パルス励起により左と右（L と R）フェノール環内のπ電子がお互いどのような相関をもって回転するかを示す[3]。L、R リング内のπ電子回転に対応する角運動量ベクトル(lL, lR)とそれ( ) ( )0.5 e
0.5
b1b 2 e b1b 2

e(b1b) 2e(b1b) 2

0.0
0.0
-0.5
-0.5
-1.0
-1.0

00

50

100
150
200
250
300
100
200
300
t/fs

図３．

スイッチできることが示された。分子振動効果や最適制御法による環電流・誘起磁場の制御・評価などは今後の課題である。
[参考文献 ]1.I. Barth, J. Manz, Y. Shigeta, J. Am. Chem. Soc.,128, 7043(2006).
2. M. Kanno, H. Kono, Y. Fujimura, Angew. Chem. Int. Ed., 45, 7995(2006); M. Kanno, H.
Kono, Y. Fujimura, S.H. Lin, Phys. Rev. Lett., 104, 108302 (2010).
3. H. Mineo et al., to be published.

3E07
線形紫外レーザーパルス誘起非平面キラル分子における π 電子回転ダイナミク
スの理論的研究： コヒーレント角運動量と環電流（国立台湾大學1，国立交通大学2，東北大・理3）峯尾 浩文1，林 聖賢2，藤村 勇
一2,3
Coherent π-electronic-state dynamics of chiral molecule induced byultrashort linearly polarized UV pulses: angular mometum and ring current( 1Institute of Applied Mechanics, National Taiwan University,
Taipei 106 Taiwan, 2Department of Applied Chemistry, Institute of Molecular Science
and Center for Interdisciplinary Molecular Science, National Chiao-Tung University,
Hsin-Chu 300 Taiwan, 3Department of Chemistry, Graduate School, Tohoku University,
Sendai 980-857)
Hirobumi Mineo1, Sheng H. Lin2, and Yuichi Fujimura2,3
In this work we theoretically show behaviors of ring currents in a non-planar
Chiral aromatic molecule with two phenol rings combined with a single chemical bond.
We remark that this talk is a continuous one of previous talk (3E06):
Quantum switching of –electron dynamics in optically induced nonplanar chiralaromatic molecule (P)-2,2’-biphenol, and use the same molecule.
The electrodynamics induced by a pulse laser E (t ) can be calculated bysolving the equations of motions of the density matrix:d  (t )
dt



i
 V (t )  (t )   (t )V (t )    i     (t )
 


 
where Vˆ (t )     E (t ) ,  is the transition dipole moment operator.  is angularfrequency difference between two electronic states  and. The angular momentum
and ring current are given as the expectation values in terms of density matrix elements:
O (t )  n  d   d 3 r1 d 3 rn  r  r1   Tr  (t )O (r1 ) ˆOˆ (r )
where n is the number of electrons, O (r )   O (r )   ,is the singleelectron operator. The angular momentum is calculated with the operatorOˆ (r )  i   x K  xK  y K  xK  n K with d   d 3 r , on the other hand the ringK  L, Re  
Oˆ (r ) 
current is given by2mei with d   dS where n is a unit vector
K

which is perpendicular to a ring K, dS represents an integral over a surface S.
The electronic angular momentum and the resultant ring current are created bya coherent excitation of two electronic excited states using the pulse laser. Theseelectronic states belong to C2 point group which have A and B representations. As issimilar to a previous talk, in this molecular system we have three electronic excited statea(A), b1(B) and b2(B), and by changing the laser frequency we obtain two coherently
excited states a,b1 (b1,b2) which belong to different (same) symmetries. The createdring current at each chemical bond oscillates back and force by a periodicity which is
given by the inverse of the energy gap between two coherently excited states.
()In Fig. 1 we show the directions of initially created currents by shining e ab1 ore (b1b) 2 where two coherently excited states are “in-phase”. From this figure we find thatcurrent or angular momentum transfer through the bridge chemical bond C1-C7 only
occurs if two coherently excited states belong to a different representation as is shownin Fig. 1(a). In this case the directions of ring currents and angular momenta on bothsides of phenol rings are toward opposite direction. On the other hand if two coherently
excited states belong to a same representation B shown in Fig. 1(b), the directions ofring currents and angular momenta on both sides of phenol rings are toward same
directions i.e., the created ring currents and angular momenta on both sides behave in asame manner.
(b) b1,b2 excitation(a) a, b1 excitationxx

Fig. 1 Behaviors of the ring currents are plotted where the two coherently excited statesare “in-phase”.

反応経路自動探索法 GRRM の特性評価と機能拡充（豊田理研 , 国立情報研 2, 東北大院理 3） ○大野公一 1、佐藤寛子 2、岩本武明 3

3E08

1

【序】ポテンシャル表面上に存在する平衡構造とそれらを結ぶ反応経路は、平衡点の周囲のポテンシャルの非調和下方歪み (ADD) に着目することによって芋づる式に探索することができ、GRRM 法として開発し実用に供してきた[1]。GRRM 法には他の方法にはない特色があり、未知の構造をも含む化学情報を個々の化学式について系統的に探索することが可能であるが、その基本特性について、これまで十分明確にしてこなかったので、今回、報告する。
さらに、GRRM プログラムをより高機能に改善する試みについてもその現状を報告する。
【方法】量子化学計算には Gaussian プログラム（g03, g09）を

表１ GRRM と Molgen で自使用し、ADD に着目する超球面探索を並列化して行うことが動探索された化学構造の数できる GRRM11 [2]を用いて、反応経路網の全面探索および限定探索（Large ADD Following 法のパラメータ LADD を使用）を行った。比較のため、原子価理論に従う構造異性体を Molgen
4.0 [3]を用いて求めた。
【結果・考察】表１に、GRRM(B3LYP/6-31G*)と Molgen によって求めた化学構造の数を示す。Molgen は、各原子が通常の原子価に従う構造異性体を迅速に探索することができるが、HCN の異性体としてよく知られている HNC や、CH3CN の異性体である CH3NC は、-NC 基が原子価を満たさないので、探索の対象にならない。また、Molgen では、回転異性体のような立体配座の異なるものも区別されない。
ニトロメタン H3CNO2 について、GRRM では 150 種類の異性体が得られるのに対し、Molgen による構造異性体は 15 種類である。ベンゼン C6H6 の異性体は、GRRM ではまだ探索中であるが既に 1400 種類
以上見出されているのに対し、Molgen よる構造異性体は 217 種類である。
GRRM で探索される構造の特徴を明らかにするめに、図１に H3CNO で探索された 30 種類の構造のエ
ネルギー分布を示す。
H3CNO について Molgen
で見出された構造は、

図１ GRRM で探索された H3CNO の異性体のエネルギー分布

NH2-CHO、HO－CH＝NH、HO－N＝CH2、CH3－NO、N(H)C(H2)O 三員環構の５種類であり、それぞれ、図１の GRRM によるエネルギー分布における構造番号、1、(2, 3, 4, 5)、(8, 9)、10、

12 に対応する。そのうち、最初の２つは、150 kJ/mol (1.6 eV) 以下の熱エネルギー領域（熱領域）にあり、(8, 9), 10, 12 は、150-300 kJ/mol (1.6-3.2 eV)の範囲の可視光に相当する光エネ
ルギー領域（光領域）にある。(2, 3,4, 5)は、いずれもビニルアルコールの＝CH2 基をイミノ
基＝NH で置換した構造に相当し、C=N 結合に対する OH 基と NH 基の向きが異なるだけである。8 と 9 は、OH 基の向きが異なるだけである。したがって、熱領域と光領域の 12 種類の構造のうち、Molgen でまったく予想できないのは(6, 7)と 11 だけである。(6, 7)はヒドロキ
シアミノカルベン HOCNH2 であり、C 原子が原子価を満たさない。また、11 は電荷分離したイオン構造 O=N＋H－C－H2 をもつ。このため、これらは、Molgen の探索対象からはずれている。したがって、GRRM では、Molgen で見出される通常の原子価を満たす構造を、熱領域や光領域などの低エネルギー領域において、多様な立体異性体を含めて見つけ出すとともに、Molgen では探索できないカルベンなどの不安定化学種や光領域を超える高エネルギー領域の化学種を多数みつけることができる。このことは、GRRM では、従来からの原子価理論では予想のつかない化学種を見つけ出すことができることを意味している。また、GRRM は、エネルギーの貯蔵に有用な高エネルギー物質の探索に優れていることがわかる。
GRRM では、反応経路網の全面探索ができるという重要な利点があるが、それには多大な計算時間を要するという問題点が付随する。しかし、低エネルギー領域を優先する限定探索オプションの LADD=n を指定することによって大幅に探索時間を短縮することができる[4, 5]。
図２に、LADD を用いない全面探索の結果と比べて、
LADD=3 として H3CNO の探索を行った結果、見つけられた異性体（赤色：Found）と見つけられなかった異性

図２ GRRM 法の LADD=3 による限定探索で見つけられた H3CNO の異性体(赤色：Found)のエネルギー分
布（白色部分は全面探索による）体（白色）のエネルギー分布(縦軸の単位は kJ/mol)を示す。
LADD=3 の場合、全面探索と比べて計算時間は約 12 分の 1 に短縮されているにもかかわらず、400 kJ/mol 以下の低エネルギー領域（光領域と熱領域を含む）の異性体はすべて見つけられていることがわかる。
【結論】GRRM は、原子価理論に従う経験的な予想では見逃されるような未知の化学種やエネルギーの貯蔵に有用な高エネルギー化学種の探索に優れている。比較的低エネルギーの光領域や熱領域の構造に限定すれば、全面探索と比べて 10 分の 1 以下に探索時間を短縮することができる。GRRM の機能拡充に関する新しい取り組みについては当日報告する。
[1] K. Ohno, S. Maeda, Chem. Phys. Lett. 384, 277 (2004); S. Maeda, K. Ohno, J. Phys. Chem. A 109,5742 (2005); K. Ohno, S. Maeda, J. Phys. Chem. A 110, 8933 (2006).
[2] 大野公一, 長田有人, 前田理, 諸熊奎治, 第 14 回理論化学討論会, 岡山 (2011), 2D1b.
[3] Molgen 4.0, http://molgen.de/?src=documents/molgenonline[4] 大野公一、第 5 回分子科学討論会、札幌 (2011), 1P122.
[5] 大野公一、豊田報告 No65, 21 (2012)、http://www.toyotariken.jp/activities/65/F03_Ohno.pdf

3E09

超球面探索法に SCC-DFTB 法を用いた
炭素原子クラスター及び炭化水素の異性体の探索（和歌山大院・システム工 1，和歌山大・システム工 2，北大院・理 3，豊田理研 4）
時子山 宏明 1，山門 英雄 2，前田 理 3， 大野 公一 4
Searching for isomers of carbon atom clusters and hydrocarbons by SCC-DFTB Methodwith the Scaled Hypersphere Search Method(Graduate School of Systems Engineering, Wakayama Univ.1; Faculty of Systems
Engineering, Wakayama Univ.2; Faculty of Science, Hokkaido Univ.3; Toyota Physical &Chemical Research Institute4)
Hiroaki Tokoyama1; Hideo Yamakado2; Satoshi Maeda3; Koichi Ohno4【序】
筆者らは 2004 年に大野・前田により開発された超球面探索法(SHS 法)1 と分子軌道法を結晶構造の予測に適用してきた。SHS 法を用いた異性体探索において、半経験的方法である
Self-Consistent-Charge Density-Functional Tight-Binding（SCC-DFTB）法 2 を用いている。
以前、ランダムに配置した炭素原子 20 個を初期構造とした C20 異性体の自動探索を行った。
Ab initio 法の代わりに SCC-DFTB 法を用いる事で高速に計算が可能となり、最安定構造であると予想されるコラヌレン様構造を含む約 2000 個の異性体が自動的に得られ、エネルギー順に並べると滑らかな曲線が得られた。SHS 法とは、非調和下方歪(ADD)を辿ることで、固有反応座標(IRC)上の平衡構造(EQ)と遷移構造(TS)を自動的に探索する方法である。
今回、炭素クラスター(Cx)及びベンゼン（C6H6）を含む炭化水素(CxHy)について報告する。
【方法】
ランダムに配置した炭素原子及び炭化水素を初期構造とし、各異性体の自動探索を行った。
DFTB 計算では Self Consistent Charge（SCC）オプションを使用した(SCC-DFTB)。SHS 計
算では ADD の大きい経路を優先して探索する l-ADDf(large-ADD-following)法（LADD=5）を用い、SCC-DFTB 計算は dftb+プログラム 3 を用いてそれぞれ行った。
【結果と考察】
ランダムな初期構造から探索を開始し、炭素クラスターCx、炭化水素 CxHy（y=x, 2x, 2x+2）の各々の見つかった異性体の数を（縦軸を指数表示にして）図 1 に示す。Cx では炭素原子の増加に伴い指数関数的に増加している。超球面探索に非常に多くの勾配計算を行うため、ab initio
計算と組み合わせた場合は莫大な計算時間がかかってしまうが、ab initio 法の代わりに
SCC-DFTB 法を用いると大幅に高速化され、C4H4 の場合約 12 時間で計算が終了した。
C6H6 では最安定構造と予想されるベンゼン(EQ139)やデュワーベンゼン(EQ276)、フルベン(EQ11)及び 4 員環を有する分子等を含む TS との反応経路について図２(c)に示す。C20 の場合同様に連続的な曲線が得られ（図２(a)参照）、急に EQ 数が増加している 400kJ/mol より下に存在する EQ の数は探索された全体の 7%程度であることが分かった。
謝辞：本計算で、自然科学研究機構
利用しており、感謝する。

岡崎共通研究施設

計算科学研究センターの電子計算機を

図１ 炭素クラスター、炭化水素の異性体数（横軸:分子内炭素原子数、縦軸:異性体数の指数表示）(a) Cx(x=3～10, 20)の異性体数、(b) CxHx(x=4～6)の異性体数、(c) CxH2x(x=2～4)の異性体数、(d) CxH2x+2(x=1～4)の異性体数

図 2 C6H6 の異性体探索(a) EQ 点をトータルエネルギー順に並び替えた EQ色枠）の拡大図(b) 左上図の 400kJ / mol 以下の部分（赤(c) 最安定構造(EQ139: ベンゼン)の周囲の反応経路

参考文献：
1) K. Ohno and S. Maeda, Chem. Phys. Lett., 2004, 384, 277; S. Maeda and K. Ohno, J. Phys. Chem.
A, 2005, 109, 5742; K. Ohno and S. Maeda, J. Phys. Chem. A, 2006, 110, 89334
2) M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, Th. Frauenheim, S. Suhai and G. Seifert,
Phys. Rev. B, 1998, 58, 7260
3) B. Aradi, B. Hourahine and Th. Frauenheim, J. Phys. Chem. A, 2007, 111(26), 5678