4P101

紫外領域における、サンスクリーン分子光安定性の理論的研究(上智大学 1 , 個人 2 ) ○島田嶺太 1 , 石田俊正 2, 南部伸孝 1

Theoretical evaluation of photostability for Sun Protect molecule in UV energy region(Sophia University1 , Individual2 ) ○Ryota Shimada1 , Toshimasa Ishida2 , Shinkoh Nanbu1【序】 地球に降り注ぐ紫外線は，三つの領域 UVA(400～315 nm)，UVB(315～280 nm)，UVC(280～10nm) に分類される。特に化粧品科学では，UVA と UVB を対象に人体に有害な紫外光を散乱あるいは吸収することにより皮膚組織を保護する機能分子を，サンスクリーン分子と呼び日焼け止めなどの主成分として用いられている．後者は紫外線吸収剤とも呼ばれており，厚生労働省が認定する幾つかの分子のみ
産業にて応用される．その中で，代表的な分子が，octyl-metoxy-cinnamate (OMC)1)であり，特に市販の日焼け止めには，広く使われている．さらに，その分子の紫外線吸収性能を，化粧品科学では SP 因子（Sun Protect Factor）で評価している．しかし，評価時に質量を基にした密度で規定した分子数を基にその因子を判定している．つまり，紫外線防御の効Fig. 1 c -MMC

果を厳密に検証していくために，最小単位である１分子あ

るいは１モル単位での値を比較すべきと思われる．本研究では，サンスクリーン分子の紫外線防御効果を厳密に解明するために１分子当たりの光吸収スペクトルを第一原理に基づく理論計算を行い，光異性化反応およびフォトクロミック過程を解明する．特に，オクチル基をメチル基とした
cis-methyl-metoxy- cinnamate(c-MMC)をモデル分子(Fig. 1)として取り上げる．c -MMC の光吸収後の遷移過程を明らかにすることを目的とする．
【 計 算 方 法 】 光 吸 収 後 の 無 輻 射 過 程 を 明 らかに す る た め ， Zhu-Nakamura （ ZN ） 理 論 に 基 づ く
Trajectory Surface Hopping（TSH）法による非断熱分子動力学計算を実施した．この時，ポテンシャルエ
ネルギー，その傾き，遷移双極子モーメントは，古典軌道を時間発展させる際，毎時間ステップ毎に電子状態計算を実施した．用いた電子状態計算法は，完全活性 SCF（CASSCF）法であり，基底関数には
Dunning らの cc-pVDZ を用いた。
8（S1, S2）上で最安定構造を取るポテンシャルエ
ネルギー値（S0 opt, S1 opt, S2 opt）を Fig. 2 に示
す．S0 opt からの垂直励起において，S1 と S2 は近接しており，遷移双極子モーメント d も大きな差はなかった．よって， S1←S0 および S2←S0 の

Potential Energy / eV

【結果と考察】 基底状態（S0 ）および励起状態

7

S2

6

S1

S2 optS1 opt5
ΔE= 6.7 eV
d2←0 = 1.8 a.u.

432

ΔE= 5.7 eV
d1←0 = 2.6 a.u.

1
0

S0 opt二つ光遷移が許容遷移であると考えられる．そ-1

こでそれぞれの場合において非断熱分子動力Fig. 2 S0 opt, S1 opt,グラフ
S2 optのエネルギー値とそれぞれ
タイトルの垂直遷移エネルギー値

学計算を行い，光吸収励起後の無輻射過程の理論的解明を行った．

0

1

2

3

4

5

6

7

8

9

S 1 ← S0 励 起 後 ， S 2 ← S 0 励 起 後 のtranectory 例を，それぞれ Fig. 3-1, Fig.
4-1 に示す．両者共に，励起後 200 fs
付近までは，S1 と S0 のエネルギー差が5 eV 程と大きい．しかし 200 fs 付近を過ぎた当たりから，徐々に S0 エネルギーが増加し始め，S1 エネルギーと近接していった．その結果，300 fs 付近で
S1→S0 非断熱遷移を観測した．
この無輻射過程における，分子の幾何構造変化の影響を調べるために
c-MMC 分子中のエチレン基における二面角θ(Fig. 5)の時間変化に着目した(Fig. 3-2, Fig4-2)．
Fig. 5 二面角θ

θ

S0 とθの時間変化を比較すると，両者の増減が，非常に同調していることが分かった．200 fs 付近までの S0 エネル
ギーが低い領域においては，分子はθが小さい，つまり cis 型の構造をとっている．200 fs 付近から始まる S0 エネ
ルギーの増加に対して，θも同調して増加し，分子は trans 型に近い構造をとるようになる．このことから，c-MMC 分子の無輻射過程は，分子の cis-trans
相互幾何異性化変化により，制御されていることが示唆された．
【まとめ】 以上，モデル分子である c-MMC の光吸収後の遷移過程を TSH 法により解明した．S1 と S2 への励起を考慮したが，両者ともに S0状態への無輻射過程を起こした．また，分子の幾何構造の変化が，同無輻射遷移過程に大きく影響を与えていることが示唆された．なお今回は，MMC の cis 異性体に着目したが，trans 異性体にも同様の計算を行い，当日の発表で両者の比較を行う．
1) M.M Jiménez , International Journal of Pharmaceutics , 272 , 1–2 ,19 March 2004 , 45-55

4P102グアニンの光励起状態緩和過程に関する理論的研究(上智大学大学院理工学専攻 1，個人 2）○伊藤 亮一 1，石田 俊正 2，南部 伸孝 1
Theoritical study of photostability for guanine in UV（Graduate School of science and technology,Sophia Univ.¹,Individual²）○Ito Ryoichi¹,Isida Toshimasa²,Nanbu Shinko¹【目的】DNA や RNA を構成している核酸塩基は、200～300 nm の紫外光に対して吸収帯を持っている。
しかし、これらが紫外光にさらされても、突然変異
やガンを誘発する光化学反応は起こりにくい。つまり、核酸塩基は光に対して安定である。この様な性質を持つ核酸塩基に対して多くの実験・理論研究がなされている。その一つとして、DNA 塩基の一つであるグアニンでは、共鳴二光子イオン化分光法

Fig. 1. 7H-keto-amoino-guanine(R2PI 法)により、エノール体のスペクトルピークは最安定構造観測されるがケト体のスペクトルピークがほとんど観測されないという結果が得られている[1]。このことから、グアニンはケト体にのみ超高速の励起状態緩和過程が存在しているのではないかと考えられている。また、グアニンの励起状態緩和時間はポンプ-プローブ紫外可視時間分解分光によって 148 ～360 fs と求められている[2]。 しかしながら、この超高速励起状態緩和メカニズムの詳細については未だ明らかでない。そこで、本研究では理論計算によるケト体グアニンの超高速励起状態緩和メカニズムの解明を目的とした。
【理論計算】グアニン互変異性体のうち、最も安定であるケト体(7H-keto-amoino-guanine;G7K)を対象として理論計算を実行した。計算方法には二次の多配置参照摂動展開法(CASPT2)を、基底関数には cc-pVDZ を用いた。
はじめに、基底(S₀)状態の構造最適化

S₂を行い、最安定構造を求めた。次に、得られた最安定構造の座標において調和振動子解析を行い、乱数発生を利S₁用して分子の初期座標および運動量を与えた。さらに、ハミルトンの運動方程式を数値的に解くことで古典軌
道の時間発展を実行した。その際に、

S₀

Zhu-Nakamura 公式に基づく古典軌道
ホップ法(ZN-TSH 法)を用いて、非断熱遷移の有無を考慮した。本研究では、
構造最適化の際に算出した遷移双極子モーメントの値を考慮し G7K を第Fig. 2. S₂状態へ励起後ポテンシャルエネルギーの時間変化

二励起(S₂)状態および第一励起(S₁)
状態に垂直励起させた後の無輻射過程を追跡した。

S₁【結果と考察】G7K の基底状態の最安定構造を Fig. 1 に示す。また、Fig. 2と Fig. 3 に ZN-TSH 法により得た 7
個の古典軌道に対するポテンシャル
エネルギーの時間変化を示した。
Fig. 2 は S₂状態に垂直励起させた後S₀のポテンシャルエネルギーの時間変化、Fig. 3 は S₁状態に垂直励起させた
後のポテンシャルエネルギーの時間

Fig. 3. S₁状態へ励起後ポテンシャルエネルギーの時間変化① S₂状態へ垂直励起後の挙動
時間変化
S₂→S₁非断熱遷移は、S₂状態への励起後、約 50fs までに起こった。Fig.2
からわかるよ変化を示している。

うに、S₂状態と S₁状態のポテンシャルエネルギーの差が約 1 eV と近接しているため、容易に非断熱遷移を起こすと考えられる。また、S₁→S₀非断熱遷移は約 120～300 fs までに起こった。特に光励起後から時間の経過に伴い、S₀状態のポテンシャルエネルギーが不安定化し、S₁状態ポテンシャルに近接し、非断熱遷移が起こりやすくなる傾向がみられた。
② S₁状態へ垂直励起後の挙動
これまでの計算において、最も時間が経過した計算結果は 93 fs までである。しかし、まだ S₁→S₀非断熱遷移は起こっていない。その一方、S₂状態へ励起後の挙動と同じように時間が経過するにつれ S₁状態と S₂状態のポテンシャルエネルギーの差が徐々に縮まっていることがわかる。したがって、S₁状態へ遷移後の緩和過程も S₂状態へ遷移後の緩和過程と類似するものと予想される。つまり，数百 fs 程度で内部転換が起こり、S₀状態へ遷移するものと思われる。しかしながら。二つの一重項状態間に三重項状態が存在し、三重項状態へ遷移後、S₀状態へ遷移する過程も考えられるため，スピン・軌道相互作用が引き起こす交間交差の可能性を探る必要がある。
以上の①、②から G7K が第二励起状態にとどまる時間は約 50 fs 以内と予想される。この結果はケト体である G7K が共鳴二光子イオン化分光法(R2PI 法)によって観測されにくいことに一致している。また、S₂状態へ励起後は約 300 fs 以内に S₀まで失活することが分かった。
S₁状態へ光励起については，S₀状態へ遷移まで計算が進んでいない。今後，計算を進めるとともに，スピン・軌道相互作用が引き起こす交間交差の可能性を探る予定である。
【参考文献】[1]

H. Saigusa, J. Photochem. Photobiol. C: Photochem. Rev., 7, 197(2006) .

[2]

C. Canuel, M. Mons, F. Piuzzi, B. Tardivel, L. Dimicoli, and M. Elhaninea J. Chem. Phys. 122,
074316 (2005)

4P103
実時間発展 TDHF/TDDFT 計算の効率化（上智大理工）○赤間 知子、佐藤 紀穂、南部 伸孝
Real-time TDHF/TDDFT calculation with efficient time evolution(Sophia Univ.) ○Tomoko Akama, Kiho Sato, Shinkoh Nanbu【緒言】光によって誘起される電子ダイナミクスは、最先端の実験で観測されつつある超高速現象であり注目を集めている。特に、分子系において光誘起電子ダイナミクスが生じる場
合には、分子デバイス材料への応用の可能性が期待されている。このような光誘起機能性分子を設計するためには、まず光誘起電子ダイナミクスのメカニズムを解明する必要がある。
このための有用な手段の一つとして、非経験的な理論計算を用いた取り組みが挙げられる。
最近では、時間依存 Hartree-Fock (TDHF)法や時間依存密度汎関数理論(TDDFT)の実時間発展(RT)計算を用いた、混合状態の波動関数による電子ダイナミクスの記述が可能になりつつある[1-3]。しかし、時間発展演算子の計算には高いコストを要することが多く、これまで適用が限られていた。RT-TDHF/TDDFT における時間発展演算子の計算には Taylor 展開法、Runge-Kutta 法、split-operator 法等が用いられてきたが、細かい時間刻みで多くのステップを計算するため、高いコストを要することが多い。これらの方法よりも低い計算コストで高精度の結果を得られる方法として、Chebyshev 展開法がある。これをさらに改良した実空間
Chebyshev 展開法[4,5]は、核波束の効率的な時間発展法として開発された。この方法により核波束ダイナミクスの計算コストが大幅に減少し、長時間のダイナミクスの追跡が可能になっている。本研究では、電子のダイナミクスを記述する RT-TDHF/TDDFT 計算に実空間
Chebyshev 展開法を適用することにより、計算を効率化することを目指す。
【RT-TDHF/TDDFT】TDHF または TD Kohn-Sham (KS)方程式は次のように表される。
i i (r, t )  Fˆ (r, t )  Vˆ ext (r, t )  i (r, t )
t





(1)

ここで、 i (r, t ) は分子軌道または KS 軌道、 F̂ は(KS) Fock 演算子、 Vˆ ext は外場との相互作用である。式(1)を変形すると、密度行列 D(t ) に関する運動方程式が得られる。
iD(t )  [F(t )  V ext (t ), D(t )]t(2)Fij (t )  i Fˆ  j
1ZA i  2  
j2
A r  RA



(3)



HF
DFT XCFij (t )  Dlk (t ) i j | kl   c Dlk (t ) i l | k j   ckl

Vijext (t )  i Vˆ ext  j   (t )eext  i r  j(4)

ここでは規格直交な基底{ϕi}を用いており、F(t ) は Fock 行列、[F(t )  V ext (t ), D(t )] は F (t )  V ext (t )と D(t ) の交換関係を表す。 V ext (t ) は系と外場との相互作用を表す行列であり、双極子近似を用いている。分極ベクトル P(t ) は密度行列 D(t ) を用いて定義され、これを Fourier 変換することにより、系の励起スペクトル Im[ P ( )] が得られる。
P(t )   i rˆ  j D ji (t )

(5)

ij

【時間発展演算子の計算】RT-TDHF/TDDFT では、式(1)または式(2)の微分方程式を数値的に積分する、すなわち実時間発展させることにより解く。本研究では、式(2)に対して実空間
Chebyshev 展開法を適用した。 V ext (t )  0 となる t の範囲に対して、演算子 L̂ を用いて密度行
列に関する運動方程式（式(2)）を次のように書き直す。
i
D(t )   F(t ), D(t )   Lˆ D(t )t(6)

演算子 L̂ の時間依存性が無視できると仮定すると、式(6)の形式解として、次式が得られる。
ˆ )D(0)D(t )  exp(iLt(7)

Chebyshev 展開法では、時間発展演算子を Chebyshev 多項式 Tk(x)で展開して計算する。
ˆ )   (2   )(i ) k J (t )T ( Lˆ )exp(iLt
k0
k
k

(8)

ˆ ( Lˆ )  T ( Lˆ ) , T ( Lˆ )  1, T ( Lˆ )  Lˆ
Tk 1 ( Lˆ )  2 LT
k
k 1
0
1

(9)

k 0

ここで、Jk(t)は第 1 種 Bessel 関数である。一方、実空間 Chebyshev 展開法では、新しい演算子
f ( Lˆ ) を定義し、これを含む密度行列の運動方程式に対して時間発展の計算を行う。
i1
D(t )  f ( Lˆ )D(t )  cos 1 (as Lˆ  bs )D(t )
t
t





Dk 1  2 as Lˆ  bs Dk  Dk 1

(10)
(12)

ここで、as 及び bs はスケーリング係数であり、 1  as Lˆ  bs  1 となるように決める。
本研究では RT-TDHF 法および RT-TDDFT を量子化学計算パッケージ GAMESS に実装した。
時間発展の方法として実空間 Chebyshev 展開法を用いた場合の計算コストと計算精度について検証し、4 次の Runge-Kutta 法や従来の Chebyshev 展開法との比較を行う。
[1] K.Yabana and G. F.Bertsch, Phys. Rev. B 54, 4484 (1996); K. Yabana, T. Nakatsukasa, J.-I.
Iwata, and G. F. Bertsch, Phys. Status Solidi B 243, 1121 (2006).
[2] C. Y.Yam, S.Yokojima, G. H. Chen, J. Chem. Phys. 119, 8794 (2003); F. Wang, C. Y. Yam, G.
H. Chen, and K. Fan, J. Phys. Chem. 126, 134104 (2007).
[3] T. Akama and H. Nakai, J. Chem. Phys. 132, 054144 (2010).
[4] S. K.Gray and G. G. Balint-Kurti, J. Chem. Phys. 108, 950 (1998).
[5] G. G. Balint-Kurti, Theor. Chem. Acc. 127, 1 (2010).

4P104周期境界条件の ONIOM 法への適用（上智大院・理工）○小林 理，南部 伸孝

Application of Ewald Sum for ONIOM method(Graduate school of Sophia Univ.) ○Osamu Kobayashi, Shinkoh Nanbu【序】諸熊らにより考案された ONIOM (our Own N-layer integrated molecular Orbitals andmolecular Mechanics)法は量子化学計算と分子力場法を統合し、巨大系に対する計算を可能にする手法として知られている[1]。ONIOM 法は系の全エネルギーを以下のように分割する。

EONIOM = EMM,real + EQM,model − E MM,model
ここで EMM,real 、 EQM,model 、 EMM, model はそれぞれ系全体に対する molecular mechasnics (MM)法によるエネルギー、model 系に対する quantum mechanics (QM)法によるエネルギー、model 系
に対する MM 法によるエネルギーである。この方法は mechanical embedding 法とよばれる。また同じく諸熊らにより QM 層と MM 層の相互作用を適切に取り入れるために、MM 層の電荷を点電荷として QM 層の計算に取り込む electronic embedding 法が開発された。この手法では MM層の原子と QM 層の電子の相互作用は、1 電子演算子で表現される。
MM が十分に大きな系であるとき、 EMM, model の計算を有限の大きさで近似して行うと、表面が生じることに起因し問題が生じる。このような問題は例えば溶液や結晶においてみられる。これは周期境界条件(PBC)の導入により解決される。Coulomb エネルギーは原子間距離に対して収束が非常に遅いことが知られており、PBC においては 2 つの収束の速い級数に分割し合計する
Ewald sum 法や計算を効率化した particle mesh Ewald sum 法(PME)が広く用いられている。
我々はこれまで非断熱現象を取り扱うために ONIOM 法を Tully’s surface hopping (TSH)法と組み合わせた QM/MM-TSH 法による分子動力学(MD)計算を行ってきた[2]。しかし、従来の方法では有限系の取扱いに限られており、系の周囲の環境が適切に取り込まれていない可能性があった。そこで
我々は ONIOM 法に Ewald 法を組み合わせた ONIOM-Ewald 法を開

NH2

非断熱
現象

発し、さらにそれを QM/MM-TSH 法に適用した。本発表では非断熱現象により励起状態における cis 型から trans 型への異性化が知られている 2,4-pentadiene-iminium cation (protonated Schiff base 3,
PSB3)のメタノール中の MD 計算を行い、本手法の効果を調べる。
【方法】シミュレーションを行う系の構造と分子力場は以下の手順に

NH2

図 1. PSB3 における非断熱的な構造変化より作成した。はじめに PSB3 を Gaussian 09 を用いて構造最適化を行った。構造最適化は密度汎関数法により行い、密度汎関数は B3LYP、基底関数は cc-pVDZ を用いた。引き続き AMBER9プログラムパッケージに付属する ANTECHAMBER を用い restrained electrostatic potential(RESP)電荷を決定した後、メタノールを PSB3 の周囲に配置した。得られた系の分子力場はgeneral AMBER force field (GAFF)を用いた。作成した系は AMBER9 の Sander モジュールを

使い 1 atm, 300K での NPT アンサンブルで 1 ns 時間発展させ、次の計算に用いた。
構造変化を調べるために、QM/MM-TSH 法を用い励起状態 MD シミュレーションを行った。
シミュレーションの開始時点において、各原子の初期速度はボルツマン分布に従いランダムに発生させ、PSB3 の電子状態は S1 にあると仮定した。シミュレーションは NVE アンサンブルで行い、タイムステップは 0.50 fs とした。遷移確率の計算には朱-中村理論を用いた。各時刻においてエネルギーと各原子にはたらく力は Gaussian 09 を用い ONIOM-Ewald 法に従い求めた。
PSB3 を QM 層に、その他の分子を MM 層に割り当てた。
8

ONIOM 法は単位セルに関して electronic embedding 法を6用い、外側の像に関しては mechanical embedding 法を用いた。QM 層の計算は Molpro 2012 を、MM 層の計算は当研究室で開発したオリジナルプログラムを用いた。MM 層の計算で用いた分子力場は先に初期構造を得るために用Potential Energy / eV

この割り当てはシミュレーション中変更していない。

(a)

4

S1

2
0

いた分子力場と同じものを用いた。また QM 層の計算は-2

CAS(6e,6o)SCF 法を用い、基底関数には MIDI-bang を用

-4
-4

S0

00

いた。比較のために周期境界条件を適用しない electronic

200
200

400
400

embedding 法による MD 計算もあわせて行った。この計
8

ONIOM-Ewald 法を用いたときと同様である。いずれの場

6

対して行い、S0 状態への非断熱遷移が確認された後はシミュレーションを停止した。
【結果と考察】はじめに分子の各層への割り当てを固定したことの妥当性は、このシミュレーションにおいては単位セルから外部への溶媒分子の拡散はほとんど観察されな

Potential Energy / eV

算の詳細は外場の取り扱いを変更した以外は合もシミュレーションは現在 100 本のトラジェクトリに

600
600

800
800

1000
1000

time / fs

(b)

4

S1

2

S0

0
-2
-4
0

かったことより支持される。これはシミュレーション時間

200

400

600

800

1000

time / fs

が溶媒の拡散速度に対して十分に小さかったためである。
図 2 にシミュレーションに用いた全てのトラジェクトリに図 2. PSB3 の S0 状態及び S1

ついて PSB3 のポテンシャルエネルギーの時間変化を示し

状態のポテンシャルエネルギーた。図 2 において(a)は 1000 fs までグラフが続き非断熱遷(a) ONIOM-Ewald 法による移せずに S1 状態を保つトラジェクトリがあることを示す(b) cut-off による一方、(b)は 400 fs 付近においてグラフが途切れ全てのト(実線が実際の電子状態を表す)

ラジェクトリにおいてこの時刻までに S0 状態に非断熱遷移したことを表している。これは ONIOM-Ewald 法を用いた場合、非断熱遷移が遅くなることを意味している。Muñoz-Losa[3]らは、非断熱遷移が無輻射遷移に要する時間は溶媒の運動に影響されると考察しており、我々の結果は Muñoz-Losa の結果を支持している。
[1] T. Vreven et. al. J. Chem. Theory Comput., 2006, 2, 815-26[2] W. C. Chung et. al. J. Phys. Chem. B, 2012, 116, 8009–8023[3] A. Muñoz-Losa et. al. J. Chem. Theory Comput., 2011, 7, 4050-4059

4P105Stieltjes モーメント法による核酸塩基の光イオン化立体ダイナミクスの研究(京大院理，理研) ○鈴木喜一，鈴木俊法

Photoionization dynamics of nucleobases studied by Stieltjes moment theory(Kyoto Univ, RIKEN) ○Yoshi-Ichi Suzuki and Toshinori Suzuki【序論】核酸塩基などの大きな分子の光電子スペクトルでは，多数の分子軌道がエネルギー的に接近しているため，スペクトルの帰属は困難である。そのような場合，おおよその帰属には量子化学計算を利用し，対称性が異なる近接した軌道には角度分解光電子光法 1, 2 が有効であるとされる。しかし，我々はイオン化確率（積分断面積の偏光依存性）に軌道対称性の違いが最も明瞭に表れることを指摘した。1, 3 た
だし，その際に用いた多重散乱法では，イオンの電子状態の違いが計算に反映されにくく，交換相互作
用も局所近似でしか取り込めない欠点があった。そこで，本研究では新たに Stieltjes モーメント法を用いて，軌道対称性による違いを検討した。
【理論】 光イオン化過程で究極的な観測量は，分子固定系での二重微分断面積 d 2σ / dRˆ / dk である。4Rˆ = {ϕ , ϑ , χ } は分子の向きを表すオイラー角，k は光電子の運動量である。直線偏光による一光子イオン化において，微分断面積を実験室系での分子の空間配向について積分すると，下記のルジャンドル多項式で表現できる。

dσ σ 0 ( E )[1 + β ( E ) P2 (cosθ )], E = = k 2
=
dk
4π
2me

(1)

θ は偏光と k ベクトルのなす角，E は光電子運動エネルギーである。β(E) ( − 1 ≤ β ( E ) ≤ 2 )は非対称パラメータと呼ばれる。図(a)に，核酸塩基の価電子軌道のイオン化における β(E) の実験値を示した。σ軌道
よりπ軌道のβ(E)が大きい傾向が認められるが，違いは小さい。エネルギーを変えてもσとπの差が大きくはならないことが，我々の多重散乱法による計算から予測されている。一方，微分断面積を光電子の運動量方向について積分し，かつ分子の Z 軸周りの角度(χ)で積分すると

dσ σ 0 ( E )[1 + (3ρ ( E ) − 1) P2 (cosϑ )]
=
4π
dRˆ

(2)

となる。 ϑ は偏光と Z 軸の成す角，ρ(E) ( 0 ≤ ρ ( E ) ≤ 1 )がイオン化確率の異方性を表す量である。 [実験
1.0

___
__
_ __
____
__
_
__
____
_
_
_
_
_
σ orbitals

0.0

π orbitals

_
_
_
_
_
__
__
_ __
__
ura
thym cil
cyto inep sineadeurinen
gua inenine

ρ(E) at PKE = 5 eV

(b)

ura
thym cil
cyto ines
pu ineade rinen
gua inenine

nine

π orbitals

ade

nine

thym

-1.0

σ orbitals

ade

0.0

ine
cyto
sine

__ __
1.0
_ _ __
_
_

ine
cyto
sine

β(E) at hν = 40 eV(a)__ __ _
_ _ _
__

thym

2.0

E
図：(a)核酸塩基の光電子非対称性パラメータの実験値(Trofimov
et al。 E2006)。
(b)多重散乱法による核酸塩基および，プリンの原子価軌道の光イオン化確率異方性の計算値。

的には，時間分解光電子分光法 5 あるいは光イオン蛍光偏光法などで ρ(E) を決定できる。] 多重散乱法による計算では，ρ(E) は明瞭に軌道対称性に依存することを示した[図(b)]。

β(E) の計算には連続状態の波動関数が必須である。一方， ρ(E) の計算では部分波間の干渉は寄与しないため， Stieltjes モーメント法 6 で振動子強度を計算するだけで済む。本研究では，HF/IVO 法 7
で振動子強度を求めた。 ρ(E) は，直線偏光による分子固定系での光イオン化積分断面積で次のように表現できる。

ρ (E) =σ Z (E)σ X (E) + σ Y (E) + σ Z (E)(3)

ここで， σ μ は偏光軸が分子内座標のμ軸と平行である場合の光イオン化積分断面積である。例えば，

ρ ( E ) = 0 は Z 軸に平行な双極子では全くイオン化しない場合， ρ ( E ) = 1 は Z 軸に平行な双極子のみでイオン化する場合である。
【結果】 DNA および RNA 塩基の HOMO から HOMO-7 までのに Rydberg 基底を 10s10p10d10f だけ加えたものである。
E > 3 eV の高エネルギー側では多重散乱法の結果同様にρ(E)σ orbitals

がほぼ定数であることが確かめられた。5 eV における計算結果を図(c)に示した。多重散乱法と比べると，π軌道の ρ(E) が全体的に 0.1 程度低下しているが，σ軌道の ρ(E) とは明らかに区別できる。一方，E < 3 eV の領域では π*型の形状共鳴状態があり，図(b)のようなσとπ軌道の ρ(E) における明瞭な違いはない。この点も多重散乱法と同じである。

0.0

π orbitals

_ __
_
_
_____
ura
thym cil
cyto ine
sine
ade
nine
gua
nine

軸とした。HF/IVO 計算の基底関数は uncontracted aug-ccpVTZ

_
_
_
_
_
__
____
_
__
____
_

ura
thym cili
cyto nesine
adengua inenine

もとで最適化[MP2/6-311++G(dp)法]し，対称面に垂直な座標を Z

ρ(E) at PKE = 5 eV

ρ(E) を 0 < E < 13 eV の範囲で計算した。分子構造は Cs 対称性の1.0E

図(c) HF/IVO/Stieltjes モーメント法
による，核酸塩基の原子価軌道の光イオン化確率異方性の計算値。

【議論】 Cs 点群では，一電子散乱波の対称性を a’と a”に区別できる。軌道対称性による ρ(E) の違いはσ(a’)，π(a”)いずれの軌道も a’状態へのイオン化が優勢であることが原因である。選択則からσ軌道は X, Y 方向の双極子で，π軌道はZ 方向の双極子でイオン化するので ρ(E) の違いが説明できる。a’状態が優勢な
理由は，形状共鳴が考えられる[図(d)]。多重散乱法の計算結果は，多くのσ*状態が 3 < E < 10 eV にあることを示しており，今回の計算でも確かめられた。σ*型
形状共鳴の位置と選択則から許される遷移は図(d)のようにまとめられる。以上より，
多重散乱法の計算結果が Sieltjes モーメント法によっても支持された。
1

Y.-I. Suzuki and T. Suzuki, J. Phys. Chem. A 112, 402 (2008).

2

A. Trofimov, et al., J. Phys. B: At. Mol. Opt. Phys. 39, 305 (2006).

3

Y.-I. Suzuki, J. Phys. B: At. Mol. Opt. Phys. 41, 215204 (2008).

4

Y.-I. Suzuki and T. Suzuki, Mol. Phys. 105, 1675 (2007).

5

T. Suzuki, Annu. Rev. Phys. Chem. 57, 555 (2006).

6

P. W. Langhoff, et al., Phys. Rev. A 14, 1042 (1976).

7

W. J. Hunt and W. A. G.oddard III, Chem. Phys. Lett. 3, 414 (1969).

図 (d) 核 酸 塩 基 分 子 の軌道エネルギーとイオン化の模式図。

4P106

ボーズ粒子を含む分子に対する時間依存多配置波動関数理論(東大院理)○ 加藤 毅，山内 薫

Time-dependent multiconﬁguration theory for molecules
including Bose particles(The University of Tokyo)○ Tsuyoshi Kato，Kaoru Yamanouchi[序] メタノールやアレンなどの炭化水素分子を，超短パルス・強レーザー場 (∼ 800 nm，∼ 50
fs，∼ 1014 W/cm2 ) と相互作用させ，単一の親分子から生成する複数のフラグメントイオン種の運動量を同時的に観測する (CMI) 実験から，分子（イオン）内超高速水素原子（あるいはプロトン，以下では水素と記す）運動の存在が報告された [1,2]．分子内水素移動過程は，入射パルスの時間幅より短い（< 60 fs）時間内で進行していることが実験的から明らかにされている [1]．すな
わち，分子内水素移動は，分子が強高度入射パルスと相互作用している間に，分子骨格を形成している結合の解離反応に先行する形で進行する．
最近のポンプ-プローブ時間分解 CMI 実験結果によれば，例えば，メタノールの１価カチオン状態における水素移動を二体解離反応によって追跡した場合，分子内水素移動過程はその時間ス
ケールに応じて２種類に分けられることが明らかとなっている [3]．一つは，分子とレーザー場との相互作用時間内 (∼ 38 fs) に進行するもので，他は，分子とレーザー場との相互作用が終了した
後の単分子解離反応 (∼ 150 fs) である．
本研究では，分子とレーザー場との相互作用時間内に進行する，初期水素移動過程に着目する．
ここで，我々は，レーザー場と相互作用している炭化水素分子における (i) 水素移動反応の時間ス
ケールを見積もること，(ii) プロトン運動は電子運動とどのような相関を持つのかを明らかにすること，(iii) プロトンの持つ量子性（波動性）がどの程度顕在化するのか，を明らかにすることを最終的な目的としている．
以上のことを理論的に解析するために，我々はこれまでに，電子に対する時間依存多配置波動関数理論（MCTDHF）[4] を拡張し，レーザー場中での超高速分子内水素移動過程を，電子-プロ
トン non-Born-Oppenherimer 分子波動関数のダイナミクスとして記述する多成分理論を提案した[5]．本研究では特に，CH3 OH, C2 H2 , C2 H4 , H5 O+
2 のような，(電子)+(プロトン)+(２つの重原子核) から構成される分子を考察の対象とする．これらの分子は，電子-プロトン構造に対して軌
道関数を導入するとき，あたかも二原子分子のように取り扱うことができるので “二原子様分子”と呼ぶ．二原子様分子は，強レーザー場中におけるプロトンダイナミクスを含む分子ダイナミク
スを考える上で最も単純なモデル分子系と見なすことができる．
我々は，拡張された MCTDHF 理論を使って，メタノール分子の電子-プロトン基底状態波動関数を数値的に求めることに成功しており，電子-プロトン基底状態波動関数から計算されるプロトン構造に対する条件付き確率密度関数の解析から，二原子様分子描像から導かれる分子構造が，BO近似に基づく分子構造と矛盾の無いものであることを示してきた [6]．
[目的] 本研究では，実験的に見出されている強光子場下におけるデューテロン置換メタノールの分子内水素移動過程における同位体効果 [1] を量子論的な観点から調べるために，ボーズ粒子であるデューテロンを軌道関数を使って取り扱うことのできる，ボーズ粒子を含む分子に対する時間依存多配置波動関数理論を提案する．ここでは簡単のため，プロトンとデューテロンの双方が同一分子に含まれている構成は考えない．

[理論]

次式の多配置展開された分子波動関数を考える．
XXΨ(re , rd , R, t) =CIB (R, t)ΦI (re , t)ΦB (rd , t)
I

(1)

B

ここで，re と rd はそれぞれ電子とデューテロンの座標を表し，R は２つの重原子核間距離である．ΦI と ΦB はそれぞれ，電子構造，デューテロン構造に対する determinant と permanent である．CI-係数に相当する CIB (R, t) は２つの重原子核間の量子論的な運動を表現する．式 (1) はいわゆる展開定理に基づいて構成されているが，ここでの特徴は，各スレーター行列式が R に露には依存しない点である．
電子，デューテロンの運動を記述するスピン軌道関数，および CIB (R, t) に対する運動方程式は，時間依存の変分原理に従って導出できる．導出方法は [7] で発表したものと同一である．文献 [8]等では，第二量子化法を使った定式化によって，フェルミ粒子系とボーズ粒子系を区別なく統一的に取り扱うことが出来るという点が強調されているが，本研究では逆に，両者の相違点を明確化する．デューテロンのスピン軌道の運動方程式を求める際，一体相互作用 ĥ，および Coulomb
相互作用 Ŵ に対する行列要素が必要となる．フェルミ粒子系と比較して質的に異なる行列要素を与えるのは主として Ŵ である．以下に，数表現された permanent |ΦB i = |n1 , n2 , n3 , · · ·i を使った場合の，０-，１-，２-粒子置換配置に対する行列要素を記す．


³
´
X
X
1
hΦB |Ŵ |ΦB i =nk (nk − 1) hk|Jˆk |ki +nl hk| Jˆl + K̂l |ki
(2)
2
|
{z
}
k
l(6=k)


³
´
X
p
ˆ
ˆ

hΦI |Ŵ |Φm
nj hm| Jˆj + K̂j |ki
kI i = nk (nm + 1) (nk − 1) hm|Jk |ki + nm hm|Jm |ki +|{z}j(6=m,k)(3)
p1
hΦI |Ŵ |Φlm
knI i = δ̄km δkn δlm (nm + 2)(nm + 1)nk (nk − 1) [mk| mk]2p+ δ̄km δkn δ̄l(k,m) (nl + 1)(nm + 1)nk (nk − 1) [lk| mk]p+ δ̄km δlm δ̄n(k,m) (nm + 1)(nm + 2)nk nn [mn| mk]nop+ δ̄km δ̄l(k,m,n) δ̄n(k,m) (nl + 1)(nm + 1)nk nn [ln| mk] + [lk| mn] .
|
{z
}

(4)

ここで，右辺に現れる添え字は各スピン軌道を区別し，1 = δij +δ̄ij なる記号を使っている．hi|Jˆl |ji =[ij|ll]，hi|K̂l |ji = [il|lj] （化学者の記法）である．それぞれの式で下括弧を付けた項は，フェルミ粒子系において対応する項を持つ．他の項は，ボーズ粒子系に特別な項である．例えば，平均場近似に対応する式 (2) から読み取れるように，交換積分の符号はフェルミ粒子の場合とは逆になる．さらに，一つのスピン軌道に複数の粒子が配置される (nk ≥ 2) とフェルミ粒子系では見られなかった第一項が出現する．また，２-粒子置換配置に対する行列要素を表す式 (4) の第一〜第三項は，パウリの禁制則が成立していなことにその起源をもつ．講演では，運動方程式の導出過程の詳細と運動方程式の性質とについて報告する．

[参考文献][1] T. Okino et al., Chem. Phys. Lett. 423, 220 (2006). [2] H. Xu, T. Okino, K. Yamanouch,

Chem. Phys. Lett. 469, 255 (2009). [3] H. Xu et al., J. Chem. Phys. 133, 071103 (2010). [4] T. Katoand H. Kono, J. Chem. Phys. 128, 184102 (2008). [5] T. Kato and K. Yamanouchi, J. Chem. Phys.
131, 164118 (2009). [6] T. Kato and K. Yamanouchi, Phys. Rev. A 85, 034504 (2012). 加藤 毅，山内薫，第 4 回分子科学討論会 2010, 2E08．[7] 加藤 毅，山内 薫，第 3 回分子科学討論会 2009, 1E16．[8] O.E.

Alon et al., Phys. Rev. A 76, 062501 (2007).

4P107
ガウシアン基底を用いた電子-核波動関数による H2 分子の量子動力学（東大院理） ○市川 雄一，加藤 毅，山内 薫

Quantum dynamics of H2 molecule by electron-nuclear wave function
using Gaussian basis(School of Science, the University of Tokyo)○Yuichi Ichikawa, Tsuyoshi Kato, Kaoru Yamanouchi【序】
水素分子の運動を原子核を含めて量子力学的に扱うために、電子と核の状態をそれぞれ結合コ
ヒーレント状態として記述し、系全体の波動関数はその掛け合わせとして扱った。ここで結合コ
ヒーレント状態とは、ガウシアン基底で表現される中心位置、運動量をパラメータに持つ状態であり、巨大な自由度を持つ量子系を計算するために用いられる手法である[1]。この方法を用いれば、波動関数のダイナミクスをガウシアン基底の中心位置と位相の時間変化として計算することができる。ガウシアン基底を用いることで、従来法であるグリッド法と比較してパラメータの数が少なくなり、その数は波動関数の空間的な広がりに依存せずに時間発展させることが可能となるため、計算コストを削減することが可能となる。本研究では、H2 を例に、電子と核の運動が少量のパラメータで効率よく求められることを示す。この手法は、強光子場中でのプロトンの超高速分子内移動や、イオン化など、波動関数が空間的に広がる場合の理論計算に適したものである。
【理論】
運動方程式は、変分原理から導かれる[2]。非相対論的量子力学においては、ラグランジアン𝐿と作用𝑆は
𝐿 = ⟨Ψ|𝑖𝜕̂ |Ψ⟩ ,
−𝐻
𝜕𝑡

(1)

𝑡2

𝑆 = ∫ 𝐿𝑑𝑡 ,

(2)

𝑡1

と表され、運動方程式は最小作用の原理から次のように導かれる。
𝛿𝑆 = 0 ⟹

∂𝐻
∂Ψ ̇
∂Ψ
= 𝑖 (⟨⟩ − ⟨Ψ̇|⟩) = ∑ 𝜂𝑖𝑗 𝑥𝑗̇
∂𝑥𝑖 |Ψ
∂𝑥𝑖
∂𝑥𝑖

(3)

𝑗

̂|𝛹⟩, 𝑥𝑗̇ = ∂𝑥𝑗 ⁄∂𝑡であり、{𝑥𝑗 }はガウシアンの中心位置や位相など、波動関数中のここで、𝐻 = ⟨𝛹|𝐻
パラメータの組を意味する。
分子の波動関数|Ψ⟩は、電子と核の波動関数の掛け合わせで、次のように表される。
|Ψ⟩ = |𝛼⟩|𝛽⟩ ,|𝛼⟩, |𝛽⟩: 電子と核のスレーター行列式(4)

これらのスレーター行列式は、ガウシアン基底|𝑧⟩
1

𝛾 4
𝛾
𝑖
𝑖𝑝(𝑡)𝑞(𝑡)⟨𝑥|𝑧(𝑡)⟩ = ( ) 𝑒𝑥𝑝 (− (𝑥 − 𝑞(𝑡))2 + 𝑝(𝑡)(𝑥 − 𝑞(𝑡)) +
)
𝜋
2
ℏ
2ℏ

(5)

を用いて構成される。ここで、𝛾はガウシアンの幅、𝑥は座標を表し、𝑝(𝑡), 𝑞(𝑡)は軌道の運動量と中心位置を表すパラメータである。

本手法では、運動方程式(3)に表れる∂𝐻⁄∂𝑥𝑖 及び𝜂𝑖𝑗 を解析的に求めることができるため、最急降
下法及び虚時間発展法によってエネルギーを最適化することが可能である。
【結果と考察】本研究では、電子系、プロトン系はそれぞれ 8, 2 基底関数を用い、電子、原子核は各々一重項状態として(3)式を評価する。以上で求められた運動方程式を H2 に適用し、そのダイナミクスをシミュレーションした。まず初めに、エネルギーの最適化を行った。具体的には適切な初期位置に置かれたガウシアン基底の位置を最急降下法で最適化する ことによって、𝐻: −1.0346 →−1.0649 a. u.となった。更に、この手続きで最適化された波動関数を虚時間発展法で最適化した。
図 1 はその時の各ガウシアン基底の中心位置の虚時間に対する時間変化を示している。実線は原子核の、点線は電子のガウシアン基底の中心位置を示している。
ガウシアン基底は x 軸
上に原点対称で配置されている。図 2 はその際の全エネルギーの変化を示している。
全エネルギーは−1.0665 a. u.に向かって収束しているが、これは NOMO 法による先行研究から得られた−1.069 a. u.[3]とほぼ一致している。ここで、図 2 に示した虚時間 50
a.u.までの最適化では完全には収束しきっておらず、わずかに平衡位置からずれたも

図 1: 虚時間発展におけるガウシアン軌道の中心位置の変化。実線 (–)は核の、点線 (--) は電子の基底の中心を示す。横軸の 1 a.u.は 24.19 アト秒である。

のとなっている。
この最適化された状態を初期位置として、実時間の運動を計算した結果を図 3 に示す。図 3 は、各ガウシアン基底の中心位置の変化を示している。核は 299 a.u. (7.2 fs)の周期で振動していて、これは、初期状態が先に述べた真の基底状態とわずかにずれていることに起因する。この周期は、水素分子の振動周期の実験値 313 a.u. (7.6 fs)に近い

̂ |𝛹⟩
図 2: 虚時間発展における全エネルギー𝐻 = ⟨𝛹|𝐻の変化。

ものとなっている。電子は核に追従して動いているが、実際には 50 as (2 a.u.)程度の周期で振動している。また、この計算においては、系全体のエネルギーは~3 × 10−6 %以内
で正しく保存した。
【参考文献】[1] A. Kirrander and D. V. Shalashilin, Phys.

Rev. A 84, 033406 (2011).
[2] P. Kramer, J. Phys.: Conf. Ser. 99, 012009

図 3: 実時間発展におけるガウシアン軌道の中心位置の変化。実線及び点線の定義は図 1 と同様である。核と電子は x 軸上を x=0 を中心として対称に運動している。

(2008).
[3] H. Nakai, Int. J. Quantum Chem. 86, 511 (2001).

4P108外場により誘起された非断熱過程からの発光の非線形ダイナミクス（東京大院・総合文化）○水野 雄太，新崎 康樹，高塚 和夫
Nonlinear dynamics of photoemission from the non-adiabatic process
induced by the external field(Tokyo University) ○Y.Mizuno,Y.Arasaki and K.Takatsuka【序】 LiF 分子はイオン結合性、共有結合性の２状態が相互作用する電子状態を持つ典型的な非断熱系である。この分子をレーザー電場中におくと、二つの透熱ポテンシャルの擬交差点は、その位置が時間的に振動することで、その上を運動している核波束に何度も遭遇することになる。核波束は擬交差点を通過するたびに分岐するので、複雑に分岐した波束同士が混ざることで、この系の挙動は大変複雑になると予想される。また、擬交差点周辺で電子状態がイオン結合性と共有結合性の間を瞬時に遷移するため、電気双極子モーメントが急激に変化し、これにより LiF 分子が発光するものと予想される。本研究では、量子核波束動力学計算の結果をもとに、外場中の LiF 分子の挙動、特に各透熱電子状態のポピュレーションや電気双極子モーメントの時
間変化とそれによる発光現象を量子系における非線形応答という観点から調べる。
【計算方法】 LiF 分子の核の波動関数は、以下のSchrö dinger方程式に従う:𝑖ℏ𝜕 Ψ1 (𝑅, 𝑡)
Ψ (𝑅, 𝑡)
̂( 1
(
)= 𝐻
)
𝜕𝑡 Ψ2 (𝑅, 𝑡)
Ψ2 (𝑅, 𝑡)

̂は、透熱表示のもとでは以下のように与えられる:
ここで、Hamiltonian 𝐻
̂ = 𝑇̂𝑁 1 + 𝑉(𝑅) + 𝑉𝐹 (𝑅, 𝑡)
𝐻
𝑇̂𝑁 は核の運動エネルギー演算子であり、𝑉(𝑅)は透熱ポテンシャル
𝑉 (𝑅) 𝑉12 (𝑅)𝑉(𝑅) = ( 11)𝑉12 (𝑅) 𝑉22 (𝑅)である。本研究では、LiF 分子の透熱ポテンシャル曲線は Giese ら[1]によって与えられたもの（図１）を使った。
𝑉𝐹 (𝑅, 𝑡)は、外場の影響を表す項であり、双極子近似を適用し、直線偏光外場に対して
𝜇 (𝑅) 𝜇12 (𝑅)
𝑉𝐹 (𝑅, 𝑡) = ( 11
) 𝐸(𝑡)
𝜇12 (𝑅) 𝜇22 (𝑅)と表されるものとした。本研究では、以上のSchrö dinger方程式を、Fourier 変換によって運動エネルギー演
算子を対角化する split-operator 法によって数値的に解いた。

図１：LiF の透熱ポテンシャル曲線。赤:V11(イオン結合性),
青:V22(共有結合性),紫:5×V12。

【計算結果】 振動基底状態にある核波束を周波数 ℏ𝜔𝑝 = 6.94 eV、FHWM 20fs の Gauss 型の包絡線を持つパルスレーザーでポンプしたあと、制御連続外場をかけた。ポンプ光の強度や、制御光の周波数・強度をいろいろに変えて計算したところ、図に示すような興味深い現象が見られた。

図２：比較的制御外場の周波数が低い場合の双極子モーメントの時間変化。下のパネルのように拡大すると、
約 12eV 相当の細かい振動が見られ、この振動の振動数は制御外場周波数に依存していないように見える。

図３：ポンプ強度を 1.0×1014W/cm2、制御外場周波数を 1.0eV に固定し、制御外場強度を変えた時の双極子モーメントの時間変化（左）とその Fourier 変換（右）。Fourier 変換は 219 as の幅のデータを用いて行い、使用データ左端を 0,100,200,300fs にしたものを上下にずらして図示した。双極子モーメントの時
間変化の制御外場強度に対する強い依存性が見られ、またその Fourier 変換には変換に用いるデータの位置に対する依存性が見られる。
【結語】 外場に駆動された LiF 分子は興味深い挙動を示すことが分かった。上記の現象の解析については、発表当日に言及したい。
【参考文献】[1]T.J.Giese and D.M.York,J.Chem.Phys.120,7939-7948(2004)[2]S.Scheit, Y.Arasaki, and K.Takatsuka, J. Phys. Chem. A116, 2644–2653(2012)[3]Y.Arasaki,S.Scheit, and K.Takatsuka,J.Chem.Phys.138,161103(2013)

4P109
集光デンドリマーをモデルするハミルトニアンの直鎖分割（東大院総合文化）○甲田 信一， 高塚 和夫

Straight-chain decomposition of model Hamiltonian for light-harvesting
dendrimers(The Univ. of Tokyo) ○ Shin-ichi Koda, Kazuo Takatsuka【序】光合成系において，光は多数の色素分子で構成される光捕集アンテナ系で吸収される．そして吸収された光のエネルギーは励起エネルギー移動によって色素分子間を渡り歩き，最終的に反応中心まで運ばれることになる．このエネルギー移動は，100 ％に近い量子収率で進むといわれている．近年では，X 線構造解析等の技術の進歩により，さまざまな光捕集アンテナ系の構造が解明されたが，その中で，色素分子が特異な配置を取っていることが明らかになった．そして現在では，色素分子の配置（大域的構造）とエネルギー移動の収率の高さの関連性を明らかにすべく，さまざまな研究が行われている．
一方で，光合成系を人工的に構築する研究も数多く行われているが，その中で，光捕集アンテ
ナ系を模した分子として，集光デンドリマーと呼ばれる分子種が提案されている．デンドリマーとは多くの枝分かれを持つ樹状の分子であり，集光デンドリマーでは，光を吸収する色素がデンドリマー中に組み込まれている．末端部で吸収された光のエネルギーは中心部に結合した分子に集められるが，このエネルギー移動の収率がデンドリマーの大域的な構造に大きく依存することが報告されている [1]．この例も，分子中の個々の色素間の相互作用だけではなく，分子全体の構造の重要性を示唆している．
励起エネルギー移動の収率を高める要因は主にふたつある．ひとつは移動速度の速さである．これは，蛍光放射など，他の過程により励起状態が失活する前に，目的の分子へエネルギーを送り届けるために必要となる．ふたつ目の要因は，移動の一方向性である．多数の色素間でランダムウォークの様にエネルギーが移動するとなると，寄り道をしてしまう分，目的の分子へエネルギーが到達するまで時間がかかることになる．したがって，光を吸収した色素から目的の分子まで一方向的にエネルギーを移動させることも，エネルギー移動の収率を高めるための要因となる．よって，これらの収率を高める要因とアンテナ系の分子の大域的構造の関連性を明らかにすることは重要な課題である．
本研究では，分子の大域的構造がエネルギー移動の一方向性を作り出し得ることを示す．具体的には，一例として集光デンドリマーを模したモデルハミルトニアンの性質を詳細に調べ，大域的構造によりエネルギーが末端から中心に集められることを確かめる．また，その原因も明らかにする．
【設定と計算】以下では，次のフレンケル励起子モデルハミルトニアンを用いる．

Ĥ =

∑
i

Ei |i⟩⟨i| +∑Jij |i⟩⟨j|

(1)

i̸=j

ここで |i⟩ は図 1 の i の色素が電子励起した状態とする．また，Ei は全て等しいものとし，簡単のためここでは Ei = 0 とする．Jij は図 1 で線分で結ばれている色素間で Jij = J ，それ以外で Jij = 0とする．本研究では，以下に定義される，熱平衡状態 ρ̂(T ) = e−Ĥ/kT /Tr(e−Ĥ/kT ) の各世代 (図 1 の∑中で同一の横並びの色素．例えば色素 2,3 は第２世代とする) の存在確率 Pn (T ) = i∈Gn ⟨i|ρ̂(T )|i⟩を計算する．Gn は第 n 世代に属する色素のラベルの集合とする．また，Pn (T ) と，エネルギー移

8

9

10

4

11

12

5

13

14

6

15

7

2

3

8’

9’

10’

11’

4’

5’

6’

7’

2’

3’

13’

14’

15’

1’

1

図 1: 色素間の相互作用0.6

図 2: 新たな状態間の相互作用
3

n=1
n=2
n=3
n=4

0.5

Pn(T)/PnR

0.4
Pn(T)

12’

0.3
0.2

n=1
n=2
n=3
n=4

2

1

0.1
0

0
0

0.5

1

1.5

2

2.5

kT/J

図 3: 世代ごとの存在確率

0

0.5

1

1.5

2

2.5

kT/J

図 4: 存在確率のランダムな場合との比動がランダムウォークで進むと見なした場合の存在確率 PnR (この場合，各色素における存在確率は全て等しい) の比も計算し，実際の熱平衡状態とランダムなケースの差異を示す．
【直鎖分割】本研究では，上記の設定のハミルトニアンが，適当な基底の変換によって単純化されることを見出した．この基底の変換は，同一世代に属する励起状態のみの重ね合わせで実現できる (変換の詳細は当日説明する)．新たに構成された状態間の相互作用は，直鎖状に連結したいく
つかのグラフで表される (図 2)．これは，一見複雑に見えるモデルハミルトニアンが，いくつかの単純な運動モードに分割できることを意味している．この新たな基底による表現を，ここでは直鎖分割と呼ぶことにする．
【結果と分析】ここでは４世代で構成されるモデルハミルトニアンを用いて計算した存在確率を示
す (図 3,4)．際立った特徴は，低温または相互作用が大きくなるにつれて，第４世代の存在確率が急激に減少することである (図 3)．この傾向は kT /J < 1 の領域で顕著である．一方，その他の世代の存在確率は増加することになるが，これとランダムにエネルギーが移動する場合との比を見ると，第１世代と第２世代において，顕著な増加が見られる (図 4)．これは，エネルギーが中心側の世代に集まりやすくなっていることを意味している．
この存在確率の偏りは，上記の直鎖分割を用いると簡明に解釈できる．まず，直鎖状の各運動モードに関し，鎖の長さが長いモードほど，状態の非局在化によりエネルギー準位の安定化が生じる．また，デンドリマー特有の樹状構造により外側ほど色素が多く存在するので，鎖が短いく
エネルギー準位が高いモードが外側に多く形成される．その結果，中心側ほどエネルギー的に安定なモードが存在することになり，熱平衡状態における存在確率は中心側に偏ることになる．
以上の結果は，たとえ同種の色素同士の複合体であっても，その大域的な構造によって，エネ
ルギー移動の一方向性が生じ得ることを示すものである．
【参考文献】 [1] W. Li and T. Aida, Chem. Rev. 109, 6047 (2009).

4P110ポリヒドロキシフラーレンの近赤外レーザーによる振動励起と反応に関する理論的研究(東北大院理 1 ，東北大理 2,

和歌山大院システム工 3，和歌山大システム工 4, 量子化学探索研究所 1,5)

池田旭伸 1，北條和彩 2，菅野学 1，時子山宏明 3，山門英雄 4, 大野公一 1,5, 河野裕彦 1

Theoretical study on the vibrational excitation and reaction processes ofpolyhydroxy fullerenes induced by near-infrared lasers(Tohoku Univ1, Wakayama Univ2, Institute for Quantum Chemical Exploration3)
Akinobu Ikeda1, Kazusa Hojo1, Manabu Kanno1, Hiroaki Tokoyama 2, Hideo Yamakado2,
Koichi Ohno3, Hirohiko Kono1【序】水酸化フラーレンC60(OH)24に定常レーザー(785 nm)を照射すると、低強度(<103W/cm2)でも高熱を持ち衝突や連続的な崩壊を経て、カーボンナノチューブや多層のフラーレンになることが知られている(Fig. 1)[1]。
カーボンナノチューブは優れた強度や弾性をもち、構造により伝導性が変化するため、高強度材料や半導体などへの応用が期待されており、この反応はナノチューブの新しい合成法として興味深い。フラーレンC60や他の官能基を持つ誘導体では同様の反応は見られず、詳しい反応機構は解明されていない。
この反応の初期段階を調べるため、定常レー
ザーによって加熱された水酸化フラーレンの挙動を理論計算により評価した。

Fig. 1 水酸化フラーレン C60(OH)24 から
ナノチューブへの反応過程
エーテル【手法】C60(OH)24の簡単なモデルとして
C60(OH)2を用いて密度汎関数緊密結合(Density Functional based Tight Binding;
DFTB)法により動力学計算を行った。温度
1700~3000 Kの範囲でC60(OH)2の動的挙動を調べた。そこで見られた反応を解析するために密度汎関数(Density Functional

エポキシド

Theory; DFT)法により遷移状態計算を行った。
【結果】一方の水酸基から隣の水酸基へ水素原子が移動し、水分子が移動する様子が複数のトラジェクトリで見られた。その結果、Fig. 2 脱水反応の生成物: エーテルとエポキシド

生成物として酸素が6員環と6員環を結ぶ単結合上に倒れこむエーテルと、酸素が5員環を結ぶ二重結合上に倒れこむエポキシドが得られた(Fig. 2)。それらの生成比は1700~3000 Kのトラジェ
クトリ全体でおおよそ5：3であった。
Fig. 3は脱水に要する時間をtとしたときのアレニウスプロットである。傾きから、エーテルが生成する反応の活性化エネルギーは2.0 eV、エポキシドの場合は1.4 eVであることが読みとれた。

Fig. 3 脱水反応のアレニウスプロット◆:生成物がエーテル ■:生成物がエポキシド

また、B3LYP/6-31G**レベルの遷移状態計算を行ったところ、前者の活性化エネルギー

中間体 A

は2.32 eV、後者は2.21 eVであり、競合して起こる反応であると予想される。
脱水後のC60Oの挙動をDFTB法で計算すると、エポキシドからエーテルへの転位やその逆の転位が見られた。B3LYP/6-31G**レベルの遷移状態計算の結果、活性化エネルギーはそれぞれ2.28 eVと2.19 eVであると見積もられた。温度2800 Kにおいてこれらの転位は1~15 psの間隔で見られ、脱水反応と同程度の時
間スケールで起こる反応である。

中間体 B

さらに2800 K以上のトラジェクトリにおいて、C60OからCOが脱離する反応が見られた。
B3LYP/6-31G**レベルの構造最適化により、CO脱離の中間体として2つの構造が得られた (Fig. 4)。これ
らの構造のエネルギーをFig.5に示した。中間体BからCOが脱離する反応の遷移状態であるTS 2のエネルギ
ーは5.44 eVと高いが、C60からC2 が脱離する反応の活性化エネルギーは10 eV程度であり[2]、水酸化Fig. 4 CO 脱離の 2 種類の中間体

フラーレンはC60より容易に崩壊する

TS2と考えられる。DFTB計算のトラジェ
クトリを解析すると、中間体Aから中間体Bを経由せずにCOが脱離する経路や、反応物から中間体Bに直接転位

TS 1
中間体 A

する経路が確認された。これらの反応

中間体 B

経路の詳細は当日発表する。
反応物(エポキシド)[1] V. Krishna, N. Stevens, B. Koopman, and B. Moudgil,

Fig. 5 CO 脱離のエネルギーダイアグラム

Nature Nanotechnology 5, 330 (2010).
[2] S.Matt, O. Echt, P Scheier, and T. D. Märk, Chem. Phys. Lett. 345, 194(2001).

4P111ガウス基底多配置核波束動力学法のプロトントンネリングへの適用(東北大院理) ○荒井 雄太, 山崎 達人, 菅野 学, 河野 裕彦

Application of the Gaussian-based multiconfiguration
wavepacket method to proton tunneling(Tohoku Univ.)○Yuta Arai, Tatsuto Yamazaki, Manabu Kanno, Hirohiko Kono
【序論】分子の量子ダイナミクスを理論的に追跡するには時間依存 Schrödinger 方程式(TDSE)を解く必要がある。TDSE の厳密解法として実空間グリッド法があるが、この手法は系の自由度の増加に伴い計算コストが指数関数的に増大する。多自由度系の量子ダイナミクスを扱うための近似解法として多配置時間依存 Hartree (MCTDH) [1]法がある。MCTDH 法は原子核の波動関数を時間依存 1 粒子関数の Hartree 積で表して多配置理論を適用することで計算コストを軽減できる。しかし、Dirac-Frenkel 変分原理から得られる波動関数の運動方程式を解くにはあらかじめ断熱ポテンシャル曲面を用意しておく必要があり、現実的な分子系の on-the-fly 動力学計算には適さない。
近年、時間依存 Gauss 関数の Hartree 積を利用する新しい多配置理論が提唱されている。Gaussian-based MCTDH (G-MCTDH)法 [2]は MCTDH 法の 1 粒子関数を Gauss
関数で置き換えたものであり、その時間発展は変分原理に従う。最近提案された Basis
Expansion Leaping Multiconfiguration Gaussian (BEL MCG)法 [3]では、時間発展に伴って新しい Gauss 基底を生成して波動関数を再展開する。これらの手法はいずれも波動関数の時間発展に必要な種々の積分を Gauss 積分にして取扱いを容易にするため、多自由度系の on-the-fly 計算が可能になると期待されている。
本研究の目的は、これらの Gauss 関数を基底とする新しい多配置波動関数理論と量子化学計算を組み合わせて分子系に適用することである。そのためには、G-MCTDH
法や BEL MCG 法がトンネル効果などの量子効果を適切に評価できるか検証する必要がある。まずは TDSE を厳密に解くことができるモデル系を用いて、これらの手法のトンネル効果に対する有効性を調べた。本発表では主に G-MCTDH 法について議論する。
【理論】G-MCTDH 法における波動関数 Ψ は時間依存 Gauss 関数{gjk(Ri, t)}の Hartree
積 (配置関数)の線形結合

R1 , R 2 ,  R N , t    Aj1 , jN t g j1 R1 , t g j2 R 2 , t  g jN R N , t (1)j1 , jN

g jk R k , t   exp R k  α jk t   R k  β jk t   R k   jk t (2)で表される。{Rk}は原子核の自由度、{Aj …j }は展開係数、{αjk, βjk, γjk}は Gauss 関数の 2次、1 次、0 次のパラメータを表している。この波動関数に Dirac-Frenkel 変分原理を適用すると展開係数 A の運動方程式
1N  S1 H  iτ A
iA(3)と Gauss パラメータ Λ={αjk, βjk, γjk}の運動方程式
  C1YiΛ

(4)が得られる。S は配置間の重なり積分、H は Hamiltonian 行列である。行列 τ, C やベ
クトル Y は配置間の位置{Rk}に関する行列や Hamiltonian 行列などを含んでいる。これらの要素は Gauss 積分を用いて高速に計算できる。
Gauss 関数には幅{αjk}が時間に依存する Thawed Gaussian (TG)と{αjk}を定数とする
Frozen Gaussian (FG)の 2 種類がある。TG を用いると柔軟に波動関数を表現できるのに対し、FG は TG よりも数値安定性が良い。本研究では FG を採用する。
【結果と考察】本研究では Fig.1 に示した 1 次元 2 重井戸型ポテンシャルに G-MCTDH法を適用した。配置関数の数は 2 配置

x, t   A1 t g1 ( x, t )  A2 t g2 ( x, t )(5)

ポテンシャルエネルギー [cm-1]とした。初期時刻で Gauss 関数の中心を-0.479 Å, 0.479 Å の位置に配置し、展開係数を A1(0) = 1, A2(0) = 0 として左側の井戸に波束を局在化させた。
核波束の時間発展を Fig.1 に示した。波束は左側の井戸から右側の井戸に移動して再び左側の井戸に戻った。この周期は 15.22 ps であった。同じ初期条件から実空間グリッド法で計算すると、波束の振動周期は 15.10 ps となった。また、2 重井戸型ポテ
ンシャルの基底状態と第一励起状態のエネルギー差 ΔE =3.8 ps0 ps-12.340 cm から算出したトンネ2420

リング周期は 15.10 ps である。
G-MCTDH 法のトンネリング0
周期はこれらの値と良く一致
7.6 ps15.1 ps
している。
G-MCTDH 法において波束
2420のトンネリング周期はHamiltonian 行列の非対角項
0
-0.5
0
0.5
1
-1
0
-0.5
0.5
1-1
H12 = H21 によって決定される。
核座標 [Å]Fig.1 波束の時間発展
今回の計算では 2H12 = -1.774(赤) G-MCTDH, (緑) 実空間グリッド法
cm-1 となり、その絶対値はエネ(青) 基底状態と第一励起状態の固有関数の重ね合わせ
ルギー差 ΔE = 2.340 cm-1 と近い値になったので精度良くトンネル効果を記述できた。
[1] H. –D. Meyer, U. Manthe, and L. S. Cederbaum, Chem. Phys. Lett. 165, 73 (1990).
[2] I. Burghardt, H. –D. Meyer, and L. S. Cederbaum, J. Chem. Phys. 111, 2927 (1999).
[3] W. Koch and T. J. Frankcombe, Phys. Rev. Lett. 110, 263202 (2013).

4P112強レーザー場誘起の多電子ダイナミクスと複数活性分子軌道のエネルギー解析(東北大院理 1, 横市大院生命ナノ 2, 東大院理 3, 大阪府大院理 4)○大村周 1, 小山田隆行 2, 加藤毅 3, 河野裕彦 1, 小関史朗 4

Multielectron dynamics induced by intense laser fields:
Energy analysis of activated molecular orbitals(Tohoku Univ.1, Yokohama City Univ.2, The Univ. of Tokyo3, Osaka Prefecture Univ.4)○Shu Ohmura1, Takayuki Oyamada2, Tsuyoshi Kato3, Hirohiko Kono1, Shiro Koseki4【序】 分子が近赤外強レーザー場 (≈ 800 nm, I > 1012 W/cm2) にさらされると, トンネルイオン化から再結合過程を経て, 基本波の整数倍の振動数を持つ高調波が生成される高次高調波発生 (high-order harmonic generation, HHG) が起こる. HHG スペクトルはイオン化していっ
た電子が属していた分子軌道の情報を反映しており, そこから分子構造[1]や電子ダイナミク
ス[2]を読み取ることができる. 従来 HHG に関する理論研究は, 単一軌道からのイオン化のみを考慮する single-active-electron モデルが中心であった. しかし近年になり, 強レーザー場において複数の分子軌道からのイオン化が観測され[3], 多電子の相関を考慮した理論解明が求められている. 我々は電子相関を考慮した多電子ダイナミクスを求めるために, 多配置時間依存 Hartree-Fock (MCTDHF) 法[4]を開発してきた. 本研究では内殻軌道が多電子ダイナミク
スに及ぼす影響を調べるために, 内殻を持つ最小系である LiH 分子に MCTDHF 法を適用する.
シミュレーションから得られた双極子モーメントを複数の分子軌道が関与する経路に分解し
HHG のスペクトルを計算する. さらに, 分子軌道のエネルギーを定義することで軌道間のエ
ネルギー交換を定量化し, HHG の複数経路が生じる要因を明らかにする.
【理論】

MCTDHF 法では電子波動関数 (t ) を(t )   C I (t ) I (t )(1)Iと展開する. ここで C I (t ) と  I (t ) はそれぞれ電子配置 I の時間依存係数とスレーター行列式である. (1) 式を Dirac-Frenkel の時間依存変分原理に代入すると, 係数と分子軌道両方の運動方程式を得ることができる. これにより効率よく系の時間発展を記述することができる.
MCTDHF 法で得られた分子軌道は各時間ステップで 1 次の密度行列を対角化する自然軌道| j(t), |k(t), …に変換することができる. したがって 1 体の演算子である電気双極子モーメン
ト演算子 d̂の期待値 d (t ) は自然軌道ごとに分解することができ, これをさらに初期時刻の自然軌道| j(0), |k(0), …に射影する.

d (t )   Re[n j (t ) k (0) | dˆ | j (t ) j (t ) | k (0) ]
j

(2)

k

nj(t)は| j(t)の占有数である. (2)式右辺の各項は, |k(0)から| j(t)に遷移し HHG 過程により|k(0)

に戻る k(0)→j(t)→k(0) という経路を表している. 各項をフーリエ変換すれば経路ごとの HHG
スペクトルを求めることができる.
(2)式のような複数経路が生成される原理を調べるために, 自然軌道のエネルギーとして

 j (t ) E (t )
n j (t )(3)で定義される時間依存軌道化学ポテンシャル[5]を導入する. ここで E(t)は全電子エネルギーであり, j(t)を用いて E (t )   j n j (t )  j (t ) と表すことができる. 軌道ごとの化学ポテンシャルの時間変化j(t)を比べることで, 軌道間のエネルギー交換を定量化することができる.
【結果と考察】 MCTDHF 法を用いて LiHの近赤外強レーザー場 (= 1,522 nm, I =3.0  1012 W/cm2) におけるイオン化シミュ
レーションを行った. 分子軸偏光の 2 サイ
クル電場を仮定し, 3 軌道 (1, 2, 3) 9 配
置 (基底電子配置 : 12) を考慮して計算した.
図 1 に全 HHG スペクトルと各経路に分解したスペクトルを示す. これを見ると, 今の場合基底配置の HOMO である 2軌道のみを遷移する 2→2 →2という経路が主

図 1: 複数経路及びトータルの HHG スペク

要 (> 99 %) であることがわかる. 我々の定トル. 2→2 →2 が主要な経路である

式化では, 異なる軌道を経由する経路のスが,異なる軌道を経由する経路も存在する.

ペクトルも得ることができる (図には2→1 →2 の み プ ロ ッ ト ). 図 2 は{j(t)} と E(t) の時間変化E(t)を表している. {j(t)} はE(t)よりも大きく変化しており, 価電子軌道と内殻軌道間でエネルギー交換していることがわかる. 異なる軌
道を遷移する HHG の経路はこの相互作用に起因していると考えられる.
本研究の MCTDHF 法を用いたシミュレーションと化学ポテンシャルを用いたエネ

図 2: { j(t)}と E(t)の時間変化E(t). E(t)とルギー解析により, 電子相関を含んだ多電比べて j(t)は大きく変化しており, 自然軌

子ダイナミクスの正しい理解へつながるこ

道間でエネルギー交換している.

とが期待される.
【参考文献】[1] J. Itatani et al., Nature 432 (2004) 867 [2] H. Niikura et al., Phys. Rev. Lett. 107(2011) 093004 [3] A. E. Boguslavskiy et al., Science 335 (2012) 1336 [4] T. Kato and H. Kono, J.
Chem. Phys. 128 (2008) 184102 [5] T. Kato et. al., Prog. Theor. Phys. Suppl. 196 (2012) 16

4P113
光格子中に捕捉された KCs 分子量子ビットを操作するマイクロ波パルスの数値設計(東北大院理) ○新井 健太、大槻 幸義、河野 裕彦

Numerical Synthesis of Microwave Pulses for Manipulating KCs Molecular Qubits
Trapped in an Optical Lattice(Tohoku Univ.) ○K. Arai, Y. Ohtsuki, and H. Kono
【序】量子コンピュータが実現すれば、量子重ね合わせを利用する超並列計算により現在のコンピュータ技術の延長では取り扱うことのできない問題を解くことができる。現在までに光子や核スピン(NMR)などを用いた数ビットの量子演算については実証実験が報告されている。しかし、実用的な問題を解くことができるサイズの量子コンピュータを実現するには多くの課題が残っている。このような中、光格子中に捕捉された冷却 KCs 分子を用いれば 104 個程度の量子ビットへ拡張できるのではないかという予測が近年報告された[1]。
このモデルでは 1 次元に配列された各量子ビット(qubit)を区別するために不均一な静電場を印加する。そのため分子の固有状態は単純な回転状態ではなく球面調和関数の重ね合わせであるペンデュラー状態で表される[2]。そこにマイクロ波パルスを照射し各量子ビットを選択的に操作する。すべての量子アルゴリズムは基本的な 1 および 2 量子ビット演算で表されることが知られている。そこで本研究では 1 量子ビット演算に着目し，最適制御法[3]を用いて高精度かつ短時間で演算するパルスを数値設計しそのスペックを提示する。
【理論】剛体回転子でモデル化した冷却 KCs 分子を考える。捕捉レーザーパルスとの分極相互作
用部分を除いて表せば，光格子中の分子のハミルトニアンはH  BJ 2  μ  E dc  μ  ε(t )(1)で与えられる。B は回転定数、J は角運動量演算子、 μ  Edc は静電場との相互作用項、 μ  ε (t ) は演算マイクロ波パルスとの相互作用項である。静電場との相互作用によって回転状態が混ざり合
いペンデュラー状態となる。このモデルでは偏光条件により M=0 の球面調和関数の重ね合わせで表された各ペンデュラー状態| 0    a j Y j ,0 ( ,  ) ,j| 1    b j Y j ,0 ( ,  ) .

(2)jを量子ビットの論理基底として用いる。本研究では最適制御法に基づく以下の手順で演算パルスを数値設計する。目的の演算子を W、演算パルス ε(t)の下での時間発展演算子を U (t f ,0,  (t )) とし、(3)式を用いて演算パルスを評価する。( | n  は量子ビットの基底) (t )
1 N 1F[ (t )]   | W | n  U (t f ,0,  (t )) | n  | 2   dt
N n 0A(t )
2(3)

右辺第 1 項はパルスの精度を評価する項、第 2 項はパルスのエネルギーを最小化するために導入したペナルティ項である。最適制御法ではシュレーディンガー方程式を拘束条件とし、変分法を用いることで F の最小化条件から演算パルスの設計方程式を導く。演算パルスの設計方程式は

N 1

 (t )   A(t ) Im  n | W †U (t f , t ,  (t ))U (t ,0,  (t )) | n 

(4)

n 0

と導出できる。(4)式をシュレーディンガー方程式と連立して解くことで最適な演算パルスを求める。なお解法には単調収束保証の繰り返しアルゴリズムを用いる [3]。
【結果】今回は量子情報処理において重要な 1 量子ビット演算の一つであるアダマール変換

WH | 0   (| 0   | 1 ) / 2

WH | 1   (| 0   | 1 ) / 2(5)を実行する演算パルスを設計した。はじめに 1 量子ビット系に対してアダマール変換を実行するパルスを設計した。(3)
式の右辺第一項を用いてこのパルスの精度を評価すると99.999%であった。
次に図 1 に示すような 3 量子ビットシステムを仮定し、中央の量子ビットのみを選択的にアダマール変換することを考える。まず先に設計したパルスをこの 3 量子ビット系に照射すると精度は 99.937%であった。なお、両端の量子ビット
は分布の変化のみを精度評価に考慮している。1 量子ビット

図 1：3 量子ビット系の概略図

系の場合より精度が低下してしまったのは、静電場による遷移周波数の分離が不十分であり両端の量子ビットの遷移が誘起されてしまったためであると考えられる。量子ビット間の静電場強度の差を大きくする、あるいは演算系への拡張および演算時間の短縮という観点からは望ましくない。
そこでこの 3 量子ビット系でパルスの再設計を行った。

Electric Field (V/m)

時間を長くすれば演算精度は改善するが、多量子ビット

(a)

100
50
0
-50

その結果、同じ静電場強度、演算時間でも演算精度を-100

99.999%と大きく改善することができた。図 2(a)に設計

1.0

(b)
|0>
|1>

ットの分布の時間変化を示す。最適パルスの下では中央のビットの演算精度を保ったまま、両端の量子ビットのPopulation

した演算パルス、図 2(b)にそのパルスの下での各量子ビ
0.5qubit2

qubit1,qubit3

分布がほぼ完全に元に戻っていることが確認できた。このように最適制御法で整形された演算パルスを用いれば短時間で高精度な演算が可能になる。

0.0
0

250

500

Time (ns)

図 2：選択的なアダマール変換【参考文献】(a)：最適パルス[1] D. DeMille, Phys. Rev. Lett. 88, 067901 (2002).

(b)：各量子ビット基底の分布の時間変化[2]Qi Wei et al., J. Chem. Phys. 134, 124107 (2011).
[3] Y. Ohtsuki, New J. Phys. 12, 045002 (2010).

4P114非断熱ダイナミクスにおける電子配置のエネルギー解析（東北大院理）◯佐藤悠介、大村周、菅野学、河野裕彦
Energy analysis of electron configurations in nonadiabatic dynamics(Graduate School of Science, Tohoku Univ.)◯Yusuke Sato, Shu Ohmura, Manabu Kanno, Hirohiko Kono
【序】イオン結合性分子の解離反応やドナー・アクセプター間での電子移動反応など多くの化学反応において、
ポテンシャルエネルギー曲面 の擬交差付近で起きる非断熱遷移が重要な役割を果たしている。
しかしながら、
非断熱遷移の全貌を明らかにするには、原子核のダイナミクスだけ ではなく，非断熱遷移の前後で急激に変化する電子状態のダイナミクスの詳細な解析が必要である。
我々は、これまで電磁場と相互作用をしている分子の動的な電子状態を評価する手法として、
電子エネルギー期待値を各分子軌道に分配し、軌道間のエネルギー交換を評価する手法を開発した[1]。水素分子の光イオン化に適用し、各分子軌道からのイオン化確率と、その軌道に分配されたエネルギーとの間に相関が見られることを見出した。電場と相互作用する分子の電子ダイナミ
クスを考える上でこの手法が有用であることを示している。
本研究ではこの手法を核ダイナミクスを含む化学反応に応用し、各電子配置に分配される擬似的なエネルギー（配置エネルギー）を導入する電子配置の分布変化に加え、電子配置同士のエネ
ルギー交換の観点から擬交差における非断熱ダイナミクスの解析を試みた。擬交差を経由する具体的な反応として LiF 分子の光解離反応を考え、電子状態及び配置エネルギーの時間変化を定量化した。
【理論】原子核の運動は古典的に扱い、核座標を R(t)で表す。核座標 R(t)における電子波動関数
(t;R(t))をスレーター行列式{i (R(t))}で多配置展開する。

(t; R(t ))   Ci (t ) i ( R(t ))(1)i

Ci(t)は展開係数である。各電子配置の配置エネルギー{Ei(t)}は、電子ハミルトニアン Ĥ e と電子波動関数のエネルギー期待値 E (t )   Hˆ e  を用いて次の式で定義される。
*
 C j (t )C j (t )E (t )ˆˆEi (t ) 


H



H

j  C (t ) i e j C * (t ) j e i 
 | Ci (t ) | 2
i
 i


(2)

この配置エネルギーを用いると全エネルギー期待値は次式で書かれる。

E (t )   Ci (t ) Ei (t )
2

(3)

i

定常状態では全ての配置エネルギー{Ei(t)}が縮退するが、非断熱的な電子状態変化により縮退が解け、各電子配置間のエネルギー交換を評価することができる。また、(1)式を電子ハミルトニア
ンに関する時間依存 Schrödinger 方程式に代入して積分すると次式を得る。
tt
C j (t ' )1 t1Ci (t )  exp   Ei (t ' )dt '  i (t ' )dt '  i  j dt '  Ci (0)0Ci (t ' )0 j
 i 0


(4)

ただし

i (t )  Im j

C j (t )Ci (t )

 i Hˆ e  j(5)(4)式の第一項は Ci(t)の位相変化に寄与する。定常状態における波動関数の位相変化との類推か
ら、Ei(t)が電子配置の瞬間的なエネルギーを表すと言える。 (5)式で定義されるi(t)は分布|Ci(t)|2の瞬間的な変化速度に寄与する。
今回の計算では核座標 R(t)と展開係数 Ci(t)の時間発展は半古典 Ehrenfest 法[2]を用いて計算した。
半古典 Ehrenfest 法では、核座標 R(t)は状態平均された Hellman-Feynman 力を受けて Newton の運動方程式に従い運動する。展開係数 Ci(t)は時間依存 Schrödinger 方程式から導かれる運動方程式に従って時間発展し、非断熱カップリング項により非断熱遷移が起こり複数の電子状態が混ざり合う。

底関数系は SA-CASSCF(2,5)/6-31G*を用いた。
活性軌道としてxy 軌道を選んだ。
得られた LiF 分子の断熱ポテンシャル曲線及び非断熱カップリング項の核間距離 R 依存性を図 1 に示した。
平衡核間距離は約 1.6 Å で、擬交差が存在する 3.1 Å 付近で非断熱カップ

断熱ポテンシャル / eV

一励起状態を考えた。電子状態の計算法と基

B1+

X1+

リング項が大きくなっていることがわかる。
始時刻 t = 0 で分子は電子基底状態の平衡核

非断熱カップリング項 / Å-1【結果】時間発展は LiF 分子の基底状態と第間距離において第一励起状態に垂直励起した核間距離 R / Å
図 1 SA-CASSCF 計算により得られた、LiF 分と仮定した。また十分な非断熱遷移を起こす子の電子基底状態・第一励起状態の断熱ポテン

ために余剰の初期運動エネルギー200 eV をシャル曲線（実線）及び非断熱カップリング項もつと仮定して時間発展を計算した。

（破線）

大きな分布を持つ、主要な４つの電子配置
ギー期待値の差を図 2 に示した。t = 0 では全ての配置エネルギーが縮退するが、時間と共に縮退が解け、それらは互いにエネルギーを交換しながら振動した。さらに、分子が擬交差を通過する t ~ 2.2 fs において全ての配置エEi(t)-Etotal(t) / eVについて配置エネルギーと系の全電子エネル

ネルギーが再び縮退することがわかった。この縮退は初期条件によらず擬交差を通過する瞬間に必ず起こる。このことから、配置エネ
ルギーは、擬交差を経由する化学反応における電子状態変化と相関を持つことがわかる。
解析の詳細は当日述べる。

time / fs
図 2 主要な 4 つの電子配置についての配置エ
ネルギーEi(t)の時間変化[1] T. Kato and H. Kono, Chem. Phys. 366, 46 (2009).
[2] H. –D. Meyer and W. H. Miller, J. Chem. Phys. 70, 3214 (1979).

4P115
誘電率に対する水の酸化還元電位依存性(金沢大・自然1 , 金沢大・大教センター2 )○ 杉森哲也1 , 杉森公一2 , 齋藤大明1 , 川口一朋1 , 長尾秀実1
Dependency of the water redox potential on the dielectric constant(Graduate School of Science, Kanazawa Univ.1 , Research Center for Higher Education Kanazawa Univ.2 )○ T.Sugimori1 , K.Sugimori2 , H.Saito1 , K.Kawaguchi1 , H.Nagao1【序論】水が関与する生体内化学反応の一つである光化学系 II では、初期過程において水分子が電子供与体 (還元剤) として (1) 式のように水素イオンと酸素に分解される。

2H2 O → 4H+ + 4e− + O2(1)

この反応は、様々な酸化還元状態を持つマンガン酸素クラスターによって触媒されており、近年構造決定された [1]。生体膜中には、光化学系 II の触媒中心以外の領域でも水分子が含まれており、膜分子の疎水性あるいは親水性の領域で電子供与体としての役割を果たす可能性がある。本研究では、生体膜内の誘電率が親水性と疎水性の領域で変化すると仮定し、水分子の酸化還元電位がどのように変化するかを予測する。
そこで、疎水性タンパク質内部をあらわす非極性環境 (誘電率  = 2) 下での水の酸化還元電位を密度汎関数
理論 (DFT) および分極連続体モデル (PCM) を用いて計算し、議論した。
【理論と計算方法】酸化還元電位を計算するため、様々な方法論が提案されている。熱力学的サイクル法を使ったインプリシットな連続体溶媒モデルは、原子論的な振る舞いが考慮されないが、酸化還元メカニズムの物理的意味を得るためには有用かつ効率的である。図 3 および (2) 式のように定義した酸化還元電位の式を使用する [2]。

E(O/R) = IPv

λ + ∆Ggas (O)∆Ggas (R) + ∆Gsol (O)∆Gsol (R)1
E(H+ / H2 )2(2)

ただし、IPv は零点エネルギー補正の差を含んだ垂直イオン化ポテンシャルを意味する。λ は IPv と断熱
イオン化ポテンシャル IPa との間の差つまり再配向エネルギーである。∆Ggas (X) は、還元型 (X = R) ま
たは酸化型 (X = O) の分配関数の使用によるエンタルピー補正、∆Gsol (X) は、気相中での最適化構造を用いた様々な溶媒和における自由エネルギー差、E(H+ / 21 H2 ) は標準水素電極 (SHE) 電位 4.43V を意味する。今回の計算では Gaussian03 を使用した [3]。水分子中の原子の初期配置は、O

HO

H 結合長を約 0.94 Å、

H 結合角を約 104.5 度に設定した。密度汎関数法の一つである B3LYP 汎関数と LANL2DZ 基底関数を用いて中性状態の水分子 H2 O とそのラジカル陽イオン H2 O・+ の構造最適化を行う。得られた最適化構造に対して同様のレベルで振動数計算をおこない最終的な電子エネルギーを求めた。これらの計算を
気相 (誘電率  = 1) 中と分極連続体モデルで近似した様々な溶媒 ( = 1 以外の誘電率) 中で繰り返し、誘電図 1: H2 O モデル ( = 1.0)

図 2: H2 O・+ モデル ( = 1.0)

B3LYP/LANL2DZ による構造最適化後の構造パラメータ

率を 0.1 から 100.0 まで 0.1 刻みに設定した。
【結果と考察】B3LYP/LANL2DZ レベルで構造最適化した H2 O と H2 O・+ のモデルを図 1 と図 2 にそれぞれ示す。真空中と誘電率  を持つ様々な媒質中の H2 O と H2 O・+ の電子エネルギーの誘電率依存性およびこれらの計算結果に基づく酸化電位の誘電率依存性をプロットした曲線を図 4 および図 5 にそれぞれ示す。
図 4 を見ると、水分子とそのラジカル陽イオンの電子エネルギーの値は比較的高い誘電率で緩やかに変化しているが、比較的低い誘電率の部分では急激に変化しており、誘電率が 1 で最大値をとるように見える。
この誘電率依存性を調べるため、それらの曲線を (3) 式でフィッティングした。

y() =a+c( + d)b(3)

ただし、y() は誘電率  に依存するパラメータ、a は  = 1 でのパラメータのスケールをあらわす量、b は減衰度、c は減衰後の収束値を意味しており、d はほぼ 0 と仮定した。図 5 を見ると、 = 1 のときの水の酸化電位は 12.27eV と見積もられ、これは真空下での水のイオン化ポテンシャル 15.43eV と概ね一致していると解釈できる [4]。 = 2 のときの水の酸化電位は 10.32eV(5.89V/SHE) と見積もられる。還元電位に関する詳細およびその他の分子種との比較は当日の学会で発表する。
表 1: 様々な誘電率 † で置き換えた溶媒中での水の電子エネルギー Eelec (eV)、ゼロ-2030

Oxidized Water
a=4.46097 b=0.947856 c=-2071.31
Neutral Water
a=0.494192 b=0.771755 c=-2079.69

50

-2079

45

-2079.1-2079.2-2045-2079.3-2050-2079.4-2055-2079.5-2060-2079.6-2065-2079.7

-2070
-2075

0

20

40
60
Dielectric constant:

80

50

35
30
25

Data Points
y=ax+b ; a=-3.80483,b=12.0889

45
40
35
30
25
20
15
10
5
-10

20

-8

-6
-4
-2
1-1/[Dielectric constant]:

015-2079.810-2079.91005

図 4: 電子エネルギー Eelec の誘電率依存性

Data Points
a=3.99168 b=0.973704 c=8.26259

40
Oxidation Redox Potential:(V)-2040

Electronic energy of neutral water:(eV)Electronic energy of oxidized water:(eV)-2035-2078.9

Oxidation Redox Potential:(V)

図 3: 熱力学的サイクル

点エネルギー補正 ZPE(eV)、断熱イオン化ポテンシャル IPa (eV)、酸化電位 Eox (V)の計算値
Eelec (eV)
ZPE(eV)
H2 O
H2 O・+H2 O
H2 O・+IPa (eV)Eox (V/SHE)
1.0-2079.23 -2066.89 0.56499 0.49223
12.27037.840321.4-2079.32 -2068.08 0.56493 0.48586
11.1608
6.73075
2.2-2079.42 -2069.19 0.56417 0.47598
10.1355
5.7054937.5-2079.66 -2070.61 0.55898 0.60874
8.37043.9403778.4-2079.67 -2071.24 0.55876 0.43438
8.30953.87953†  の値はそれぞれ真空 (1.0)、アルゴン (1.4)、ベンゼン (2.2)、アセトニトリル(37.5)、水 (78.4) の誘電率。

0

10

20

30

40
50
60
Dielectric constant:7080

図 5: 酸化電位 Eox の誘電率依存性【参考文献】[1]Y. U. Umena et al. Nature. 473, 55-61(2011).
[2]P. Winget et al. Phys. Chem. Chem. Phys. 2, 1231-1239(2000).
[3]M. J. Frisch et al. Gaussian03 , revision B.01, Gaussian, Inc(2004).
[4]C. G. Elles, et al. J . Chem. Phys. 125(4), 44515-44515(2006).
[5]R. Jono, M. Sumita, Y. Tateyama, K. Yamashita, J . Phys. Chem. Lett. 3, 3581-3584(2012).

90

100

4P116
粗視化シミュレーションを用いた脂質二重層膜の構造安定性に関する理論的研究(金沢大院∗ , 東京薬科大∗∗ )⃝ 田村 一真∗ , 吉田 津日門∗ , 高須 昌子∗∗ , 齋藤 大明∗ , 川口 一朋∗ , 長尾 秀実∗

Theoretical study on the structural stability of the lipid bilayer membrane
by a coarse-grained simulation(Kanazawa Univ.∗ , Tokyo Univ. of Pharm. and Life Sci.∗∗ ) ⃝Kazuma Tamura∗ ,Tsuhito Yoshida∗ , Masako Takasu∗∗ , Hiroaki Saito∗ , Kazutomo Kawaguchi∗ , Hidemi Nagao∗1

研究背景
脂質二重層膜は細胞膜を構成していることで知られており、代表的な構成分子であるジオレオイルホス

ファチジルコリン (DOPC) やジパルミトイルホスファチジルコリン (DPPC) は古くから盛んに実験・研究されている。
この膜は親水部と疎水部から成る脂質分子が集まり、上下逆向きの二層になって形成している。脂質膜は常温常圧で巨視的には 1 枚の膜に見えるものの、微視的には脂質分子１つ１つが膜面方向に拡散によって流動している。また、脂質平面膜が閉じて球状になったベシクルはエネルギー的に安定であり、そのサイズ
は閉じる前の平面膜の大きさによることが知られている。脂質膜の全原子シミュレーションはこれまでも広く行われているが、全原子シミュレーションの時間スケールは膜で起こる多くの現象の時間スケールよりも短く、それらの現象は観察されていない。
長時間のシミュレーションを可能にし、時間スケールの大きな現象を観察するために開発された手法が粗視化モデルである。粗視化モデルを用いること(a) 粗視化モデル

により時間スケールのみならず距離スケールも大きな系を扱うことができる
ため、AFM をはじめとした直接観察された実験結果と比較することができる[1]。脂質膜の粗視化のシミュレーションは、脂質 1 分子を 10〜数 10 個の粒子で粗視化し水をあらわに計算する方法がよく使われる。しかしながらこの粗視化の程度だと時間スケールや距離スケールが中小規模しか扱えない。そこで、
脂質 1 分子を 1 個または 3 個の粒子で粗視化し水を陰溶媒で扱うモデルで計算が行われた [2]。それに対し、膜に対して水は大きな影響を与えるため、水(b) 全原子モデル

もあらわに計算する方法を Kawamoto らが行った [3]。Kawamoto らは脂質とペプチドをそれぞれ 3 粒子と 4 粒子で置き換えたモデルを用い、ミセルの形
成やペプチドの膜透過の計算を行った。
我々は孤立巨大平面膜及びベシクルの構造安定性の研究のために、Kawamotoらが行った脂質と水の粗視化モデル (図 1(a)) を用いたシミュレーションを実行できるプログラムを作成した。図 1(b) には比較のためそれぞれに対応した
全原子モデルを掲載してある。本研究では先行研究との比較のため、周期境
界のかかった平面膜を用いてシミュレーションを行い、分子間相互作用のパラ
メータと膜の構造安定性の関係について議論する。

2

図 1: DOPC と水の粗視化モデルと全原子モデル

計算手法
本研究では、Kawamoto らが行った計算と同じモデルを用いた [3]。Kawamoto らの用いたモデルは 1 つの親水性粒子と 2 つの疎水性粒子で脂質を表現し、複数の水分子を 1 つの粒子で表現している。ポテンシャ
ルは次の式 (1) で与えた。
∑
∑
∑ { σ
σ }
Utotal =Kbond (r − σ0 )2 +Kangle (θ − θ0 )2 +4ϵ ( )12 − ( )6
R
R
bond

angle

LJ

(1)

粒子質量はすべての粒子で等しく計算中での単位質量とした。分子内相互作用は結合長ポテンシャルに
式 (1) 第一項と結合角ポテンシャルに式 (1) 第二項を用いた。

σ0 は平衡結合長を表し本計算の単位長とした。r は分子内隣接粒子間距離を表している。式 (1) 第一項及び第二項の係数 K はそれぞれ単位エネルギー ϵ0 を用いて Kbond = 200ϵ0 /σ02 , Kangle = 1.0ϵ0 とした。
分子間相互作用は式 (1) 第三項の LJ ポテンシャルを各粒子に対して用いた。R は他分子粒子間距離を表している。カットオフ半径は 2.5σ0 とした。式 (1) のパラメータ ϵ, σ は親水性粒子と疎水性粒子で異なる値をとるため、この値を変化させて分子間相互作用と構造安定性の関係性について議論する。当日は様々なパラメータについての計算結果を発表予定である。
項目
アンサンブル
境界条件
水分子数
脂質分子数
時間
時間刻み
設定温度
カットオフ

値
NPT
周期境界
512
3200
2000τ (200,000steps)
0.01τ
0.5ϵ0
2.5σ0

ϵ
Lipid Head
Lipid Tail
Waterσ
Lipid Head
Lipid Tail
Water

Lipid Tail
0.2
1.0

Lipid Head
1.0

Lipid Tail
1.2
1.0

Water
1.0
0.6
0.2
Water
1.0
1.2
1.0

表 2: LJ ポテンシャルパラメータ

表 1: 計算条件

3

Lipid Head1.0

結果
温度 0.5ϵ0 におけるシミュレーションを行った。図２ (a) にスナップショットを掲載した。赤, 灰, 青の粒

子はそれぞれ脂質親水部, 脂質疎水部, 水を表している。図２（b) には粒子種ごとの分布をプロットした。
赤線は脂質親水部の分布であり、赤線のピークを境にピンクの線で描かれた水の分布が極めて小さくなることから、脂質親水部のピーク間距離が脂質分子の層、即ち膜の厚さに相当すると考えられる。図２ (b) の脂質親水性粒子のピーク点の差から膜厚は 4.3σ0 と計算した。図２ (c) には膜面積の時間変化をプロットした。膜面積は 1.5σ02 で一定になっている。現在は平衡状態での計算を行なっているところであり、その結果を用いて構造安定性の解析を行う。温度以外の条件を統一し、kB T = 0.1ϵ0 〜1.0ϵ0 での計算も同時に行なっており、当日はその結果とパラメータを変化させて行った計算の両方の結果から、膜の構造安定性と分子間相互作用の関係について議論する予定である。
0.09

2
Lipid Head
Lipid Tail
Lipid Tail
Water

0.08

membrane surface area / lipid

0.07
1.5

Distribution

2

surface area (σ0 )

0.06
0.05
0.04

1

0.03
0.5

0.02
0.01
0
0

2

4

6

8

10

12

14

Z axis(σ0)(a) スナップショット(t = 2000τ )

0
0

20000

(b) 粒子の分布 (t = 2000τ )

40000

参考文献
[1] S. J. Marrink, et al., Biochimica et Biophysica Acta, 1778, 149-168(2009).

[3] S. Kawamoto, et al., J. Chem Phys, 134, 095103-095108(2011).

80000 100000 120000 140000 160000 180000 200000 220000
steps(c) 膜面積の時間変化

図 2: 膜の構造[2] G. Brannigan, F. L. H. Brown., J. Chem. Phys., 120, 1059-1071(2009).

60000

4P117
電解質中での OH-ダイナミクスと
局所水和構造に関する理論解析(九大稲盛セ)○石元孝佳・古山通久【序論】
バイオ燃料を直接利用した燃料電池は水素製造を必要とせず、バイオマスを用いて得られるグル
コースやエタノールなどの燃料をそのまま使用することが出来るという利点がある。しかし、これらの燃料は水素に比べて反応性に乏しく、燃料極に適した触媒の開発や反応機構の解明が求められている。近年の研究でアルカリ溶液を電解質として使用したアルカリ形燃料電池(AFC)で高いエタノー
ル酸化能を有することが実験により示された[1]。また、酸性条件下で動作する固体高分子形燃料電池(PEFC)とは異なり、アルカリ環境下にある金属触媒は OH-により被覆されている可能性も示唆されている[2]。これらの結果はアルカリ溶液が反応機構に大きな影響を及ぼしていることを意味している。以上の観点から、我々は金電極触媒上でのグルコースの酸化反応を解析した結果、アルカリ溶
液中の OH-が直接反応に関与することによりグルコン酸がスムーズに生成することを見出した[3]。
AFC の高性能化のためには電極反応はもちろんのこと、アルカリ環境下での OH-自身の拡散挙動
や OH-を含んだ局所水和構造を正しく理解することが不可欠である。本研究では、量子化学計算を用いて OH-の拡散挙動や局所水和構造の解析を行った。
【方法】
本研究では、図 1 に示すような(H2O)99 と OH-からなるクラスターを取り上げた。水の拡散挙動および動径分布関数は HF/STO-3G レベルの非経験的分子動力学計算より解析した。
また、OH-を含んだ局所水和構造の詳細を解析するために、プロトンの波動性を考慮した量子化
学計算を実行した。全ての計算は B3LYP レベルで行い、電子の基底関数には 6-31G**、プロトンの基底関数には[1s]GTF を用いて PCM 法によりクラスター周辺の溶媒効果を考慮した。

(H2O)OH-(H2O)99OH-(H2O)4OH-(H2O)2OH-(H2O)3OH-(H2O)5OH-

図 1 OH-含有水クラスター((H2O)99OH-)と局所水和構造モデル((H2O)nOH-, n = 1 ~ 5)

【結果】まず、(H2O)99OH-クラスター中での拡散挙動を解析した。図 2 には得られた最小二乗変位(MSD)
をプロットしたグラフを示す。ここでは比較のために(H2O)100 クラスターと重水モデル((D2O)99OD-、(D2O)100)の計算結果も合わせて示した。水のみからなる(H2O)100 クラスターよりも(H2O)99OH-の方が大きな MSD をとっており、水クラスター中に OH-が加わることで拡散係数が大きくなった。また水素(H)の同位体である重水素(D)からなる D2O モデルでは H2O に比べて拡散係数が小さくなった。これは質量の違いに起因するものである。酸素間距離に関する動径分布関数を解析すると、(H2O)100 ク
ラスターよりも(H2O)99OH-ではブロードなピークが得られ、乱雑な水素結合ネットワーク構造をとっていた(図 3(a))。一方、(H2O)OH-、(D2O)OD-モデルではピーク構造はほぼ一致しており、本計算手法
では水素結合強度の違いによる構造変化が十分に記述できていないことがわかる(図 3(b))。
0.3

4

4

0.1

0

(b)

(H2O)99OH(H2O)1003

0.2

gO-O(r)MSD (Å2)

(a)

3

2

2

1

1

0
0

1

2

3

4

5

Time (ps)

系列1(H2O)99OH系列2(D2O)99OD-002

図 2 OH-含有水クラスターの MSD

4

6

8

0

2

4

6

8

r (Å)
r (Å)
図 3 酸素間距離に関する動径分布関数

次に H/D の波動性の違いが安定化エネルギーや水素結合構造に与える影響を解析するために、
局所水和構造モデルを解析対象とした。OH-と水からなるモデル構造の安定化エネルギーを表 1 に示す。ここでは比較のため、通常の量子化学計算で得られた値も示している。H/D の波動性を考慮することで水素結合強度の違いが適切に記述され、従来法より大きな安定化エネルギーが得られた。
水分子の数が増加することで安定化エネルギーは大きくなっている。水 4 分子までは OH-(OD-)に直接相互作用した構造が最も安定な構造となっていた(図 1)。
これら 2 種類の解析手法により得られた詳細な解析結果については当日報告する。
【謝辞】

九州大学稲盛フロンティア研究センターの研究活
動は京セラ(株)の支援により行われた。また本研究の一部は「産学官地域連携による水素実証研究」
および科研費(25709012)の助成により行われた。
関係各位に感謝する。

表 1 局所水和構造モデルの安定化エネル
ギー (kJ/mol)(X2O)nOX1
2
3
4
5

X=H
-118.5-193.6-255.7-331.3-393.5

X=D-114.5-188.1-249.7-323.1-382.0

X = H(conv.)-105.0-175.3-235.9-306.2-357.8【参考文献】[1] Y. Kwon, S. C. S. Lai, P. Rodriguez, and M. T. M. Koper, J. Am. Chem. Soc., 133, 6914 (2011).
[2] P. Rodriguez, Y. Kwon, and M. T. M. Koper, Nat. Chem., 4, 177 (2011).
[3] T. Ishimoto, H. Kazuno, T. Kishida, and M. Koyama, Solid State Ionics, submitted.

4P118ポリヒドロキシスチレンにおける
パラトルエンスルホン酸拡散挙動の分子動力学的解析(名大院・情報科学 1、名大院・VBL2、京大・ESICB3、富士フイルム・ATC4)○野口高史 1,4、高柳昌芳 1,2、長岡正隆 1,3

Molecular dynamical study of diffusion behavior of
p-toluenesulfonic acid in polyhydroxystyrene(Graduate School of Information Science, Nagoya Univ.1,
VBL, Nagoya Univ.2, ESICB, Kyoto Univ.3, ATC, Fujifilm Corp.4)○Takafumi Noguchi1,4, Masayoshi Takayanagi1,2, Masataka Nagaoka1,3【序】半導体デバイス製造において、化学増幅型レジストは微細化の鍵を握る重要な材料である。化学増幅型レジストは、高分子、光酸発生剤と酸失活剤からなる材料で、光照射により発生した酸が高分子樹脂中を拡散することにより微細パターンを実現する。従って、酸の拡散係数は、レジストの性能 (解像性、Line Width Roughness) を左右する重要な物性であるが、拡散係数を規定する要因が必ずしも明確ではない。このため、理論計算による酸拡散挙動の解析が望まれているが、拡散係数の値が極めて小さいため、これまで殆ど解析されてこなかった。本研究では、化学増幅型レジスト
における代表的な高分子であるポリヒドロキシスチレン (PHS) 中のパラトルエンス
ルホン酸 (PTS) の挙動を例に、分子動力学 (MD) 計算を用いて、実測酸拡散係数の再現と酸拡散機構の解明を試みる。
【理論】拡散係数 D を求める方法として、Einstein 公式(1)、Green-Kubo 公式(2)、第二種揺動散逸定理を前提とした公式等が知られている。

<| x (t ) − x (0) |2 >D = limt →∞6t< v ( s )v (0) >dst →∞ 03

D = lim ∫

(1)

t

(2)

ここに、x は位置、v は速度、<A>は物理量 A のアンサンブル平均を表す。式(1)と式(2)は理論的には等価であると考えられているが、式中に含まれる極限操作を有限項で打切る数値計算においては等価では無い。例えば、公式適用の前提条件は、前者では平均二乗変位 (<|x(t)−x(0)|2>) が時間に対して線形になることだが、後者では速度自己相関関数 (<v(s)v(0)>) にエイリアシング (サンプリングに伴う折り返し誤差) が発生しないことが必要である。本研究では、MD 計算により得られたトラジェクトリに対して、式(1)と式(2)を適用し、両者の比較検討も行った。
【計算方法】化学増幅型レジストのモデル対象系として、8 量体 PHS 鎖 9 本と PTS 1分子を用いた。PTS は強酸のため解離していると想定し、H3O+と PTS−として取り扱った。PHS と PTS−には GAFF 力場と AM1-BCC 電荷を、H3O+には文献[1]の力場と電荷を用いた。初期構造は、低密度で配置した対象系を NPT アンサンブルで圧縮 (2ns@1 atm、2 ns@1000 atm、2 ns@1 atm の合計 6 ns) した後、PHS ガラス転移点より高い温度 (443 K) で緩和させ作製した。MD 計算は、Amber10 の sander モジュールを用い、NVT アンサンブル、周期境界条件、カットオフ半径 9Å、SHAKE 法により水素原子を含む結合を拘束し、時間刻み 1 fs にて行った。初期構造は 5 種、温度は PHS
ガラス転移点以下の温度で 5 種 (表 1) 用意し、それぞれ 52.4 ns の計算を行った。原
子速度は速度自己相関関数にエイリアシングが発生しない 10 fs の間隔で出力した。
1

【結果と考察】
1) H3O+実測拡散係数の再現
H3O+ の平均二乗変位は、劣拡散挙動(時間に対して上に凸) を示し、Einstein
公式(1)適用の前提条件が成立しなかった。このため、式(1)を適用して求めた拡散係数は、実測値に対して 4 桁程度大きな値を示した (表 1)。
これに対して、Green-Kubo 公式(2)では実測値と同程度の拡散係数が算出された (表 1)。この一致は、緩和時間が長く現実に実行可能な MD 計算時間では式(1)が適用できない系でも、式(2)による拡散係数算出が可能であることを示している。
2) 酸拡散機構の解明
式(2)を用いて H3O+と PTS−の拡散係数を計算した。何れも昇温に伴い拡散係数が増大したが、H3O+の拡散係数が常に大きかった (図 1)。また、両者の速度自己相関関数は異なる形状を示した (図 2)。
以上から、両イオンはイオン対を為さな
いこと、異なる拡散機構に基いていることが明確になった。なかでも H3O+は、PHS 側鎖 (−C6H4OH) の水酸基に囲まれて、速度自己相関関数の形状もケージ効
果を示していた。
一般に、正負イオンはイオン対として拡散すると考えられている[4]が、PHS側鎖配向やベンゼン環間の相互作用により、PHS 中での H3O+と PTS−は、イオン対を為さずに拡散していると言える。

表1

H3O+拡散係数の比較 (nm2/s)

温度(K)

Einstein
公式(1)*

Green-Kubo
公式(2)*

363

1.7×104

6.9

373

実測値

1.8-3.6[2]1.7×104

11.8

10.7×104

16.0 1.0-13.0[3]403

13.5×10

4

20.4

423

55.3×104

31.6

383
393

*5 本の MD 計算結果の平均値

図1【結論】MD 計算と Green-Kubo 公式を併用することにより、化学増幅型レジス
トの酸拡散挙動を解析する目途を得た。
従来、酸拡散係数の低減にはアニオンサ
イズの増大が有効であると考えられてきた[2,3]が、今回の解析から、アニオン
サイズの増大だけでは限界があることと、PHS 側鎖の配向制御が有効であることが示唆された。

図2

酸拡散係数の温度依存性

酸の速度自己相関関数 (393 K)【参考文献】[1] M. Baaden, M. Burgard, G. Wipff, J. Phys. Chem. B, 105, 11131 (2001)[2] T. Itani, H. Yoshino, M. Fujimoto, K. Kasama, J. Vac. Sci. Technol. B, 13, 3026 (1995)[3] H. Tsubaki, T. Tsuchihashi, T. Tsuchimura, Proc. SPIE7273, 72731K (2009)[4] E. L. Cussler, Diffusion 2nd Ed., p142, Cambridge UP (1997)2

4P119動力学計算に向けた鉄錯体の励起状態ポテンシャル関数の開発	
 （名大院・情報科学）○井内 哲、古賀 伸明
Development of model electronic Hamiltonian for dynamics studies of iron(II) complexes(Graduate School of Information Science, Nagoya University)○Satoru Iuchi, Nobuaki Koga
	  鉄(II)錯体が持つ多様な励起スピン状態の間で起こるダイナミクスは、光照射によるスピン
クロスオーバー現象などで注目されている。しかし、遷移金属錯体の励起状態では、内部変換、項間交差、振動緩和が競合的に起こるため、スピン状態変化や構造緩和を特徴づける要因はよくわかっていないことが多い 1)。例えば、水溶液中の [Fe(bpy)3]2+ (図１)では、光励起後の電荷移動励起状態から中間状態を経ずに 5 重項 d-d 状態(5T2)に緩和するという実験的知見が得られている 2)。しかし、複数の d-d 励起状態（図 2 参照）が存在する中で、そのような緩和が起こる理由はよくわかっていない。	
 

25

1

3

-1

Energy (10 cm )

20
15

T1
3

3

10

T1

T2

5

5
1

0
1.8

図１：トリスビピリジン鉄(II)錯体(bpy: 2,2′-bipyridine, 赤: Fe, 青: N, 黒: C, 白: H）

	
 

5

1.9

A1

2

2.1
Fe-N (Å)

E

T2
2.2

2.3

2.4

図２： 1A1 と 5T2 状態の安定構造を内挿した経路に対するポテンシャルエネルギー曲線。横軸は構造変化の指標としての Fe-N の長さ。本研究で改良したモデルを用いた結果。

	  このような実験研究を基に、発表者らは動力学シミュレーションによって励起状態ダイナ
ミクスを研究することを目指している。そのため、分子動力学シミュレーションで広く用いられるポテンシャル関数と同等の簡便性を持ちながらも、複数の励起スピン状態を同時に記述する方法の開発を進めている。これまでに	 [Fe(bpy)3]2+に対して簡便な有効ハミルトニアン
モデルを開発し、古典分子動力学シミュレーションを用いてモデルが水溶液中の 1 重項低ス
ピン状態(1A1)と 5 重項高スピン状態(5T2)を妥当に記述することを示した 4)。しかし、1, 5 重項
状態の交差領域を調べたところ、3 重項状態とのエネルギー関係も同時にバランスよく記述
するには、さらなる改良が必要であることもわかった 5)。	
 
	  そこで本研究では、光励起後のダイナミクスを扱う動力学計算に使用することを目的に、

広い領域で図２に示す 1, 3, 5 重項状態(1A1, 5T2, 3T1, 3T2, 1T1)をバランスよく記述できる有効ハ
ミルトニアンモデルに改良することを目指した。そして、[Fe(bpy)3]2+錯体に対する分子動力学シミュレーションを用いて、3 つのスピン状態の記述の妥当性を検討した。例えば、図３
は 1, 5 重項間の seam 上の構造をサンプルし、モデルと DFT 計算のエネルギーを比較したものである。ここでは、光励起後に生成する 5 重項状態(5T2)での熱平衡(300 K)を仮定し、seam
上に拘束した分子動力学シミュレーションを行うことで、様々な構造をサンプルしている。
図３から、定量性は改良の余地があるが、モデルが 1, 3, 5 重項状態(1A1, 5T2, 3T1)を同時にバランスよく記述していることがわかる 6)。
	  ダイナミクス研究には、さらにスピン状態変化の起因であるスピン軌道相互作用を、できる限り正確に記述することも重要である。そこで複数の構造においてモデルの枠組みでスピン軌道相互作用の大きさを計算し、CAS-CI 計算の結果と比較した。	
 図４は、3 つのスピン状態(1A1, 5T2, 3T1)の安定構造における比較であり、モデルが CAS-CI 計算に近い精度を持つことがわかる。	
 
600

35

400

-1

Model (cm )3-1Model (10 cm )

30

500

singletquintettriplet

300

25

200
100

20
20

25

3

-1

DFT (10 cm )3035

図 ３ : 	  50 点 の 構 造 に お け る モ デ ル と DFT(TPSSh/TZVP) 計算から得られたエネルギーの比較。分子動力学シミュレーションを用いて、1A1 と5
T2 状態間の seam 上の構造をサンプルした。エネル
ギーの基準は基底状態の最安定構造。

0
0

100

200

300
400
-1
CAS (cm )

500

600

図４: モデルと CAS-CI 計算によるスピン軌
道相互作用の大きさの比較。1A1, 5T2, 3T1 の安定構造における 1A1, 5T2, 3T1, 3T2, 1T1 状態が関係するもの。計算条件や内容の詳細は当日に
発表する。

1)

例えば M. Chergui, Dalton Trans. 41, 13022 (2012).

2)

例えば Ch. Bressler et al., Science 323, 489 (2009).

3)

例えば C. de Graaf and C. Sousa, Int. J. Quantum Chem. 111, 3385 (2011).

4)

S. Iuchi, J. Chem. Phys. 136, 064519 (2012).

5)

第 6 回分子科学討論会 2012 東京.

6)

DFT 計算では、3 重の擬縮退状態となる 5T2 と 3T1 の最低エネルギー状態のみ計算可能である。有効ハミルトニアンモデルでは、その性質上、3 重の擬縮退状態が同時に求められる。

4P-120
フッ素化アルキルスルホン酸ダイマーの赤外スペクトルの AIMD シミュレーション（技術研究組合 FC-Cubic1, 産総研 2）山口真 1、大平昭博 1,2

AIMD simulation of infrared spectra of fluorinated alkyl sulfonic acid dimers(FC-Cubic1, AIST2）Makoto Yamaguchi1, Akihiro Ohira1,2【背景と目的】固体高分子型燃料電池の実用化に向けた高温低加湿条件での性能向上には、湿度および温度変化に伴う電解質膜中の水やスルホン酸基の状態変化を把握し、プロトン伝導性との関連を明らかにすることが不可欠である。IR スペクトル測定は最も有効な分析手法であるが、湿度や温度に伴うスペクトルの変化は複雑であり、ピークの帰属には議論の余地が残されている。
Nafion などのフッ素化スルホン酸ポリマーの側鎖モデル化合物の IR スペクトルの 600～
1500cm-1 の領域は DFT 計算で非常に良く再現できたが 1)、水分量の変化に伴う酸の解離とプロトンの水和がもたらす IR スペクトルの変化はむしろ1500cm-1 以上の領域に顕著に現れる。
トリフルオロメタンスルホン酸 CF3SO3H は水素結合により環状ダイマーを形成し、図 1（実線）に示すよ
うに 3000cm-1 付近に OH 非対称伸縮振動によるブロー
ドな強い吸収帯を示す 5)。また、有機溶媒中の一水和物のスペクトルは、図１（点線）に示すように、ブロー
ドなピークが 2800cm-1 付近にシフトし、2200 および
1700cm-1 付近にこれよりも小さく幅の狭いピークが現れる 6)。これらのスペクトルは、高温での真空排気処理で完全に脱水した Nafion 膜、および室温で真空排気処

図 1 トリフルオロメタンスルホン酸液体(赤・実線)および一水和物溶液(青・点線)の実測 IR スペクトル

理した Nafion 膜の IR スペクトル 5)とそれぞれよく一致している。
こうした実験結果から、昨年の討論会では、H 型 Nafion や側鎖モデル化合物が低湿度下で2800cm-1 付近に示すブロードなピークがスルホン酸ダイマー水和物の酸解離構造に起因する可能性について、トリフルオロメタンスルホン酸の基準振動解析に基づき指摘したが、ピークの位
置や強度はよく再現できたものの、ブロードな形状などは本質的に再現できていなかった 2)。今
回の発表では、この吸収帯の特徴的な形状であるブロードな形状を第一原理分子動力学計算に基づき再現することで、吸収帯の帰属を確実なものとすることを試みた。
【計算方法】AIMD 計算には CP2K(Ver.2.4)を用いた 3)。立方格子の周期境界条件で構造最適化
および MD 計算 (NPT アンサンブル、ステップ幅 0.5fs) を行い、Travis コード 4)を用いて各原子および Wannier 関数の軌跡から IR スペクトルを求めた。
【計算結果】水素結合した OH 基のブロードな吸収帯の形状が MD 計算に基づく IR スペクトルでシミュレーション可能なことを確認するため、PBE 汎関数と DZVP 基底関数および GTH 型有効ポテンシャルを用いてトリフルオロメタンスルホン酸ダイマーの MD 計算を行い、IR スペクトルのシミュレーションを行った。得られた IR スペクトル（図 2（実線））には、3000cm-1 にブロードな吸収、2500cm-1 付近に OH 変角振動の結合音による小さいピークが現れており、この方法
により非調和性の強いブロードな赤外吸収帯がよく再現できることを確認した。

次に一水和物ダイマーの計算を行った。DFTによる構造最適化計算では、図 3 に示す酸解離型と非解離型の最適化構造のエネルギー差が非常に小さかったが、上述のダイマーと同じ汎関数と基底関数の組み合わせ、および PBE0 ハイブリッ
ド型汎関数と 6-31G*基底関数の組み合わせで、初期構造を解離型として MD 計算を行ったとこ
ろ、いずれの場合も最終的に非解離型構造へと移行して長時間存在する結果となった。このときのIR スペクトル（図 2（点線））には、2200cm-1付近に水素結合したスルホン酸の OH 基のブロ3700cm-1

ードな伸縮振動ピーク、3400 および付近に水の OH 伸縮振動ピークが現れ、図 1 に示図 2 トリフルオロメタンスルホン酸(赤･実線)
および一水和物(青・点線)ダイマーのIR スペクトルシミュレーション

したトリフルオロメタンスルホン酸一水和物の溶液の IR スペクトルとは明らかに異なっていた。
そこで、ヒドロニウムイオンに由来する IR ス
ペクトルの特徴的なピーク形状について確認するため、ヒドロニウムイオンを介したダイマー構
造をとることが知られているフッ化ホウ素酸水
和物ダイマー7) (BF4-･H3O+)2 の MD 計算を PBE0
汎関数+6-31G*基底で行った。図 4 に分子構造の図 3 トリフルオロメタンスルホン酸一水和物
ダイマーの解離および非解離最適化構造

スナップショットを示すが、2 個のヒドロニウム
イオンは MD 計算の全時間領域にわたり安定に存在し、2 個のフッ化ホウ素アニオンとの間の水素結合を次々と組み替えていた。得られた IR ス
ペクトルは、図 1 に示したトリフルオロメタンス
ルホン酸一水和物のスペクトルと比較すると高波数側にシフトしているが、約 1000cm-1 にわたりブロードな吸収が出現し、その低波数側の端に細いピークを伴い、また 1700cm-1 付近にもピークが現れており、全般的に特徴が再現されている。
これらの MD 計算の結果と実験結果の比較から、トリフルオロメタンスルホン酸一水和物のIR スペクトルに現れるブロードな吸収帯は、水素結合を組み替えながら運動しているヒドロニ

図 4 フッ化ホウ素酸一水和物ダイマーの構造のスナップショットおよびIR スペクトルシミュレーション

ウムイオンの OH 伸縮振動に帰属することが妥当であると現時点では考えている。今後は汎関数や基底関数の選択および計算時のパラメータ設定を見直し、トリフルオロメタンスルホン酸一水和物ダイマーでヒドロニウムイオンが安定に存在しうる条件の探索を継続する予定である。
【文献】1)Yamaguchi and Ohira, JPCA, 116, 10850 (2012). 2)第 6 回分子科学討論会、2P094 (2012).
3)www.cp2k.org. 4)www.travis-analyzer.de. 5)Buzzoni et al., JPC, 99, 11937 (1995). 6)Stoyanov et al.,JACS, 128, 1948 (2006). 7)Fărcaşiu et al., J.Mol.Model., 6, 171 (2000).
【謝辞】本発表は経済産業省資源エネルギー庁および独立行政法人新エネルギー・産業技術総合開発機構の「燃料電池先端科学研究事業」によるものである。計算の一部は東北大学サイバーサイエンスセンターの大規模科学計算システムで実施した。