平成 1３年度

卒業論文

液晶相転移の分子動力学シミュレーション

高知工科大学工学部
知能機械システム工学科
知能流体力学研究室

高橋

裕一郎

目次
第１章

第２章

第３章

第４章

緒言
1.1 はじめに1

1.2 分子動力学

1

1.3 液晶3

計算方法
2.1 verlet 法

6

2.2 周期境界条件

8

2.3 Lennard-Jones potential

10

2.4 Gay-Berne potential

12

計算結果
分子数 100 個，温度 5k，t*=0〜t*=2，t*=12 の結果

14

分子数 100 個，温度 100k，t*=0〜t*=1，t*=12 の結果

21

分子数 100 個，温度 90k〜1k，t*=1 の結果

25

結言
考察

30

参考文献

31

謝辞

32

第１章

1.1

緒言

はじめに

液晶にはさまざまな種類があり，その挙動も多種多用である．基本的な液晶分子の形は棒状や，円盤状をしている．その特異な分子形状が，重心位置の3次元性を失わせ分子の配向を残し，液体・固体とは違う4つ目の状態としての挙動をしめす．
温度変化で液晶状態を形成するとき，それをサーモトロピック液晶（温度転移型液晶）とよび，溶媒に溶解させたとき形成される多成分系液晶をリオトロピック液晶（濃度転移型液晶）と呼ぶ．また，高分子液晶は主に主鎖型高分子液晶を用いた高強度・高弾性率繊維(ケ
ブラー繊維など)や精密機械部品などに利用されており，低分子液晶は，液晶性に直接関係する流動性や外部刺激に対する高い応答性などの機能を利用してディスプレイ素子などに利用されている．
液晶についての従来の研究は液晶材料を連続体的に扱っておりそれによって出される結果は限界に来ている．そこで本研究では利用価値の高い液晶材料を連続体的ではなく多体系として扱い，
分子レベルでの液晶の挙動を把握し連続体としてでは見出せなかった現象について検証し明らかにしようというものである．

1.2

分子動力学法

分子動力学法は近年多くの分野でさかんに利用されている．それは最近の計測および実験
技術と計算機の進歩により，物理現象を分子レベルで理解してその利用技術のより高度な発展を考えたり，ナノテクノロジーにおける必然的な分子原子レベルでの現象の制御等，分子
原子レベルでの現象を観察理解する必要性が急増してきているからである．これらの現象を実験的に観察し理解するには，非常に高価な実験装置および高度な実験技術が必要であり，容易でない場合が多い．これに対して計算機を用いて分子原子の状態・運動を分子動力学的に観察する場合は，その計算条件さえはっきりしていれば容易にかつ詳細な情報を得ることが出来るようになる．
しかし，分子動力学による物理化学現象の解析において，その現象が分子動力学的に計算可能であり，かつ計算の前提が明確であるかどうかという問題がある．すなわち，その現象の解析に分子動力学が用いられ得るかということである．
分子動力学が成り立つかどうかは分子動力学の成立条件を慎重に検討する以外に方法が無く，結果から帰納法的に判断することは難しい．
また，分子動力学により計算された分子，原子レベルの現象とマクロな現象とが必ずしも1線形加算的なつながりを持たないという問題がある．現在の計算機では取り扱える原子，分子数が制限され，計算結果は限定された原子分子数の範囲の分子原子レベルでの現象ということになる．問題はこの結果がより大きなマクロな現象の一部（ミクロ）を表すかということである．マクロな現象は確かに分子原子レベルでの現象の何らかの平均過程の結果であるが，その平均過程においてはミクロ現象の単なる線形加算平均ではなく，非線形的な結合であったり非一様性のみが強調されたりすることが多い．問題はこうしたミクロとマクロの繋がりを分子動力学にどう取り入れるかということである．
マクロ性が特に強調される熱流体の現象では，これらのミクロ，マクロのつながりの問題は特に重要であり分子動力学の活用は慎重な考慮が必要である．しかし，その活用は計り知れない効用を持つ為，今後ますます盛んに研究されると考えられる．

1.3

液晶分子

一般的に物質は常温，常圧のもとで固体・液体・気体のいずれかの状態で存在している．

2

そして物質の状態は温度や圧力に依存して変化する．
しかし液晶状態を示す物質は他の物質とは異なり，結晶から液体には直接に転移しない．
液晶は結晶と液体の両方の性質を示す中間の状態を経て液体になる．
液晶状態を示す物質は温度Tmになると結晶から融解して粘りのある濁った液体になり，さらに温度がTi に上昇す
ると透明な流動性の高い液体になる．このTmとTiの間の状態が液晶状態（物質の第四の状態）
である．
この時，濁った液体状態（液晶状態）は工学的に異方性であり，透明な状態（液体あるいは等方相と呼れる）は光学的に等方相である．

気化
液体

気体
凝縮
昇華

昇華

固
凝
解

(c)

(d)

融
(b)
結晶

液晶
(a)

図1‑1

物質の状態
(a)：結晶から液晶への転移（昇温過程），Tmに対応(b)：液晶からの液体への転移（昇温過程），Tiに対応(c)：液体から液晶へ転移（降温過程）(d)：液晶から結晶への転移（降温過程）

Tm
結晶

Ti
液晶

温度
液体

(等方性流体) (等方性)3

（１）ネマチック液晶，Tm＝21℃，Ti＝47℃

CH 3O

CH = NC4 H 9（２）スメクチック液晶，Tm＝114℃，Ti＝120℃

N =N

C 2 H 6 OCO

図1‑2

COOC 2 H 5

主な液晶物質の相図

液晶分子は主に棒状・円盤状であり，その特異な分子形状によって，重心位置の3次元性は失っているが分子の配向が残り，液体・固体とは違う4つ目の状態としての挙動をしめす．
温度転移型液晶相を形成する液晶分子の中で，棒状有機低分子はもっとも有名な分子であり，例えば図1‑3に示すように(a)p‑アゾキシアニソール(PAA)，(b)N‑(p‑メトキシベンジリデン)‑p‑ブチルアニリン(MBBA)，(c)4‑シアノ‑4 ‑ペンチルビフェニール(5CB)などがある．
液晶相を形成する棒状分子はかなりの種類が存在し，またさまざまな分子が合成されている．液晶分子の構造は主に剛直コア部である2個以上のベンゼン環で構成される部位と，その両端にそれぞれ炭化水素鎖，あるいは極性基をもつ．剛直コア部は液晶秩序を形成する為に大きな役割をはたしており，工業的機能性にかかわる屈折率，誘電率の異方性，磁気異方性など重要な性質に関係していると考えられる．またベンゼン環のコア部だけでは液晶相とはならずに結晶化してしまうことから，両サイドに付く炭化水素鎖も液晶相の安定化に重要な役割を担っていると考えられる．

(a)p‑アゾキシアニソール(PAA)

CH 3O

N =N

OCH 3

O
4

(b)N‑(p‑メトキシベンジリデン)‑p‑ブチルアニリン(MBBA)

CH 3O

CH = NC4 H 9(c)4‑シアノ‑4 ‑ペンチルビフェニール(5CB)C 5 H 11NC

図1‑3温度転移型液晶を形成する棒状有機低分子の例

側鎖部の長さにより，融点・沸点と分子量の関係から液晶相転移温度は上昇することがわ
かる．そして，スメクチック相などの対象性の低い液晶相は側鎖の長い方がが安定に現れるということが一般的である．
分子内に不斉が存在するような場合には，巨視的に光学活性を持つ液晶相が発現する．例
えばネマチック相を形成する場合は，コレステリック相が現れ，スメクチックC相の場合は，強誘電性液晶であるキラルスメクチックC相となる．分子内の局所的な双極子の存在も，液晶相の性質に大きく左右し，特に双極子相互作用は，応用的な見地からもネマチック相の誘電異方性，強誘電性液晶相の自発分極の大小を制御するという意味で注目されている．分子形状に不飽和結合や分岐など，分子全体の形状を左右するような性質を側鎖に導入すると，
回転楕円体で表されるような棒状分子の1軸対称性が失われる．このようなミクロな分子形状の2軸性や分子頭尾の区別の存在により巨視的な液晶相の構造・性質に大きな影響が現れ
る．

5

第2章

計算方法

分子動力学法は多粒子系のモデルのそれぞれの粒子一つ一つについて相互作用している力を計算し，運動方程式に代入し数値的に解析し次の瞬間の粒子の配列を求めるというものである．この計算の中で最も困難であるのが力の計算であり，この力の計算を行い個々の粒子の軌跡を求める為に数値積分法を導入する必要がある．そこで分子動力学シミュレーションを適用する時間に渡ってエネルギーが保存される為には，少なくとも２次のアルゴリズムが必要となる．以下に記す verlet 法は代表的な２次のアルゴリズムである．

2.1 verlet 法

verlet 法は非常によく利用される基本的なアルゴリズムの１つであり，verlet 法は方向性の無い系と方向性のある剛体系の粒子の並進運動に用いられる．
周期境界条件を用いて任意の時刻ｔにおける基本セル内の粒子 i の位置ベクトルを Ri(t)とし粒子に働く力を Fi(t)とすると，Newton の運動の第２法則より&& (t )Fi (t ) = mi Ri
（3.1）
一方微小な時間幅⊿t 後の粒子の位置 Ri(t+⊿t)は，Taylor 展開により&& (t )(∆t ) 2 &R&& (t )(∆t ) 3
R&
Ri (t + ∆t ) = Ri (t ) + Ri (t )∆t + i
+ i
+•••
2!
3!
（3.2）
になる．同様に&& (t )(∆t ) 2 &R&& (t )(∆t ) 3
R&
Ri (t − ∆t ) = Ri (t ) − Ri (t )∆t + i
− i
+•••
2!
3!
（3.3）
これらの両辺を加え，⊿t を十分に小さい値をとれば

6

&& (t − ∆t ) + (∆t ) Fi (t )Ri (t + ∆t ) = 2 Ri (t ) − R
i
mi
2

（3.4）
が導かれる．つまり，時刻 t-⊿t と t における位置，および t における i に働く力 Fi(t)を知れば t+⊿t における位置を知ることが出来る．
これを次々に繰り返せば，i の位置の時間経過すなわち軌跡を知ることができる．この方法は verlet 法と言われる．(3.4)式から分かるように多くの計算を必要とする部分は，力 Fi(t)を求める部分である．

∂φ ij (r ) Ri − R j
)(
)
∂r
r
j ≠i

Fi (t ) = − Σ gradφ ij = − Σ (
j ≠i

(3.5)
の関係を用いて 2 体ポテンシャルφij(r)から求められる．ここでｒは粒子 i,j 間の距離である．

r = Ri − R j(3.6)

φij(r)は粒子 i,j 間の距離の関数として表した 2 体力ポテンシャルで MD でもっとも重要な入力情報である．
対象とする系や研究の目的に応じてどのようなポテンシャルを用いるかを十分に考慮しなければならない．
(3.5)式において∂φ ij (r ) / ∂r * (1 / r )(3.7)の部分は２体力ポテンシャルφij の形が分かっていれば，あらかじめ計算しておくとよい．
また，速度はvi (t ) ={Ri (t + ∆t} − R j (t − ∆t )
2∆t
(3.8)

から求めることができる．

7

2.2 周期境界条件

周期境界条件は，現象が空間的に周期性を持つときはもちろんであるが，分子数の少ない系を見かけ上無限の分子数の系に見立てるときに仮定される．
すなわち，周期的な大きさ，あるいは考えている分子系の空間の大きさを L とするとき，基本位置 Xi(i=1,2,・・)に有る分子をxi ± nL(n=0,1,2,・・・)(4.1)に配置する（図 2-1）．したがって，分子 Xi のポテンシャルはU ( x i ) = Σ Σ U ( ( x j ± mL) − ( x i ± mL)) − ( x i ± nL))
j n,m

(4.2)

U ( xi ) = Σ {U ( x j − xi ) + U ( x j − xi + L ) + U ( x j − xi − L )j+ U ( x j − xi + 2 L ) + U ( x j − xi − 2 L ) + • • •}(4.3)
となる． L が十分大きな系では，ほとんどの場合最初の 3 項で決まり．この中でも特に分子間距離の最小なものが支配的である．すなわち，

xij ,min = min{ xi − x j , x j − xi + L , x j − xi − L }
(4.4)
とすると

U ( xi ) = ΣU ( xij , min )j(4.5)と表される．L の大きさはその値によって結果に影響が現れないような十分大きな値を選ばなければならない．

8

(Xi,Yi+Ly)(Xj,Yj+Ly)(Xi+Lx,Yj)

Ly(Xi,Yi)(Xj,Yj)(Xj+Lx,Yj)Lx

図 2-1 周期境界条件：Xi→Xi±nL

92.3 Lennard-Jones potential
Lennare-Jones potential は分子動力学においてよく使われる potetial の一つである．分
子間ポテンシャルとは粒子間相互作用のことであり，厳密に記述しようとすれば，電気的に中性な分子についてのエネルギーは，原理的に量子力学的な計算により第一原理から構成することが出来るがその計算は非常に困難であるので，通常は簡単な現象論的なものを採用する事が一般的である．
単純液体について最も重要な特徴は r が小さいところでの強い斥力と r が大きいところでの弱い引力である．小さな r での斥力はパウリ(Pauli)の排他原理によるものであり，つまり，2 つの分子の電子雲は重なりを避ける為に変形しなければならず，結果として電子のい
くつかは異なる量子状態に入ることになる．これらの効果により運動エネルギーが増大して，電子間にはコアの(芯の)斥力と呼ばれる実効的な斥力が生じる．r が大きい時に支配的となる弱い引力は，各分子がお互いに分極を起こすことによるものであり，この引力はファン・デル・ワールス（van del waals）力と呼ばれている．
こ れ ら の 現 象 を よ く 表 す 現 象 論 的 な u(r) の 形 の 1 つ が レ ナ ー ド ・ ジ ョ ー ン ズ（Lennard-Jones）ポテンシャルである．

σσu (r ) = 4ε [( ) n − ( ) m ]
r
r
(5.1)
が，一般的なレナードジョーンズポテンシャルであり特に n=12，m=6 のものを 12-6 型のレナードジョーンズポテンシャルとしてよく用いられる．
ポテンシャルの谷の深さをεによってしめし，σはその時の原子・分子間位置を与える
パラメータである．

10

potential U

L-J (Lennard-Jones)
U=4ε{(σ/R)12 -(σ/R)6}

atom and molecular distance
図 3-1 Lennard-Jones potential

この potential を無次元化して使用する．無次元化することによって 2〜3 種類の分子で用いるときには，一つの分子に対して基礎式を無次元化することにより他の分子の影響をそれと相対的に考えることが出来非常に使いやすくなる．

11

2.4 Gay-Berne potential

Gay‑Berne potentialは回転楕円体の形をした1軸性の分子について，その複雑な分子間の相互作用を，分子の重心の位置と，対称軸の方位のみを関数として，解析的に書き表したものである．
回転対称軸をもつ1軸性の分子を考えるとき，i番目の分子の対称軸の向きを表すベクト
ルと重心の位置を表す座標をそれぞれ uˆ i , rˆ とする．また，ｊ番目の分子とi番目の分子の重心間の相対位置を表すベクトルの絶対値を rij = ri − r j で表し，その方向の単位ベクトルをr̂ij で表す．このとき，これらの2分子間のGay‑Berne potentialは次式で表される．
基礎式

U (uˆ i , uˆ j , rˆ) = 4ε (uˆ i , uˆ j , rˆ) * [{σ0r − σ (uˆ i , uˆ j , rˆ) + σ 0}12 − {σ0r − σ (uˆ i , uˆ j , rˆ) + σ 0

}6 ]

(6.1)

ここで Uˆ 1 , Uˆ 2 は 2 つの集団の配向を示す単位ベクトルであり， r は分子の位置を表すベ
クトルである．

σ (uˆ i , uˆ j , rˆ) は
12

σ (uˆ i , uˆ j , rˆ) = σ 0 [1 − χ {( rˆ ⋅ uˆ i + rˆ ⋅ uˆ j ) 2
1 + χ (uˆ i ⋅ uˆ j )+( rˆ ⋅ uˆ i − rˆ ⋅ uˆ j ) 2
1 − χ (uˆ i ⋅ uˆ j )}] −1 / 2

(6.2)
となり，非等方性パラメータχはχ = {(σ e / σ s ) 2 − 1} /{(σ e / σ s ) 2 + 1}

（χ>0）
(6.3)

である．
χは分子間斥力の異方性を調整するパラメータで形状の異方性を決定する．
すなわち，χ>1 のときは棒状の分子を表す，またχ<1 のときは円盤状の分子を表す．

σ e は end-to-end の配置の時のパラメータ，σ s は side-by-side の配置の時のパラメー
タを表す．

χ ′ = {1 − (ε e / ε s )1 / µ }{1 + (ε e / ε s )1 / µ }12（χ´>0）

(6.4)
χ´は引力ポテンシャルの深さの異方性を調整するものであり，μ，νは調整用の指数である．
ポテンシャル谷の深さはε (uˆ1 , uˆ 2 , rˆ) = ε 0 ε ν (uˆ1 , uˆ 2 )ε µ (uˆ1 , uˆ 2 , rˆ)(6.5)
で表すことができる．
また，
−

1

ε 0 ε (uˆ1 , uˆ 2 ) = ε 0 {1 − χ 2 (uˆ1 ⋅ uˆ 2 ) 2 } 2
(6.6)
である．
μ，ν を適 当な値 を入 れるこ とで 性質を 変え ること がで き，目 的に あった 性 質 のGay-Berne

potential を使用しなければならない．Gay-Berne モデルでは，等方相，ネマチック相，スメクチック B 相などを表すことが出来る．

13

第3章

計算結果

以下に計算結果を示す．以下の図において，点は系における分子を表し，黒い部分は系の範囲を示す．初期値は，速度・位置ともにランダムによって与えられ，図 1-A〜図 1-L は
分子数 100 個，温度 5K，時間 t*=0.2 刻みで t*=0〜t*=2 までと，t*=12 の結果を示す．図1-A での初期配置から図 1-B に到る t*=0.2 の間にほぼ均一に系全体に分子が広がっているのが分かる．その後，結晶化構造特有の六角格子状の配列に到るのが図 1-L で見える．図1-L で結晶化に到ったことから，その過程において相転移が起こっていたと考えられる．
図 2-A〜図 2-L は，分子数 100 個，温度 100K，時間 t*=0.2 刻みで t*=0〜t*=１までと，
t*=12 の結果を示す．図 2-A は先と同様に初期配置であり，図 2-B に至る t*=0.2 の間に系
全体に分子が広がっているのが分かる．そして図 2-B 以降，図 2-G に至るまで先の結果に見られたような結晶化の兆候は見られず常に流動していた．よって 100K においては液体状態で安定していると考えられ，相転移を確認することは出来なかった．
図 3-A〜図 3-J は分子数 100 個，温度 90K〜1K，時間 t*=1 の結果を示す．図 3-A から
図 3-H までは分子が一様に散らばって流動しており，結晶化の兆候は見られなかった．しかし 10K まで温度を下げたときの図 3-I で図の一部に結晶化の兆候が見られ，さらに 1Kにまで下げた図 3-J においてはっきりと結晶化構造特有の六角格子状の配列が，欠陥が多数見られるが確認できた．これより，10K 辺りから相転移が起こり始めるのではないかと考
えられる．

14

分子数 100 個，温度 5K，時間 t*=0〜t*=2 で（t*=0.2 間隔），t*=12 の結果を図 1-A〜図
1-L に示す．

図 1-A

分子数 100 個，温度 5K，時間 t*=0

図 1-B

分子数 100 個，温度 5K，時間 t*=0.215図 1-C

分子数 100 個，温度 5K，時間 t*=0.4

図 1-D

分子数 100 個，温度 5K，時間 t*=0.6

16図 1-E

分子数 100 個，温度 5K，時間 t*=0.8

図 1-F 分子数 100 個，温度 5K，時間 t*=1.017図 1-G

分子数 100 個，温度 5K，時間 t*=1.2

図 1-H 分子数 100 個，温度 5K，時間 t*=1.418図 1-I

分子数 100 個，温度 5K，時間 t*=1.6

図 1-J

分子数 100 個，温度 5K，時間 t*=1.819図 1-K

分子数 100 個，温度 5K，時間 t*=2.0

図 1-L

分子数 100 個，温度 5K，時間 t*=1220分子数 100 個，温度 100K，時間 t*=0〜t*=1（t*=0.2 間隔），t*=12 の計算結果を図 2-A〜図 2-L に示す

図 2-A

分子数 100 個，温度 100K，時間 t*=0

図 2-B

分子数 100 個，温度 100K，時間 t*=0.221図 2-C

分子数 100 個，温度 100K，時間 t*=0.4

図 2-D

分子数 100 個，温度 100K，時間 t*=0.622図 2-E

分子数 100 個，温度 100K，時間 t*=0.8

図 2-F 分子数 100 個，温度 100K，時間 t*=1.023図 2-G

分子数 100 個，温度 100K，時間 t*=1224分子数 100 個，温度 90K〜10K（10K 間隔），1K まで，時間 t*=1 の結果を図 3-A〜図
3-J に示す．

図 3-A

分子数 100 個，温度 90K，時間 t*=1

図 3-B

分子数 100 個，温度 80K，時間 t*=125

図 3-C

分子数 100 個，温度 70K，時間 t*=1

図 3-D

分子数 100 個，温度 60K，時間 t*=126

図 3-E

分子数 100 個，温度 50K，時間 t*=1

図 3-F 分子数 100 個，温度 40K，時間 t*=127

図 3-G

分子数 100 個，温度 30K，時間 t*=1

図 3-H 分子数 100 個，温度 20K，時間 t*=128

図 3-I

分子数 100 個，温度 10K，時間 t*=1

図 3-J

分子数 100 個，温度 1K，時間 t*=1

29

第4章

結言

考察
本 シ ミ ュ レ ー シ ョ ン で は 温 度 の 変 更 と ， 時 間 の 推 移 に よ る 変 化 を Lennard-Jones
potential を利用した分子モデルを作成し，相転移を観察した．最初に 5K で t*=0 から t*=12まで観察すると，分子は初期状態から分子はゆっくりと相転移状態を見せて結晶状態に向かい安定する．同様に 100K 状態で t*=0 から t*=12 まで観察すると，分子は液体状態のままであり，結晶化する様子は確認できなかった．次に 90K から 10K づつ順に温度を下げそれぞれの t*=1 の計算結果を見た．90K から 20K 辺りまでは結晶化の兆候は見られなかったが，10K の辺りで結晶化してきているように分子が配置してきているように見える．このことから，およそ 10K 前後でこのモデルは相転移を起こし結晶化が始まるのではないかと考えられる．また，相転移は一様に始まるのではなく，一部から全体に広がっていくよ
うに見える．よってモデル中の相転移し易い配置になっている場所から順に相転移が起こ
り，先の現象が観察できたと考えられる．これはモデル中の分子が一様に分布していないことに起因していると思われる．
今回は図による観察のみであった為どの時点で相転移が始まり，どの時点で相転移が完了したのかがはっきりわからなかった．それは平衡分布関数や，動径分布関数などの相関関数を使用しなかった為と考えられ，相転移時の相の変化を明確にする為にこれらの方法を取り入れる必要がある．
本研究では Lennard-Jones モデルを使用したが，液晶分子の挙動をシミュレーションする為には Gay-Berne potential を使用したモデルが必要である．今回使用したシミュレーションモデルは球状分子であり分子に方向を持たない．そのため液晶分子の挙動をシミュレ
ーションすることは出来ない．液晶分子の挙動をシミュレーションするには，方向性を持
つ分子形状を使用する必要があり，その為の potential として Gay-Berne モデルは適している．Gay-Berne モデルでは回転楕円体を分子形状に利用するため分子に方向が発生し液晶に近い挙動を示す．このモデルを使用し，今後，液晶の相転移を観察する必要がある．

30

参考文献
（１） 小竹進，熱流体の分子動力学，丸善株式会社，（１９９８）（２） 液晶便覧編集委員会編, 液晶便覧,丸善株式会社,（２０００）（３） ハーベイ･ゴールド，ジャン・トボチニク，計算物理学入門，（株）ピアソン・
エデュケーション，（２０００）（４） Luckhurst GR, Stephens RA, Phippen RW, Liquid Crystals, vol8, no4, 451‑464,(1990)

31

謝辞
本研究を行うにあたって蝶野重臣教授，辻知宏助手には終始丁寧な御指導を頂き深く感謝しております．また，知能流体力学研究室の皆さんからは，多大な援助を頂き，合わせて感謝します．

32