Dynamics of Doubly Ionized OCS Molecules and Hidden Reaction Pathways

Abstract

The global reaction route mapping (GRRM) strategy enables an exhaustive search of static chemical reaction pathways, providing deep insight into chemical reactions. On the other hand, ab initio molecular dynamics (AIMD) simulations explicitly account for the momenta of atoms, revealing the realistic dynamical motion of molecules. In this study, we constructed the reaction path network for the ground-state OCS²⁺ species, which has attracted significant interest from spectroscopists. Recently, a new dissociation mechanism of S⁺via the COS²⁺ isomer has been proposed. By referencing the results of AIMD simulations and focusing on the isomerization process, we emphasize the crucial role of dynamical effects and the concept of dynamically hidden reaction paths.

Translated Abstract

The global reaction route mapping (GRRM) strategy enables an exhaustive search of static chemical reaction pathways, providing deep insight into chemical reactions. On the other hand, ab initio molecular dynamics (AIMD) simulations explicitly account for the momenta of atoms, revealing the realistic dynamical motion of molecules. In this study, we constructed the reaction path network for the ground-state OCS²⁺ species, which has attracted significant interest from spectroscopists. Recently, a new dissociation mechanism of S⁺via the COS²⁺ isomer has been proposed. By referencing the results of AIMD simulations and focusing on the isomerization process, we emphasize the crucial role of dynamical effects and the concept of dynamically hidden reaction paths.

1 背景

OCS分子は準安定なジカチオン状態を持ち，光2重イオン化後のクーロン爆発解離過程を調べるためのモデルとして用いられてきた [1]．光イオン–光イオンコインシデンス（PIPICO）分光法は光照射によって気相分子を多重イオン化させ，そこから得られる解離生成カチオンを同時観測することにより，進行する解離反応の種類とエネルギー閾値を入射光エネルギーの関数として網羅的に測定する．1998年，MilliéらはPIPICO実験によりOCS²⁺からのCO⁺ + S⁺解離（S⁺解離），CS⁺ + O⁺解離と2種類の3体解離反応に関するエネルギー閾値（AE）を報告し，これらの反応への励起状態の関与を示唆した [2]．Endoらは近年，2倍波を混ぜることによって電場構造を非対称化したレーザー光を用いてPIPICO実験を行い，解離生成物に分配された並進エネルギーである運動エネルギー放出（KER）が2つのピークを持つことを見出して，異なる励起状態を始状態とする解離反応として帰属した [3]．これを受けてElandらは追証実験を行い，高感度装置を用いたコインシデンス実験により，Milliéらが報告したS⁺解離の閾値（AE^ES = 33.5 eV ± 0.5 eV）に加えてより低いエネルギー領域にS⁺解離の低強度なシグナルを発見した [4]．Elandらはこの新たに見出された反応様式はEndoらが報告したKERのピークに対応するものであるとしつつも，強レーザー場を用いた実験においてOCS²⁺からOS⁺ + C⁺の解離反応の進行が微弱ながら観測された結果 [5]を踏まえ，従来想定されていたOCS²⁺構造からの解離ではなくCOS²⁺への異性化の後に進行するS⁺解離機構を提案するとともに，基底状態において進行する過程も含む4つの解離機構を新たに提案し，その解離閾値（AE^GS? = 31.7 eV ± 0.4 eV）を報告した．ただしElandらは，電子状態計算によって求めた反応中間体種のエネルギーに基づきS⁺解離を議論しており，動力学効果については言及していない．主要な実験観測値をTable 1に示す．

Table 1. Comparison of observed thresholds for dissociation channels (in eV): appearance energy (AE) and kinetic energy release (KER).

	Milliéet al. [2]	Endoet al. [3]	Elandet al. [5]
Type	AE	KER	AE
Double ionization	30	–	29.5
CO+ + S+	33.5	4.0 & 5.2	31.7 & 34.1
CS+ + O+	40.2	5.6	36.9
O + C+ + S+	40.2	–	–
O+ + C + S+	47	–	–

大野，前田らは反応経路自動探索法を開発してGlobal Reaction Route Mapping（GRRM）プログラムを整備し，反応素過程を表す経路（IRC）からなる「反応経路ネットワーク」概念を創出した [6]．OCS²⁺の三重項基底状態について反応経路ネットワークを求めれば基底状態におけるS⁺解離反応の可能性を検証できる．ただし反応経路ネットワークは静的なものであり，動力学効果は考慮されない．近年当グループでは，CF₃⁺ + COに対し反応経路ネットワークが与える活性化障壁がイオンガイドビーム実験で観測された反応閾値に整合しない現象を見出し，反応動力学の観点から議論を行った [7]．当該研究において我々は，この不一致の原因は反応が迅速に進行することで平衡構造（EQ）において熱平衡が成立しないことが原因であると結論付け，IRCは存在するものの動力学効果により進行が阻害される過程を「動力学的に隠される反応経路」と名付けた．本稿ではこの描像を念頭に置き，Elandらが新たに提唱した4つのS⁺解離機構のうち電子基底状態におけるCOS²⁺異性化を経由する過程について，反応経路ネットワークと第一原理分子動力学（AIMD）シミュレーションからその妥当性を検討し，動力学的効果の観点から議論を行う．

2 計算手法

OCS²⁺の3重項基底状態（X̃ ³ Σ⁻）における反応素過程に対し，GRRM17プログラムに実装された非調和下方歪み追跡（ADDF）法を用いて反応経路ネットワークを決定した．ADDF計算では解離反応も調べるためにDown-DC=13とした．遷移状態構造（TS）の探索にはUB3LYPを用い，得られたTSをCASPT2(14e,12o) //CASSCFで再最適化した後にIRCを求めた．基底関数はすべてaug-cc-pVTZとした．さらに，所属研究室で開発しているSPPR [8]を用いてCASSCF(10e,10o)/cc-pVTZレベルでAIMDシミュレーションを行った．電子状態計算にはすべてMolpro2012を用いた．

3 結果

ADDF計算により，2個のEQと5個のTSを含む反応経路ネットワークが得られた．Figure 1に反応経路に沿ったエネルギー変化図とEQ, TSの構造，中性状態からのFranck-Condon（FC）構造におけるX̃ ³ Σ⁻状態とÃ ³ Δ, ã ¹Δ，ã ⁵ Σ⁺状態のエネルギー準位および実験で観測されたS⁺解離反応の閾値（青線）を示す．EQ0，EQ1はそれぞれOCS²⁺，COS²⁺に対応しており，この間の異性化反応のTS1のエネルギー31.1 eV（活性化障壁3.06 eV）は，Elandらが見出した解離閾値E^GS?を下回っていることから，IRC描像ではElandらの報告通りCOS²⁺を経由したS⁺解離が進行しうることが示唆された．なお，光イオン化の場合はAuger過程やcharge transferイオン化の場合とは異なり，ジカチオンの3重項状態に優先的に遷移することが指摘されており [2]，また実験により観測された二重イオン化エネルギーよりもFC構造における5重項最低状態のエネルギーは高いことから，これ以降では3重項状態のみについて考慮することとする．

Figure 1.

 (a) The energy profile for the reaction route map of OCS²⁺ in the ground state (X̃ ³ Σ⁻) calculated at the CASPT2 level of theory, including the energy levels of the à ³Δ, ã ¹Δ, and ã ⁵ Σ⁺ states at the Franck-Condon (FC) structure. The blue lines represent the experimental threshold for CO⁺ + S⁺ dissociation. The given energy is relative to the global minimum of neutral OCS; (b) The geometric structures of FC structure and each stationary point.

S⁺解離反応における反応動力学の影響を議論するために，FC構造からOCS²⁺のX̃ ³ Σ⁻状態で反応が進行する様子をCASSCFレベルのAIMDシミュレーションで調べた．CASSCFレベルでは，OCS²⁺構造からCOS²⁺構造に至る活性化障壁は3.53 eVであり，CASPT2よりも高めに見積もられた．はじめに300 Kでボルツマン分布に従うエネルギーを基準振動に基づく構造変位・運動量として与えて50個の初期条件を生成し，古典軌道を500 fsまで時間発展させたが，COS²⁺への異性化やS⁺解離は起こらなかった．これは，与えた初期エネルギーでは TS0, TS1, TS2を乗り越えることができなかったことを意味しており，OCS²⁺がマイクロ秒程度の寿命を持つ準安定なジカチオンである [4]ことを裏付けている．

Elandらが提案したCOS²⁺異性化の関与について調べるために，電子基底状態の振動励起状態に遷移した場合を想定して300 Kのボルツマン分布からサンプリングしたエネルギーに加え，FC構造におけるÃ ³Δ状態とX̃ ³ Σ⁻状態のエネルギーの差分4.423 eVを振動エネルギーとして各振動モードにランダムに加えたAIMDシミュレーションを行った^注1．50本の古典軌道のうち10本についてはS⁺解離に至らずEQ0近傍で激しい振動運動を続けたが，これらの古典軌道では初期条件でC-O伸縮や2つの変角振動モードに由来する分子軸回りの回転運動に大きなエネルギーが分配されていた．一方50本のうち40本がS⁺解離に至ったが，そのうち5本は直線形を保持してTS0近傍を経由したものであり，35本はOCS⁺の変角振動にエネルギーが多く分配され，CO⁺がS⁺に対して回転しながらS⁺解離が起こる挙動を示した．また、それらのうち1本についてはOCS²⁺ (EQ0) → TS1 → COS²⁺ (EQ1)と異性化の挙動を示したことから，詳しい解析を行った．

該当する古典軌道に沿った構造変化とエネルギー，内部座標の変化をFigure 2に示す．初期条件では，2つの伸縮振動に合計0.69 eVが分配されたのに対しO-C-S変角振動モードには3.79 eVが分配されており，r(C-O)があまり変化しないままCO⁺がS⁺に対して回転して50 fsでCOS²⁺に到達した．しかしその後さらに回転が続き，110 fsでOCS²⁺に戻った後にS⁺解離が起こった．一般に反応素過程において分子系は遷移状態を越えた後に分子内振動エネルギー緩和（IVR）により反応方向に集中したエネルギーが他の自由度へ分配されることで熱平衡に至るが，本系のようにCOS²⁺（EQ1）が準安定中間体でポテンシャル井戸が浅い場合には慣性で分子運動が進行し，S⁺に対するCO⁺のさらなる回転運動によってOCS²⁺（EQ0）へと構造が戻ってからS⁺解離に至るという動的描像が得られた．

Figure 2.

 (a) Variations in potential energy (upper) and structural parameters (lower); (b) Snapshot of the trajectory of OCS²⁺ → COS²⁺ → OCS²⁺ → CO⁺ + S⁺.

Figure 3.

The changes in interstructural distances between the structures along the trajectories and the EQ0 (OCS²⁺). structure. The dotted line indicates 0.2 Å.

本計算では原子核座標と速度の時間発展に速度ベルレ法を用いているため古典軌道は時間反転対称性を満たす．このことを用いて，OCS²⁺からCOS²⁺に到達できる初期条件について検討するために，両構造をつなぐTS1より50本の古典軌道を走らせた．なお，初期条件は各振動モードに対しボルツマン分布に従って生成し，反応座標方向についてはOCS²⁺方向に進行するように運動量を与えた．Figure 3に，古典軌道に沿った分子構造とOCS²⁺の平衡構造の間の構造間距離 [9] ^注2の時間変化を示す．これより，古典軌道はOCS²⁺平衡構造近傍に到達した後に単調に離れていく様子が可視化された．OCS²⁺平衡構造に対して最近接距離が0.2 Å未満となった29本の古典軌道に対し，OCS²⁺の基準振動を用いて最近接構造におけるC-S伸縮振動，C-O伸縮振動，O-C-S変角振動の振動量子数を計算し，Figure 4にヒストグラムの形で示した．この図は，OCS²⁺平衡構造付近からTS1に到達するために各自由度が持っている必要のあるエネルギーの程度を表している．特徴的なのはO-C-S変角振動モードの振動量子数が300を上回り，約0.1 hartreeほどの大きなエネルギーを持たないとTS1に到達するのは難しい点である．これはOCSが二重イオン化によりCOS²⁺へと構造変化することがそもそも難しいことを示している．すなわち，Elandらが提案したCOS異性体を経由するS⁺解離経路は，少なくとも電子基底状態では好まれないことが示唆された．

Figure 4.

 The distributions of vibrational quantum numbers for C-S/C-O stretching and O-C-S bending vibrations, estimated for molecular structures closest to the equilibrium structure of OCS²⁺ in the ground state (X̃ ³Σ⁻) along each classical trajectory.

4 結論

本研究では，GRRMによって計算される反応経路ネットワークとAIMDシミュレーションを用いてOCS²⁺から生じるCO⁺ + S⁺解離反応について，COS²⁺異性体の関与を考慮しながら議論した．静的な反応経路であるIRCでは，COS²⁺へと異性化した後にS⁺が解離するような反応経路は存在するが，分子運動を考慮したAIMD計算では，一度COS²⁺構造をとった分子は慣性の影響を受けてさらにOCS²⁺へと構造変化しS⁺解離が進行した．Elandらは，このような結合組み換え過程は多重イオン化された分子種の解離反応において主要な機構となりうると考えたが，この機構は動力学的には好まれないものであり，主要な機構とはなり得ないことが示唆される．

OS⁺ + C⁺の解離反応は，IRCが存在しているのにも拘らず実験的には非常に微弱なシグナルしか観測されていない．今後は，この反応についても動力学的に隠される反応経路という観点から解析を行っていく．また，Elandらは多配置・多参照理論を用いてポテンシャルエネルギー曲線を計算し，励起状態におけるCOS²⁺を経由するS⁺解離機構も提案している．今後の研究では，OCS²⁺の三重項励起状態における解離反応について非断熱遷移AIMDシミュレーションを駆使しながら検討する．

謝辞

本研究で行った計算の一部は，分子科学研究所の計算科学研究センター（RCCS）を利用して行った（Project: 23-IMS-C016, 24-IMS-C017）．

注釈

注1：振動エネルギー分布の偏りの実験的な観測はなされていない．一方でFC構造とEQ0の構造を比較するとC-O結合長はほとんど変化しないのに対しC-S結合長は約0.24 Å増加することから，三重項状態における初期運動はC-S結合伸縮振動が励起された状態であることが示唆される．

注2：古典軌道に沿った分子構造とEQ0の2つの分子構造に対してそれぞれの重心を座標原点に移動し， Kabschアルゴリズムにより座標軸をできるだけ重ね合わせた後，3N次元（Nは原子数）デカルト座標空間における直線距離を算出した．

参考文献

• [1]   P. Lablanquie, I. Nenner, J. H. D. Eland, J. Delwiche, M. J. Hubin-Franskin, P. Morin, Photophysics and Photochemistry above 6 eV, edited by F. Lahmani (Elsevier, Amsterdam, 1985), p. 53.
• [2]   P. Millié, I. Nenner, P. Archirel, P. Lablanquie, P. Fournier, J. H. D. Eland, J. Chem. Phys., 84, 1259 (1986). doi:10.1063/1.450518
• [3]   T. Endo, K. M. Ziems, M. Richter, F. G. Fröbel, A. Hishikawa, S. Gräfe, et al., Front Chem., 10, 859750 (2022). doi:10.3389/fchem.2022.859750 PMID:35464205
• [4]   M. Jarraya, M. Wallner, S. Ben Yaghlane, E. Olsson, V. Ideböhn, R. J. Squibb, et al., Phys. Chem. Chem. Phys., 25, 19435 (2023). doi:10.1039/D3CP01688A PMID:37435640
• [5]   S. Zhao, B. Jochim, P. Feizollah, J. Rajput, F. Ziaee, K. Raju, et al., Phys. Rev. Accel. Beams, 99, 053412 (2019).
• [6]   S. Maeda, Y. Harabuchi, Wiley Interdiscip. Rev. Comput. Mol. Sci., 11, e1538 (2021). doi:10.1002/wcms.1538
• [7]   K. Oda, T. Tsutsumi, S. Keshavamurthy, K. Furuya, P. B. Armentrout, T. Taketsugu, ACS Phys. Chem. Au, 2, 388 (2022). doi:10.1021/acsphyschemau.2c00012 PMID:36193292
• [8]   T. Tsutsumi, K. Oda, L. Qu, J. Sakugawa, Y. Minegishi, A. Fujita, et al., SPPR (a developmental version), 2022, Hokkaido University: Sapporo, Japan.
• [9]   T. Tsutsumi, Y. Ono, Z. Arai, T. Taketsugu, J. Chem. Theory Comput., 14, 4263 (2018). doi:10.1021/acs.jctc.8b00176