高分子材料シミュレーションの概要と事例

高分子材料は、空間スケールで言えばnmからμmまでのマルチスケール特性を有しているため、現象や物性の発現がどのスケールに起因するか注意して考える必要がある。実験・計測技術の発展に支えられて多くのメカニズムが明らかにされてきているが、それだけでは捉えきれない部分に関して、シミュレーション技術が重要な役割を果たしている。また、最近ではデータサイエンスを用いた物性予測なども進んできているが、シミュレーション技術はそこでも相補的な役割を担っている。
以下では、高分子のシミュレーション技術について代表的ないくつかの手法の概要、それぞれの連携、ソフトウェアなどを紹介する。各手法がターゲットとする空間スケールに注目して見ていくとわかりやすいため、小さなスケールから順番に説明する。なお、各手法の文献を挙げていくときりがないため、基本的には参考文献[1]を参照頂きたい。

主に数nm以下をターゲットとする手法で、分子や結晶内の電子状態をSchrodinger方程式に基づいて計算する。特に励起や分極、分子間力、化学反応など、物質の電子的な性質を評価する際に用いられる。分子軌道法（Molecular Orbital method : MO）や、電子密度を計算対象とすることで計算負荷が低減できる密度汎関数理論（Density Functional Theory : DFT）があるが、電子相関の取り入れ方や、基底関数の選び方などによって手法が細分化される[2]。

図1 DFTによるCu表面への分子の吸着エネルギーの評価
（左：propane/butane/pentane/hexane分子の吸着エネルギー曲線、右：計算モデル）

図1はCuの結晶表面に分子が吸着した際の、表面からの距離とエネルギーの関係をvan der Waals力を考慮したDFTで計算した結果（分子種ごとに点で表示）を示しており、後述の分子動力学で用いるためのLennard-Jones（LJ）ポテンシャルにフィッティングした曲線も併せて表示している。このような表面吸着エネルギーは界面の設計において重要となるが、妥当なLJポテンシャルパラメータが存在しないケースもあり、DFTなどによる評価が必要となる。

全原子分子動力学（Full Atomistic Molecular Dynamics : FAMD）は、一つの粒子（質点）が一つの原子を表し（図2左）、各質点のダイナミクスをNewtonの運動方程式に基づき計算する。各原子に働く力は関数とパラメータによって与えられ、その中には前述のLJポテンシャルなどが含まれる。パラメータはデータベースとして公開されていることが多いが、より精度を上げるため、前述のDFTのような一つ下のスケールの手法と組み合わせる。FAMDを用いることで、官能基など分子構造を変更した際の分子形状やダイナミクスへの影響を評価することができる。例えばバルク状態の高分子の密度や弾性率、自由体積の分布やその中の気体分子の拡散、配向（結晶）構造などが挙げられる。様々な用途に適用できるため普及しているが、一般的なコンピュータを用いる限り扱える空間スケールが10nm程度、時間スケールは最近のソフトウェアとハードウェア（GPUなど）を用いる場合でも100nsのオーダーであるため、例えば高分子の長時間の緩和現象を扱うには不足する。

図2　左：FAMDによるPolyisopreneの分子構造、右：CGMDのイメージ

粗視化分子動力学（Coarse Grained Molecular Dynamics : CGMD）では、一つの粒子（質点）がいくつかの原子の集まりを表現する（図2右）。これによって、FAMDと比べて時間空間ともに大きなスケールの現象を扱うことが可能となる（10〜100nm、μsなど）例えばゴム材料の引っ張り特性を評価するのであれば、ゴム領域を扱えるような時間（変形速度）の計算が必要となるが、FAMDでは難しくCGMDを用いることで可能となる。粗視化ユニット間の相互作用（粗視化ポテンシャル）はFAMDをベースに決定する場合[3]と、簡易的なパラメータを与えて計算後に意味づけをする場合がある。
図3ではStyrene-Isoprene-Styrene（SIS）のようなトリブロック共重合体が形成する相分離構造（熱可塑性エラストマー）の一軸伸長時の挙動を評価している[4]。初期構造の構築には後述の平均場法の計算結果を用いている。伸長時のドメインや分子鎖の変形の様子が捉えられている。

図3　CGMDによるトリブロック共重合体の相分離構造の一軸伸長計算
（左：初期構造、右：変形中の構造（上側は一本の鎖、下側はドメイン構造））

散逸粒子動力学（Dissipative Particle Dynamics：DPD）[5]はCGMDの一種であるが、粒子間相互作用にすり抜けを許すソフトポテンシャルが用いられることで、相分離やフィラー分散構造の評価などに適している。ダイナミクスには流体効果も含まれている。平均場法（Mean Field : MF）では、高分子の各成分を濃度（体積分率）場で表現することで相分離のダイナミクスや平衡状態までを評価する。各成分の拡散および流体力学とのカップリングが可能である。中でも自己無撞着場理論（Self-Consistent Field Theory : SCFT）を用いることで高分子鎖の形状も考慮できる。これらの手法では成分（粒子）間の相互作用としてFlory-Hugginsのχパラメータが用いられるが、前述のFAMDや、量子化学計算を用いて推算する手法[6]も提案されている。
図4（左）はDPDを用いて計算した高分子電解質と水の相分離構造、図4（右）は平均場法を用いて計算した高分子3成分系の相分離（コア・シェル）構造である。数100nmまでの空間スケールがターゲットとなる。

図4　左：DPDによる高分子電解質の相分離構造、
右：平均場法による高分子3成分の相分離構造

高分子鎖を絡み合い点間分子量の単位で粗視化することで、高分子鎖の絡み合い構造に基づくレプテーション（reptation）ダイナミクスを計算することができる。Slip LinkモデルやPrimitive Chain Network（PCN）モデル、前述のDPDに高分子鎖の絡み合い効果を加えることが可能なSlip Springモデル[7]などがある。
図5は、Polyisoprene（分子量約48,800）のG’（貯蔵弾性率）とG”（損失弾性率）を評価したものである。スナップショット中の折れ曲がりが絡み合い点に相当する。

図5　PCNによる動的粘弾性（左：スナップショット、右：マスターカーブ）

前述のDPDや平均場法で得られた相分子構造を用いて、連続体モデルに基づいた有限要素法（Finite Element Method : FEM）の計算を実施することで、平均弾性率や熱伝導率などの物性評価が可能となる。各成分の材料物性はパラメータとして入力される。図6（左）は図4（右）で得られた相分離構造に対して変形を加えた際の界面のひずみエネルギー分布を示しており、図6（右）はDPDで得られた高分子中のナノフィラーの分散構造に対して熱伝導計算を実施した際の熱流束の分布である。図6の計算では何かしらのシミュレーションで得られた構造を元にしているが、その代わりに散乱試験やCTなどのデータを用いることや、仮想的な構造を作り出すことで物性と構造の関係を数値実験するような活用法もある。

図6　有限要素法による機械特性の評価
（左：相分離構造中のひずみエネルギー分布、右：ナノフィラー分散構造中の熱流束分布）

量子化学・密度汎関数理論についてはSIESTAとABINIT-MPが同梱されており、Gaussianなどとのインターフェイスが提供されている。全原子・粗視化分子動力学についてはCOGNACとVSOP、GENESISが同梱されており、LAMMPS、GROMACS、HOOMD-blueとのインターフェイスが提供されている。散逸粒子動力学（DPD）については分子動力学と同様であるがGENESISとGROMACSには含まれない。平均場法についてはSUSHIが同梱されている。レプテーションダイナミクスについてはCOGNACおよびVSOPのSlip-Springモデルを適用したDPD機能、PASTA、NAPLESが該当する。連続体モデルについてはMUFFINが該当するが、その他商用ソフトウェアとのインターフェイスも提供されている。