SENAC Vol. 52, No. （1 2019. 1）― 3 ―［共同研究成果］

複雑系に有効な非断熱分子動力学計算法の開発と分子モーター光異性化反応への適用
菅野 学,1 高橋みなみ,2 花崎 浩太,3 原田 宣之,4 河野裕彦 1
1
東北大学大学院理学研究科化学専攻
2 名古屋大学大学院理学研究科物質理学専攻
3
京都大学福井謙一記念研究センター
4 東北大学多元物質科学研究所

2016 年ノーベル化学賞の対象となった分子モーターは, 光や熱による異性化を利用して一方向の回転運動を達成した. このような光化学反応の機構を解明するためには, 異なる断熱電子状態
間の非断熱遷移を考慮することが不可欠である. 非断熱分子動力学の代表的な計算法である
surface hopping 法は, 断熱ポテンシャル曲面を走る古典トラジェクトリに遷移確率を与えて他のポテンシャル曲面へのホップ（遷移）を起こすことによって, 量子力学的な波動関数の分岐（生成物の分岐比）を求める. 近年, 1 次元非断熱問題の完全解である Zhu–Nakamura 理論を多次元に拡張して非断熱遷移確率を評価する新しい surface hopping 法が提案された. この方法は, 計算コストの高い電子状態間非断熱結合が不要で大きな分子の非断熱動力学計算にも適用可能であるが,ポテンシャル形状によっては遷移確率を誤る場合があった. そこで, 我々はこの問題を克服した改良法を開発した. 2 つの比較的小さな分子系に適用してその実用性を確認した後に, 分子モータ
ーの光異性化の経路を探索した.

1. 序論
非断熱遷移とは, 分子が電子の運動と原子核の運動（分子振動）の間に働く振電相互作用に起因して, 異なる断熱電子状態の間を乗り移ることを指す. 断熱電子状態は核座標をパラメーターとして扱った電子ハミルトニアンの固有状態であり, その固有値が形成する多次元曲面は断熱ポテ
ンシャル曲面と呼ばれる. 非断熱遷移は主に断熱ポテンシャル曲面の交差（1 次元では同じ対称性を持つ断熱電子状態のポテンシャル曲面は交差できないので擬交差, 2 次元以上の場合には円錐交差と呼ばれる）点の近傍で起こる. 複数の電子状態が関与する光化学反応において, 非断熱遷移は反応の分岐比を左右する重要な量子効果である.
非断熱遷移を伴う光化学反応の例として, 2016 年ノ bcーベル化学賞の受賞テーマとなった分子モーターが挙熱げられる. 図 1 はノーベル賞受賞者の Feringa と本稿の著者の 1 人である原田らが 1999 年に世界で初めて報告した一方向に回転する分子モーター[1]である. この分子はキラリティーを有しており, それによって回転光
光方向が一意に決まる. まず安定なトランス体 a に紫外d
光を照射すると電子励起が起こり, そこから電子基底 a状態に非断熱遷移することで反応障壁を越え, シス体
b へ異性化する. 次に熱を加えると, ナフタレン環同士
熱の立体障害を乗り越えてシス体 c になる. また光を照射するとトランス体 d になり, 加熱によってメチル基
図 1 世界で初めて合成された一方向とナフタレン環の反発を越えて初期構造 a に戻る. こに回転する分子モーター[1].
のように, キラル分子モーターは光異性化と熱異性化

― 4 ―

SENAC Vol. 52, No. （1 2019. 1）の組み合わせにより, 一方向の連続回転を実現させた.
一般に, このような光化学反応の機構（反応経路・時間など）を探る上で非断熱遷移の考慮は欠かせない. 分子モーターのような大きな系に適用可能で非断熱遷移を取り入れた分子動力学計算法が求められている. 非断熱分子動力学計算法の例として surface hopping 法が挙げられる. この方法では, Newton 方程式に従い断熱ポテンシャル曲面を走る古典トラジェクトリに遷移確率を与えて他の断熱ポテンシャル曲面へのホップ（遷移）を起こす. 適切な遷移確率の評価によって, 反応経路の分岐に対応する量子力学的な波動関数の分岐を再現する手法である. 代表的なものにTully の fewest switches アルゴリズム[2]があるが, これは時間発展のステップ毎に計算コストの高い電子状態間非断熱結合行列の評価が必要であり, 大きな系への適用が困難という欠点を抱えていた.
近年になって, Zhu らが大規模系への適用を可能とする非断熱結合計算の不要なアルゴリズムを提案した[3]. これは, 1 次元非断熱問題に対する完全解を与える Zhu–Nakamura 理論[4–7]を多次元に拡張して取り入れたものである. Zhu–Nakamura 理論では, 断熱電子状態の重ね合わせとして表される透熱電子状態の情報から非断熱遷移確率を得ることができる. この重ね合わせの係数はポテンシャル交差（擬交差や円錐交差）点において劇的に変化し, 透熱電子状態は交差点の前後で異なる断熱電子状態に漸近する. Zhu らが提案した確率計算アルゴリズムでは, 交差点における透熱電子状態のポテンシャル曲面（透熱ポテンシャル曲面）の勾配を断熱ポテンシャル曲面の勾配のみから推定し, 遷移確率を求める. しかし, この方法では透熱ポテンシャル勾配の推定に線形補間を利用しており, ポテンシャル曲面の形状によっては遷移確率の評価を誤ってしまう場合がある.
そこで, 我々は断熱ポテンシャル曲面の曲率を利用して透熱勾配を推定することによって上記の問題を克服した改良アルゴリズムを開発した[8]. 本稿ではまず, Zhu らの従来法が実験データを再現することが報告されているアゾベンゼン光異性化反応[3], および従来法が実験データを再現できないピラジン超高速無輻射失活[9]に 2 つの手法を適用した結果を比較し, 我々の改良法の実用性を確認する. そして, より複雑な系である分子モーター, 特に原田らが図 1 の二重結合で繋がった 6 員環を 5 員環に置き換えた新型の分子モーター（立体障害が小さく回転頻度が上昇することが報告されている[10,11]）に改良法を適用した結果から, 回転の第一段階であるトランス–シス
光異性化反応の機構の理解を試みる.

2. 非断熱遷移確率の定式化
最初に, Zhu–Nakamura 理論から導かれる非断熱遷移確率の表式とその多次元系への拡張について概説する. この遷移確率の評価および遷移方向の決定には, 多次元の透熱勾配の情報を要する.
そこで, Zhu らの従来法と我々の改良法における透熱勾配推定の違いを述べ, 従来法が失敗する場
合でも改良法は透熱勾配を正しく見積もれることを示す.

2.1

Zhu–Nakamura 理論の多次元系への拡張

Zhu–Nakamura 理論は 1 次元 2 状態のポテンシャル交差問題を扱う理論であり, 擬交差の近傍で透熱ポテンシャルを座標の 1 次関数, 透熱結合を定数とみなしたときの Schrödinger 方程式の解析解を与える[4–7]. 擬交差を 1 回通過したときの非断熱遷移確率 P は, 厳密にはエネルギーや結合強度の大きさによって場合分けされた複雑な式となるが, 分子が高いエネルギーを持って擬交差を通過する場合には以下の簡単な式で近似することができる.

P exp  
 4 a2
パラメーターa と b はそれぞれ

4b 1 


2
b2 

(1)

を通過する場合には以下の簡単な式で近似することができる.
複雑系に有効な非断熱分子動力学計算法の開発と
分子モーター光異性化反応への適用


― 5 ―




4
b 1 
2


P exp  4 a2


b2 

(1)

パラメーターa と b はそれぞれ
a2 

2
2

F1 F2 F1  F2

 2V12 

F1  F2b 2   Et  E X ,

3

F1 F2  2V12 (2)と定義される. (2)式において, μ は換算質量, F1 と F2 は透熱ポテンシャルの傾き（透熱勾配）, V12は透熱結合, Et は系の全エネルギー, EX は交差点におけるポテンシャルエネルギーである. (1)式の符号は F1F2 > 0（F1 と F2 が同符号）ならば+, F1F2 < 0（F1 と F2 が異符号）ならば−を意味する. (1)式や(2)式に非断熱結合が含まれていないことが Zhu–Nakamura 理論の優れた特徴の 1 つである.
2 原子分子より大きな多次元系では, 各原子の質量を考慮すべきであり, 透熱勾配はベクトルとなる. そこで, Zhu らは(2)式に現れる質量と透熱勾配を含む 1 次元パラメーターを次のように置き換えることで(1)式を多次元系に適用できる形式へと拡張した[3].

F1  F2





1

F1 F2

2

 m F   F  ,
i

1

i

2

i



i



1

 m F   F 
i

1

i

2

i

(3)

i

mi は i 番目の原子の質量であり, F1  と F2  は i 番目の原子の透熱勾配ベクトルである. 多次元系でもう 1 つ重要な点は, 遷移方向の決定である. Zhu らは, (3)式から遷移方向を次式で定義した.
i

i

si 

F1   F2 
i

i

(4)

mi

surface hopping 法では, ホップの前後でトラジェクトリのポテンシャルエネルギーが不連続変化するので, それを相殺するように運動エネルギー（運動量）を設定しなくてはならない. si に垂直な運動量成分は不変とし, 平行な運動量成分は全エネルギーEt の保存則を満たすように調整する.
ii
拡張した Zhu–Nakamura 理論を導入した surface hopping 法を実装する上では, F1  , F2  , V12, EXを量子化学計算から得られる断熱ポテンシャルとその勾配から見積もる必要がある. V12 と EX に関しては, 交差点において以下の関係式が厳密に成立する.

V12 

U  U,
2

EX 

U  U
2(5)

U+と U−は断熱ポテンシャル（U+ > U−）を表す. 残るのは F1  と F2  の推定である.
i

2.2

i

透熱ポテンシャル勾配の推定

Zhu らの従来法では, 次のように断熱勾配を線形補間することで交差点での透熱勾配を評価す
る.

FJi  t   U   t3  iU   t1  i1 R  t   Ri  t1   R  t   Ri  t3   

i
i
i
R
R  t3   R  t1   R

i

(6)

Riα は i 番目の原子の核座標（α = x, y, z）である. J = 1 のとき複号とはそれぞれ+と−に対応し, J= 2 に対しては逆になるとする. 次頁の図 2 は(6)式を模式的に示した概念図である. 図 2(a)のように時刻 t1 と t3 が十分に離れている（時間ステップが大きい）場合には, 透熱勾配を精度良く近似できる. しかし, 図 2(b)のように時間ステップが小さく, 断熱ポテンシャル曲面が滑らかな場合には, 交差点（t = t2）における透熱勾配を誤って等しく見積もり, (1)式の遷移確率が 0 になってしまう.
一方, 我々が開発した改良法では, 透熱勾配を次式で推定する（導出の詳細は文献[8]を参照されたい）.

― 6 ―

SENAC Vol. 52, No. （1 2019. 1）

 U  U

U




FJ   
 G
 R
 

2
4 G







(7)

∆U ≡ U+ − U− および ∆G ≡ U+ − U− と定義しており, (7)式はトラジェクトリに沿った断熱ポテ
 （実際の計算では数値差分）, 各原子の速度を使っ
ンシャル, 断熱勾配, 断熱勾配の時間微分 Gて求められる. (7)式の右辺第 2 項の分母は以下のように変形できる.
22
   U   U
G
R
t 2
t 2

(8)

(8)式には断熱ポテンシャルの時間に関する 2 階微分が現れている. つまり, 改良法における透熱勾配の推定には, 断熱ポテンシャルの傾きのみでなく曲率の情報も含まれている. 故に, 図 2(b)のように線形補間に基づく従来法が失敗する場合でも遷移確率を適切に計算することができる.

図 2 Zhu らの従来法における透熱勾配推定（線形補間）の概念図. 実線は断熱ポテンシャル, 破線は透熱ポテンシャル（赤破線が J = 1, 緑破線が J = 2）, 矢印はポテンシャル勾配を表す. (a) 成
功する場合. (b) 失敗する場合.

3. 従来法と改良法の比較
従来法と改良法を 2 つの異なる分子系に適用した数値計算の結果を比較する. 1 つ目の系はアゾベンゼンの光異性化反応, 2 つ目の系はピラジンの超高速無輻射失活である.

3.1

アゾベンゼン光異性化反応

アゾベンゼンは 2 つのベンゼン環がアゾ基（N = N）で結合した分子である（図 3）. シス体は環同士の立体反発に
より平面構造から歪むため, トランス体と比べて共役による安定化が小さい. アゾベンゼンの光異性化反応は古くから盛んに研究され, 溶媒中でシス体からトランス体への異図 3 アゾベンゼンの異性化反応.
性化（シス–トランス異性化）はフェムト秒, その逆反応（トランス–シス異性化）はピコ秒の時間スケールで進行することが知られている[12]. 両方向の異性化反応に関して, Zhu らは従来法が量子収率の実験値を再現すると報告した[3]. 我々は改良法も同様にアゾベンゼン光異性化反応の実験データを再現できるかを確認するための計算を実行した.
本稿では, シス–トランス異性化に焦点を絞ってその結果を示す.
Zhu らは surface hopping 法の代表的な計算プログラム Newton-X[13–15]に従来法を組み込んでおり, 我々はこのプログラムの提供を受けて改良法を追加実装した. 文献[3]に倣って, アゾベンゼ

複雑系に有効な非断熱分子動力学計算法の開発と分子モーター光異性化反応への適用― 7 ―ンの断熱ポテンシャルおよび断熱勾配の評価には量子化学計算パッケージ Molpro[16,17]を使い, 6
電子 6 分子軌道から成る活性空間を設定した状態平均多配置自己無撞着場（SA-CASSCF）
法[18,19]を採用した. 分子軌道の記述にはガウス基底関数系 6-31G[20]を用いた. Newton-X には複数の量子化学計算パッケージに対応したインターフェースが搭載されているが, 現時点で Molpro には非対応であるため, Molpro の出力を読み込めるように Newton-X のソースコードを修正した. 電子基底状態 S0 のシス体に相当する振動基底状態を反映した Wigner 分布[21]を仮定し, アゾベンゼンを構成する各原子の初期位置と初期速度をランダムに生成して第 1 電子励起状態 S1 に垂直励起した状
態を時刻 t = 0 とした. 古典トラジェクトリを走らせながら, 断熱ポテンシャルの差（∆U）の最小点を検出した. その点における透熱勾配を(6)式または(7)式に従って推定して(1)式の遷移確率を算出し, 0 から 1 の一様乱数との大小関係からホップの判断を行った. 遷移確率が乱数より大きけ
れば S1 から S0 へトラジェクトリをホップさせて(4)式で与えられた遷移方向の運動量成分を調整
し, 遷移確率が乱数より小さければそのまま S1 のポテンシャル曲面上を運動させた. 時間ステッ
プは 0.5 fs と設定した.
図 4 は従来法と改良法のそれぞれから得られた S1 の分布変化である. 両者は良く一致している
ことが分かる. 表 1 ではトランス体の量子収率と S1 の励起寿命の計算値を実験値と比較している.
これらの結果から, アゾベンゼン光異性化反応に対する従来法と改良法の計算精度はほぼ同じと言える. しかし, 量子収率の計算値は実験値と一致しているが, 励起寿命の再現性は比較的低い.
この反応に関しては, 電子状態計算の手法や基底関数系が概ね同一であれば, 従来法は fewest
switches アルゴリズムを用いた場合と近い結果を与えることが示されている[3]. そこで, fewest
switches アルゴリズムで算出された励起寿命の文献値を表 2 に並べた. 電子状態計算の手法を SACASSCF 法から時間依存密度汎関数（TDDFT）法に変えると実験値と合致している. このことから, 励起寿命の再現性が低い要因の 1 つは, 非断熱遷移アルゴリズムではなく電子状態計算の精度にあるとみなせる. 他の要因として, 実験が液相中で行われたのに対して計算は気相を仮定したことが影響している可能性も考えられる.
表 1 量子収率と励起寿命の計算値および実験値.
実験値[12]
（溶媒中）

従来法

改良法

トラジェクトリ数

167

170

量子収率

0.48

0.51

0.41–0.56

励起寿命（fs）

40.05

40.76

100–170

表 2 fewest switches アルゴリズムで計算された励起寿命. SA3-CASSCF(10,8)は 10 電子 8 軌道の活性空間で第 2電子励起状態 S2 も含む 3 状態を平均した SA-CASSCF 法.
図 4 S1 の分布の時間変化 .
赤線は従来法 , 青線は改良法
の結果.

3.2

励起寿命（fs）

SA3-CASSCF(10,8) [22]TDDFT [23]

67

121

ピラジン超高速無輻射失活

ピラジン C4H4N2 はヘテロ芳香環を有する生体分子（核酸塩基など）のプロトタイプとみなせる.
紫外光を吸収したピラジンは, 光学許容 S2 状態と光学禁制 S1 状態の間に存在する円錐交差を経由して超高速の無輻射失活を起こすことが知られている. 時間分解ポンプ・プローブ光電子分光実験から, 気相中における S2 の励起寿命は 22  3 fs と見積もられた[24].

― 8 ―

SENAC Vol. 52, No. （1 2019. 1）

図 5 ピラジン超高速無輻射失活の有効 2 自由度モデル. (a) 基準振動モード Q6a と Q10a. 白球は水素 H, 黒球は炭素 C, 青球は窒素 N を表す. (b) S2（青）と S1（赤）の 2 次元断熱ポテンシャル曲面.

Transition probability

これまでの我々の研究から, ピラジンの S2 から S1 への非断熱遷移は図 5 の有効 2 自由度モデルを使って良く記述できることが分かっている[9]. このモデルは, 図 5(a)に描いた面内骨格振動モード Q6a と CH 面外変角振動モード Q10a の 2 つの自由度を採用したものである. 図 5(b)は Molpro
を使用して多参照配置間相互作用（MRCI）法[25–27]と基底関数系 6-311++G**[20]から得られたS2 と S1 の 2 次元断熱ポテンシャル曲面である. これらのポテンシャル曲面上の量子核波束動力学
計算を実行したところ, S2 の励起寿命は 19.4 fs と算出され,1
上記の実験値を再現した[9]. そこで, 我々は図 5(b)の 2 次元
断熱ポテンシャル曲面の数値データを読み込んで surface
0.9
hopping 計算を行う Fortran プログラムを作成した. アゾベン
0.8
ゼンの場合と同様に, S0 の振動基底状態に対応する Wigner 分
布から初期条件のアンサンブルを作成して S2 へ垂直励起し,0.7
従来法と改良法の計算結果を比較した.
図 6 は横軸を時間ステップの対数として, 縦軸に遷移確率
0.6（円錐交差を 1 回通過して S2 から S1 へ移った分布の量）を0.5
プロットしたものである. 改良法は時間ステップの大きさ
0.01
0.1
1
をある程度変えても遷移確率がほぼ同じ値になった. S2 の励Time step / fs
起寿命を見積もると約 18 fs であり, 上記の実験値や量子核
波束計算値と一致した. 対照的に, 従来法は時間ステップを
図 6 遷移確率の時間ステップ小さくすると, 2.2 節で述べたように遷移確率を著しく過小依存性 . □は従来法 , ○は改良評価した. これらの結果から, 改良法は従来法より安定で適法の結果 . いずれの場合もトラ用範囲が広いことを確認できた.
ジェクトリの総数は 8000 本.

4. 分子モーター光異性化反応への適用
3 章の議論から改良法の実用性を確認できたので, 光駆動分子モーターの非断熱分子動力学計算に取り組む準備が整った. 図 1 の 6 員環が二重結合で繋がった分子モーターは熱異性化が回転の律速段階であり, b → c の過程
ではナフタレン環同士, d → a ではメチル基とナフタレン環の立体障害が大きいために回転速度が遅かった. これを改善するために, 原田らは立体障害が小さくなるように上記の 6 員環を 5 員環で置き換えた分子モーター（図7）を設計・合成し, 回転頻度を上昇させることに成功した[10,11]. 回転の第一段階であるトランス体 a からシス体ab

図 7 5 員環が二重結合で繋がった
新型分子モーターのトランス–シス
光異性化反応[10,11].

複雑系に有効な非断熱分子動力学計算法の開発と分子モーター光異性化反応への適用― 9 ―

b への光異性化は, メタノール溶媒中において光照射開始から 30 秒ほどで光平衡（a と b の比が
6 : 94）に達することが分かっている. しかし, トランス体 a 以外の過渡的な構造は観測されておらず, 詳細な回転機構は未だ明らかになっていない. 改良法をこの反応に適用し, 異性化機構の解明を目指した.

4.1

2 状態モデル

分子サイズを考慮し, 電子状態計算には時間依存密度汎関数強束縛（TDDFTB）法[28]を用いた.
TDDFTB 法は TDDFT 法のハミルトニアンをパラメーター化した半経験的手法であり, TDDFT 法に近い精度でより高速に電子励起状態を求めることができる. その中でも Kohn-Sham エネルギー汎関数を電子密度揺らぎに関して 3 次の項まで展開する TDDFTB3 法[29]を選択した. パラメーターセットとして 3ob-3-1[30], 計算プログラムには DFTB+[31]を使った. それに伴い, Newton-X のインターフェースを修正して DFTB+の TDDFTB ルーチンにも対応させた. トランス体 a に相当す
る振動基底状態の Wigner 分布を基に用意したトラジェクトリを S0 から S1 へ垂直励起し, その後の反応動力学を追跡した.
図 8 はトランス体 a からシス体 b への異性化を起こした典型的なトラジェクトリの結果である.
図 8(a)にはトラジェクトリに沿った S1 と S0 の断熱ポテンシャル, 図 8(b)には回転軸（C = C 結合軸）周りの二面角の時間変化を示した. 励起から 3.62 ps 後に S1 から S0 へ非断熱遷移が起こり, 追随するように二面角が著しく減少して異性化する. 図 8(c)は非断熱遷移に伴う構造変化を描いた
スナップショットである. 遷移直前の 3.50 ps の時点では二面角の値が 155.3°でありトランス体 aを形成しているが, 遷移を終えた 4.76 ps には 2 つのナフタレン環が近付いてシス体 b（二面角は−7.3°）になった. このように, 数 ps という極めて短い時間で異性化が起こる経路が見出された.

図 8 2 状態モデルにおいてトランス体 a からシス体 b への異性化を起こした典型的なトラジェク
トリ. (a) S1（紫線）と S0（黄線）の断熱ポテンシャルの時間変化. 緑線はトラジェクトリの時間発展に使われたポテンシャル. (b) 回転軸（C = C 結合軸）周りの二面角の時間変化. (c) 励起から 3.50
ps および 4.76 ps 後の分子構造. 紫色の 4 原子で回転軸周りの二面角を定義した.

4.2 5 状態モデル
前節で提示した計算結果は S1 と S0 の 2 状態を考慮したものであった. 続いて, 第 4 電子励起状
態 S4 まで含めた 5 状態モデルの下で同様の非断熱分子動力学計算を行ったところ, 得られたトラジェクトリは 2 状態モデルの場合とは異なる挙動を示した. 次頁の図 9(a)は 5 状態モデルの典型的なトラジェクトリにおける S1 から S4 までの断熱ポテンシャルの時間変化である. 4 つの電子励起状態が接近し合っており, それらの間で頻繁に非断熱遷移が起こっている. S0 には遷移しにくくなり, 回転軸周りの二面角は 220°から 170°の範囲に留まって異性化する様子は見られなかった.

― 10 ―

SENAC Vol. 52, No. （1 2019. 1）

図 9 (a) 5 状態モデルの典型的なトラジェクトリにおける S1 から S4 までの断熱ポテンシャルの時間変化. ○はトラジェクトリの時間発展に使われたポテンシャル. (b) 異性化の有効ポテンシャ
ル曲線.
図 9(b)は回転軸周りの二面角をパラメーターとして他の自由度を S1 状態上で最適化することで求めた有効ポテンシャル曲線である. この部分最適化には量子化学計算パッケージ Gaussian
09[32]を使い, TDDFT 法（汎関数は ωB97X-D[33], 基底関数系は 6-31G）を採用した. S1 の有効ポ
テンシャル曲線の傾きから, 励起直後はシス体へ向かって二面角が小さくなると予想される. しかし, 図 9(b)からも電子励起状態が互いに接近していることが分かる. S0 と同じく有効ポテンシャ
ル曲線の極小点が 200°近傍に存在している S2 や S3 へ遷移することで, 異性化が抑制されると考えられる. 以上の有効ポテンシャル曲線に基づく考察は, 分子動力学計算の結果と合致する.
異性化が起こりにくければ, トランス体 a に戻る経路が存在するはずである. この場合, 分子モ
ーターは再び光を吸収して異性化を試みることになる. 原田らの実験はメタノール溶媒中で行われているため, S1 より高い電子励起状態は溶媒効果でクエンチされる（急速に失活する）可能性が考えられる. 実際には 2 状態モデルの方が現実に近いのかもしれない. 今後は, 溶媒分子に囲まれた分子モーターのシミュレーションなども実施し, 検討を進めていきたい.

5. まとめと展望
動的機能分子の代表例として有望視される光駆動分子モーターの回転動力学を解明するために,複雑系に有効な新しい非断熱分子動力学計算法を開発した. これは Zhu–Nakamura 非断熱遷移確率の評価に必要な透熱ポテンシャル勾配を断熱ポテンシャル曲面の曲率を利用して推定することにより実現できた. 本稿で紹介した計算結果の解析から得られた成果を以下にまとめる.
＜アゾベンゼン光異性化反応＞
Zhu らの従来法と我々の改良法は同じ計算精度で実験データを再現した.
＜ピラジン超高速無輻射失活＞
改良法は従来法より時間ステップの大小に対して安定である.
＜分子モーター光異性化反応＞
2 状態モデルでは, 数 ps という極めて短い時間でトランス体 a からシス体 b へ異性化する経路が見つかった. 状態数を増やすと, 近接した電子励起状態の存在によって S0 への非断熱遷移及び異性化が起こりにくくなる.

複雑系に有効な非断熱分子動力学計算法の開発と分子モーター光異性化反応への適用― 11 ―

分子モーターを始めとするナノマシンの研究分野に対する社会の期待は大きく, 実用的な機能の発現が望まれている. 例えば, 原田らによって DNA に結合させた分子モーターの回転を利用して二重らせん構造を制御する試みがなされているが, 成功には至っていない[34]. 化学修飾による回転挙動の変化を予測できれば, 分子モーターに力学的な仕事をさせることが可能になると見込まれる. これはナノテクノロジーが掲げる究極の目標の 1 つであり, 産業・医療など多方面への応用が期待される.
図 7 の新型分子モーターでも光平衡に要する時間は秒の単位であり, 4 章の計算結果から示唆されるように様々な緩和経路が含まれると考えられる. 量子収率なども未知であり, 機構解明と高速回転可能な分子モーターの設計が喫緊の課題である. 大型計算機を活用したシミュレーション
研究は, 化学修飾や溶媒の存在によって複雑さを増す分子モーターの回転動力学の全貌を解き明
かす上で重要な鍵となる. 計算から得られる回転挙動の変化を系統的に整理して分子モーターの設計指針を確立し, ナノマシンとしての実用化に向けて一方向の高速回転に起因する機能発現の予測に繋げていきたい.
アゾベンゼンおよび分子モーターの光異性化シミュレーションに用いた Newton-X は, 世界中
で広く使われている分子動力学計算プログラムである. 通常の断熱的な（単一の断熱ポテンシャ
ル曲面上を走る）古典トラジェクトリの時間発展はもちろん, fewest switches アルゴリズムに基づいて非断熱遷移を考慮した surface hopping 計算も容易に実行可能である. 非商用・非営利目的であれば, ユーザー登録をすることで最新版である version 2.0, build 17（2017 年 5 月リリース）を無償で入手できる[35]. Newton-X のインターフェースは, 化学の分野で最も頻繁に利用されている
Gaussian（本来は有償だが, 東北大学サイバーサイエンスセンターがアプリケーションサービスとして提供している）を始め, 無償ソフトウェアである GAMESS[36]や COLUMBUS[37]など, 複数の量子化学計算パッケージに対応している. 我々が分子モーターの非断熱分子動力学計算に使った DFTB+も, 電子基底状態上の断熱的なトラジェクトリ追跡であれば対応済みである. 非対応のソフトウェアであっても, 我々が行ったようにその出力を読み込めるよう Newton-X のソースコードを書き換えてコンパイルすれば利用できる. 分子動力学のみでなく, 吸収・発光スペクトルや
光電子スペクトルを評価する機能も備えており, 利用価値の高いプログラムと言える.

謝辞
本研究の計算の一部は, 東北大学サイバーサイエンスセンターの並列コンピュータシステム
NEC LX 406Re-2 を利用して実行された. 従来法を組み込んだ Newton-X のソースコードを提供して下さった台湾国立交通大学の L. Yue 博士と C. Zhu 教授, および TDDFTB 法に関して貴重な助言を下さった Université Claude Bernard Lyon 1（フランス）の T. A. Niehaus 教授に感謝の意を表したい. なお, 本研究は JSPS 科研費 JP16H04091, JP26810002, JP18K05022 の助成を受けて行われた.

参考文献[1] N. Koumura, R. W. J. Zijlstra, R. A. van Delden, N. Harada, and B. L. Feringa, Light-driven
monodirectional molecular rotor, Nature 401, 152 (1999).
[2] J. C. Tully, Molecular dynamics with electronic transitions, J. Chem. Phys. 93, 1061 (1990).
[3] L. Zhu, C. Xu, Y. Lei, C. Zhu, and Z. Wen, Trajectory-based nonadiabatic molecular dynamics without
calculating nonadiabatic coupling in the avoided crossing case: trans ↔ cis photoisomerization inazobenzene, Phys. Chem. Chem. Phys. 16, 25883 (2014).
[4] C. Zhu, H. Nakamura, N. Re, and V. Aquilanti, The two-state linear curve crossing problems revisited.
I. Analysis of Stokes phenomenon and expressions for scattering matrices, J. Chem. Phys. 97, 1892
(1992).
[5] C. Zhu and H. Nakamura, The two-state linear curve crossing problems revisited. II. Analytical

― 12 ―

SENAC Vol. 52, No. （1 2019. 1）

approximations for the Stokes constant and scattering matrix: The Landau–Zener case, J. Chem. Phys.
97, 8497 (1992).
[6] C. Zhu and H. Nakamura, The two-state linear curve crossing problems revisited. III. Analytical
approximations for Stokes constant and scattering matrix: Nonadiabatic tunneling case, J. Chem. Phys.
98, 6208 (1993).
[7] C. Zhu and H. Nakamura, Two-state linear curve crossing problems revisited. IV. The best analyticalformulas for scattering matrices, J. Chem. Phys. 101, 4855 (1994).
[8] K. Hanasaki, M. Kanno, T. A. Niehaus, and H. Kono, An efficient approximate algorithm for
nonadiabatic molecular dynamics, J. Chem. Phys. 149, 244117 (2018).
[9] M. Kanno, Y. Ito, N. Shimakura, S. Koseki, H. Kono, and Y. Fujimura, Ab initio quantum dynamical
analysis of ultrafast nonradiative transitions via conical intersections in pyrazine, Phys. Chem. Chem.
Phys. 17, 2012 (2015).
[10] T. Fujita, S. Kuwahara, and N. Harada, A New Model of Light-Powered Chiral Molecular Motor with
Higher Speed of Rotation, Part 1 – Synthesis and Absolute Stereostructure, Eur. J. Org. Chem., 4533
(2005).
[11] S. Kuwahara, T. Fujita, and N. Harada, A New Model of Light-Powered Chiral Molecular Motor with
Higher Speed of Rotation, Part 2 – Dynamics of Motor Rotation, Eur. J. Org. Chem., 4544 (2005).
[12] H. M. Bandara and S. C. Burdette, Photoisomerization in different classes of azobenzene, Chem. Soc.
Rev. 41, 1809 (2012).
[13] M. Barbatti, G. Granucci, M. Persico, M. Ruckenbauer, M. Vazdar, M. Eckert-Maksić and H. Lischka,
The on-the-fly surface-hopping program system Newton-X: Application to ab initio simulation of thenonadiabatic photodynamics of benchmark systems, J. Photochem. Photobio. A 190, 228 (2007).
[14] M. Barbatti, M. Ruckenbauer, F. Plasser, J. Pittner, G. Granucci, M. Persico, and H. Lischka, NewtonX: a surface-hopping program for nonadiabatic molecular dynamics, WIREs: Comput. Mol. Sci. 4, 26(2014).
[15] M. Barbatti, G. Granucci, M. Ruckenbauer, F. Plasser, R. Crespo-Otero, J. Pittner, M. Persico, and H.
Lischka, NEWTON-X: a package for Newtonian dynamics close to the crossing seam, version 1.4, 2013.
[16] H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby and M. Schütz, Molpro: a general-purpose quantum
chemistry program package, WIREs: Comput. Mol. Sci. 2, 242 (2012).
[17] H.-J. Werner et al., MOLPRO, version 2012.1, a package of ab initio programs, http://www.molpro.net,
2012.
[18] H.-J. Werner and P. J. Knowles, A second order multiconfiguration SCF procedure with optimum
convergence, J. Chem. Phys. 82, 5053 (1985).
[19] P. J. Knowles and H.-J. Werner, An Efficient Second-Order MCSCF Method for Long Configuration
Expansions, Chem. Phys. Lett. 115, 259 (1985).
[20] I. N. Levine, Quantum Chemistry, 6th ed. (Prentice Hall, New Jersey, 2009), pp. 471–635.
[21] E. Wigner, On the Quantum Correction For Thermodynamic Equilibrium, Phys. Rev. 40, 749 (1932).
[22] M. Pederzoli, J. Pittner, M. Barbatti, and H. Lischka, Nonadiabatic Molecular Dynamics Study of thecis–trans Photoisomerization of Azobenzene Excited to the S1 State, J. Phys. Chem. A 115, 11136 (2011).
[23] A. J. Neukirch, L. C. Shamberger, E. Abad, B. J. Haycock, H. Wang, J. Ortega, O. V. Prezhdo and J. P.
Lewis, Nonadiabatic Ensemble Simulations of cis-Stilbene and cis-Azobenzene Photoisomerization, J.
Chem. Theory Comput. 10, 14 (2014).
[24] Y. Suzuki, T. Fuji, T. Horio, and T. Suzuki, Time-resolved photoelectron imaging of ultrafast S2 → S1
internal conversion through conical intersection in pyrazine, J. Chem. Phys. 132, 174302 (2010).
[25] H.-J. Werner and P. J. Knowles, An efficient internally contracted multiconfiguration-reference
configuration interaction method, J. Chem. Phys. 89, 5803 (1988).
[26] P. J. Knowles and H.-J. Werner, An Efficient Method for the Evaluation of Coupling Coefficients in複雑系に有効な非断熱分子動力学計算法の開発と分子モーター光異性化反応への適用― 13 ―

Configuration Interaction Calculations, Chem. Phys. Lett. 145, 514 (1988).
[27] P. J. Knowles and H.-J. Werner, Internally contracted multiconfiguration-reference configurationinteraction calculations for excited states, Theor. Chim. Acta 84, 95 (1992).
[28] T. A. Niehaus, S. Suhai, F. Della Sala, P Lugli, M. Elstner, G. Seifert, and Th. Frauenheim, Tight-binding
approach to time-dependent density-functional response theory, Phys. Rev. B 63, 085108 (2001).
[29] Y. Yang, H. Zhu, D. York, Q. Cui, and M. Elstner, Extension of the self-consistent-charge densityfunctional tight-binding method: Third-order expansion of the density functional theory total energy andintroduction of a modified effective coulomb interaction, J. Phys. Chem. A 111, 10861 (2007).
[30] M. Gaus, A. Goez, and M. Elstner, Parametrization and Benchmark of DFTB3 for Organic Molecules,
J. Chem. Theory Comput. 9, 338 (2013).
[31] B. Aradi, B. Hourahine, and Th. Frauenheim, DFTB+, a Sparse Matrix-Based Implementation of the
DFTB Method, J. Phys. Chem. A 111, 5678 (2007).
[32] M. J. Frisch et al., GAUSSIAN 09, Revision A.02, Gaussian, Inc., Wallingford, CT, 2009.
[33] J.-D. Chai and M. Head-Gordon, Long-range corrected hybrid density functionals with damped atomatom dispersion corrections, Phys. Chem. Chem. Phys. 10, 6615 (2008).
[34] F. Nagatsugi, Y. Takahashi, M. Kobayashi, S. Kuwahara, S. Kusano, T. Chikuni, S. Hagihara, and N.
Harada, Synthesis of peptide-conjugated light-driven molecular motors and evaluation of their DNAbinding properties, Mol. BioSyst. 9, 969 (2013).
[35] http://www.newtonx.org から入手可能.
[36] http://www.msg.chem.iastate.edu/gamess から入手可能.
[37] https://www.univie.ac.at/columbus から入手可能.