博士論文
A Reactive Force Field Molecular-Dynamics Simulation Study on High-Density Ice
2021 年 9 月
足立 優司
岡山大学大学院
自然科学研究科
1
目次
1. Chapter 1 序論 I. はじめに
II. 高圧氷
III. Plastic 結晶 IV. 反応力場 : ReaxFF
2. Chapter 2 ReaxFF を用いた分子動力学シミュレーションによる高圧氷の構造と相挙 動
I. 研究の背景
II. シミュレーションモデルと計算手法
III. 結果と考察
IV. 結言
V. 補助資料
3. Chapter 3 高圧氷の構造と相挙動への ReaxFF ポテンシャルパラメーターの影響 I. 研究の背景
II. シミュレーションモデルと計算手法 III. 結果と考察
IV. 結言
4. Chapter 4 ReaxFF を用いた分子動力学シミュレーションによるカーボンナノチュー ブ内の高圧氷の構造
I. 研究の背景
2 II. シミュレーションモデルと計算手法 III. 結果と考察
IV. 結言
5. Chapter 5 結言 I. 本論文の結言 II. 引用文献
III. 論文発表および学会発表 IV. 謝辞
3
Chapter 1
序論
4
1-1 はじめに
氷 VII は 2 GPa 以上の高圧下において存在する高密度氷であり，2 つの低圧立方晶氷格 子の相互貫入により構成された体心立方構造を有する．近年，氷 VII の温度圧力領域にお いて様々な異状性を示す実験結果が報告されている．氷 VII の温度圧力領域における高密 度氷は従来想定されてきた以上に複雑な構造・ダイナミクス・相挙動を示すことが明らか になりつつあるが，その全容は未解明である．一方，氷 VII を含む高圧氷の理論的研究で は，古典分子動力学（MD）計算や第一原理（ab initio）MD 計算により Plastic 結晶の可 能性が報告されている．Plastic 結晶は，結晶中で分子がその格子点に固定されていても分 子の配向が定まっておらず，いわゆる回転拡散の状態にある．TIP4P/2005 等の剛体水モ デルを用いる古典 MD 計算は，大規模・長時間計算を行うことができるため，未知の相転 移の探索に有効である．ただし，剛体水モデルでは共有結合長・結合角が固定されている ため，それを用いる古典 MD 計算では，高圧下における分子構造変化や共有結合の架け替 えを伴う相挙動が正しく表現できない可能性がある．一方，ab initio MD 計算では，計算 コストが高いため，古典 MD に比べ原子数・計算時間を大幅に抑える必要があり，巨視的 な系の相転移を正しく評価できない場合がある．そこで本研究では，結合の生成・開裂を も取り扱うことができる反応力場である Reactive Force Field（ReaxFF）を用いて氷 VII の安定領域を含む温度・圧力領域において MD 計算を行い，高密度氷の構造・ダイナミク ス・相挙動を調べた．Chapter 1 では，本研究に関わる概要を述べる．Chapter 2 では ReaxFF を使用した MD 計算によって得られた水の高圧相挙動を報告する．Chapter 3 で は，Chapter 2 で明らかにした高圧相挙動への ReaxFF ポテンシャルパラメーターの影響 を報告する．Chapter 4 では，ReaxFF を用いた分子動力学シミュレーションによるカーボ ンナノチューブ内の高圧氷の構造と相挙動について報告する．Chapter 5 では，本研究の 結論を述べる．
5
1-2 高圧氷
水は水素結合を有する身近な物質の一つで，17 種類以上の結晶構造の存在が確認されて いる[1]．通常の氷とは異なる結晶構造は最初に Tammann らによって発見され[2]，それ以 降新たに発見された氷多形にはローマ数字が振り分けられた．これまで発見された水の固 相は多彩な構造を有しているため，完全な全体像は未だに明らかにされていない可能性が 高い．例えば氷中の水素結合のトポロジーとゼオライトのトポロジーとの類似性から，計算 上 74,963 種類もの可能な構造が提案されている[3]．これらすべての構造が安定相として 存在する可能性は低いが，一方ですべての構造が明らかにされたわけではない可能性を示 唆しており，将来的にはさらなる結晶構造が発見される可能性がある．水の固相は，水分子 の配向の乱れによって以下の 3 つに分類される．
1. 水分子の配向が無秩序である hydrogen-disordered 氷相 2. 水分子の配向が秩序化した hydrogen-ordered 氷相 3. H2O が分子特性を失う対称相
多くの hydrogen-disordered 氷相は，冷却すると水素原子位置が秩序化され hydrogen- ordered 氷相へと転移する[1]．つまり原則として水素無秩序相と水素秩序化された対応相 からなるペアが存在する．また超高圧下においては 3 番目の対称相が出現し，そこでは水 素結合と共有結合の区別がなくなり分子特性を失うため，水素の秩序相/無秩序相のペアは 存在しない．
氷の結晶構造は驚くべき多様性をもつが，一方で 2 GPa 以上の圧力となると氷 VII，氷 VIII，氷 X の 3 相のみが確認されており，その多様性は影をひそめる．図 1.1 に GPa 圧力 領域における H2O 相図を示す．氷 VII は，2 GPa 以上の高圧下において存在し，2 つの低圧 立方晶氷格子の相互貫入により構成された体心立方（bcc）構造を有している．氷 VIII は氷 VII の水分子の配向が秩序化した hydrogen-ordered 氷相である．氷 X は氷 VII の水素結合 が対称化した bcc 構造を有している．
6 図 1.1：GPa 圧力領域における高圧水の相図
氷を圧縮することにより水素結合が対称化するとの予測は 1960 年代にはすでになされ ていた[4]．加圧による水素結合変化は O–H⋯O 結合距離の減少をもたらし，いずれある圧力 で水素原子が O–O 原子間の中点に位置すると予測されている．Goncharov らは，対称化す る圧力，すなわち氷 VII から氷 X への相転移圧力が H2O で 60 GPa，D2O で 70 GPa 程度と 報告している[5]．これは氷 VII から氷 X への相転移を観察した最初の報告であり，その後 多くの研究が続くことになる．Benoit らの ab initio 計算によると[6]，O–O 原子間距離が かなり長い 2.85 Å では各プロトンは所属する酸素原子に結合したままであるが，圧縮され 2.43 Å になるとプロトン密度分布が二峰性を示し，さらに圧縮されることで対称化に至る．
プロトン密度分布が二峰性を示す状態は両端の酸素原子付近で最小値を持つ二重井戸型ポ テンシャルで模式的に表現される．
氷 VII から氷 X への相転移過程の詳細は，いまだに議論の中にある．Sugimura らは 40 GPa での圧縮率の急激な変化を報告しており，これはプロトンが非局在化した状態である
7
氷 VII’相であると予測した[7]．また Hernandez と Caracas は，氷 VII から氷 X に至る過程 で超イオン状態である Superionic 氷 VII 相が現れると ab initio 計算から予測した[8]．
Superionic 氷 VII 相は，共有結合の架け替えが発生する状態にありながら H2O 分子自身が 回転することによりプロトンが結晶中を拡散している状態にある．一方で Kuriakose らは，
氷中の音速と単結晶弾性率に関する研究の結果から氷 X への相転移圧力は少なくとも 82 GPa 以上であると報告しており[9]，氷 VII から氷 X への相転移過程のさらなる研究が求め られている．
氷 VII は 2 GPa から 60 GPa までの広い圧力領域で安定であり，この圧力領域の中では 構造を保っているように見える．しかし，およそ 10 GPa を境に低圧側と高圧側で異なるト レンドを示すことが多くの実験により報告されている．例えば，Okada らは 2 GPa から 40 GPa までの電気伝導度について調査した結果およそ 10 GPa にピークが存在すること報告 している[10]．また，Noguchi らは水素の拡散速度が 10 GPa で最大に達することを報告し ている[11]．最近では Komatsu らが氷 VII から氷 VIII への転移速度が 10 GPa で最小にな ることを報告している[12]．以上のように氷 VII の安定領域と思われていた 10 GPa 付近に おいて異常な動的挙動が存在している．
ここまで見てきたように，氷 VII 相の物性は「氷 X への相転移過程」と「10 GPa 付近の 異常な動的挙動」の大きく分けて 2 つの圧力領域において解明されておらず，その詳細に ついてさらなる研究が求められている．
8
1-3 Plastic 結晶
四面体配位のような三次元対称性が高い分子構造を持ち，比較的軽い元素からなる分子 が高圧条件下で固体を形成するとき，柔粘性結晶（Plastic 結晶）となることが知られてい る．一般的に結晶とは分子の重心位置に規則性がある状態をいい，一方で液体とは分子運動 により分子の配向や重心位置の規則性が失われた状態をいう．Plastic 結晶とは固体と液体 の中間の性質を持ち，結晶中で分子がその格子点に固定されていながらその配向が定まっ ておらず，いわゆる回転拡散の状態にある固相のことを指す．
Plastic 結晶の分類基準として水素結合の有無がある．水素結合が存在しない Plastic 結 晶の例として，四塩化炭素（CCl4）などの CX4系分子が挙げられる．CX4系分子はそれ自身 が四面体構造をとるため，圧力誘起構造相転移により面心立方格子の Plastic 結晶となるこ とが知られている[1]．水素結合が存在する Plastic 結晶の例としてとして硫化水素（H2S）
が挙げられる．H2S は H2O と比べ分子がより円形に近いため回転しやすいと考えられてい る．H2S は大気圧下で温度を下げると，液体相 - (187.6 K) - 固体 I 相 - (126.2 K) - 固体 II 相 - (103.5 K) - 固体 III 相の 3 つの固体相を示す[2]．I 相は面心立方格子の Plastic 相と呼 ばれ，回転拡散の状態にある．以上のように水素結合の有無に関わらず現在実験で観察され ている Plastic 結晶は面心立方格子をとることが知られている．
H2O の Plastic 結晶はまだ実験で発見されていない．また，面心立方格子構造の高圧氷は 最近発見された超高圧氷のみであり[3]，これらの違いについて詳細な研究が求められてい る．
9
1-4 反応力場 : ReaxFF
原子スケールの分子動力学シミュレーションは，新しい理論の発見や新規の材料開発の ための有用な手法を提供する．量子力学的手法（QM）に基づくシミュレーション手法は，
近年のソフトウェアの発達により非常に人気が高まっているが，一方で計算コストにより シミュレーションスケールが大幅に制限される問題を抱えている．QM と比較して大幅に少 ない計算リソースによって化学反応のダイナミクスを原子スケールで記述するために，半 経験的な方法による反応力場が開発されている[1-3]．
Reactive Force Field（ReaxFF）は，計算コストを下げるために精度をトレードオフにし て，QM で扱いにくい規模を超えるシミュレーションスケールの計算を可能にする[1]．
ReaxFF は，炭化水素に関する最初の報告[1]以来，高エネルギー材料，シリコン/シリコン 酸化物，ポリマー分解，不均一触媒，および界面を含む複雑な化学反応モデリングに適用さ れている[2]．ReaxFF は半経験的な原子間ポテンシャルを利用して，原子の位置関数として システムのエネルギーを計算する．調和振動で表現される角度ひずみやファンデルワール ス力，クーロン相互作用などは古典的な近似を適応し，原子の結合性の変化をモデル化する ための結合次数を介して反応性イベントを記述する．ReaxFF における系のポテンシャルエ ネルギーは，以下のエネルギー要素の合計で表現される．
𝐸
𝑠𝑦𝑠𝑡𝑒𝑚= 𝐸
𝑏𝑜𝑛𝑑+ 𝐸
𝑜𝑣𝑒𝑟+ 𝐸
𝑎𝑛𝑔𝑙𝑒+ 𝐸
𝑡𝑜𝑟𝑠+ 𝐸
𝑣𝑑𝑊𝑎𝑎𝑙𝑠+ 𝐸
𝐶𝑜𝑢𝑙𝑜𝑚𝑏+ 𝐸
Specific（式 1.１）
10
ReaxFF は，分極可能な電荷とともに結合次数を使用して，原子間の反応性相互作用と非 反応性相互作用の両方を記述する．これにより ReaxFF は，さまざまな分子の共有結合と 静電相互作用の両方を正確にモデル化できる．ReaxFF ポテンシャルへのエネルギーの寄 与は，次のように要約される．𝐸𝑎𝑛𝑔𝑙𝑒と𝐸𝑡𝑜𝑟𝑠は，3 体結合角ひずみと 4 体ねじれ角ひずみ に関連するエネルギーである．𝐸𝑜𝑣𝑒𝑟は，過剰な結合を防止するための原子価規則に基づい たエネルギーペナルティである．例えば，炭素原子が 4 つを超える結合を形成する場合，
エネルギーペナルティが適用される．𝐸𝑣𝑑𝑊𝑎𝑎𝑙𝑠と𝐸𝐶𝑜𝑢𝑙𝑜𝑚𝑏は，すべての原子間で計算される 静電力および分散力の関連するエネルギーである．𝐸Specificは，非共有電子対，共役，水素 結合などのエネルギーを含む．
ポテンシャルは，結合次数に依存するものと独立しているものに分けられる．結合次数 は，以下の式 1.2 を使用して原子間距離から直接計算される．
𝐵𝑂𝑖𝑗 = 𝐵𝑂𝑖𝑗𝜎 + 𝐵𝑂𝑖𝑗𝜋+ 𝐵𝑂𝑖𝑗𝜋𝜋
= exp [ 𝑝𝑏𝑜 1 (𝑟𝑖𝑗 𝑟𝑜𝜎)
𝑝𝑏𝑜 2
] + exp [ 𝑝𝑏𝑜 3 (𝑟𝑖𝑗 𝑟𝑜𝜋)
𝑝𝑏𝑜 4
] + exp [ 𝑝𝑏𝑜 5 (𝑟𝑖𝑗 𝑟𝑜𝜋𝜋)
𝑝𝑏𝑜 6
]
（式 1.２）
ここで，𝐵𝑂𝑖𝑗は原子 i と原子 j の間の結合次数，𝑟𝑖𝑗 項は原子間距離，𝑟𝑜 項は平衡結合長，
𝑝𝑏𝑜 項は経験的パラメーターである．式 1.2 は連続的で，σ，π，および ππ 結合特性 間の遷移による不連続性を含まない．これにより，原子間力の計算に必要な微分可能なポ テンシャルエネルギー曲面が生成される．この結合次数式は，遷移状態にある分子構造の 相互作用にも対応可能であり，反応障壁を正確に予測できる形式になっている．結合エネ ルギーや角度ひずみなど，結合次数に依存するポテンシャルの項は，修正された結合次数 から直接計算される．最後に，電荷平衡化処理が適用され，原子の部分電荷が計算され る．これを使用して，クーロン相互作用を計算する．ReaxFF は，古典的手法に結合次数
11
形式を組み合わせることにより，高価な QM 計算を行わずに化学結合を暗黙的に記述する ことができる．そのため ReaxFF は，QM と従来の古典的方法のシミュレーションスケー ルのギャップを埋めるのに役立つ可能性がある．
水は単純な分子構造を持っているにもかかわらず，非常に複雑な性質を持ち合わせてい る．水の微視的構造を巨視的特性に結びつけるために多くの H2O ポテンシャルが提唱され てきた．フレキシブル H2O モデルのほとんどは，調和振動による結合ポテンシャルと角度 ポテンシャルを採用しており，共有結合の生成と切断を表現することができない[4]．一方 で ReaxFF は，様々な結合をモデル化するために必要な項目が含まれているため，共有結 合の生成と切断が表現可能である．Duin らの ReaxFF によるバルク水のシミュレーション 結果によると[5]，ReaxFF による H2O 分子の振る舞いは，QM 計算による H2O 分子の振 る舞いだけでなく実験結果をも非常によく再現している．例えば，バルク水での O–O 動径 分布関数を ReaxFF による計算と実験の結果で比較すると，第一ピーク（2.8 Å），第二ピ ーク（4.4 Å），第三ピーク（6.7 Å）の溶媒和殻距離は極めて正確に再現されることが報告 されている．一方で定圧条件下でのバルク水密度について，実験結果を再現せず 90%程度 になることが報告されている．
ReaxFF の H2O ポテンシャルは水と他分子との相互作用を考慮することができるので，
水相における広範囲の複雑な化学プロセスの反応性分子動力学シミュレーションが可能に なる．ReaxFF を用いたバルク水シミュレーションの研究は，化学的プロセスの解明に広 く用いられている．例えば Rahaman らは銅錯体の特性を理解するための ReaxFF を用い た研究を報告している[6]．また Jeon らは金属酸化物と水溶液の界面で発生する化学反応 を理解するための ReaxFF を用いた研究を報告している[7]．金属酸化物と水性溶媒の界面 で発生する化学的相互作用は複雑であり，これらの現象を原子スケールから理解すること 重要な課題である．Sengul らは超臨界条件での酢酸-水混合物の凝集について ReaxFF を 用いた研究を報告している[8]．以上のように ReaxFF を用いたバルク水シミュレーション
12
は，多くの水溶液中のプロセスの理解を深めることにも大きく貢献する可能性を有する計 算手法である．
13
Chapter 2
ReaxFF を用いた分子動力学シミュレーション
による高圧氷の構造と相挙動
14
2-1 研究の背景
実験によって確認されている氷結晶構造の共通の特徴は，各分子に 4 つの水素結合した 隣接分子があることである．17 種類以上の結晶構造の存在が実験によって確認されている 一方で，分子動力学（MD）シミュレーション研究からは，水分子が格子サイトで自由に または頻繁に回転する Plastic 氷相の可能性が報告されている[1-5]．シミュレーションが 予想する Plastic 氷相は，氷 VII 領域の圧力と温度で安定して存在していることがわかる．
Takii ら[1]，Aragones ら[2]，および Aragones と Vega[3]は，水の剛体モデルを使用した MD 計算において，圧力一定条件下の温度上昇に伴い氷 VII が体心立方（bcc）格子の Plastic 相に相変化し，その後液体に転移することを報告している．また，面心立方
（fcc）格子の Plastic 氷相も，モデルの選択によっては安定相として存在する可能性が報 告されている．Aragones と Vega[3]は，水の一般的な 3 つのモデルの高圧相図を報告し た．これらの相図にはすべて，安定相として bcc と fcc の両方の Plastic 氷がある．bcc と fcc の Plastic 氷の構造的および力学的特性が詳細に研究されている．
実験ではまだ Plastic 氷相は確認されていないが，いくつかの実験的研究で氷 VII の安定 領域においていくつかの物性に特異な挙動が報告されている．ラマン分光法の研究によ り，𝜈1 OH 伸縮モードの圧力依存は 12 GPa で最小を示し[7，8]，氷 VII のプロトンの拡 散係数は最大で 10 GPa となり[9]，圧力依存のラマンモードは，13 GPa〜15GPa で異常 な変化を起こすこと[10]が分かった．X 線回折の研究では，11 から 14 GPa の間で圧縮す ると立方格子から正方晶構造への構造変化が示唆されている[11，12]．中性子回折研究で は、重陽子が酸素副格子の格子間サイトを占める構造への 13 GPa から始まる変化が示唆 されたが[13]，一方で最近の研究では侵入型モデルを裏付ける証拠は示されていない [14]．新しい NMR 手法に基づく最近の研究では，室温で 20 GPa と 75 GPa で不連続に変 化することを報告している[15]．これは最初の遷移が高障壁水素結合と低障壁水素結合を 持つ氷 VII の 2 つのドメインの間にあることを示唆している．以上のように氷 VII の異常
15
な構造的および動的挙動は，VII-X 相転移よりもはるかに低い圧力で観察されており，異 常の原因は解明されていない．
ab initio MD 計算において，プロトン密度分布が隣接する酸素原子のペア間で単峰性を 示す氷 VII から構造変化を起こし，二峰性を示す氷 VII’と呼ばれる相に至ることが報告さ れている[17-22]．最近の高圧 NMR 研究は，氷 VII 中の陽子の量子力学的トンネリングの 証拠を報告している[15]．最近になって Hernandez と Caracas は，ab initio 計算によ り，初めて bcc Plastic 氷の出現を報告した[23]．
この章では Reactive Force Field（ReaxFF）を使用した MD 計算によって得られた水の 高圧相挙動を報告する．ReaxFF は，その発明以来主要な力場に発展し，さまざまな反応 システムに適用されてきた[25]．反応性 MD は，分子があらかじめ定義されているのでは なく，構成原子から分子形成される．したがって，ReaxFF は，分子構造が変形する，あ るいは異なる分子種が形成される可能性があるより広い範囲の圧力および温度に潜在的に 適用可能である．さらに，数千の原子で構成されるシステムを妥当な計算コストでシミュ レーションすることができるため，巨視的システムの相転移を調べるときに ab initio MD 計算よりも計算コストの面で格段に優れている．本章では，まずモデルシステムとシミュ レーションの方法について説明する．次にこの水モデルのシミュレーション結果，構造解 析，動力学解析，および相図を示す．最後に研究の結論と今後の研究への展望を示す．
16
2-2 シミュレーションモデルと計算手法
ReaxFF を使用して，NpTおよびNVT MD 計算を実行した．この方法では，与えられた 位置にある原子で構成されるシステムのポテンシャルエネルギーは，結合エネルギー𝐸𝑏𝑜𝑛𝑑 とエネルギーペナルティ𝐸𝑜𝑣𝑒𝑟，3 体（結合角）エネルギー𝐸𝑎𝑛𝑔𝑙𝑒，4 体（ねじれ角）エネル ギー𝐸𝑡𝑜𝑟𝑠，ファンデルワールス𝐸𝑣𝑑𝑊𝑎𝑎𝑙𝑠およびクーロンエネルギー𝐸𝐶𝑜𝑢𝑙𝑜𝑚𝑏，および水素結 合などの相互作用を説明するために必要なエネルギー𝐸Specificの合計で表現される．本研究 における H および O 原子の ReaxFF パラメーターは，参考文献[26]の補足資料に記載され ているものを使用した．ReaxFF ポテンシャルは，バルク，界面，および閉じ込められた水 の MD 計算に使用されている[27–30]．バルク水中の H2O の分子構造は実験で決定された ものに非常に近く，バルク水中の構造的および動的特性の両方が実験結果と一致している [27]．以前のシミュレーション研究で発見され[31–33]，後に実験で確認された[34，35]，
カーボンナノチューブに閉じ込められた水とグラフェン層間の水の構造と相挙動 は，
ReaxFF ポテンシャルを使用した MD 計算でも十分に再現されている[28–30]．極端な条件 下で水をシミュレートするための ReaxFF の有効性を鑑みると，ReaxFF が水の剛体モデル が信頼できなくなる高圧および高温条件下での高密度の氷および液体の水に適応可能であ ると予想できる．
システムは，1024 個の酸素原子と 2048 個の水素原子で構成されている．すべての MD 計 算 は ， Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS, 17 November 2016 version)によって実施された．温度Tと圧力pは，Nosé–Hoover 温度制御 法と圧力制御法によって一定に保たれた．NpT MD 計算におけるシミュレーションセルは 周期境界条件下の直方体であり，3 次元の比率は固定されていない．NVT MD 計算におい て，周期境界条件下の立方体シミュレーションセルを使用した．初期構造として，酸素原子 が格子定数 a = b =3.1 Å および c =3.244 Å の氷 VIII を形成するように配置され，水素原子 は氷 VII の各副格子の酸素原子の中間点に配置された．平衡化を 8 GPa および 100 K で 10
17
ps 実行し，続いて 300 K で終了する等圧加熱プロセスを実行したところ，その過程でプロ トン無秩序構成を持つ氷 VII （bcc 格子）が得られた．また氷 X の構造から始めて，氷 VII の同じ構造が得られることも確認した．本研究で調べた熱力学的状態は，300 K から 600 K の温度範囲，4 GPa から 42 GPa までの圧力範囲である．温度と圧力は，それぞれ 5 K〜10 K と 0.01 GPa〜0.1 GPa ずつ段階的に変化させた．以下に示すように，この範囲の熱力学 的状態において，結晶-Plastic，結晶-結晶，および固液の自発的な相変化が観察された．
18
2-3 結果と考察
Plastic 氷の最大の特徴は分子の回転運動にある．ReaxFF ポテンシャルによってモデル 化された H2O の Plastic 氷相の存在を調べるために，最初にいくつかの圧力と温度条件で の水分子の再配向動力学を評価した．図 2.1 は，水分子の双極子モーメントの単位ベクトル u(𝑡) の再配向相関関数 C(𝑡) を示している．
C(𝑡) = 〈 u(𝑡) ∙ u(0) 〉
（式 2.1）
各分子の双極子モーメントは，分子内の電荷の中性と 2 つの H 原子の電荷対称性を仮定し て，H 原子と O 原子の瞬間的な位置から計算される．
図 2.1（a）に示されているのは，3 つの異なる温度（300 K，400 K，500 K，すべて 8 GPa 条件下）での C(𝑡) である．400 K および 500 K では，C(𝑡) が急速に減衰し 3ps 以内に消 失した．これは，水分子が 3 ps より短い時間で自由に回転することを意味する．一方 300 K では，相関関数は 3 ps で約 0.8 に減衰した．これは，水分子が頻繁に回転しないことを 意味する．図 2.1（b）に示されているのは，3 つの異なる圧力（8 GPa，27 GPa，40 GPa，
すべて 400 K 条件下）の C(𝑡) である．C(𝑡) は 8 GPa で非常に急速に減衰し，3 ps 以内に 消失した．40 GPa は 3 ps で[ C(𝑡) 〜0.2]であるが，27 GPa で相関が消えることはなかっ た．したがって C(𝑡) の圧力依存性は非単調であり，水分子は低圧域と高圧域で多かれ少な かれ自由に回転するが，中圧域では回転しない．この結果は，Plastic 氷相のリエントラン トな存在形態を意味する．図 2.1（c）の C(𝑡) の片対数プロットから，8 GPa と 40 GPa で の C(𝑡) が固有時間 𝜏：C(𝑡) = exp(−𝑡 / 𝜏 ) で指数関数的に減衰することが確認された．一 方，27 GPa, 400 K での C(𝑡) は，片対数プロットと両対数プロットからわかるようにサブ ピコ秒のべき乗則の減衰：C(𝑡) ~ 𝑡−𝛼 を示している．水分子の C(𝑡) におけるべき乗則の
19
減衰は，水素結合の再配列が短期間でたまにしか発生せず，そうでなければ水分子が水素結 合パートナーを変化させないことを示している[37]．
20 図 2.1：高圧氷中の H2O 分子の再配向相関関数 C(𝑡)
(a) 300 K，400 K，および 500 K（すべて 8 GPa）の C(𝑡) (b) 8 GPa，27 GPa，および 40 GPa（すべて 400K）の C(𝑡)
(c) （b）のC(𝑡)の片対数（上の行）および両対数（2 番目の行）プロット
21
表 2.1 に，C(𝑡) の指数関数減衰時間 𝜏 を示す．𝜏 の範囲は 0.35 ps〜9 ps であった．室 温での液体水中の水素結合の平均寿命が 2 ps 〜 3 ps という結果[38] を考慮すると，𝜏 の値が有限である氷は回転運動に関して「液体」，つまり Plastic 相と見なすことができ る．表 2.1 の空白のセルは，C(𝑡) が 3 ps 以内に感知できるほど減衰しないか，べき乗則 で減衰する場合である．
高圧での ReaxFF 液体 H2O の動的特性を確認するために，400 K で 0.1 GPa，1 GPa，2 GPa，および 3 GPa での水分子の並進拡散係数 𝐷𝑇 および回転拡散係数 𝐷𝑅 [39]を計算し た．𝐷𝑇 および 𝐷𝑅 の値は補足資料に記載されている．𝐷𝑇 は圧力の増加とともに減少する が，𝐷𝑅 はほぼ一定である．この傾向は，Bove ら[40] によって得られた実験結果と定性的 に同じである．ただし，ReaxFF ポテンシャルから得られた拡散係数は，対応する実験デー タよりも一貫して大きくなっている．
22
300 K 400K 500 K
40 GPa 8.78 2.02 1.00
27 GPa - - 0.58
8 GPa - 0.66 0.35
表 2.1：0.2 ps〜1.2 ps の間隔でフィッティングして得られた再配向相関関数 C(𝑡) の指数 関数減衰時間 𝜏（ps）．空白のセルは，C(𝑡)が 3 ps で 0.7 未満の値に減衰しない非 Plastic 状態に対応する．
23
ReaxFF は，従来の古典力場とは異なり各分子を相互作用する原子として扱うため，分 子は共有結合を解離して再結合する可能性がある．したがって H2O の回転運動の可能性に 加えて，2 つの隣接する酸素原子間の水素原子が 1 つのポテンシャル井戸から別のポテン シャル井戸に移動し，共有結合および水素結合酸素原子を交換する可能性がある．このよ うなプロトンホッピング運動は，ab initio MD 計算で観察される[23]．
H2O の回転運動とプロトンホッピング運動の 2 つの可能性を想定して，基準時間 𝑡0 か
らの ∆𝑡 = 3 ps の間の水素原子相対座標の分布を調べた．図 2.2（a）に示すのは座標系で
ある．原点は焦点を当てる水素原子に最も近い酸素原子の位置であり，x 軸は 2 番目に近 い酸素原子の位置によって定義される．座標系は 𝑡 = 𝑡0 で定義され，2 つの酸素原子の位 置が変動する間，∆𝑡 にわたって固定されたままになる．水素原子座標 𝑥𝐻，𝑦𝐻，𝑧𝐻は，∆𝑡
の間に時間発展する．図 2.2（b）は，∆𝑡 中の 𝑥𝐻 の分布 𝑓𝑥 を示している．これは，
0.1 ps の等間隔で 𝑡0 を選択し，システム内のすべてにわたって平均化されている．関数 𝑓𝑥(𝑥) は，共有結合 O–H が無傷で，x 軸方向に沿った水素結合が ∆𝑡 の間維持されている 場合，約 1 Å に単一のピークを持つ．水素原子が O–H⋯O と O⋯H–O に対応する 2 つのポテ ンシャル井戸間をホップする場合，𝑓𝑥(𝑥)は 2 つのピークを持ち，1 つは約 1 Å，もう 1 つ は他の酸素原子から約 1 Å になる．𝑓𝑥(𝑥)は，分子が ∆𝑡 の間に頻繁に回転する場合，𝑥 <
0 の範囲に裾野が広がる分布を示す．図 2.2（b）から，9 つのいずれの条件においても，
𝑓𝑥(𝑥)は 𝑥 > 0 の範囲で二重ピークがなく，∆𝑡 中にプロトンホッピング運動が発生しない ことが確認された．（300 K, 8 GPa），（300 K, 27 GPa），および（400 K，27 GPa）の 3 つ の条件では，𝑓𝑥(𝑥) に鋭い単一のピークがあり，これらの固体が結晶氷であることを示し ている．残りの 6 つの条件のうち 4 つ，（400 K, 8 GPa），（500 K, 8 GPa），（500 K, 27 GPa），および（500 K, 40 GPa）では𝑓𝑥(𝑥)は −1 < 𝑥 < 1 の範囲に広がり，他の 2 つの圧 力温度条件下（300 K，40 GPa）と（400 K，40 GPa）では分布の裾野は 𝑥 < 0 に広がる が，幅は前者よりも狭くなる．上記の 6 つの条件では，図 2.1 に示すように， C(𝑡) が指数
24
関数的に急速に減衰することが確認され，これらの固体が Plastic 氷であることを示して いる．
図 2.2（c）は，y，z 平面に投影された水素原子の分布𝑓𝑦𝑧を示している．（300 K, 8 GPa），
（300 K, 27 GPa），（400 K, 27 GPa）の 3 つの条件では，分布が他の場合よりもシャープ で，中央にピークがある．（400 K, 8 GPa）と（500 K, 8 GPa）の場合，分布は非常に広く，
分子がほぼ自由に回転していることを示している．最後に，40 GPa の 3 つの条件では，8 GPa での広い分布とは異なる，広い環状の分布が見られる．図 2.2（b）の対応する分布 𝑓𝑥(𝑥) とともに 𝑓𝑦𝑧 のリング状の分布は，O–H⋯O 原子が共線的である可能性が低く，水分子が歳 差運動を起こし，時折回転して水素を変化させることを示している．
25 図 2.2：隣接する酸素原子間の水素原子分布 (a) 𝑡0 で定義された座標系
(b) 𝑥𝐻 の分布𝑓𝑥
(c) 𝑦𝐻 と 𝑧𝐻の分布𝑓𝑦𝑧，トーンが明るいほど，分布は大きくなる
26
水の古典 MD 計算では，2 つの安定した Plastic 氷相が存在することが示されている．1 つは体心立方（bcc）格子で，もう 1 つは面心立方（fcc）格子構造である．bcc 格子と fcc 格子は，図 2.3（a）の z 座標をスケーリングするとマルテンサイト変態を介して相互 に変態する可能性がある．そこで Plastic 氷の格子構造を調べた．図 2.3 は，図 2.2 と同じ 9 つの圧力温度条件下における各酸素原子の最近傍酸素原子の空間分布を示している．図 2.3（a）で bcc 格子および fcc 格子の xyz 座標系が定義されている．図 2.3（b）は，x，y 平面での分布の投影を示し，図 2.3（c）は，x，z 平面での分布の投影を示している．9 つ の条件うち，3 つの条件（300, 8 GPa），（300 K, 27 GPa），および（400 K, 27 GPa）は，
氷 VII 相に対応し，他の 6 つの状態は Plastic 氷相に対応する．氷 VII の 3 つの条件のそれ ぞれで，x，y および x，z 平面上の投影に 4 つのスポットが表示され，8 つの最近傍とい う bcc 格子の特性を反映している．Plastic 氷（500 K, 27 GPa），（300 K, 40 GPa），（400 K, 40 GPa），（500 K, 40 GPa）の 4 つの条件では，氷 VII の場合と同様に 4 つのスポット がある．ただし，Plastic 氷の他の 2 つの条件（400 K, 8 GPa）と（500 K, 8 GPa）では，
x，y 平面の投影に 8 つのスポットがあり，x，z 平面の投影に 7 つのスポットがあった．
これは，それらの fcc 格子構造を反映している．2 つの条件でのシミュレーションボック スの相対的な寸法が fcc 格子および bcc 格子と一致していることも確認した．
補足資料には氷 VII，bcc Plastic，fcc Plastic の外観図を示している．S2.4（a）は氷 VII，S2.4（a）は bcc Plastic 氷，S2.4（a）は fcc Plastic 氷を示す．赤色が酸素原子，白 色が水素原子を表す．いずれの条件においても酸素原子が格子状に整列していることが確 認できる．S2.4（a）の水素原子が O-O 直線状に存在している確率が高いのに対して S2.4
（b）および S2.4（c）の水素原子は酸素原子周辺に一様に存在していることが確認でき る．
27
図 2.3：高圧氷の各中心酸素原子からの最も近い隣接酸素原子の密度分布 (a) 酸素原子の bcc および fcc 格子
(b) x，y 平面 (c) x，z 平面
28
補足資料には，（500 K，27 GPa）および（500 K，8 GPa）での O-O 動径分布関数 𝑔(𝑟) が示されている．27 GPa および 8 GPa で得られた 𝑔(𝑟) は，それぞれ，TIP4P/2005 ポテ ンシャルで得られた bcc および fcc Plastic 氷の 𝑔(𝑟) と類似していた．
fcc Plastic 氷は，氷 VII および bcc Plastic 氷よりも低い圧力（例えば，8 GPa）にあるこ とに注意が必要である．つまり fcc Plastic 氷は，同じ温度の氷 VII および bcc Plastic 氷よ りも密度が低くなる．一方で TIP4P/2005 ポテンシャルの場合，fcc Plastic 氷は bcc Plastic 氷よりも密度が高くなる．
fcc と bcc の両方の Plastic 氷が，ReaxFF のそれぞれの圧力，温度条件下で自発的に形 成されることを確認したので，次に 2 つの Plastic 氷相の間の相変化を調べた．図 2.4
（a）に，500 K での 8 GPa と 20 GPa の間の圧縮および減圧時の密度変化を示す．8 GPa から圧縮を開始すると，密度は 18 GPa ~ 19 GPa の間の圧力で不連続ジャンプが発生する まで連続的に増加し，bcc Plastic 氷が fcc Plastic 相から自発的に形成された．減圧する と密度は連続的に減少し 9 GPa ~ 10 GPa の圧力までヒステリシスループを形成した．こ こで，bcc から fcc への Plastic 相変化が発生し，その後密度は圧縮と同じ経路をたどっ た．
図 2.4（b）は，低圧，中圧，高圧の範囲での密度，エンタルピー，ポテンシャルエネルギ ーの変化を示している．ヒステリシスループの圧力範囲におけるポテンシャルエネルギー
∆𝑈，圧力-体積項 𝑝∆𝑉，および fcc と bcc の Plastic 氷間のエンタルピー ∆𝐻 の差を表 2.2 に示す．ここで，∆𝑋 は 𝑋𝑓𝑐𝑐− 𝑋𝑏𝑐𝑐 を意味する．fcc Plastic 氷のポテンシャルエネルギー は，3 つの圧力のいずれにおいても bcc Plastic 氷のポテンシャルエネルギーよりも低く，
圧力が高いほどその差は小さくなった．すでに述べた通り，各圧力での fcc Plastic 氷のモ ル体積は bcc Plastic 氷のモル体積よりも大きくなる．さらにもう 1 つの直感に反する事実 は，モル体積差 ∆𝑉 が圧力の増加とともに増加し，その結果圧力-体積項 𝑝∆𝑉が圧力に対し て線形よりも急速に増加することである．エンタルピー差
∆𝐻 = ∆𝑈 + 𝑝∆𝑉
は，この29
範囲で圧力が増加するとその符号が（負から正に）変化する．つまりこれは，ヒステリシス 範囲の平衡圧力に応じて，fcc と bcc の Plastic 氷の間の相境界の勾配 𝑑𝑝/𝑑𝑇 が正または 負になる可能性があることを意味する．∆𝑈 と 𝑝∆𝑉 を使用すると，相境界の無次元勾配
(𝑇/𝑝)(𝑑𝑝/𝑑𝑇) はクラペイロン方程式で与えられる．
𝑇 𝑝
𝑑𝑝
𝑑𝑇 = ∆𝐻
𝑝∆𝑉 = 1 + ∆𝑈 𝑝∆𝑉
（式 2.2）
30
図 2.4：500 K で 8 GPa から 20 GPa の間の圧力の増減に伴う密度，エンタルピー，および ポテンシャルエネルギーの等温変化．ヒステリシス間隔内の任意の圧力での低密度相と高 密度相は，それぞれ bcc Plastic 氷および fcc Plastic 氷である．青と赤の色は圧力を増減し たときに得られるデータを示す．
(a) 密度の変化
(b) 低圧，中圧，および高圧範囲での密度，エンタルピー，およびポテンシャルエネルギー の変化
31
𝑃 ∆𝑈 ∆𝑉 𝑝∆𝑉 ∆𝐻 𝑇
𝑃 𝑑𝑝 𝑑𝑇
9.9 GPa −1.04 0.046 0.46 −0.59 −1.29
14.0 GPa −1.17 0.110 1.54 0.37 0.24
18.2 GPa −0.64 0.107 1.95 1.31 0.67
表 2.2：fcc と bcc の Plastic 氷間のポテンシャルエネルギー ∆𝑈（kJ / mol），モル体積∆𝑉
（cm3 / mol），圧力-体積項𝑝∆𝑉（kJ / mol），およびエンタルピー∆𝐻（kJ / mol）の差 温度は 500K に固定されている．
32
図 2.5 は，等圧および等温密度の変化としての結晶性氷-Plastic 氷相転移を示している．
図 2.5（a）は，8 GPa で 200 K から 400 K に昇温したときの等圧密度の変化を示している．
これには，2 つの相転移に関連する密度ジャンプが含まれた．1 つは，氷 VIII から VII への 移行であり，密度のわずかな増加を伴った．もう 1 つは，密度の大幅な低下を伴う氷 VII か ら fcc への Plastic 相転移であった．図 2.5（b）は，300 K での圧縮と減圧による密度変化 を示している．低圧側の固相は氷 VII であり，高圧側の固相は bcc Plastic 氷である．
すでに見たように，ReaxFF 水モデルでは，fcc Plastic 相は，それぞれの相境界で氷 VII および bcc Plastic 相よりも密度が低くなっている．次に，隣接する酸素原子間の O–O 距離 𝑑𝑜𝑜 は，bcc 相（氷 VII または Plastic）から fcc Plastic 相への遷移時に増加する必要があ る．fcc と bcc の Plastic 氷相の間の遷移，および氷 VII と fcc の Plastic 氷相の間の遷移時 に，𝑑𝑜𝑜 がどのように変化するかを確認するのは興味深いことである．平衡圧力を 500 K で 14 GPa と仮定すると（図 5 を参照），bcc Plastic 性結晶の 𝑑𝑜𝑜 は 2.73 Å であるのに対し，
fcc Plastic 性結晶の 𝑑𝑜𝑜 は 2.82 Å であることが確認された．320 K および 8 GPa での氷 VII と fcc Plastic 氷の間の遷移の場合，氷 VII の 𝑑𝑜𝑜 は 2.79 Å であり，fcc Plastic 氷の 𝑑𝑜𝑜 は 2.89 Å であった．したがって，上記のケースは，fcc Plastic 氷の 𝑑𝑜𝑜 が氷 VII また は bcc Plastic 氷の 𝑑𝑜𝑜 よりも 3〜4％長いことを示している．
33
図 2.5：VII–VII，VII–fcc Plastic，および VII–bcc Plastic の相転移に関連する密度変化．
(a) 8 GPa での密度と温度
(b) 300 K での密度と圧力．青と赤の色は，減少および増加したときに得られたデータを 示す．
34
表 2.3 は，シミュレーションの等圧および等温プロセス中に観察された相変化をまとめた ものである．500 K および 600 K での等温経路での fcc–bcc Plastic 相変化，および 300 K での等温経路での VII - bcc Plastic 相変化でヒステリシスが観察された．
35
Fixed T or p Range of T or p Phase change 300 K 6 GPa–3 GPa fcc → liquid at 3.8 GPa 300 K 32 GPa–42 GPa VII → bcc at 38 GPa 300 K 42 GPa–32 GPa bcc → VII at 35 GPa 400 K 8 GPa–4 GPa fcc → liquid at 4.4 GPa 450 K 8 GPa–20 GPa fcc → bcc at 17 GPa 500 K 8 GPa–20 GPa fcc → bcc at 18.5 GPa 500 K 20 GPa–8 GPa bcc → fcc at 9.5 GPa 600 K 8 GPa–20 GPa fcc → bcc at 19.2 GPa 600 K 20 GPa–8 GPa bcc → fcc at 9.5 GPa 600 K 8 GPa–4 GPa fcc → liquid at 6.0 GPa 8 GPa 300 K–400 K VII → fcc at 320 K 8 GPa 700 K–800 K fcc → liquid at 750 K 14 GPa 350 K–450 K VII → bcc at 410 K 18 GPa 400 K–500 K VII → bcc at 440 K 27 GPa 400 K–500 K VII → bcc at 430 K 30 GPa 400 K–500 K VII → bcc at 420 K
表 2.3：ReaxFF MD シミュレーションで観察された相変化
36
2-3 結言
図 2.6 に，ReaxFF MD 計算から得られた H2O の相図を示す．この温度範囲で固液相境界 を持つ氷相は fcc Plastic 相である．融解曲線は正の勾配を持ち，fcc Plastic 氷が平衡状態 の液体よりも密度が高いことを示している．fcc Plastic 相は，約 410 K 未満の温度で氷 VII と相境界を持つ．また，約 410 K を超える温度で bcc Plastic 氷との相境界がある．bcc Plastic 相は，fcc Plastic 相との相境界の高圧側にある．fcc-bcc 相境界は，図 2.6 の圧力温 度範囲で温度軸にほぼ平行である．相境界の推定勾配を表 2.2 に示す．氷 VII，fcc Plastic，
および bcc Plastic 相の三重点は約 10 GPa および 410K にある．VII-bcc 相境界は，三重点 から開始して高圧領域に向かって伸び，その相境界には次の点がある．最高温度以下の圧力 では，氷 VII と bcc Plastic 氷の密度は同じである．最高温度を超える圧力では，bcc Plastic 氷は氷 VII よりも密度が高くなる．
図 2.6： ReaxFF を用いた MD 計算から得られた高圧水の相図
37
ReaxFF MD 計算から得られた H2O の相図は，剛体分子モデルの相図と比較して，類似 点と相違点の両方がある[3，5]． ReaxFF モデルにおいて，TIP4P/2005，TIP4P，および SPC / E モデルと同様に，fcc と bcc の両方の Plastic 氷相の存在を確認した[3]．ReaxFF モデルによる Plastic 氷は，剛体分子モデルの結果と一致して，実験によりに氷 VII に割り 当てられた領域の一部にあり，融解線を持つ安定相として存在した．TIP4P/2005 および SPC / E モデルの相図では，fcc Plastic 相は氷 VII の領域の高圧，高温側の広い領域を占 め，bcc Plastic 氷相は，氷 VII，fcc Plastic，および液体水との 3 つの境界で囲まれた低 圧領域にわずかに存在する[3]．顕著な違いは，fcc Plastic 相が低圧領域に存在し，bcc Plastic 相が高圧領域に存在することである．すでに述べたように，fcc Plastic 相は，氷 VII および bcc Plastic 相よりも密度が低くなっている．H2O の分子間相互作用は単純液体 の相互作用よりもはるかに複雑であるから，fcc 格子を持つ固相が相境界に bcc 格子を持 つ固相よりも密度が低いことは不可能ではない．実際 ab initio 計算から作成された相図に は，2 つの超イオン氷相の間に相境界があり，それに沿って fcc 相は bcc 相よりも密度が 低くなっている[41]．剛体水モデル，ReaxFF モデルおよび ab initio 計算から得られた相 図は，fcc 氷に対する bcc 氷の相対的な安定性が非常に小さいため，力場パラメーター，
計算手法またはその他条件の選択によって逆転する可能性を示唆している．ただし，これ らの不一致の原因メカニズムを特定することは，本研究の範囲を超えている．
Plastic 氷相の存在を調べる ab initio MD 計算研究はほとんどない．2018 年に Hernandez と Caracas は，NVTアンサンブル ab initio 計算から得られた Plastic 氷の最初の証拠を報 告し，Plastic 相を含む相図を提案した[23]．Plastic 氷相は，現在の低圧部分に存在するこ とがわかり，これは従来の古典的な MD 計算結果[3]および本研究結果と一致している．し かし ab initio 計算では，より高い圧力において Plastic 相は観察されず，代わりに 2 つの隣 接する酸素原子間のプロトン密度分布が二峰性を示す氷 VII’と呼ばれる固相が見つかった．
つまり，プロトンは O–H⋯O と O⋯H–O の間で位置が変化する．ReaxFF は，分子の解離を
38
含む化学反応を計算できるように設計されているため，氷 VII’などの相が観察される可能性 がある．しかし，本研究で調べた温度圧力領域では，bcc Plastic 氷相は，氷 VII の領域の高 圧側で安定していることがわかり，プロトンホッピング運動は観察されなかった．ReaxFF ポテンシャルを使用した予備的なNpT MD 計算では，氷 VII’に関連する 2 つの隣接する酸 素原子間のプロトンホッピング運動または並進は，ここで調べたよりもはるかに高い圧力 と温度，例えば 1000 K, 170 GPa と 850 K, 190 GPa で観察され観察され，実験と ab initio MD 計算によって非局在化が見られるものよりも高くなっている．この不一致は，現在の形 式の H2O の ReaxFF モデルが，超イオン相を含む相図を正確に再現していないことを示唆 している可能性がある．
次にシミュレーション結果を実験結果と比較する．まず本研究で調べた圧力範囲におい て fcc 構造を持つ氷を報告している実験研究はないことに注意が必要である．衝撃圧縮 X 線 回折に基づく最近の研究では，100 GPa を超える fcc 氷が報告されている[42]．Pruzan ら は， 𝜈1 OH 伸縮モードの幅が室温で約 12 GPa で最小であると報告した[7，8]．非単調な振 る舞いは，氷の水素結合がその圧力で最も強く，その点よりも高い圧力と低い圧力で著しく 弱いことを意味する．別の実験的研究は，氷 VII が異なる長さの水素結合を含むという証拠 を示しており，氷 VII の構造が以前に考えられていたよりも複雑であることを示している [43]．氷 VII の水素結合長の分布は，予測された Plastic 氷の分子配向の分布によるもので あるという説明が考えられる．最近の時間分解 X 線回折研究は，液体の水が急速な圧縮で 氷 VII に変化することを示しているが，その氷が結晶性であるか Plastic であるかは未定で ある[44]．リエントラントな挙動自体は ReaxFF ポテンシャルに固有のものではなく，
TIP4P/2005 モデルで観察される[3]．ただし，ab initio MD 計算では観察されない[17–22]．
最後に本章の結果をまとめる．ReaxFF ポテンシャルを使用した NpT アンサンブル MD 計算を実行して，H2O の高圧相挙動を調査した．fcc と bcc の格子を持つ 2 つの安定した Plastic 性結晶相は，H2O 分子の再配向相関関数とプロトンの空間分布によって確認された．
39
ReaxFF 水は，8 GPa の固定圧力で温度が上昇すると，VIII–VII および VII–fcc の Plastic 相 変化を示した．VII - fcc の相変化は，14 GPa 以下の圧力で温度を上げると観察され，300 K で圧力を上げると観察された．450 K 以上の温度で圧力を変えると，fcc – bcc の相変化 が観察された．fcc Plastic 相は氷 VII の低圧側に，bcc Plastic 性相は高圧側に存在した．
この特徴は，水の剛体モデルから得られた結果とは異なる．fcc Plastic 氷と bcc Plastic 氷 のプロトンの空間分布は互いに異なり 2 つの Plastic 相における H2O 分子の異なる再配向 運動を反映している．fcc Plastic 氷では，H2O 分子は自由に回転し，それによって酸素原 子の周りにほぼ均一なプロトンの分布を持つ．一方，bcc Plastic 性相では，各プロトンの 分布は，2 つの酸素原子を結ぶ軸を中心に環状であり，その軸に投影すると長い尾を持つ．
これは，水分子が歳差運動を起こし時折回転して水素結合パートナーを変化させることを 意味する．H2O 分子の 3 種類の特徴的な回転運動，低圧 fcc での「自由」回転，中圧氷 VII での「無」回転，および高圧 bcc Plastic 相での「歳差運動」回転は密接に関連している．
40
2-4 補助資料
並進および回転拡散係数 𝐷𝑇 および 𝐷𝑅 の数値結果を以下に示す．𝐷𝑇 は平均二乗変位の 傾きから得られ，𝐷𝑅は参考文献に記載されている式から得られた．実験データ 2 と比較す ると，ReaxFF ポテンシャルから計算された 𝐷𝑇 と 𝐷𝑅 はどちらも実験データよりも一貫し て大きいことが確認されたが，𝐷𝑇 と 𝐷𝑅の計算された圧力依存性は実験結果と定性的に同 じである．
𝑃 (GPa) 0.1 1 2 3
𝐷𝑇
(10−5 cm2 /s)
12 8.6 6.4 4.5
𝐷𝑅 (ps−1 ) 1.1 1.0 0.93 1.1
𝑃 (GPa) 0.1 1 2 3
𝐷𝑇𝑒𝑥𝑝 (10−5 cm2 /s)
9.2 5.2 3.7 2.95
𝐷𝑅𝑒𝑥𝑝 (ps−1 ) 0.44 0.47 0.48 0.49
表 S2.1：ReaxFF MD 計算から得られた 400 K での液体水の拡散係数
41 図 S2.1：O–O 距離の動径分布関数
(a) bcc Plastic 氷（500 K, 27 GPa）
(b) fcc Plastic 氷（500 K, 8 GPa）
42
図.S2.2： 600 K での 8 GPa から 20 GPa の間の圧力の増減に伴う密度，エンタルピー，お よびポテンシャルエネルギーの等温変化．ヒステリシス間隔内の任意の圧力での低密度相 と高密度相は，それぞれ bcc および fcc Plastic 氷である．
(a) 密度の変化
(b) 低圧，中圧，および高圧範囲での密度，エンタルピー，およびポテンシャルエネルギー の変化．
青と赤の色は，圧力を増減したときに得られるデータを示す．図 2.4 と比較するとヒステ リシス間隔内圧力領域がほとんど一致した．この結果は fcc と bcc の Plastic 氷の間の相境
界勾配 (𝑇/𝑝)(𝑑𝑝/𝑑𝑇) が極めて小さい値である結果と矛盾しない．
43
図.S2.3：VII–bcc Plastic の相転移に関連する密度変化（NVTアンサンブル）
図 2.5（b）と比較すると，ヒステリシス間隔内圧力領域がNVTとNpTでほとんど一致 し，相挙動がアンサンブルに依存したものではないことが確認された．
44
図.S2.4：氷 VII，bcc Plastic，fcc Plastic のスナップショット 赤：酸素原子，白：水素原子
(a) 氷 VII（300 K, 8 GPa）
(b) bcc Plastic 氷（500 K, 40 GPa）
(c) fcc Plastic 氷（500 K, 8 GPa）
45
Chapter 3
高圧氷の構造と相挙動への
ReaxFF ポテンシャルパラメーターの影響
46
3-1 研究の背景
ab initio MD 計算において，プロトン密度分布が隣接酸素原子間で単峰性を示す氷 VII から構造変化を起こし，二峰性を示す氷 VII’と呼ばれる相に至ることが報告されている[1- 6]．最近になって Hernandez と Caracas は，ab initio MD 計算よる研究において bcc Plastic 氷や氷 VII’，氷 X を報告した[7]．一方，従来の古典力場による MD 計算で氷 VII’
および氷 X を再現した報告はない．
ReaxFF には，Water-2010[8]と Water-2017[9]と呼ばれる異なる H2O ポテンシャルの系 統が存在する．Water-2017 は，Water-2010 の修正版でありプロトンホッピング運動のメ カニズムを説明するために提案された[9-11]．本章では Water-2017 ReaxFF ポテンシャル を使用した高圧氷のシミュレーション結果を報告する．Water-2017 はプロトンホッピング 運動のメカニズムを説明するために提案されたモデルであるため，プロトンホッピング運 動が Water-2010 で観察された圧力よりも低い圧力領域で観察される可能性がある．
本章では，まずモデルシステムとシミュレーションの方法について説明する．次に Water- 2010 と Water-2017 の双方でのシミュレーション結果，構造解析，動力学解析，および相 図を示す．最後に結論と今後の研究への展望を示す．
47
3-2 シミュレーションモデルと計算手法
ReaxFF を使用して，NpTシミュレーションを実行した．この方法では，与えられた位置 にある原子で構成されるシステムのポテンシャルエネルギーは，結合エネルギー𝐸𝑏𝑜𝑛𝑑とエ ネルギーペナルティ𝐸𝑜𝑣𝑒𝑟，3 体（結合角）エネルギー𝐸𝑎𝑛𝑔𝑙𝑒，4 体（ねじれ角）エネルギー 𝐸𝑡𝑜𝑟𝑠，ファンデルワールス𝐸𝑣𝑑𝑊𝑎𝑎𝑙𝑠およびクーロンエネルギー𝐸𝐶𝑜𝑢𝑙𝑜𝑚𝑏，および水素結合な どの特定の相互作用を説明するために必要なエネルギー𝐸Specificの合計で表現される．
モデルシステムは，1024 個の酸素原子と 2048 個の水素原子で構成されている．すべて の MD 計算は，Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS, 17 November 2016 version)によって実施された．温度Tと圧力pは，Nosé–Hoover 温度制御 方と圧力制御方によって一定に保たれた．シミュレーションセルは周期境界条件下の直方 体であり，3 次元の比率は固定されていない．酸素原子および水素原子の初期配置は Chapter 2 の計算結果から得られたものを採用した．本研究で調べた熱力学的状態は，300 K から 600 K の温度範囲，4 GPa から 70 GPa までの圧力範囲である．温度と圧力は，それぞれ 5 K〜10 K と 0.01 GPa〜0.1 GPa ずつ段階的に変化させた．以下に示すように，この範囲の 熱力学的状態において，結晶-Plastic および固液の自発的な相変化が観察された．
48
3-3 結果と考察
Plastic 氷の最大の特徴は分子の回転運動である．ReaxFF ポテンシャルによってモデル 化された H2O において Plastic 氷相の存在を調べるために，最初にいくつかの圧力と温度 条件での水分子の再配向動力学を評価した．図 3.1 は，水分子の双極子モーメントの単位ベ クトル u(𝑡) の再配向相関関数を示している．各分子の双極子モーメントは，分子内の電荷 の中性と 2 つの H 原子の電荷対称性を仮定して，H 原子と O 原子の瞬間的な位置から計算 される．
図 3.1（a）に示されているのは，3 つの異なる温度（300 K，400 K，500 K，すべて 8 GPa 条件下）での C(𝑡) である．400 K および 500 K では，C(𝑡) は急速に減衰し 3 ps 以 内に消失した．これは，水分子が 3ps より短い時間で自由に回転することを意味する．一 方 300 K では，相関関数は 3 ps で約 0.9 に減衰した．これは，水分子が頻繁に回転しない ことを意味する．図 3.1（b）にプロットされているのは，3 つの異なる圧力（8 GPa，27 GPa，40 GPa，すべて 400 K 条件下）の C(𝑡) である．C(𝑡) は 8 GPa で非常に急速に減 衰し，3 ps 以内に消失したが，27 GPa と 40 GPa で相関が消えることはなかった．この 結果は Water-2010 の結果と一致せず，モデル依存が存在することを示唆している．C(𝑡) の指数関数減衰時間 𝜏 を計算すると，8 GPa，400 K において 0.63 と Water-2010 と非常 に近い値を示した．この結果から（8 GPa, 400 K）と（8 GPa, 500 K）は Plastic 相である ことが示唆された．
49 図 3.1 高圧氷中の H2O 分子の再配向相関関数 C(𝑡)
(d) 300 K，400 K，および 500 K（すべて 8 GPa）の C(𝑡) (e) 8 GPa，27 GPa，および 40 GPa（すべて 400K）の C(𝑡)
50
Water-2017 は Water-2010 に比べ脱プロトン化しやすいモデルとして開発されている．
したがって Water-2017 は Water-2010 よりも低い圧力領域で共有結合および水素結合酸素 原子を交換する可能性がある．そこで H2O の回転運動とプロトンホッピング運動の 2 つの 可能性を想定して，基準時間 𝑡0 からの ∆𝑡 = 3 ps の間の水素原子相対座標の分布を調べ た．座標系は Chapter 2 の図 2.2（a）と同様で，原点は焦点を当てる水素原子に最も近い 酸素原子の位置であり，x 軸は 2 番目に近い酸素原子の位置によって定義される．座標系は
𝑡 = 𝑡0 で定義され，2 つの酸素原子の位置が変動する間，∆𝑡 にわたって固定されたままに
なる．水素原子座標 𝑥𝐻，𝑦𝐻，𝑧𝐻は，∆𝑡 の間に時間発展する．図 3.2（a）は，∆𝑡 中の 𝑥𝐻
の分布 𝑓𝑥(𝑥) を示している．これは，0.1 ps の等間隔で 𝑡0 を選択し，システム内のすべて にわたって平均化されている．関数 𝑓𝑥(𝑥) は，共有結合 O–H が無傷で，x 方向に沿った水 素結合が ∆𝑡 の間維持されている場合，約 1 Å に単一のピークを持つ．陽子が O–H⋯O と O
⋯H–O に対応する 2 つのポテンシャル井戸間をホップする場合，𝑓𝑥(𝑥) は 2 つのピークを持 ち，1 つは約 1 Å，もう 1 つは他の酸素原子から約 1 Å になる．𝑓𝑥(𝑥) は，分子が ∆𝑡 の間 に頻繁に回転する場合，𝑥 < 0 の範囲にテールが広がる広い分布を示す．図 3.2（a）から，
𝑓𝑥(𝑥)は 𝑥 > 0 の範囲で二重ピークがなく，9 つのいずれの条件においても ∆𝑡 中にプロト ンホッピング運動が発生しないことが確認された．（300 K，8 GPa），（400 K，27 GPa），
（500 K，27 GPa），（400 K，40 GPa），（500 K，40 GPa）の 5 つの条件では，𝑓𝑥(𝑥) に鋭 い単一のピークがあり，これらの固体が結晶氷であることを示している．残りの 4 つの条 件のうち 2 つ，（400 K，8 GPa），（500 K，8 GPa）では，𝑓𝑥(𝑥) は −1 < 𝑥 < 1 の範囲に広 がり，他の 2 つの条件（600 K，27 GPa）と（600 K，40 GPa）では，分布のテールは𝑥 < 0 に広がるが，前者よりも幅は狭くなる．上記の 4 つの条件では，図 3.1 に示すように，再配 向相関関数 C(𝑡) が指数関数的に急速に減衰することが確認され，これらの固体が Plastic 氷であることが示唆された．
図 3.2（b）は，y，z 平面に投影された水素原子の分布 𝑓𝑦𝑧 を示している．いずれの条件
51
においても中心にピークがあるが，一方で分布の広がり方は条件によって異なった．また Water-2010 で観察された広い環状の分布は見られなかった．以上のことから Water-2017 は，水分子が歳差運動を引き起こさず，Water-2010 とは異なる挙動が示唆された．
52 図 3.2：隣接する酸素原子間の水素原子分布
8 GPa は左から 300K, 400K, 500K，27 GPa と 40 GPa は左から 400K, 500K, 600K (a) 𝑥𝐻 の分布𝑓𝑥
(b) 𝑦𝐻 と 𝑧𝐻の分布𝑓𝑦𝑧，トーンが明るいほど，分布は大きくなる
53
Manzano らによる Water-2017 モデルの ReaxFF ポテンシャルを用いた超臨界水シミュ レーションによると，過剰なプロトン添加によってプロトンホッピング運動が観察された [10]．このことは，高圧氷にも適応される可能性がある．そこで水素原子を添加した系を調 べる．
図 3.3 は水素原子を 1 原子添加した系（500 K，40 GPa）における基準時間 𝑡0 から始ま る時間間隔 ∆𝑡 = 3 ps における 2 つの隣接する酸素原子間の水素原子分布を示している．
図 3.3（a）は，∆𝑡 中の 𝑥𝐻 の分布𝑓𝑥(𝑥)を示している．これは，0.1 ps の等間隔で𝑡0 を選 択し，システム内のすべてにわたって平均化されている．図 3.4（a）の左上には y 軸拡大 図を示す．図 3.3（a）から，𝑓𝑥(𝑥) には 𝑥 > 0 の範囲で二重ピークが存在し，プロトンホ ッピング運動が発生していることが確認された．図 3.3（b）は，y，z 平面に投影された陽 子の分布𝑓𝑦𝑧を示している．この結果は図 3.2（b）の結果と同様であり，水分子の回転運動 は違いがないことがわかる．
ここでプロトンホッピング運動が発生する確率について考える．この時 3 ps 以内に少な くとも一回はホッピングが発生する確率が 0.85%であった．この確率は以前報告された ab initio MD 計算によるホッピング確率と非常に近い値を示した．また Water-2010 モデル ReaxFF ポテンシャルでも同様にプロトン添加した系を調査した．その結果，（500 K，40 GPa）ではプロトンホッピング運動が観察されなかった．このことから，Water-2017 と Water-2010 のプロトンホッピング運動には差異が存在することが示唆された．
54
図 3.3：プロトン添加した系における隣接する酸素原子間の水素原子分布．（500 K，40 GPa）
(a) 𝑥𝐻 の分布𝑓𝑥 (b) 𝑦𝐻 と 𝑧𝐻の分布𝑓𝑦𝑧
