(57)【要約】【課題】層状化合物の挙動について効率的にシミュレーションする。【解決手段】層状化合物について各層を力学的単位となる剛体シートとして粗視化し、粗視化した各層の剛体シートについて、隣接層の剛体シートとの層間相互作用、内部振動に起因する内部作用を決定する（Ｓ１〜Ｓ３）。その後に、複数層の剛体シートからなる層状化合物の挙動についての計算を行う（Ｓ５）。【選択図】図１８
2011058947.tif 000002
【特許請求の範囲】
【請求項１】
層状化合物について各層を力学的単位となる剛体シートとして粗視化し、
粗視化した各層の剛体シートについて、隣接層の剛体シートとの層間相互作用、内部振動に起因する内部作用を決定し、
その後に、複数層の剛体シートからなる層状化合物の挙動についての計算を行うことを特徴とするシミュレーションシステム。
【請求項２】
請求項１に記載のシミュレーションシステムにおいて、
前記層間相互作用は、相互ポテンシャルを含み、相互ポテンシャルは電子レベルおよび原子レベルの微視的なシミュレーションにより決定することを特徴とするシミュレーションシステム。
【請求項３】
請求項２に記載のシミュレーションシステムにおいて、
相互作用ポテンシャルが、二次元の周期性および末端効果の少なくとも１つを有するものであることを特徴とするシミュレーションシステム。
【請求項４】
請求項１に記載のシミュレーションシステムにおいて、
前記内部作用は、原子レベルの微視的なレベルにおいて決定することを特徴とするシミュレーションシステム。
【請求項５】
請求項１〜４のいずれか１つに記載のシミュレーションシステムにおいて、
さらに、２つの固体原子層に挟まれた層状化合物に対して圧力およびせん断を印加し、せん断応力を含む力学的応答を得ることにより、層状化合物の挙動をシミュレーションすることを特徴とするシミュレーションシステム。
【請求項６】
コンピュータに、
層状化合物について各層を力学的単位となる剛体シートとして粗視化させ、
粗視化した各層の剛体シートについて、隣接層の剛体シートとの層間相互作用、内部振動に起因する内部作用を決定させ、
その後に、複数層の剛体シートからなる層状化合物の挙動についての計算を行わせる、
ことを特徴とするシミュレーションプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、グラファイトや雲母といった層状化合物についての挙動をシミュレーションするシミュレーションシステムに関する。
【背景技術】
【０００２】
グラファイトや雲母といった層状化合物の中のある種の物質は、低摩擦特性を示すことがある。雲母はそれほど低摩擦ではないが、グラファイトや、二硫化モリブデンといった固体は摩擦係数が低い。こうした化合物は、宇宙空間などオイルの使えない人工衛星の部品などにおいて固体潤滑剤として用いられたり、結晶のもととなる化合物を自動車のエンジンオイルに添加剤として混ぜてしゅう動面で結晶化させて低摩擦にする、というように使われる。
【０００３】
層状化合物において、層間相互作用がvan der Waals力によるグラファイトのような結晶であっても、二硫化モリブデンのようなイオン結晶であっても、摩擦力が低い理由は、層間ですべりが起こるためであると説明されている。
【０００４】
このような層状化合物について、詳細な解析が可能であれば、用途に応じてどのような層状化合物を使用すべきなどの構造設計などが行える。
【０００５】
これに関連して、特許文献１には、流体中の高分子（紐状ミセル）を粗視化したブラウン動力学シミュレーションが示されている。また、高分子の粗視化手法については、特許文献２、３などに示されている。さらに、トライボ関係の分子シミュレーションについては、特許文献４などがある。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特許０３９８０６００号公報
【特許文献２】特開２００２−８０４０６号公報
【特許文献３】特開２００６−３３８４４９号公報
【特許文献４】特開平８−３５８３１号公報
【発明の概要】
【発明が解決しようとする課題】
【０００７】
特許文献１〜３では、二次元的な周期性を持つ相互作用を層間相互作用として規定していないため、層状化合物の運動を現実的に扱うことが出来ない。また、液体に対応する手法であるため、拡散係数などの各粗視化粒子の内部振動に起因する作用が等方的である系を取り扱い、異方性の作用が大きい層状化合物の特性を再現できない。従って、特許文献１に挙げられた方法では、層状化合物の挙動を解析できない。
【０００８】
特許文献４では、原理的には全原子モデルにより層状化合物の摩擦挙動を扱うことができる。しかし、層状化合物についての現実的な挙動を解析しようとすると、対称となる原子数が非常に大きく、従って全原子モデルを利用すると膨大な計算時間を必要となる。
【課題を解決するための手段】
【０００９】
本発明は、層状化合物について各層を力学的単位となる剛体シートとして粗視化し、粗視化した各層の剛体シートについて、隣接層の剛体シートとの層間相互作用、内部振動に起因する内部作用を決定し、その後に、複数層の剛体シートからなる層状化合物の挙動についての計算を行うことを特徴とする。
【００１０】
また、前記層間相互作用は、相互ポテンシャルを含み、相互ポテンシャルは電子レベルおよび原子レベルの微視的なシミュレーションにより決定することが好適である。
【００１１】
また、相互作用ポテンシャルが、二次元の周期性および末端効果の少なくとも１つを有するものであることが好適である。
【００１２】
また、前記内部作用は、原子レベルの微視的なレベルにおいて決定することが好適である。
【００１３】
さらに、２つの固体原子層に挟まれた層状化合物に対して圧力およびせん断を印加し、せん断応力を含む力学的応答を得ることにより、層状化合物の挙動をシミュレーションすることが好適である。
【００１４】
また、本発明は、このようなシミュレーションを行うためのプログラムである。
【発明の効果】
【００１５】
本発明によれば、層状化合物の摩擦力のような力学的特性、熱安定性などの挙動を効率よく計算できる。これによって、所望の分子特性を示す層状化合物の分子設計を行うことが可能となる。
【図面の簡単な説明】
【００１６】
【図１】層状化合物の摩擦メカニズムを説明する図である。
【図２】グラフェンシートの粗視化と６自由度を説明する図である。
【図３】層状化合物のマルチスケール解析スキームを示す図である。
【図４】グラファイトの分子構造を示す図である。
【図５】炭素原子間のLennard-Jonesポテンシャルを示す図である。
【図６】炭素原子−グラフェンシート間のポテンシャル（Ca,Cbはa,b間上の格子点）を示す図である。
【図７】炭素原子−グラフェンシート間のポテンシャル（ｘ，ｙ空間）を示す図である。
【図８】炭素原子−グラフェンシート間のポテンシャル（ｘ，ｙ空間）を示す図である。
【図９】シミュレーションのモデルを示す図である。
【図１０】熱平衡状態におけるシートのｘ方向の変位を示す図である。
【図１１】熱平衡状態におけるシートのｙ方向の変位を示す図である。
【図１２】熱平衡状態におけるシートの全エネルギーＵのゆらぎを示す図である。
【図１３】ＨＯＰＧ（高配向性グラファイト）上の”移着片”モデルを示す図である。
【図１４】非対称系におけるすべり運動下のシートのｘ方向の変位を示す図である。
【図１５】非対称系におけるすべり運動下のシートのｙ方向の変位を示す図である。
【図１６】非対称系におけるすべり運動下のシートの全エネルギーＵの経時変化を示す図である。
【図１７】非対称系におけるすべり運動下のシートのｘ，ｙ方向の変位を示す図である。
【図１８】計算のフローチャートである。
【発明を実施するための形態】
【００１７】
以下、本発明の実施形態について、図面に基づいて説明する。
【００１８】
「シミュレーションシステムの構成」
本実施形態は、層状化合物について、モデルを用いたシミュレーションを行うシミュレーションシステムおよびそのプログラムに関し、ＣＤ、ＤＶＤ、通信などによって、コンピュータの記憶部に前記モデルを利用したシミュレーションを行うプログラムをインストールし、プログラムを実行し、各種条件、データを入力して、摩擦現象や、摩擦係数などの物性値を求めるためのシミュレーションを行う。なお、シミュレーションプログラム自体は、すでに販売されているものを用い、そのシミュレーションプログラム内のモデルを本実施形態のプログラムによって設定することも好適である。
【００１９】
そして、シミュレーションは、グラファイトや雲母といった現実の層状化合物について行われ、層状化合物からなる機能性分子の構造設計が行われる。
【００２０】
「層状化合物の挙動解析」
例えば、無限に長い板（グラフェン[1]）が積層しているグラファイトの下地の上を、半径数ｍｍ程度の鋼球が接触して相対すべり運動を行っている状態を考える。このとき、グラフェン層間ですべりが起こることは有り得ない（単位面積当たりの層間摩擦力が有限であれば層間すべりに要する力は無限になる）。また実際の物質では、高配向性のグラファイトであってもμｍオーダーの湾曲があるので、下地の固体内部で剪断が起こることは幾何学的に困難である。
【００２１】
したがって、現実的に考えると、図１に示すように、下地のグラファイトから鋼球側に、グラファイト断片の移着が生じ、その移着片と下地のグラファイトとの間ですべり摩擦が生じていると考えるべきであろう。機械試験に使う鋼球にはサブμｍオーダーの表面粗さがあるので、その先端に移着片があると考えられ、実際に、この存在を３次元表面構造解析顕微鏡や透過型電子顕微鏡（Transmission Electron Microscope;TEM）[2]により観察した研究もある。それでは、移着片と下地グラファイト面との摩擦の機構が疑問である。
【００２２】
本実施形態では、移着片と下地との間ですべりが起きているのか、それとも移着片の層間ですべりが起こっているのか、また移着片の形状、大きさとの関係は、といった問題を扱う。
【００２３】
グラファイトの低摩擦機構については、二層間の挙動については全原子分子動力学（Molecular Dynamics;MD）法により解析した例があるが[3]、原子間力顕微鏡（Atomic Force Microscope; AFM）の先端に移着したグラファイト片でさえ実際は数十万原子からなるため、そのような系全体のダイナミクスを全原子分子動力学法で扱うことは困難である。そこで、本実施形態では移着片を粗視化して、計算の簡略化を図る。
【００２４】
ところで、移着片を構成するグラフェンシートの内部振動には多数のモードがあるが、各グラフェンシートをそれぞれ一つの剛体とみなすと、この内部振動に起因する電子雲のゆらぎは、van der Waals力のゆらぎとして隣接する剛体の質点に及ぼす揺動力と見なすことができる。また、この剛体には並進３、回転３の６自由度があり、この中で重要な運動は並進の３個と、層方向の軸に対する回転であるｙａｗ角の４個である。移着片のすべり運動は、この４自由度におけるブラウン運動を有するグラフェンシート剛体が積層したものの運動として記述できる（図２参照）。
【００２５】
本実施形態では、これをMonte Carlo Brownian Dynamics（MCBD）法[4]で計算し、層状化合物の摩擦挙動をシミュレートする。MCBD法はコロイド溶液の計算手法として1991年に菊地らによって提案された手法であるが、後述するようにLangevin動力学と等価であり、熱ゆらぎが粒子間ポテンシャルと同等のオーダーとなるような現象においては有効な粗視化シミュレーション手法となる。
【００２６】
「シミュレーションの手法」
「ＭＣＢＤ法」
多粒子系を扱うための計算機シミュレーションの代表的な方法として、分子動力学（MD）法と、Monte Carlo（MC）法が用いられている。MD法は、各粒子についての運動方程式を連立させ、これを解くことにより、粒子の軌跡を追跡する。一方、MC法は統計力学における配置積分を、乱数を用いて計算する。それぞれの方法は、粗視化の度合いや求める物理量によって多くのバリエーションがあるが、アルゴリズムの基礎を、前者はNewton力学に、後者は統計力学に置いているという違いがある。
【００２７】
Brown運動は巨視的にはFokker-Planckの拡散方程式によって記述されるが、微視的にはLangevin方程式から議論を出発する理論が多い[5]。Langevin方程式は、運動方程式に確率的に加わる、溶媒分子の熱運動に起因する揺動力による項が入った形をとる。シミュレーションの方法として、このLangevin方程式の揺動項に乱数を用いるMD法が開発されている[6]。一方、MC法の代表的な方法にMetropolisらの方法[7]がある。これは、時点ｎにおける遷移確率が時点ｎ−１の状態のみによって決まるという遷移を繰り返す、Markov鎖の概念に基づいている。通常のMetropolis MC法においては、状態が遷移する間隔と物理的な”時間”との明確な対応はない。
【００２８】
ところで、Brown運動においては、粒子の移動をMarkov的であると仮定するので、Metropolis MC法において、状態が遷移する過程において”時間”との対応をとることができれば、同法もまたBrown運動のシミュレーションの方法となるはずである。
【００２９】
Metropolis MC法は平衡系を仮定しており、遷移の過程において生成されるMarkov鎖は、熱平衡状態における配置積分を得るために導入される。したがって、過渡過程をも含めた非平衡系におけるBrown運動のシミュレーションに直接適用することはできない。しかし、Brown運動における遷移の過程もMarkov的であるため、遷移確率の与え方に一定の条件を加えれば、Metropolis MC法で実現されるMarkov鎖は、物理的な”時間"に対応するBrown粒子の軌跡となるはずである[4]。
【００３０】
MCステップにおいて、系の状態を識別するラベルとして”時間”ｔを導入し、”時間”ｔにおける配置Ｘについての確率分布関数をＰ（Ｘ，ｔ）とする。このときサンプリングはMarkov過程を仮定しているので、Ｐ（Ｘ，ｔ）の”時間”変化は、次のマスター方程式により決定される[5]。
【数１】
2011058947.tif 000003
【００３１】
ここで、Ｗ（Ｘ｜Ｘ’）は配置Ｘ’からＸへの遷移についての単位”時間”あたりの遷移確率である。マスター方程式は系の大きさを表わすパラメータで展開することにより、近似的に解くことができる。ただし、Ｗ（Ｘ｜Ｘ’）は次の正準形をしているものとする。
【数２】
2011058947.tif 000004
【００３２】
ここでｘ’＝Ｘ’／Ω，ｒ＝Ｘ−Ｘ’である。このとき、ジャンプ能率_αν，λ（ｘ）を次式で定義する。
【数３】
2011058947.tif 000005
【００３３】
α１，０のとき、この展開は非線形のFokker-Planck方程式の第一近似として次式を与える。
【数４】
2011058947.tif 000006
【００３４】
ただし、τ＝ω−２ｔである。Metropolis MCの場合、ωは”時間”間隔Δｔにおける最大変位幅の逆数ととることができる。したがって、式（2）の右辺第一項は、
【数５】
2011058947.tif 000007
と計算される。そして、たとえば系のポテンシャルエネルギーＵ（ｘ）のｘについての導関数Ｕ’（ｘ）が負の場合は、十分大きいΩに対して
【数６】
2011058947.tif 000008
である。したがって、ジャンプ能率は、
【数７】
2011058947.tif 000009
となり、式（4）は次の拡散方程式に帰着する。
【数８】
2011058947.tif 000010
ここで、拡散定数Ｄは次式で定義される。
【数９】
2011058947.tif 000011
【００３５】
すなわち、Metropolis MC法における最大変位幅Ω−１を、粒子に働く力がMC"時間”間隔Δｔの間一定とみなせるほど十分小さくとることにより、Metropolis MC法はBrown運動のシミュレーションの方法となることが示された。この際、MC”時間”は物理時間と式（7）を介して対応する。
【００３６】
拡散係数は、実験値を用いるか、あるいはより微視的（原子レベル）なシミュレーション手法により決定することができるが、実験値を用いた例としては高分子電解質溶液中の低分子イオンの例がある[8]。
【００３７】
MD法を基礎とした方法に対するこの方法の利点は、粒子に働く力を計算する必要がなく、ポテンシャルエネルギーＵ（ｘ）の変化のみを計算すればよいということである。また、一般のMC法と比較して熱平衡値だけでなく、目的とする物理量の時間変化も与えるが、そのためには多数のトラジェクトリにわたる積算が必要である。
【００３８】
図３に、層状化合物のダイナミクスを扱うための本MCBD法を中心とするマルチスケール解析スキームの例を示す。上述したように、MCBD法におけるシミュレーションの物理化学的な現実性は、層間相互作用のポテンシャルエネルギーＵ（ｘ）およびシートの拡散係数Ｄを現実的なものとすることである。
【００３９】
本発明における一例として、Ｕ（ｘ）はグラフェン系（グラファイト、カーボンナノチューブ、Ｃ６０など）の炭素原子間のvan der Waalsポテンシャルとして有名なGirifalcoらの提案した文献値を用いる[9]。また、拡散定数Ｄについては、信頼性の高い実験値が存在する場合はそれを用いても良く、あるいは、よりミクロな全原子の分子動力学を行い決定することも可能である。全原子の分子動力学に必要な力場パラメータが存在しない、あるいは信頼性が問題となる場合は、さらに力場パラメータをab initio分子軌道計算などにより量子化学的に決定できる。
【００４０】
このマルチスケール解析スキームを用いることにより、グラファイトやh-BNのようなvan der Waals系のみならず、二硫化モリブデンや雲母のようなイオン系を含めた任意の層状化合物のダイナミクスを扱うことが可能となる。したがって、エンジン油添加剤および固体潤滑剤として有名なモリブデン系の層状化合物の解析・開発などにもこの手法を適用可能である。
【００４１】
「グラフェンシート間相互作用」
グラファイトの分子レベルにおける構造を図４に示す。熱平衡状態のＡＢスタック構造における炭素原子間距離ｌaは１．４２Åであり、シート間距離ｌzは３．３５Åである。グラフェンシート間の相互作用は、Girifalcoらによる炭素原子間のvan der Waalsポテンシャルの文献値[9]をもとにした。このポテンシャルでは、安定構造であるＡＢスタック構造も再現されており、グラフェンシート間の引きはがし応力なども再現されている。このポテンシャルを剛体間の相互作用として導入することによって、運動方向による層間相互作用の異方性やエネルギー地形を取り込むことができる。
【００４２】
「炭素原子間ポテンシャル」
原子間相互作用Uatom-atom（ｒa）は以下のLennard-Jones型の関数として表わされる。
【数１０】
2011058947.tif 000012
【００４３】
ここで、ｒaは原子間距離、Ａ，Ｂは定数である。この関数は、ポテンシャルの深さを表わす指標εと、平衡距離に関係するσとを使って書き直すことができる。
【数１１】
2011058947.tif 000013
【００４４】
εと平衡距離ｒa0は以下で表される。
【数１２】
2011058947.tif 000014
2011058947.tif 000015
【００４５】
グラフェン上の炭素原子間ポテンシャルとして、つぎの値を用いた。Ａ＝１５．２ｅＶ×Å６，Ｂ＝２４．１×１０−３ｅＶ×Å１２，ｒ0＝３．８３Å，ε＝２．８９ｍｅＶ。以上により表わされる原子間ポテンシャル式（８）を図５に示す。
【００４６】
「炭素原子−グラフェン間ポテンシャル」
つぎに、式（８）をもとに炭素原子−グラフェン間のポテンシャルを次式により決定する。
【数１３】
2011058947.tif 000016
【００４７】
ここで、シート直交座標系におけるｘをすべり方向、ｙをすべりと直角方向とし、右辺の和は平面ｚ＝０上に配置した一辺ｌaのハニカム構造の炭素原子群と、ｚ＝ｌzに配置した炭素原子との相互作用とする。シミュレーションにおいては、系の対称性を考慮して主軸ベクトルａ＝（（√３／２）ｌa，（１／２）ｌa），ｂ＝（０，ｌa）を有する座標系において、ａ，ｂ方向にそれぞれ１００分割した格子点上でUatom-sheetをあらかじめ計算し（図６）、これをシートの運動に際して直交座標系に戻しシート間ポテンシャルの計算を行う。各格子点におけるUatom-sheetは、原子近傍から遠方に差分ΔUatom-sheet＜１０−６ｅＶとなるまで積算した。直交座標系におけるUatom-sheet（ｘ，ｙ）を図７に示す。
【００４８】
「グラフェン−グラフェン間ポテンシャル」
二つの向かい合うハニカム構造に配置する原子は、周期性により幾何的に２種類に区別されるので、グラフェンシート間のポテンシャルUsheet-sheetはUatom-sheet（ｘ，ｙ）を用い、
【数１４】
2011058947.tif 000017
により計算できる。ここで、Δｘ，Δｙは、シート間のそれぞれｘ、ｙ方向の座標の差、Ｎｘ，Ｎ’ｘ、Ｎｙ，Ｎ’ｙは相対するシートのｘ，ｙ方向の原子数である。式（１３）により、有限の大きさのシートの相互作用における末端効果を取り入れることが可能となる。Ｎｘ＝Ｎ’ｘ＝１００，Ｎｙ＝Ｎ’ｙ＝１００の系におけるUsheet-sheetのポテンシャル面を、図８に、Δｘ，Δｙ軸上の断面をそれぞれ示す。Δｘ＝０，Δｙ＝０のときがＡＢ構造に対応し最安定であり、シート同士がずれることによりシート間相互作用が弱まることが示されている。
【００４９】
シートの片面が末端の存在しない高配向性グラファイト板（HOPG; highly oriented pyrolytic graphite）である場合のポテンシャルUsheet-plateは、
【数１５】
2011058947.tif 000018
により表わす。ここでｘ，ｙはシートの座標、Ｎｘ，Ｎｙはシートのｘ，ｙ方向の原子数である。この相互作用においては、図８に示されるポテンシャル面のベースラインの変化は存在せず、周期関数となる。
【００５０】
「シミュレーションモデル」
グラフェン移着片のＭＣＢＤシミュレーションのモデルを図９に示す。上下のスライダーの間にＮｚ層（Ｎｚ個の層）のグラフェンシートを挟む。それぞれのシートはスライダーのすべり方向ｘ、直角方向ｙに対してそれぞれＮｘ，Ｎｙ個の原子からなるとする。たとえば、（Ｎｘ，Ｎｙ，Ｎｚ）＝（１００，１００，２０）の系では移着片の原子数は２０万となり、全原子モデルによる分子動力学によってダイナミクスを扱うことは大変困難である。一方、本実施形態のＭＣＢＤ法では計算時間はＮｚに比例し、Ｎｘ，Ｎｙにはよらない。これは、上述した層間相互作用の計算において、先に末端効果のない原子−グラフェン間相互作用ポテンシャルUatom-sheetを決定し、その上で、末端効果を含めたグラフェン−グラフェン間相互作用ポテンシャルUsheet-sheetを決定するという手順をとったことによる。（Ｎｘ，Ｎｙ，Ｎｚ）＝（１００，１００，２０）の系における移着片サイズは、（ｘ，ｙ，ｚ）方向にそれぞれ１２．３ｎｍ、２１．３ｎｍ、６．４ｎｍであり、ＴＥＭによるＡＦＭチップ先端に付着した移着片サイズと同じスケールである。
【００５１】
ＭＣＢＤ法における物理的な時間ΔｔとＭＣステップ間の時間とは、拡散係数Ｄと式（４）によって関係づけられるが、Ｄの値は実験値、あるいはよりミクロスケールのシミュレーション、たとえば有限のグラフェンシートの運動についての全原子の分子動力学計算を実行することによって求めることができる。
【００５２】
最大変位幅はΩ−１＝０：０１Å、層間距離は極めて軽荷重条件として熱平衡値３．３５Å、すべり速度は１０−６Å／ＭＣｓｔｅｐとした。温度Ｔは３００Ｋの室温、以下に示す大半のシミュレーションにおいてシートの初期配置として、すべての層がスタックしているｘ＝０，ｙ＝０とした。
【００５３】
上述したように、各グラフェンシートを剛体とみなすと、それぞれ並進３、回転３の自由度を有する（図２右）。例として、ｘ，ｙの運動のみに着目した結果を示す。
【００５４】
「結果」
「平衡状態におけるグラファイトの熱ゆらぎ」
図１０、図１１、図１２に、各シートが１００×１００原子、２０層のグラフェンシートの系の平衡状態におけるｘ，ｙおよび系の全エネルギーのゆらぎをプロットする。第０（上スライダー），１，５，１０，１５，２０番目のシートについて示した。
【００５５】
ここで、図１０は熱平衡状態におけるシートのｘ方向の変位である。以下、図中凡例の数字は第０（上スライダー），１，５，１０，１５，２０番のシートについてプロットしていることを示す。Ｎｘ＝Ｎｙ＝１００、Ｎｚ＝２０、自由度：ｘ，ｙ、末端効果：上下あり、移着片形状：直方体、すべり速度：０とした。図１１は、熱平衡状態におけるシートのｙ方向の変位である。Ｎｘ＝Ｎｙ＝１００、Ｎｚ＝２０、自由度：ｘ，ｙ、末端効果：上下あり、移着片形状：直方体、すべり速度：０とした。図１２は熱平衡状態におけるシートの全エネルギーＵのゆらぎである。Ｎｘ＝Ｎｙ＝１００、Ｎｚ＝２０、自由度：ｘ，ｙ、末端効果：上下あり、移着片形状：直方体、すべり速度：０とした。
【００５６】
上スライダー、下スライダーの両者について末端効果を考慮しており、ポテンシャル面は図８に示したように各シートにおける平衡位置（原点）からのずれはエネルギー的に不利となっている。ｘ，ｙいずれの方向においてもゆらぎは小さく、原点付近にとどまっており、各シートは図８におけるポテンシャルの凸部を超えることはない。
【００５７】
詳細に観察すると、上下のスライダーに接している第１および２０番のゆらぎは他のスライダーと比較すると抑制されている。
【００５８】
「すべりのダイナミクス」
実験における移着片の状態を考える。移着片は鋼球と一体化して移動する上スライダーに付着していると考えられるため、上スライダーは有限の面積すなわち末端効果を含むシートと仮定するのが自然である。下スライダーは高い平坦性が期待される高配向性グラファイトであるから、式（１４）で表わされる周期的関数と置くことが相応しいと考えられる。この非対称な上下スライダーの関係を図１３に示す。
【００５９】
せん断シミュレーションの結果を図１４に示す。この図１４は非対称系におけるすべり運動下のシートのｘ方向の変位を示し、Ｎｘ＝Ｎｙ＝１００、Ｎｚ＝２０、自由度：ｘ，ｙ、末端効果：上あり、移着片形状：直方体、すべり：ありという条件である。上スライダーの移動とともに、移着片全体が追随している様子が見て取れる。
【００６０】
さらに、下スライダー（ＨＯＰＧ側）と接する第２０層の運動に着目すると、上スライダーに追随する運動とともに周期的に振動していることがわかる。
【００６１】
上記シミュレーションにおけるｙ方向の変位を図１５に、エネルギーの変化を図１６にそれぞれ示す。この図１５，１６もＮｘ＝Ｎｙ＝１００、Ｎｚ＝２０、自由度：ｘ，ｙ、末端効果：上あり、移着片形状：直方体、すべり：ありという条件である。ｙ方向においても、第２０層が周期的に振動していることが読み取れる。また、エネルギーは第２０層の２倍の周期で振動しており、摩擦発現が移着片の最下層とＨＯＰＧ面との間の相互作用と密接に関係していることが示唆される。
【００６２】
第２０層のｘ，ｙ平面における変位を図１７に示す。この図１７もＮｘ＝Ｎｙ＝１００、Ｎｚ＝２０、自由度：ｘ，ｙ、末端効果：上あり、移着片形状：直方体、すべり：ありという条件である。図８のポテンシャル面と比較すると、等ポテンシャル面に沿うように第２０層が運動していることがわかる。
【００６３】
以上より、鋼球−移着片−ＨＯＰＧ面を模した非対称の系において、移着片は鋼球とともに運動し、ＨＯＰＧの等ポテンシャル面に沿うような運動によりすべっているということがわかる。
【００６４】
移着片の視点に立ち考察すると、移着片は下スライダーに対してはｘ軸のどちらの方向に進んでもエネルギー的な損得はないが、上スライダーは末端効果を有するため、有限である上スライダーに追随した方がエネルギー的に有利である。スライダーが有限の大きさを有することが重要であるということがわかる。
【００６５】
以上のように、ＭＣＢＤ計算により、様々な条件における層間すべりの様子を捉えることができる。下地は周期ポテンシャルで表わされるグラフェンとし、移着片の上部は一定速度で移動するスライダーとした。一辺の大きさが１００原子（１２ｎｍ）程度の移着片においては、移着片はスライダーとともに移動し、その最下層と下地との間で摩擦が生じる。
【００６６】
層状化合物のすべり摩擦においては、このように現実的なすべり速度、温度におけるシート間のエネルギー地形が非常に重要であることがわかる。
【００６７】
本実施形態では、ｘ，ｙの自由度について挙げたが、シミュレータをｚ方向およびｙａｗ角の自由度を取り入れることも当然可能である。同時に拡散係数Ｄの全原子分子動力学による決定も行うことにより、決定することができる。
【００６８】
本実施形態のシミュレータを活用することにより、低摩擦発現に効果のある適切な移着片サイズなどが判明すると、その直接観察や、積極的に移着を促すような表面性状の改良などが考えられる。またグラファイト系にかぎらず、図３で示したマルチスケール解析スキームを用いることによって、新規低摩擦材料の探索を行うことも可能である。
【００６９】
「処理手順および効果」
図１８には、本実施形態における計算の全体的な処理の流れが示されている。まず、層状化合物（層状分子）の候補を選択する（Ｓ１）。層状化合物には、上述のように、グラファイト、雲母、二硫化モリブデンなど各種のものがあり、計算の対象となる物質を特定する。
【００７０】
対象となった物質の文献値、実験値または量子計算の結果を基に、層間相互作用を確定し入力する（Ｓ２）。すなわち、層状化合物が特定されれば、その層状化合物の原子配列など基本的な構成は文献等から得られる。従って、文献値、実験値などから層間の相互ポテンシャルなどの相互作用を決定する。また、相互ポテンシャルなどの相互作用については、そのまま文献値などを使えない場合も多く、基本的な層化合物の構成に基づき、上述のような層対層（シート対シート）の相互作用を決定することが好適である。例えば、式（１３）の値などが決定される。
【００７１】
特に、本実施形態によれば、シート対シートの相互作用（相互ポテンシャル）を計算する前に、末端効果を考慮せずに原子対シートの相互作用を計算し、その結果に基づき末端効果を含めシート対シートの相互作用を計算する。従って、末端効果を有効としたシート対シートの計算が容易になる。そして、このような計算によって、層間の相互作用を１つのシートの特性として規定することができる。
【００７２】
次に、文献値、実験値または原子レベルの計算結果を基に、各層の内部振動に起因する作用を示すパラメータ（拡散係数や揺動力など）を確定し入力する（Ｓ３）。拡散係数Ｄなどについても、上述のように基本的なデータを基づき計算によって、決定することができる。また文献などによって、直接決定できる場合には、それを利用する。例えば、式（７）の値がこれによって決定される。
【００７３】
このように、拡散係数なども、層の特性として記述できるため、各層を層間の相互ポテンシャル、拡散係数を特定した剛体シートとして粗視化することができる。
【００７４】
各層の初期配置、外場などの初期条件を入力する（Ｓ４）。「シミュレーションモデル」の項において説明したように、最大変位幅、層間距離、温度、シートの配置などが決定される。
【００７５】
そして、構造計算あるいは動力学計算を実行する（Ｓ５）。例えば、「平衡状態におけるグラファイトの熱揺らぎ」、「滑りのダイナミクス」の項で示したような計算が実行される。
【００７６】
このようなシミュレーションは、粗視化された剛体シートについての力を印加したり上部のものを移動させたりして行う。従って、計算を比較的簡単なものとして、かつ相互ポテンシャルや、拡散係数については原子レベルの影響を考慮した計算を行うことができる。
【００７７】
このような計算が終了した場合には、実行結果が出力される（Ｓ６）。すなわち、対象となった層状化合物が各種の条件でどのような挙動を示すかシミュレートされ、この結果により、対象となる層状化合物の評価や、所定の条件に適切な層状化合物の構造などを把握できる。
【００７８】
従来、層状化合物片が一層のシートにつきＮｘｙ個あり、全体がＮｚ層により構成される場合、全原子模型による分子シミュレーションではＮｘｙ×Ｎｚ個の粒子を扱う必要があり膨大な計算時間を必要とした。一例として、原子間力顕微鏡上のチップに付着したグラファイト移着片は、Ｎｘｙが数万、Ｎｚが数十であり、ひとつの物理量を計算するために大型計算機により数ヶ月の計算を必要とした。本実施形態によれば、層状化合物の各層を構成する原子の振動の影響は各層単位での相互作用および拡散係数に繰り込むことができるという特徴に基づき内部自由度を限定することで、取り扱う粒子数はＮｚ個に減少しＮｘｙ倍計算が高速化される。
【００７９】
その上で、層間相互作用などは微視的な計算に基づくため、各化合物ごとの特性を再現することができる。これにより新規層状化合物の分子設計を行うことができる。さらに本手法では、層の膜厚方向の拡散係数は小さく、これに直角方向の拡散係数は大きい、というように系の異方性を取り入れるため、層状化合物の特性を反映した計算を実現できる。
【００８０】
「参考文献」
上述の説明において、引用した文献を以下に示す。
[1] A.K. Geim and K.S. Novoselov. aaa. Nature Materials, Vol. 6, p. 183, 2007.
[2] M. Dienwiebel et al. aaa. Surface Science, Vol. 576, p. 197, 2005.
[3] M. Tsukada N. Sasaki, K. Kobayashi. Atomic-scale friction image of graphite in atomic-force microscopy. Phys. Rev. B., Vol. 54, pp. 2138{2149, 1996.
[4] K. Kikuchi, M. Yoshida, T. Maekawa, and H. Watanabe. Metropolis monte carlo method as a numerical technique to solve the fokker-planck equation. Chem. Phys. Lett., Vol. 185, pp. 335{338, 1991.
[5] N. G. van Kampen. Stochastic processes in physics and chemistry, revised and enlarged edition, chapter 11. North-Holland, Amsterdam, 1992.
[6] D. L. Ermak and J. A. McCammon. Brownian dynamics with hydrodynamic interac-tions. J. Chem. Phys., Vol. 69, pp. 1352{1360, 1978.
[7] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller. Equation of state calculations by fast computing machines. J. Chem. Phys., Vol. 21, pp. 1087{1092, 1953.
[8] M. Yoshida and K. Kikuchi. Metropolis monte carlo brownian dynamics simulation of the ion atmosphere polarization around a rodlike polyion. J. Phys. Chem., Vol. 98, pp. 10303{10306, 1994.
[9] M. Hodak L.A. Girifalco and R. S. Lee. aaa. Phys. Rev. B., Vol. 62, p. 13104, 2000.
【図５】
2011058947.tif 000019
【図１２】
2011058947.tif 000020
【図１５】
2011058947.tif 000021
【図１６】
2011058947.tif 000022
【図１７】
2011058947.tif 000023
【図１８】
2011058947.tif 000024
【図１】
2011058947.tif 000025
【図２】
2011058947.tif 000026
【図３】
2011058947.tif 000027
【図４】
2011058947.tif 000028
【図６】
2011058947.tif 000029
【図７】
2011058947.tif 000030
【図８】
2011058947.tif 000031
【図９】
2011058947.tif 000032
【図１０】
2011058947.tif 000033
【図１１】
2011058947.tif 000034
【図１３】
2011058947.tif 000035
【図１４】
2011058947.tif 000036