修 士 学 位 論 文
題 名
ナ ノ 細 孔 に お け る 水 の 相 転 移 と 輸 送 : カ ー ボ ン ナ ノ チ ュ ー ブ を 用 い た 研 究
指 導 教 員
真 庭 豊
教 授
2 0 20
年
1月
1 0日
提 出
首都大学東京大学院
理 学 研 究 科
物 理 学
専 攻 学修番号 18844411
氏 名 小 倉 宏 斗
学位論文要旨（修士（理学） ）
論文著者名 小倉 宏斗
論文題名：ナノ細孔における水の相転移と輸送
：カーボンナノチューブを用いた研究
[研究背景と目的]
ナノサイズの空間に閉じ込められた水は、バルクとは異なる性質を示す。水を閉じ込め ることが可能なナノ空洞を提供する物質の1つに、カーボンナノチューブ（CNT）がある。
CNT はナノメートル単位の直径を持つ円筒状物質である。CNT に内包された水の構造や相 転移挙動は、空洞直径Dに著しく依存することが報告されている。例えば、D < ~1.5 nm で は、内包水は低温で筒状の氷（ice-NT）を形成し[1]、D > ~1.6 nmでは、低温でwet-dry 転 移（水がCNT外部へ排出される現象）を起こす[2]。CNT内の水の挙動は空洞のサイズだけ でなく、空洞形状にも強く影響されると予想される。しかし系統的にCNT形状を変えて内 包水の挙動を調べた研究は、我々の知る限りではまだ報告されていない。
ナノ細孔内の水は、その動的挙動もバルクとは異なることが指摘されている。例えばCNT の円筒空洞内に水を流すと、その水輸送速度はマクロな流体力学から予測される値よりも4 ケタ以上大きいという報告がある[3]。このような特異な水輸送特性が現れる要因として、
空洞壁と水との間に摩擦が生じない（もしくは非常に摩擦が小さい）とする説がある。こ の摩擦には内包水のナノ構造が大きく影響すると予想されるが、現時点では理論・実験共 に数が限られており、その詳細な理解に至っていない。
そこで本研究では、(1)内包水の挙動の空洞形状への依存性を系統的に明らかにする、(2) 水のナノ構造と摩擦に着目し、CNT空洞内での高速水輸送の機構を明らかにする、(3)変形 CNT内の水の挙動を実験により検証することを目的とした。
[実験方法と結果]
(1) 変形CNT内の水：古典分子動力学（MD）計算を用いて、一軸方向の外力で変形させ たCNTに内包した水の構造と相転移挙動を調べた（図1）。CNTはD = 1.24nmと1.51nm の2種類とし、水分子にはSPC/Eモデルを用いた。その結果、内包水の液体-固体様相転移 温度と低温構造はCNTの変形割合γ (%)に著しく依存することが明らかになった（図2，図 3）。とくに、扁平化したCNTでは、低温で新規な氷構造（“リボン状氷”）が形成されるこ とが見出された。リボン状氷では水分子のプロトンの配向が秩序化しており、CNT の長さ 方向に並ぶ水分子鎖の本数が奇数（偶数）では（反）強誘電的になることも示唆された。
また短幅d ~ 0.9nmのとき、内包水はDに依らず低温まで液体状態であることが示された。
(2) CNT内の水輸送：CNTの両端に水容器をつないだモデルで、チューブ軸方向に圧力差
Pを加えて水を流す古典MD 計算を行った。計算結果から、CNT内の水の体積流量QCNTを
求めた。図4に、ハーゲン・ポアズイユ（H-P）方程式と比べた体積流量の増大率QCNT/QH-P
の直径依存性を示す。P=200MPa では内包水はD に依らず液体状態であり、D が小さいほ どQCNT/QH-Pの値は増大した。一方P < 100MPaのとき、D = 1.15nmと1.24nmでは内包水 がice-NT構造を形成しており、QCNT/QH-Pは減少した。QCNT/QH-Pが減少する原因には、CNT
内部の ice-NT（固相）と外部のバルク水（液相）との相境界で流れが妨げられる、または
液体に比べてice-NTでは空洞壁との間の摩擦が大きい、などの可能性が考えられる。そこ で、CNTと内包水との間の摩擦を直接的に解析可能なモデルを現在試作している。
(3) CNT試料の高圧実験 :アンビルセルを用いて、CNT試料の高圧ラマン散乱実験と粉末
X線回折実験を行った。乾燥したフィルム状試料をダイヤモンドアンビルで一軸加圧した結 果、高圧下でCNT試料のバンドル構造およびCNTの円筒形状が変化することが示唆された。
今後は構造解析手法を確立し、(1)で述べたリボン状氷の実証実験を行いたいと考えている。
[結論]
CNT内包水の構造と相転移挙動は空洞形状に著しく依存し、扁平化したCNT内ではリボ ン状の氷が形成されることが見出された。このリボン状氷の実証に向けて、高圧実験に着 手している。CNT内の水輸送において、内包水がice-NT構造のときは液体構造のときに比 べて流量が減少することが示された。この要因を明らかにするため、水とCNT壁との間の 摩擦を直接的に解析可能な計算モデルを現在試作している。今後、空洞壁と水との間の摩 擦が、水のナノ構造によってどのように変化するかを明らかにしたい。
[参考文献] [1] Y. Maniwa, et al. Chem. Phys. Lett. 401, 534-538 (2005) [2] H. Kyakuno, et al. J.
Chem. Phys. 145, 064514 (2016) [3] J. K. Holt, et al. Science. 312, 1034-1037 (2006) 図1. CNTの変形（上図）と、変形CNT
内包水の計算モデル（下図）。変形の割合 はγ (%)= (D-d)/D×100と定義した。
図 3. 変形 CNT 内包水の温度-変形相図。Fused ice-NTは、2種類のice-NTが融合した構造である。
図4. H-P方程式と比べた体積流量の増大率の直径 依存性。挿入図は使用した計算モデル。
図2. 変形CNT内包水の低温構造（D = 1.24nm）。左：断面図、右：側面図。右 図ではCNTを省略。
0 10 20 30 40 50 60 70
0.8 1.2 1.6 2 2.4 2.8 3.2 QCNT/QH-P
D (nm) T=300K
200MPa 50MPa
20MPa 100MPa
目次
第
1章 序論
1.1 ナノ細孔に閉じ込められた水の構造と相転移 1.2 一次元ナノ細孔内での水輸送特性
1.2 本研究の目的
第
2章 実験原理
2.1 古典分子動力学計算 2.2 ラマン散乱分光 2.3 X線回折
第
3章 実験方法
3.1 変形SWCNT内包水の構造と相転移 3.2 SWCNT内の水輸送
3.3 SWCNT試料の加圧実験
第
4章 結果と考察
4.1 変形SWCNT内包水の構造と相転移 4.1.1 構造と相転移温度
4.1.2 二体相関関数
4.1.3 拡散係数と回転相関時間 4.1.4 リボン状氷の自発分極 4.1.5 SWCNTの変形に必要な圧力
4.1.6 不均一に変形したときの内包水の低温構造
4.2 SWCNT内部の水輸送
4.2.1 流速分布と体積流量の解析：モデル①
4.2.2 水とSWCNT壁との摩擦の解析：モデル②
4.2.3 水とSWCNT壁との摩擦の解析：モデル③
4.3 SWCNT試料の加圧実験 4.3.1 ラマン散乱分光実験 4.3.2 粉末X線回折（XRD）実験
4.3.3 ラマン散乱分光実験とXRD実験の結果比較 4.3.4 XRD計算
第
5章 結論
第
6章 付録
第
7章 参考文献
第
8章 謝辞
1
第 1 章 序論
水は私たちにとって最も身近な物質のひとつであり、あらゆる環境下に存在する。中でも、
ナノサイズの制限空間内に閉じ込められた水は、バルクとは異なる振る舞いを示すことが 知られている。このような水の性質を理解することは、生体膜のチャネル機能解明から水環 境下ナノデバイスの開発まで、幅広い分野において必要とされており、理論・実験共に盛ん な研究が行われている。以下に、その例を紹介する。
1.1 ナノ細孔に閉じ込められた水の構造と相転移
ナノ細孔に閉じ込められた水の研究の 1 つとして、構造や相転移挙動に着目した研究が 盛んに行われている。例として、単層カーボンナノチューブ（SWCNT）を用いた研究があ る。SWCNTはナノメートル単位の直径を有する円筒状物質である。その空洞内にはバルク に匹敵する密度の水を内包することができる。内包された水は、SWCNTの空洞直径Dや温 度Tに依存して様々な振る舞いを示す。D < ~1.5nmでは、液体-固体相転移を伴って低温で 筒状の氷（アイスナノチューブ、Ice-NT）を形成する [1]-[7]。このIce-NTは、最初、原子 レベルの計算機シミュレーションによって予測され[1]、その後、粉末 X 線回折（XRD）実 験によってその存在が実証された[2][3]。これにより、SWCNT内包水が液体からIce-NTに なる際は液体-固体相転移様挙動を伴い、かつその相転移温度が直径Dの増加と共に急激に 減少することが見出されている。このIce-NT形成における相転移温度については、計算機 シミュレーションによる詳細な解析も行われ、実験結果を再現する直径依存性が得られて いる[4][5]。一方、D > ~1.6nmでは内包水と空洞壁との親和性が変化してwet-dry転移が起 こることが、XRD、核磁気共鳴（NMR）実験、電気抵抗測定などを用いて明らかにされた[8]- [10]。客野らは、広範囲にわたる直径で内包水の振る舞いを系統的にまとめたグローバル相 図を提案している（図1.1）[10]。
SWCNT以外にも、2層の疎水性平板による2次元ナノ細孔や、ゼオライト鋳型炭素（ZTC）
による3次元ナノ細孔に内包された水についても研究がなされてきた。2層の疎水性平板で は、平板間距離hに応じて低温で1層から3層の二次元氷を形成し、その際に液体-固体相 転移を伴う（図1.2）[7][11]-[14]。ZTCでは、低温でガラス転移を伴って低密度水（LDL）
様のアモルファス構造を形成することが報告されている（図1.3）[15]。
以上のように、ナノ細孔に内包された水の構造や相転移挙動は、空洞のサイズや次元によ って著しく変化することが明らかにされてきた。これらに加え、ナノ空洞の形状も内包水の 物性に大きく影響を与えると考えられる。たとえばSWCNTは、その弾性ゆえに空洞形状を 系統的に変化させる上で適していると思われる。しかし、SWCNT形状を系統的に変化させ たときの内包水の振る舞いについては我々の知る限り報告されていない。
2
a b
c
図1.1.1 SWCNT内包水の温度（T）-直径（D）相図（文献[10]より引用）。
図1.1.2 2層の疎水性平板間における二次元氷。(a)1層の二次元氷（文献[12]より引用）。 (b)2層の二次元氷（文献[13]より引用）。(c)3層の二次元氷（文献[14]より引用）。
図1.1.3 ZTC内における3次元のアモルファス氷。（文献[15]より引用）。
3 1.2 一次元ナノ細孔内での水輸送特性
ナノ細孔内の水の研究では、構造や相転移挙動だけでなく、輸送特性についても注目され ている。従来のマクロな流体力学では、水が円筒状の管の中を流れるときの体積流量はハー ゲン・ポアズイユ（H-P）方程式を用いて
(1.1) と表される。このときΔPは管の両端の圧力差、rは管の半径、Lは管の長さ、ηは水の粘性 係数である。しかし、CNTのようなナノサイズの管では、このようなマクロな流体力学から 逸脱した振る舞いがみられることが理論・実験共に報告されてきた。
理論的な観点では古典 MD 計算による報告が主である。例として、Hummer らは D = 0.8nm[16]、ThomasらはD = 1.66-4.99nm[17]のSWCNTについてMD計算を行い、得られ た体積流量がH-P方程式から予測される値の102-104倍であることを明らかにした。また、
KalraらはD = 0.8nmのSWCNT両端に浸透圧を与えたとき、内部における水の伝達速度が アクアポリンの空洞内部における輸送速度に匹敵することを見出した[18]。
実験では、HoltらはD = 1.3-2nmの2層カーボンナノチューブ（DWCNT）[19]、Majumder らはD = 7nmの多層カーボンナノチューブ（MWCNT）[20]、更にWhitbyらはD = 44nmの カーボンナノパイプ[21]を用いた研究を行い、体積流量がH-P方程式の予測値の102-104倍 であると報告している。以上のようなCNT 内の高速水輸送についてまとめると、図 1.2.1a のようになる[22]。
これらの報告は、マクロな流体力学における”滑りなしの条件”が CNT のナノ空洞では破 綻することを示唆する結果であると考えられている。すなわち図 1.2.1b に示すように、巨 視的な管では管壁における流速はゼロになり、“滑りなしの条件”に則った振る舞いを示す。
しかし、CNTのようなナノサイズの管では、管壁における流速はゼロにならず、”滑りなし の条件”が破綻すると考えられている。
以上のように、CNT 内ではマクロな流体力学から逸脱した水輸送特性が得られる。その 原因として、CNT 壁が原子レベルでなめらかであるために水との間に摩擦が生じない（も しくは非常に摩擦が小さい）とする説があり、原子レベルでの詳細な議論が求められている。
例えばJosephらはD ~ 1.6nmのSWCNTで古典MD計算を行い[23]、SWCNT内の速い水輸 送は疎水性や原子レベルの滑らかさによるものだとした。また、Falk らは同様の見解の下 で、摩擦の SWCNT 直径依存性を解析し、直径が小さいほど摩擦が小さいことを見出した
[24]。こうした一方で、SWCNT と同じ原子配列をもつ窒化ホウ素ナノチューブでは水輸送
速度は、SWCNT より1ケタ以上小さいという報告[25]もある。以上より、ナノチューブ内 の流動現象については、単純に空洞壁の幾何学的構造や空洞サイズのみで特徴づけること はできないと言える。
更に、CNT 内の水輸送は、内包水のナノ構造にも強く影響を受けると予想される。例え ば、Thomasらは内包水がIce-NT様の構造のときは流れが抑制されると報告している[26]。
4
また、WangらはSWCNT動径方向の流速分布と密度分布が類似した振動構造を描き、その
流速はNavier-Stokes方程式による予測値よりもはるかに速いと主張している[27]。しかし、
内包水のナノ構造と水輸送特性との相関（水クラスターとCNT壁との間の摩擦など）につ いては、十分に議論されているとは言い難いのが現状である。
1.3 本研究の目的
以上の研究背景をふまえ本研究では、(1)内包水の挙動の空洞形状への依存性を系統的に 明らかにすること、(2)水のナノ構造と摩擦に着目し、CNT空洞内での高速水輸送の機構を 明らかにすること、(3)変形CNT内の水の挙動を実験により検証することを目的とした。
研究手法には、古典分子動力学（MD）計算、高圧ラマン散乱分光実験、高圧X線回折実 験を用いた。
a b
図1.2.1 ナノサイズの管における特異な水輸送。(a)H-P方程式と比べた体積流量の増大
率の空洞径依存性（文献[22]より引用）。(b) 巨視的な管とナノサイズの管における動径 方向の速度分布のモデル図。
5
第 2 章 実験原理
文献[28]-[35]を用いて、古典分子動力学計算、ラマン散乱分光、X線回折の原理を述べる。
2.1 古典分子動力学計算
古典分子動力（MD）計算は、原子にニュートンの運動方程式を適用することで、原子・
分子の運動を記述する手法である。系の初期状態（初期構造・速度）を決め、原子間の相互 作用ポテンシャルより得られる力を用いて数値積分を行い、系の時間発展を記述する。
図 2.1-1 に MD 計算の具体的なフローチャートを示す。このフローチャートに記された
個々の計算プロセスについて以下で説明する。
(i) 原子配置
MD計算では、まずユニットセルの中に対象の原子・分子を並べなければならない。単に セル内に原子を並べた場合、それは表面を持ったクラスターとなる（図2.1-2a） 。表面で は原子の再配列などの複雑な現象が起こる場合があるため、表面に注目しないなら周期境界 条件を用いる必要がある。
周期境界条件を適用させた場合、ユニットセルと全く同じセル(イメージセル)が周りに連 続して続いていると想定するため、疑似的に無限のサイズであると考えることができる（図
2.1-2b）。これにより、表面のないバルクでの性質を調べることができる。
周期境界条件を用いて計算をする場合、ユニットセルのサイズに注意する必要がある。各 原子間ポテンシャルにはカットオフ距離 𝑟𝑐 が設定され、その原子からポテンシャル相互作 用が及ぶ範囲が定められている。もしセルサイズがこの範囲より小さいと、イメージセル内 の同じ原子から重複してポテンシャル相互作用を受け、原子の挙動が不自然になってしまう
（図2.1-2c）。それを防ぐため、セルサイズを 𝐿𝑥、𝐿𝑦、𝐿𝑧としたとき、
図2.1-1 MD計算のフローチャート
6 𝑟𝑐<𝐿𝑥
2, 𝑟𝑐<𝐿𝑦
2, 𝑟𝑐<𝐿𝑧
2 (2.1)
となるようにセルサイズを設定しなければならない。
また、速やかに系を平衡状態にするため、できるだけ平衡状態に近い初期状態から計算を 始めるのが望ましい。この初期状態というのは人為的に決められたものであるため、解析の 際は、初期状態の影響が十分緩和されているか否かに注意する必要がある。
(ii) 原子にはたらく力を計算
原子や分子の挙動を記述する際に、電子状態について考えることは不可欠である。しか し、古典MD計算は古典力学に基づくものであり、電子に関する量子力学は考慮されない。
原子・分子間の相互作用を設定する際は、現実の実験値をよく再現するような経験的ポテン シャルが導入される。例として、非結合二体力相互作用を記述する最も基本的な関数である レナード・ジョーンズポテンシャルと、炭素原子間やシリコン原子間の共有結合など、三体 力相互作用を表現するためによく用いられる、Tersoff ポテンシャルを紹介する。
c
図2.1-2 MDセルのイメージ図。簡単のため二次元で表示している。(a) ユニットセル(b) 周期境界条件を課した時のイメージセル。(c) カットオフ距離の設定。
a b
7 (a) レナード・ジョーンズポテンシャル
レナード・ジョーンズポテンシャルは、2原子間のファンデルワールス相互作用を表現し たポテンシャル関数である。原子間距離を 𝑟 としたとき、
𝜙(𝑟) = 4𝜀 {(𝜎 𝑟)
12
− (𝜎 𝑟)
6
} (2.2)
と定義される。ここで、パラメータ𝜀、𝜎は原子の種類によって異なる値となり、右辺第 1 項は斥力、第2項は引力相互作用を表している。これをグラフにすると図2.1-3のようにな る。ポテンシャルが極小値をとる原子間距離を 𝑟0 としたとき、原子間距離が 𝑟0 より小さい ときは斥力が、大きいときは引力が働く。そのため、2 原子は原子間距離 𝑟0 を中心に振動 することになる。
(b) Tersoff ポテンシャル
Tersoff ポテンシャルは、シリコンやダイヤモンドなどの共有結合素材によく用いられて
いる。関数内に角度依存性の項を加えることで、sp2や sp3混成軌道を上手く表現できる。
その一方で、グラフェンの分散を再現できないといった課題も存在した。
そこで、Optimized Tersoffポテンシャルは、Tersoffポテンシャルのパラメータを修正す ることでこの課題を克服した[36]。これにより、上記の分散だけでなく、グラフェンや
SWCNTの熱伝導特性にも最適化された。
粒子 𝑖 と 𝑗 のポテンシャルは、
𝑉𝑖𝑗(𝑟) = 𝑓𝐶(𝑟)[𝑓𝑅(𝑟) + 𝑏𝑖𝑗𝑓𝐴(𝑟)] (2.3) とおける。ここで
𝑟 = |𝒓𝑖− 𝒓𝑗| (2.4)
である。力の到達距離を 𝑅 とすると、カットオフを決める項 𝑓𝐶 は
図2.1-3レナード・ジョーンズポテンシャルにおけるポテンシャルエネルギー（PE）の
原子間距離 r 依存性。
r r0
PE
0 σ
-ε
8 𝑓𝐶(𝑟) =
{
1, 𝑟 < 𝑅 − 𝐷 1
2−1
2𝑠𝑖𝑛 [𝜋(𝑟 − 𝑅) 2𝐷 ] 0, 𝑟 > 𝑅 + 𝐷
, 𝑅 − 𝐷 < 𝑟 < 𝑅 + 𝐷 (2.5)
と書ける。また、斥力項 𝑓𝑅(𝑟) と引力項 𝑓𝐴(𝑟) はそれぞれ
𝑓𝑅(𝑟) = 𝐴𝑖𝑗exp (−𝜆1𝑟) (2.6) 𝑓𝐴(𝑟) = −𝐵𝑖𝑗exp (−𝜆2𝑟) (2.7) と表せる。式(2.3)の右辺第2項の 𝑏𝑖𝑗 は結合角を表す項であり、
𝑏𝑖𝑗 = 𝜒𝑖𝑗(1 + 𝛽𝑖𝑛𝑖𝜁𝑖𝑗𝑛𝑖)−1/2𝑛𝑖 (2.8) である。更に、 𝜁𝑖𝑗 は
𝜁𝑖𝑗 = ∑ 𝑓𝐶(𝑟)𝜔𝑖𝑘𝑔(𝜃𝑖𝑗𝑘)𝑒𝑥𝑝[𝜆33(𝑟𝑖𝑗− 𝑟𝑖𝑘)3]
𝑘≠𝑖,𝑗
(2.9) と書かれ、 𝑒𝑥𝑝 の項は結合長依存を示す。角度依存を表す 𝑔(𝜃𝑖𝑗𝑘) は
𝑔(𝜃𝑖𝑗𝑘) = 1 +𝑐2
𝑑2− 𝑐2
[𝑑2+ (ℎ − 𝑐𝑜𝑠𝜃)2] (2.10) と定義される。Tersoff ポテンシャルがダイヤモンド構造やグラファイト構造で安定になる のは、この角度依存項の効果が大きい。
Optimized Tersoffポテンシャルのパラメータを以下の表2.1-1に示す。
表2.1-1 : Optimized Tersoffポテンシャルのパラメータ
パラメータ 𝐴𝑖𝑗(𝑒𝑉) 𝐵𝑖𝑗(𝑒𝑉) 𝜆1(1/Å) 𝜆2(1/Å) 𝛽 値 1.3936×103 4.3×102 3.4879 2.2119 1.5724×10-7
𝑛 𝑐 𝑑 ℎ 𝜆3(1/Å) 𝑅(Å) 𝐷(Å)
7.2751×10-1 3.8049×104 4.3484 -9.3×10-1 0 1.95 1.5×10-1
𝜒𝑖𝑗 𝜔𝑖𝑘
1 1
(c) ポテンシャルのカットオフ
原子間距離が、ある距離 𝑟𝑐以上離れたときに相互作用を打ち切ることをカットオフとい う。カットオフ距離である𝑟𝑐はポテンシャル関数内に与えられている場合が多いが、与えら れていない場合は自身で設定する必要がある。
カットオフ距離の設定は計算に影響が出ないように十分長い距離で設定しなければなら ない。また、単純に相互作用を打ち切った場合、カットオフ距離 𝑟𝑐 を境にエネルギーが不 連続になってしまう。それを防ぐため、エネルギー値が 𝑟𝑐 でゼロになるように、ポテンシ ャル関数をシフトする必要がある（図2.1-4）。
9
𝜙𝑠= { 𝜙(𝑟) − 𝜙(𝑟𝑐) (𝑟 ≤ 𝑟𝑐)
0 (𝑟 ≥ 𝑟𝑐) (2.11)
𝜙𝑠 をシフトポテンシャルと呼ぶ。しかし、単にシフトするだけではエネルギーの微分値（力）
が不連続になってしまうため、エネルギーの微分値が 𝑟𝑐 でゼロになるようにした、シフト フォースポテンシャルがよく用いられている。
𝜙𝑠𝑓 = { 𝜙(𝑟) − 𝜙(𝑟𝑐) −𝑑𝜙(𝑟𝑐)
𝑑𝑟 (𝑟 − 𝑟𝑐) (𝑟 ≤ 𝑟𝑐) 0 (𝑟 ≥ 𝑟𝑐)
(2.12)
(d) 計算の高速化手法
MD 計算ではカットオフ距離 𝑟𝑐 内でのみ相互作用の計算を行うため、どの原子が 𝑟𝑐 内に あるかの台帳が必要になる。しかし、1ステップの計算では台帳はほとんど変化しないため、
毎ステップごとに台帳を更新するのは明らかに無駄である。そこで、カットオフ距離 𝑟𝑐 よ り少し大きい範囲内 (𝑟𝑐+ 𝑑𝑟) に存在する原子の番号を台帳に記録しておき、台帳の外の原
子が𝑟𝑐 内に入ってくる可能性がないときには台帳の更新を行わないことにより、計算量を大
幅に削減することができる。これをブックキーピング法と呼ぶ。
(iii) 差分法による原子の速度・位置の計算
MD 計算では、多数の粒子から成る力学系の運動方程式を差分法で近似して数値積分し、
粒子の速度や位置を時間の関数として求める。この数値積分では、できるだけ精度よく、且 つ短い計算時間で積分を実行することが求められる。ここでは、分子動力学の数値積分によ く用いられる予測子 - 修正子法と、これを改良したGear法について述べる。
(a) 予測子-修正子法
図 2.1-4 ポテンシャル関数シフトの模式図。点線がシフト前で、実線がシフト後の関数
である。カットオフ距離rc でポテンシャルの値がゼロになるようにする。
10
運動方程式は、座標𝒓𝑖と速度𝒗𝑖について以下のような1階の連立微分方程式に分けられる。
𝑑𝒓𝑖 𝑑𝑡 = 𝒗𝑖 𝑑𝒗𝑖
𝑑𝑡 = 𝑭𝑖 𝑚𝑖
𝑖 = 1, 2, ⋯ , 𝑁. (2.13)
これは6N 個の連立微分方程式であるが、取り扱いは各成分で同じなので、従属変数𝑦(𝑡)に ついての1階常微分方程式
𝑑𝑦
𝑑𝑡 = 𝑓(𝑦, 𝑡) (2.14)
を考える。ここで、時間刻み幅をℎとし、𝑡𝑛 = 𝑛ℎ、 𝑦(𝑡𝑛) = 𝑦𝑛、𝑓𝑛 = 𝑓(𝑦𝑛, 𝑡𝑛)とおく。
まず、(𝑟 + 1)個の時刻𝑡𝑛, 𝑡𝑛−1, ⋯ , 𝑡𝑛−𝑟において𝑦(𝑡)が𝑦𝑛, 𝑦𝑛−1, ⋯ , 𝑦𝑛−𝑟であるとしたときの、
時刻𝑡𝑛+1における𝑦(𝑡)の近似解を求める。式(2.14)を𝑡𝑛から𝑡𝑛+1まで積分すると、
𝑦𝑛+1= 𝑦𝑛+ ∫ 𝑓(𝑦(𝑡), 𝑡)𝑑𝑡
𝑡𝑛+1 𝑡𝑛
(2.15) である。ここでの被積分関数は未知関数𝑦(𝑡)を含むので、これだけでは𝑦𝑛+1を求めることは できない。そこで、被積分関数を時刻𝑡𝑛−𝑘において値𝑓𝑛−𝑘≡ 𝑓(𝑦𝑛−𝑘, 𝑡𝑛−𝑘), 𝑘 = 0, 1, 2, ⋯ 𝑟を とる𝑡について𝑟次の多項式で近似する。時間刻み幅ℎが等間隔で与えられているとき、被積 分関数𝑓(𝑦(𝑡), 𝑡)はニュートンの後退補間公式を用いて、
𝑓(𝑦(𝑡), 𝑡) = 𝑓𝑛+(𝑡 − 𝑡𝑛)
ℎ ∇𝑓𝑛+(𝑡 − 𝑡𝑛)(𝑡 − 𝑡𝑛−1)
2! ℎ2 ∇2𝑓𝑛+ ⋯ +(𝑡 − 𝑡𝑛)(𝑡 − 𝑡𝑛−1)(𝑡 − 𝑡𝑛−2) ⋯ (𝑡 − 𝑡1)
𝑛! ℎ𝑛 ∇𝑛𝑓𝑛
= 𝑓𝑛+(𝑡 − 𝑡𝑛)
ℎ ∇𝑓𝑛+(𝑡 − 𝑡𝑛)(𝑡 − 𝑡𝑛+ ℎ)
2! ℎ2 ∇2𝑓𝑛+ ⋯
(2.16)
+(𝑡 − 𝑡𝑛)(𝑡 − 𝑡𝑛+ ℎ)(𝑡 − 𝑡𝑛+ 2ℎ) ⋯ (𝑡 − 𝑡𝑛+ (𝑛 − 1)ℎ)
𝑛! ℎ𝑛 ∇𝑛𝑓𝑛 (2.17)
と表される。これを式(2.15)に代入して積分を実行すると、
𝑦𝑛+1𝑝 = 𝑦𝑛+ ℎ (𝑓𝑛+1
2∇𝑓𝑛+ 5
12∇2𝑓𝑛+3
8∇3𝑓𝑛+251
720∇4𝑓𝑛+ ⋯ ) (2.18) となる。右辺は𝑟次の多項式の場合、∇𝑟の項で終わる。∇は後退差分演算子と呼ばれ、
∇𝑓𝑛 = 𝑓0,
∇𝑓𝑛 = 𝑓𝑛− 𝑓𝑛−1,
∇2𝑓𝑛= ∇(∇𝑓𝑛) = 𝑓𝑛− 2𝑓𝑛−1+ 𝑓𝑛−2,
⋮
(2.19)
11
∇𝑘𝑓𝑛= ∑(−1)𝑖(𝑘 𝑖)
𝑘
𝑖=0
𝑓𝑛−𝑖
である。更に、式(2.18)を𝑓𝑛, 𝑓𝑛−1, ⋯ , 𝑓𝑛−𝑟を用いた線形の式で表すと、
𝑦𝑛+1𝑝 = 𝑦𝑛+ ℎ(𝛽𝑟,0∗ 𝑓𝑛+ 𝛽𝑟,1∗ ∇𝑓𝑛−1+ 𝛽𝑟,2∗ ∇2𝑓𝑛−2+ ⋯ + 𝛽𝑟,𝑟∗ 𝑓𝑛−𝑟) (2.20) となる。この式はアダムス-バシュフォース公式(AB公式)と呼ばれ、右辺の既知量から左辺 の未知量(“予測子”)を与えるための式である。尚、𝑟 = 1~5における𝛽𝑟,𝑖∗ の値は以下の表のよ うになる。
表2.1-2 AB公式 (2.20) の係数 𝛽𝑟,𝑖∗ 𝑖
𝑟 0 1 2 3 4 5
𝑓𝑛 𝑓𝑛−1 𝑓𝑛−2 𝑓𝑛−3 𝑓𝑛−4 𝑓𝑛−5
1 2𝛽1,𝑖∗ 3 -1
2 12𝛽2,𝑖∗ 23 -16 5
3 24𝛽3,𝑖∗ 55 -59 37 -9
4 720𝛽4,𝑖∗ 1901 -2774 2616 -1274 251
5 1440𝛽5,𝑖∗ 4277 -7923 9982 -7298 2877 -475 以上のように予測子を求めたら、次は修正子を求める。そのためには、既知量𝑓𝑛, 𝑓𝑛−1, ⋯ , 𝑓𝑛−𝑟 に加えて、未知量𝑓𝑛+1= 𝑓(𝑦(𝑡𝑛+1), 𝑡𝑛+1)も用いて式(2.15)の被積分関数を(𝑟 + 1)次の多項式 で近似する必要がある。これを積分すると、
𝑦𝑛+1𝑐 = 𝑦𝑛+ ℎ (𝑓𝑛+1−1
2∇𝑓𝑛+1− 1
12∇2𝑓𝑛+1− 1
24∇3𝑓𝑛+1
− 19
720∇4𝑓𝑛+1− ⋯ )
(2.21)
となり、更に後退差分演算を実行すると、
𝑦𝑛+1𝑐 = 𝑦𝑛+ ℎ𝛽𝑟+1,−1𝑓𝑛+1+ ℎ(𝛽𝑟+1,0𝑓𝑛+ 𝛽𝑟+1,1𝑓𝑛−1+ 𝛽𝑟+1,2𝑓𝑛−2
+ ⋯ + 𝛽𝑟+1,𝑟𝑓𝑛−𝑟) (2.22) という形が得られる。尚、𝑟 + 1 = 1~5における𝛽𝑟+1,𝑖, 𝑖 = −1, 0. 1, 2, ⋯ , 𝑟の値を以下の表に 示す。
12 表2.1-3 AM公式 (2.21) の係数 𝛽𝑟,𝑖
𝑖
𝑟 + 1 -1 0 1 2 3 4
𝑓𝑛+1 𝑓𝑛 𝑓𝑛−1 𝑓𝑛−2 𝑓𝑛−3 𝑓𝑛−4
1 2𝛽1,𝑖∗ 1 1
2 12𝛽2,𝑖∗ 5 8 -1
3 24𝛽3,𝑖∗ 9 19 -5 1
4 720𝛽4,𝑖∗ 251 646 -264 106 -19
5 1440𝛽5,𝑖∗ 475 1427 -798 482 -173 27
この式は、アダムス-モールトンの公式(AM公式)と呼ばれ、右辺第2項の𝑓𝑛+1に未知量𝑦𝑛+1 が含まれている。そこに、まずはAB公式の𝑦𝑛+1𝑝 (≡ 𝑦𝑛+1[0] とおく)を初期値として代入し、こ れを修正した𝑦𝑛+1𝑐 を次の近似解𝑦𝑛+1[1] とする。これを再び𝑓𝑛+1に代入し次の近似解𝑦𝑛+1[2] を求め る。この反復を繰り返すことで、数値解の誤差を少なくしていくのが予測子-修正子法であ る。
ここで更に、予測子-修正子法の反復の構造を𝑟 = 2の予測子と𝑟 + 1 = 3の修正子を例にと って説明する。表2.1-2, 2.1-3より、予測子と修正子は
𝑦𝑛+1[0] = 𝑦𝑛+ ℎ
12(23𝑓𝑛− 16𝑓𝑛−1+ 5𝑓𝑛−2) (2.23) 𝑦𝑛+1[𝑚+1]= 𝑦𝑛+3
8ℎ𝑓𝑛+1[𝑚]+ ℎ
24(19𝑓𝑛− 5𝑓𝑛−1+ 𝑓𝑛−2), 𝑚 ≥ 0 (2.24) である。ここで、反復回数を上付き括弧[m]で表している。まず、𝑚 = 0としたときの式(2.23) と式(2.24)の差をとると、
𝑦𝑛+1[1] − 𝑦𝑛+1[0] =3
8ℎ𝑓𝑛+1[0] −3
8ℎ(3𝑓𝑛− 3𝑓𝑛−1+ 𝑓𝑛−2) (2.25) となる。ここで、
𝑓𝑛+1[−1]= 3𝑓𝑛− 3𝑓𝑛−1+ 𝑓𝑛−2 (2.26) とおくと、
𝑦𝑛+1[1] − 𝑦𝑛+1[0] =3
8ℎ(𝑓𝑛+1[0] − 𝑓𝑛+1[−1]) (2.27) が得られる。同様に、式(2.24)から
𝑦𝑛+1[𝑚+1]= 𝑦𝑛+1[𝑚] +3
8ℎ(𝑓𝑛+1[𝑚]− 𝑓𝑛+1[𝑚−1]) (2.28) という形も得られる。これが修正子の反復公式であり、𝑚 = 0の場合も含んでいる。右辺第 2項は反復に伴う修正を与え、𝑚 → ∞で𝑦𝑛+1[𝑚+1]が解に収束するならばこの修正項は0になる。
13
そして、|𝑓𝑛+1[𝑚]− 𝑓𝑛+1[𝑚−1]|の大きさは反復回数の判定に使うことができる。
(b) 予測子-修正子法の行列表示
予測子-修正子法は行列の形で表すと見通しが良くなる。まず、予測子は
( 𝑦𝑛+1[0]
ℎ𝑓𝑛+1[−1]
ℎ𝑓𝑛[−1]
ℎ𝑓𝑛−1[−1])
= (
1 23/12 −16/12 5/12
0 3 −3 1
0 1 0 0
0 0 1 0
) ( 𝑦𝑛
ℎ𝑓𝑛 ℎ𝑓𝑛−1
ℎ𝑓𝑛−2
) (2,29)
とかける。右辺の4 × 4の正方行列を𝐁とおく。このベクトル式の第1, 第2成分はそれぞれ 式(2.23)と式(2.26)を表しており、第3, 第4成分はそれぞれ𝑓𝑛[−1]= 𝑓𝑛, 𝑓𝑛−1[−1]= 𝑓𝑛−1である。
行列𝐁は差分式の形で決まるものであり、時間刻み幅ℎには関係しない。修正子の反復公式 は式(2.28)より
( 𝑦𝑛+1[𝑚+1]
ℎ𝑓𝑛+1[𝑚]
ℎ𝑓𝑛[𝑚]
ℎ𝑓𝑛−1[𝑚])
= (
𝑦𝑛+1[𝑚]
ℎ𝑓𝑛+1[𝑚−1]
ℎ𝑓𝑛[𝑚−1]
ℎ𝑓𝑛−1[𝑚−1])
+ ℎ (𝑓𝑛+1[𝑚]− 𝑓𝑛+1[𝑚−1]) ( 3/8
1 0 0
) , 𝑚 ≥ 0 (2.30)
で表される。第1成分は式(2.28)であり、第2成分は𝑓𝑛+1[𝑚−1]を𝑓𝑛+1[𝑚]に修正するためのもので ある。
これらについて、予測子が𝑟次、修正子が(𝑟 + 1)次の多項式で与えられる場合の公式でま とめてみる。このとき、(𝑟 + 2)次元のベクトル𝒀𝑛, 𝒀𝑛+1[0] , 𝒀𝑛+1[𝑚]および補正ベクトル𝑪と、(𝑟 + 2) × (𝑟 + 2)の正方行列𝑩を以下のように導入する。
𝒀𝑛= (
𝑦𝑛
ℎ𝑓𝑛 ℎ𝑓𝑛−1
⋮
⋮ ℎ𝑓𝑛−𝑟)
, 𝒀𝑛+1[0] =
( 𝑦𝑛+1[0]
ℎ𝑓𝑛+1[−1]
ℎ𝑓𝑛[−1]
⋮
⋮ ℎ𝑓𝑛−(𝑟−1)[−1] )
, 𝒀𝑛+1[𝑚] =
( 𝑦𝑛+1[𝑚]
ℎ𝑓𝑛+1[𝑚−1]
ℎ𝑓𝑛[𝑚−1]
⋮
⋮ ℎ𝑓𝑛−(𝑟−1)[𝑚−1] )
(2.31)
𝑪 = (
𝛽𝑟,−1
1 0
⋮
⋮ 0 )
, 𝑩 = (
1 𝛽𝑟,0∗ 𝛽𝑟,1∗ 𝛽𝑟,2∗ ⋯ 𝛽𝑟,𝑟−1∗ 𝛽𝑟,𝑟∗ 0 𝛿𝑟,0∗ 𝛿𝑟,1∗ 𝛿𝑟,2∗ ⋯ 𝛿𝑟,𝑟−1∗ 𝛿𝑟,𝑟∗
0 1 0 0 ⋯ 0 0
⋮ 0 1 0 ⋯ 0 0
⋮ ⋮ ⋮ ⋮ ⋯ ⋮ ⋮
0 0 0 0 ⋯ 1 0 )
(2.32)
ここで、
𝛿𝑟,𝑘∗ = 1
𝛽𝑟,−1(𝛽𝑟,𝑘∗ − 𝛽𝑟,𝑘), 𝑘 = 0, 1, ⋯ , 𝑟 (2.33) であり、
14
𝑓𝑛+1[−1]= ∑ 𝛿𝑟,𝑘∗ 𝑓𝑛−𝑘
𝑟
𝑘=0
(2.34) で表せる。よって、初期値ベクトル𝒀𝑛を与えると予測子は
𝒀𝑛+1[0] = 𝑩𝒀𝑛 (2.35)
となり、修正子は
𝒀𝑛+1[𝑚+1] = 𝒀𝑛+1[𝑚] + ℎ (𝑓𝑛+1[𝑚]− 𝑓𝑛+1[𝑚−1]) 𝑪 (2.36) と与えられる。
(c) Gear法
前述の予測子-修正子法は、過去の複数のステップから次のステップの値を求める手法で あるため、多段階法と呼ばれる。それに対し、Gear法は1つの時刻における高階微分から 計算を行う。その際、𝑦𝑛, ℎ𝑦𝑛′, (ℎ2/2!)𝑦𝑛′′, (ℎ3/3!)𝑦𝑛′′′, ⋯ のうち、必要なだけの値を指定する ため、多値法と呼ばれる。多値法は多段階法と違い、1つの時刻における量だけで数値積分 が出発できる、時間刻み幅ℎを積分の途中で自由に変更できるといった利点がある。しかし、
力の高階微分を各時刻で求めなければならないため、計算の負荷が増すというデメリットも 存在する。
多値法のベクトルを𝒀𝑛𝑁とかき次のように表すこととする。
𝒀𝑛𝑁=
( 𝑦𝑛 ℎ𝑦𝑛′ (ℎ2/2!)𝑦𝑛′′
(ℎ3/3!)𝑦𝑛′′′
(ℎ4/4!)𝑦𝑛(4)
⋯
⋯ )
(2.37)
このとき、𝑦𝑛+1, 𝑦𝑛+1′ , 𝑦𝑛+1′′ , ⋯ について𝑡𝑛 まわりのテイラー展開を実行すると、
𝒀𝑛+1𝑁 = (
𝑦𝑛+1 ℎ𝑦𝑛+1′ (ℎ2/2!)𝑦𝑛+1′′
(ℎ3/3!)𝑦𝑛+1′′′
(ℎ4/4!)𝑦𝑛+1(4)
⋯ )
=
(
𝑦𝑛+ ℎ𝑦𝑛′+ (ℎ2/2!)𝑦𝑛′′+ (ℎ3/3!)𝑦𝑛′′′+ (ℎ4/4!)𝑦𝑛(4)+ ⋯ ℎ𝑦𝑛′ + 2(ℎ2/2!)𝑦𝑛′′+ 3(ℎ3/3!)𝑦𝑛′′′+ 4(ℎ4/4!)𝑦𝑛(4)+ ⋯
(ℎ2/2!)𝑦𝑛′′+ 3(ℎ3/3!)𝑦𝑛′′′+ 6(ℎ4/4!)𝑦𝑛(4)+ ⋯ (ℎ3/3!)𝑦𝑛′′′+ 4(ℎ4/4!)𝑦𝑛(4)+ ⋯
(ℎ4/4!)𝑦𝑛+1(4) + ⋯
⋯ )
= (
1 1 1 1 1 ⋯
0 1 2 3 4 ⋯
0 0 1 3 6 ⋯
0 0 0 1 4 ⋯
0 0 0 0 1 ⋯
⋯ ⋯ ⋯ ⋯ ⋯ ⋯)(
𝑦𝑛 ℎ𝑦𝑛′ (ℎ2/2!)𝑦𝑛′′
(ℎ3/3!)𝑦𝑛′′′
(ℎ4/4!)𝑦𝑛(4)
⋯ )
(2.38)
となる。更に
15 𝐀 =
(
1 1 1 1 1 ⋯
0 1 2 3 4 ⋯
0 0 1 3 6 ⋯
0 0 0 1 4 ⋯
0 0 0 0 1 ⋯
⋯ ⋯ ⋯ ⋯ ⋯ ⋯)
(2.39)
とおくと、
𝒀𝑛+1𝑁 = 𝑨𝒀𝑛𝑁 (2.40)
と表すことができる。また、𝑓𝑛, 𝑓𝑛−1, 𝑓𝑛−2⋯ を𝑡𝑛+1 まわりで展開すると
𝒀𝑛+1 = (
𝑦𝑛+1 ℎ𝑓𝑛+1
ℎ𝑓𝑛
ℎ𝑓𝑛−1 ℎ𝑓𝑛−2
⋯ )
=
(
𝑦𝑛+1 ℎ𝑦𝑛+1′
ℎ𝑦𝑛+1′ − 2(ℎ2/2!)𝑦𝑛+1′′ + 3(ℎ3/3!)𝑦𝑛+1′′′ − 4(ℎ4/4!)𝑦𝑛+1(4) + ⋯ ℎ𝑦𝑛+1′ − 4(ℎ2/2!)𝑦𝑛+1′′ + 12(ℎ3/3!)𝑦𝑛+1′′′ − 32(ℎ4/4!)𝑦𝑛+1(4) ⋯ ℎ𝑦𝑛+1′ − 6(ℎ2/2!)𝑦𝑛+1′′ + 27(ℎ3/3!)𝑦𝑛+1′′′ − 108(ℎ4/4!)𝑦𝑛+1(4) ⋯
⋯ )
= (
1 0 0 0 0 ⋯
0 1 0 0 0 ⋯
0 1 −2 3 −4 ⋯
0 1 −4 12 −32 ⋯
0 1 −6 27 −108 ⋯
⋯ ⋯ ⋯ ⋯ ⋯ ⋯)(
𝑦𝑛+1 ℎ𝑦𝑛+1′ (ℎ2/2!)𝑦𝑛+1′′
(ℎ3/3!)𝑦𝑛+1′′′
(ℎ4/4!)𝑦𝑛+1(4)
⋯ )
(2.41)
となる。これにより、多段階法のベクトルが多値法のベクトルに変換されている。ここで、
𝐓 = (
1 0 0 0 0 ⋯
0 1 0 0 0 ⋯
0 1 −2 3 −4 ⋯
0 1 −4 12 −32 ⋯
0 1 −6 27 −108 ⋯
⋯ ⋯ ⋯ ⋯ ⋯ ⋯)
(2.42)
とおけば、
𝒀𝑛+1= 𝐓𝒀𝑛+1𝑁 と表せる。これを式(2.36)に代入すると、
𝒀𝑛+1𝑁[𝑚+1]= 𝒀𝑛+1𝑁[𝑚]+ ℎ (𝑓𝑛+1[𝑚]− 𝑓𝑛+1[𝑚−1]) 𝑪𝑁 (2.43) が得られる。ここで、補正ベクトル 𝑪𝑁は
𝑪𝑁= 𝑻−1𝑪 (2.44)
である。以上より、Gear法の予測子・修正子は
𝒀𝑛+1𝑁[0]= 𝐴𝑌𝑛𝑁 (2.45)
𝒀𝑛+1𝑁[𝑚+1]= 𝒀𝑛+1𝑁[𝑚]+ ℎ (𝑓𝑛+1[𝑚]− 𝑓𝑛+1[𝑚−1]) 𝑪𝑁 (2.46) となる。
Gear法と共によく用いられる数値積分法として、速度ベルレ法がある。この方法では、N ステップと N+1 ステップの中間で部分速度 𝑣𝑖𝑁+0.5 を求めることで精度を向上させており、
16
Gear法に比べて計算負荷が低い。しかし、時間刻みが小さい場合は、Gear法の方が高精度 であることが分かっている。
(iv) 物性値の制御
差分法で原子の動きを計算した後、設定した物性値（温度、圧力など）に即した振舞いに なるように系を制御する。例えば、温度制御においては仮想の熱浴との熱のやり取りを考え た能勢法が有名だが、最もシンプルな方法として、速度スケーリング方がよく用いられる。
以下に、速度スケーリング法の詳細を述べる。
速度スケーリング法とは、系全体の運動エネルギーが設定温度 𝑇̃ での熱エネルギーと等し くなるよう、各粒子の速度をスケールする手法である。をスケールする手法である。スケー ル前の粒子の速度を 𝑣𝑖 、スケール後の速度を 𝑣̃ 、スケーリングパラメータを𝛼とすると、 𝑖
𝑣̃ = 𝛼𝑣𝑖 𝑖 (2.47)
である。粒子𝑖の質量を 𝑚𝑖 、粒子数をNとしたとき、𝛼は
∑1 2𝑚𝑖𝑣̃𝑖2
𝑖
= 𝛼2∑1 2𝑚𝑖𝑣𝑖2
𝑖
=𝑔
2𝑁𝑘𝐵𝑇̃ (2.48)
となる。ここで 𝑔 は系の自由度、𝑘𝐵 はボルツマン定数である。これより、
α = √∑1 2𝑚𝑖𝑣̃𝑖2
𝑖
∑1 2𝑚𝑖𝑣𝑖2
𝑖
⁄ = √𝑔
2𝑁𝑘𝐵𝑇̃ ∑1 2𝑚𝑖𝑣𝑖2
𝑖
⁄ (2.49)
が得られる。
このようにして、系の物性を制御する。他にも、圧力を制御する方法としてアンデルセン 法、応力を制御する方法としてパリネロ‐ラーマン法がよく用いられている。
(v) 結果の分析
計算が終わったら結果を使って二次解析を行う。以下にその例を示す。
(a) 二体相関関数と積算配位数
二体相関関数とは、ある原子から距離 𝑟 離れた場所にどれだけの原子が存在するかの指 標であり、構造の特徴を詳細に解析する際に用いられる。MD 計算では、半径 𝑟−𝛿𝑟/2 と 𝑟+𝛿𝑟/2 の2 つの曲面に挟まれた粒子数 𝑛k(𝑟−𝛿𝑟/2,𝑟+𝛿𝑟/2) を各原子について求めて平均 することにより、二体相関関数を求めることができる。すなわち、
𝑔(𝑟) = 𝑉 𝑁(𝑁 − 1)
1
4𝜋𝑟2𝛿𝑟∑ 𝑛𝑘
𝑁
𝑘=1
(𝑟 − 𝛿𝑟/2, r + 𝛿𝑟/2) (2.50) と表せる。ここで 𝑉 はユニットセルの体積、𝑁 は原子数である。SCIGRESS MEでは上式の ように、ユニットセル内における原子の平均密度 (N-1) / V を単位として二体相関関数が出 力される。しかし本研究ではユニットセルに配置されたSWCNT 内部の水に着目している
17
ため、(N-1) / V を単位とするのは無意味である。そこで本研究では、SCIGRESS MEの出力
データに (N-1) / V をかけることで、 𝑟±𝛿𝑟/2の曲面内における原子の平均密度として二 体相関関数を求めた。
積算配位数は、注目した原子から距離 𝑟 以内にある他の原子の数として 𝑁(𝑟) =𝑁
𝑉∫ 4𝜋𝑟2𝑔(𝑟)𝑑𝑟
𝑟 0
(2.51) と計算することができる。本研究では、水分子の平均配位数を求めるために積算配位数を利 用した。すなわち、注目した水分子に対して距離0.32nm以内に存在する水分子を最隣接分 子とみなし、𝑟=0.32nmでの積算配位数を水分子の平均配位数とした。
(b) 平均二乗変位と拡散係数
平均二乗変位（MSD）は、ある初期状態を基準に時間 t 後の変位の差の二乗を平均した ものである。3次元系におけるMSDは
〈|𝑟(𝑡) − 𝑟(0)|2〉 (2.52)
と表され、拡散係数D は
𝐷 =〈|𝑟(𝑡) − 𝑟(0)|2〉
6𝑡 (2.53)
と計算できる。一次元系においては、位置z についての拡散係数を 𝐷𝑧=〈|𝑧(𝑡) − 𝑧(0)|2〉
2𝑡 (2.54)
と計算すればよい。この拡散係数は液体から固体になると値が急激に小さくなることが知ら れており、相変化の指標としてよく用いられる。
本研究では、変形SWCNT内の水の拡散係数を解析し、内包水の相変化を解析した。詳細 は4.1.3節で述べる。
(c) 回転相関関数と回転相関時間
回転相関関数（RCF）は、原子や分子の回転軸の回転運動を示した関数であり様々な定義 がある。その中で、ラマン散乱測定やNMR測定との関連付けが可能なRCFは次式で表され る。
𝐶2(t) =1
2〈3{𝒖(𝑡) ∙ 𝒖(0)}2− 1〉 =1
2〈3cos2𝜃(𝑡) − 1〉 (2.55) ここで、𝒖(𝑡)は時間 t における回転軸のベクトルであり、𝜃(𝑡) は時間 t における初期状態 からの回転軸の傾きである。また、𝐶2(t)は
𝐶2(t) ∝ exp{−(𝑡/𝜏)𝛽} (2.56) とも表され、これにより回転相関時間 𝜏 を求めることができる。この 𝜏 はNMR測定からも 求めることができ、MD計算との相互比較が可能である。パラメータ 𝛽 は1以下の値であり、
18
これが1に近いほど 𝜏 の分布が小さく、動的に均一であると考えられている。液体から固体 になると回転運動は妨げられ、𝜏 の値は急激に大きくなるため、相変化の指標とすることが できる。
本研究では、変形SWCNT内の水の回転相関時間を求め、内包水の相変化を解析した。ま た、SWCNT の変形に伴う内包水分子の回転方向の異方性を見出した。詳細は 4.1.3 節で述 べる。
2.2 ラマン散乱分光
光が物質に入射して分子と衝突すると、その一部は散乱される。この散乱光の波長を調べ ると、大部分の成分は入射光と同じ波長（レイリー散乱光）だが、極わずかな成分として、
入射光と異なった波長の光が含まれている。この入射光と異なった波長をもつ光の振動数は、
分子の固有振動数になっており、この入射光と異なった波長をもつ光をラマン散乱光と呼ぶ。
ラマン散乱光の性質を調べることにより、物質の分子構造や結晶構造などを知ることができ る。入射光の波長は、理論上ではどんな波長でも構わないが、ラマン散乱光の強度は、レイ リー散乱光の強度に比べて極めて微弱なため、レーザー光のような高強度光源を用いる必要 がある。ラマン散乱光が検出されるとき、レイリー散乱光より短波長側に検出されるラマン 散乱光をアンチストークス線、長波長側に検出されるものをストークス線と呼ぶ（図2.2.1）。 一般的にはより強度の大きいストークス線が解析に用いられる。分光された各波長の情報は 波数（1/波長）に換算し、入射光の波数との差を計算する。この波数の差を横軸に、強度 を縦軸にとったものがラマンスペクトルである。
ここで、SWCNTのラマンスペクトルを図2.2.2に示す。SWCNTではその格子振動に由来 する特徴的なスペクトルが現れ、それぞれRBM（radial breathing mode）、D-band、G-band と呼ばれる。これらのスペクトルについて詳細に述べる。
(i)RBM（100 cm-1から350cm-1）
RBM はナノチューブ特有のラマンスペクトルであり、直径が振動することに由来する。
断面が円の形を保って直径が振動すると、sp2結合のC-Cボンド長が変化する振動になり、
RBMが観測される。RBMの振動数は直径D (nm)に反比例し、248/D (cm-1)で表すことがで きる。RBMはナノチューブだけで観測されるスペクトルで、低振動数側に現れる。RBMは
図2.2.1 ラマン散乱分光のモデル図。
19
ナノチューブ合成後の試料評価などによく用いられている。
(ii)D-band（1300cm-1付近）
D-band は SWCNTやグラフェンの結晶格子における欠陥に由来するスペクトルである。
そのため、D-band は試料の品質を評価する場合によく用いられる。特に、後述の G-band との強度比の大小によって、品質の良し悪しが議論される場合が多い。
(iii)G-band（1590cm-1付近）
SWCNT やグラフェンは sp2結合で結晶格子を作る。この sp2結合における格子振動によ
ってC-Cボンド長が変化するスペクトルが1600cm-1近傍で現れる。これをG-bandと呼ぶ。
SWCNTのG-bandはG+とG-の2つのスペクトルに分裂する。これはSWCNTが円筒面を有 することに由来する。G+バンドは円筒軸方向の振動に対応するため、直径に依らず一定値
（~1590cm-1）を持つのに対し、G-バンドは円筒の赤道方向の振動に対応するため、直径が 小さくなるとスペクトルの位置が低いところに現れる。ここで、炭素以外の他の元素（Si など）の場合、炭素より原子量が大きいので、ラマンスペクトルはG-bandよりはるかに低 いところに現れる。故に、G-bandを観測できれば、sp2結合をもった炭素材料であると判断 できる。
図 2.2.2 SWCNT 試料のラマンスペクトル。アーク放電法によって作製された平均直径
1.46nmのSWCNT試料（名城ナノカーボン社製ArcSo）で測定されたデータの一例であ
る。(a)全体図。(b)RBMの拡大図。(c)D-bandの拡大図。(d)G-bandの拡大図。
a b
c
d
0 500 1000 1500 2000
Intensity
Raman Shift(cm-1)
100 120 140 160 180 200 220 240 260
Intensity
Raman Shift(cm-1) RBM
1200 1250 1300 1350 1400 1450 1500 Intensity D-band
Raman Shift(cm-1)
1500 1550 1600 1650 1700 Intensity G-band
Raman Shift(cm-1) G+ G-
![図 2.3.2 デバイの式を用いて計算された XRD 回折パターン（文献[34]より引用） 。上図](https://thumb-ap.123doks.com/thumbv2/123deta/10114922.1948921/28.892.133.393.310.729/図232デバイの式を用いて計算さ回折パターン文献より引用上図.webp)
