JP7456260B2

JP7456260B2 - 高分子材料のシミュレーション方法

Info

Publication number: JP7456260B2
Application number: JP2020074289A
Authority: JP
Inventors: 真一上野
Original assignee: Sumitomo Rubber Industries Ltd
Current assignee: Sumitomo Rubber Industries Ltd
Priority date: 2020-04-17
Filing date: 2020-04-17
Publication date: 2024-03-27
Anticipated expiration: 2040-04-17
Also published as: JP2021174021A

Description

本発明は、高分子材料のシミュレーション方法に関する。

下記特許文献１には、高分子材料のシミュレーション方法が記載されている。この方法は、高分子材料の分子構造をモデル化した分子構造モデルを設定するステップと、分子構造モデルをセルに配置して、構造緩和計算を行うステップとを含む。

特開２０１３－１９５２２０号公報

一般に、複数の分子構造モデルは、セルの内部にランダムに初期配置されるため、隣接する分子構造モデル間の距離が必要以上に近くなることがある。また、前記初期配置において、分子構造モデルが不自然な形状になる場合もある。このような初期配置を持った分子構造モデルの運動を進展させると、構造緩和計算が不安定になり、いわゆる計算落ちが多発するという問題があった。

本発明は、以上のような実状に鑑み案出されたもので、構造緩和を安定して計算することができる高分子材料のシミュレーション方法を提供することを主たる目的としている。

本発明は、高分子材料のシミュレーション方法であって、前記高分子材料の分子鎖に基づいて、数値計算用の分子鎖モデルをコンピュータに入力する工程と、前記高分子材料の一部に対応する仮想空間であるセルを、前記コンピュータに定義する工程とを含み、前記コンピュータが、前記セルの内部に、複数の前記分子鎖モデルを配置する工程と、複数の前記分子鎖モデルを対象に、分子動力学計算に基づく構造緩和を計算する緩和工程とを含み、前記緩和工程は、予め定められた第１の時間刻みで、前記分子鎖モデルの運動を進展させる第１計算工程と、前記第１計算工程の後に、前記第１の時間刻みよりも大きな第２の時間刻みで、前記分子鎖モデルの運動を進展させる第２計算工程とを含むことを特徴とする。

本発明に係る前記高分子材料のシミュレーション方法において、前記第１の時間刻みは、前記第２の時間刻みの０．１～０．５倍であってもよい。

本発明に係る前記高分子材料のシミュレーション方法において、前記第１の時間刻みは、０．０５～０．２ｆｓであり、前記第２の時間刻みは、０．３～１．２ｆｓであってもよい。

本発明に係る前記高分子材料のシミュレーション方法において、前記第１計算工程での繰り返し計算数、及び、前記第２計算工程での繰り返し計算数は、１００～１００００００回であってもよい。

本発明に係る前記高分子材料のシミュレーション方法において、前記緩和工程は、前記第１計算工程に先立ち、複数の前記分子鎖モデルを対象に、分子力学計算を行う工程をさらに含んでもよい。

本発明に係る前記高分子材料のシミュレーション方法において、前記緩和工程は、前記第２計算工程の後に、前記第２の時間刻みよりも大きな第３の時間刻みで、前記分子鎖モデルの運動を進展させる第３計算工程をさらに含んでもよい。

本発明の高分子材料のシミュレーション方法は、上記の工程を採用することにより、構造緩和を安定して計算することができる。

高分子材料のシミュレーション方法を実行するためのコンピュータの一例を示す斜視図である。高分子材料のシミュレーション方法の処理手順の一例を示すフローチャートである。分子鎖モデルの一例を示す概念図である。分子鎖モデルが配置されたセルの一例を示す概念図である。高分子材料モデルの一例を示す概念図である。緩和工程の処理手順の一例を示すフローチャートである。本発明の他の実施形態の緩和工程の処理手順の一例を示すフローチャートである。本発明のさらに他の実施形態の緩和工程の処理手順の一例を示すフローチャートである。

以下、本発明の実施の一形態が図面に基づき説明される。なお、各図面は、発明の内容の理解を高めるためのものであり、誇張された表示が含まれる他、各図面間において、縮尺等は厳密に一致していない点が予め指摘される。

本実施形態の高分子材料のシミュレーション方法（以下、単に「シミュレーション方法」ということがある。）は、複数の分子鎖モデルを対象に、分子動力学計算に基づく構造緩和が計算される。図１は、高分子材料のシミュレーション方法を実行するためのコンピュータの一例を示す斜視図である。

コンピュータ１は、本体１ａ、キーボード１ｂ、マウス１ｃ及びディスプレイ装置１ｄを含んで構成されている。この本体１ａには、例えば、演算処理装置（ＣＰＵ）、ＲＯＭ、作業用メモリ、磁気ディスクなどの記憶装置、及び、ディスクドライブ装置１ａ１、１ａ２が設けられている。記憶装置には、本実施形態のシミュレーション方法を実行するためのソフトウェア等が予め記憶されている。

本実施形態の高分子材料には、少なくとも１種類の分子鎖（ポリマー）が含まれている。なお、高分子材料には、フィラーやカップリング剤がさらに含まれてもよい。本実施形態の分子鎖は、スチレンブタジエンゴムである場合が例示されるが、特に限定されるわけではない。

図２は、高分子材料のシミュレーション方法の処理手順の一例を示すフローチャートである。本実施形態のシミュレーション方法では、先ず、高分子材料の分子鎖に基づいて、数値計算用の分子鎖モデルが、コンピュータ１に入力される（工程Ｓ１）。図３は、分子鎖モデル３の一例を示す概念図である。

本実施形態の分子鎖モデル３は、分子鎖（本例では、スチレンブタジエンゴム）の実際の構造に基づいた全原子モデルである場合が例示されるが、特に限定されない。分子鎖モデル３は、例えば、ユナイテッドアトムモデル（図示省略）や、粗視化分子モデル（Kremer-GrestモデルやDPDモデルなど）であってもよい。

本実施形態の分子鎖モデル３は、分子鎖を構成する原子をモデリングした粒子モデル４と、粒子モデル４、４間を結合するボンド（結合鎖モデル）５とを含んで構成されている。

粒子モデル４は、後述の分子動力学計算に基づいたシミュレーションにおいて、運動方程式の質点として取り扱われる。粒子モデル４には、例えば、質量、直径、電荷、又は、初期座標などのパラメータが定義される。本実施形態の粒子モデル４は、分子鎖の炭素原子をモデリングした炭素粒子モデル４Ｃと、分子鎖の水素原子をモデリングした水素粒子モデル４Ｈとが含まれる。

ボンド５は、粒子モデル４、４間を拘束するものである。本実施形態のボンド５は、炭素粒子モデル４Ｃ、４Ｃを連結する主鎖５Ａ、及び、炭素粒子モデル４Ｃと水素粒子モデル４Ｈとの間を連結する側鎖５Ｂとを含んでいる。

分子鎖モデル３には、各粒子モデル４、４間の結合長さである結合長や、ボンド５を介して連続する３つの粒子モデル４がなす角度である結合角が定義される。さらに、分子鎖モデル３には、ボンド５を介して連続する４つの粒子モデル４において、隣り合う３つの粒子モデル４が作る二面角などが定義される。これにより、分子鎖モデル３の三次元構造が定義される。このような分子鎖モデル３は、慣例に従い、外力又は内力を受けることによって、結合長、結合角及び二面角が変化する。これにより、分子鎖モデル３の三次元構造が変化する。

結合長、結合角及び二面角は、適宜定義されうる。本実施形態の結合長、結合角及び二面角を定義するためのポテンシャルの各定数は、論文１（Huai Sun等著、「COMPASS II: extended coverage for polymer and drug-like molecule databases」、J. Mol. Model. vol. 22、no. 2、article 47、2016年、pp.1-10）のCOMPASS II力場に基づいて定義される。

上記のような分子鎖モデル３は、例えば、ダッソー・システムズ社製のシミュレーションソフトウェア（Materials Studio（「Materials Studio」は登録商標））等が用いられることで、容易に作成されうる。分子鎖モデル３は、コンピュータ１に記憶される。

次に、本実施形態のシミュレーション方法では、高分子材料の一部に対応する仮想空間であるセルが、コンピュータ１に定義される（工程Ｓ２）。図４は、分子鎖モデル３が配置されたセル６の一例を示す概念図である。

本実施形態のセル６は、少なくとも互いに向き合う一対の面７、７、本実施形態では、互いに向き合う三対の面７、７を有している。本実施形態のセル６は、直方体又は立方体（本実施形態では、立方体）として定義されている。各面７、７には、周期境界条件が定義されている。セル６の一辺の各長さＬ１は、例えば、特許文献（特開２０１９－０３２２７８号公報）の記載に基づいて、適宜設定されうる。セル６は、コンピュータ１に記憶される。

次に、本実施形態のシミュレーション方法では、コンピュータ１が、セル６の内部に、複数の分子鎖モデル３を配置する（工程Ｓ３）。本実施形態の工程Ｓ３では、複数の分子鎖モデル３がランダムに配置される。

分子鎖モデル３をランダムに配置する手順は、特に限定されない。本実施形態の工程Ｓ３では、例えば、ＤＢＭＣ法（Density Biased Monte Carlo）に基づいて、分子鎖モデル３がランダムに配置される。このような分子鎖モデル３の初期配置は、例えば、上記のシミュレーションソフトウェア等が用いられることで、容易に求めることができる。

セル６の内部に配置される分子鎖モデル３の個数は、例えば、コンピュータ１の計算能力や、分子鎖モデル３の大きさ等に基づいて、適宜設定されうる。分子鎖モデル３が配置されたセル６は、コンピュータ１に記憶される。

次に、本実施形態のシミュレーション方法では、コンピュータ１が、セル６内に配置された分子鎖モデル３、３間に、力場Ｆ１を定義する（工程Ｓ４）。本実施形態の工程Ｓ４では、隣接する分子鎖モデル３、３について、それらの粒子モデル４、４間に、力場Ｆ１がそれぞれ定義される。これらの力場Ｆ１については、適宜定義することができる。本実施形態の力場Ｆ１は、例えば、上述のCOMPASS II 力場に基づいて定義される。このような力場Ｆ１の定義は、例えば、上記のシミュレーションソフトウェア等が用いられうる。力場Ｆ１は、コンピュータ１に記憶される。

次に、本実施形態のシミュレーション方法では、コンピュータ１が、複数の分子鎖モデル３を対象に、分子動力学計算（ Molecular Dynamics ：ＭＤ）に基づく構造緩和を計算する（緩和工程Ｓ５）。緩和工程Ｓ５では、例えば、セル６について所定の時間、分子鎖モデル３が古典力学に従うものとして、ニュートンの運動方程式が適用される。そして、緩和工程Ｓ５では、予め定められた時間刻み（シミュレーションの単位時間（微小時間））Δｔで、分子鎖モデル３の運動が進展される。

本実施形態の構造緩和の計算は、セル６において、圧力及び温度が一定、又は、体積及び温度が一定に保たれる。これにより、緩和工程Ｓ５では、実際の高分子材料の分子運動に近似させて、分子鎖モデル３の初期配置が精度よく緩和されうる。

緩和工程Ｓ５では、分子鎖モデル３の初期配置が十分に緩和されるまで、構造緩和が計算される。これにより、緩和工程Ｓ５では、分子鎖モデル３が平衡状態（構造が緩和した状態）となった高分子材料モデル１０が作成されうる。このような緩和計算は、例えば、上記のシミュレーションソフトウェアが用いられうる。図５は、高分子材料モデル１０の一例を示す概念図である。

ところで、図４に示されるように、分子鎖モデル３は、セル６の内部にランダムに初期配置されるため、隣接する分子鎖モデル３、３間の距離（例えば、粒子モデル４、４の距離）Ｄ１が必要以上に近くなることがある。また、初期配置において、分子鎖モデル３が不自然な形状になる場合もある。このような初期配置を持った分子鎖モデル３の運動を進展させると、例えば、隣接する粒子モデル４、４間でエネルギーが発散して、構造緩和計算が不安定になる場合がある。このため、構造緩和の計算の異常終了（いわゆる計算落ち）が多発するという問題がある。

本実施形態の緩和工程Ｓ５では、時間刻みが段階的に大きくされて、分子鎖モデル３の運動が徐々に進展される。図６は、緩和工程Ｓ５の処理手順の一例を示すフローチャートである。

本実施形態の緩和工程Ｓ５は、予め定められた第１の時間刻みで、分子鎖モデル３（図４に示す）の運動を進展させる第１計算工程Ｓ５１と、第１の時間刻みよりも大きな第２の時間刻みで、分子鎖モデル３の運動を進展させる第２計算工程Ｓ５２とを含んでいる。第２計算工程Ｓ５２は、第１計算工程Ｓ５１の後に行われる。緩和工程Ｓ５では、それぞれの第１計算工程Ｓ５１及び第２計算工程Ｓ５２において、分子鎖モデル３の運動の進展（構造緩和）が、時間刻み（第１の時間刻み、又は、第２の時間刻み）ごとに繰り返し計算される。

第１計算工程Ｓ５１では、第２の時間刻みよりも小さな第１の時間刻みで、分子鎖モデル３（図４に示す）の運動が進展される。このため、第１計算工程Ｓ５１では、例えば、隣接する粒子モデル４、４間の距離Ｄ１（図４に示す）が必要以上に近い場合や、分子鎖モデル３が不自然な形状になった場合でも、分子鎖モデル３の構造を少しずつ緩和することができる。

一方、第２計算工程Ｓ５２では、第１計算工程Ｓ５１で構造が緩和された分子鎖モデル３を対象に、構造緩和が計算されるため、第１の時間刻みよりも大きな第２の時間刻みで、分子鎖モデル３の運動が進展されたとしても、計算落ちの発生を防ぐことができる。さらに、第２計算工程Ｓ５２では、相対的に大きな第２の時間刻みで分子鎖モデル３の運動が進展されるため、例えば、第１の時間刻みのみで分子鎖モデル３の運動が進展される場合に比べて、計算コストが増大するのを防ぐことができる。

このように、本実施形態の緩和工程Ｓ５では、時間刻みを段階的に大きくしながら、分子鎖モデル３の運動を徐々に進展させることができるため、計算コストの増大を防ぎつつ、上記のようなエネルギーの発散を防ぐことができる。したがって、本実施形態の緩和工程Ｓ５では、構造緩和を安定して計算することができる。

このような作用を効果的に発揮させるために、第１の時間刻みは、第２の時間刻みの０．１～０．５倍が望ましい。第１の時間刻みが、第２の時間刻みの０．５倍以下に設定されることにより、第２計算工程Ｓ５２に比べて、第１計算工程Ｓ５１での分子鎖モデル３の運動を、緩やかに進展させることができる。これにより、緩和工程Ｓ５では、上記のようなエネルギーの発散を防ぎつつ、安定した構造緩和を計算することができる。一方、第１の時間刻みが、第２の時間刻みの０．１倍以上に設定されることにより、構造緩和の計算時間（計算コスト）が必要以上に大きくなるのを防ぐことができる。このような観点より、第１の時間刻みは、第２の時間刻みの０．４倍以下が望ましく、また、０．２倍以上が望ましい。

第１の時間刻みは、第２の時間刻みよりも小さければ、適宜設定されうる。上記のような作用を効果的に発揮させるためには、第１の時間刻みが０．０５～０．２ｆｓに設定されるのが望ましい。第１の時間刻みが０．２ｆｓ以下に設定されることにより、従来の構造緩和計算に用いられている一般的な時間刻み（例えば、１．０～２．０ｆｓ）よりも確実に小さく設定されうる。これにより、第１計算工程Ｓ５１では、上記のようなエネルギーの発散を確実に防ぎつつ、安定した構造緩和を計算することができる。一方、第１の時間刻みが０．０５ｆｓ以上に設定されることにより、構造緩和の計算時間（計算コスト）が必要以上に大きくなるのを防ぐことができる。このような観点より、第１の時間刻みは、０．１５ｆｓ以下が望ましく、また、０．１ｆｓ以上が望ましい。

第２の時間刻みは、第１の時間刻みよりも大きければ、適宜設定されうる。上記のような作用を効果的に発揮させるためには、第２の時間刻みが０．３～１．２ｆｓに設定されるのが望ましい。第２の時間刻みが１．２ｆｓ以下に設定されることにより、従来の構造緩和計算に用いられている一般的な時間刻みよりも大きくなるのを防ぐことができるため、上記のようなエネルギーの発散が発生するのを防ぐことができる。一方、第２の時間刻みが０．３ｆｓ以上に設定されることにより、第１計算工程Ｓ５１に比べて、分子鎖モデル３の運動を早く進展させることができるため、構造緩和の計算時間（計算コスト）が必要以上に大きくなるのを防ぐことができる。このような観点より、第２の時間刻みは、１．０ｆｓ以下が望ましく、また、０．５ｆｓ以上が望ましい。

第１計算工程Ｓ５１での繰り返し計算数、及び、第２計算工程Ｓ５２での繰り返し計算数は、分子鎖モデル３の運動を進展させて、構造緩和を計算することができれば、適宜設定されうる。上記のような作用を効果的に発揮させるためには、第１計算工程Ｓ５１及び第２計算工程Ｓ５２での各繰り返し計算数が１００～１００００００回に設定されるのが望ましい。それぞれの繰り返し計算数が１００回以上に設定されることにより、上記のようなエネルギーの発散を防ぎつつ、分子鎖モデル３の運動をより確実に進展させることができる。一方、それぞれの繰り返し計算数が１００００００回以下に設定されることにより、構造緩和の計算時間（計算コスト）が必要以上に大きくなるのを防ぐことができる。このような観点より、繰り返し計算数は、好ましくは５００回以上、さらに好ましくは１０００回以上であり、また、また、好ましくは５０００００回以下であり、さらに好ましくは１０００００回以下である。

なお、初期配置での分子鎖モデル３、３間の距離や、分子鎖モデル３の不自然な形状によっては、上記の繰り返し計算数に基づいて、第１計算工程Ｓ５１が実行されたとしても、次の第２計算工程Ｓ５２において計算落ちが発生する可能性が僅かに存在する。このため、第１計算工程Ｓ５１では、例えば、全ての分子鎖モデル３について、粒子モデル４、４間の距離Ｄ１が、予め定められた閾値（例えば、１～２Å）以上となるまで、分子鎖モデル３の運動の進展が計算されてもよい。これにより、第１計算工程Ｓ５１の後に実施される第２計算工程Ｓ５２では、隣接する粒子モデル４、４間において、エネルギーが発散するのを抑制することができるため、計算落ちの発生をより確実に防ぐことが可能となる。

本実施形態の第２計算工程Ｓ５２では、分子鎖モデル３の初期配置が十分に緩和されるまで、分子鎖モデル３の運動の進展が、第２の時間刻みで繰り返し計算される。これにより、緩和工程Ｓ５では、構造緩和の計算の異常終了（計算落ち）を防ぎつつ、分子鎖モデル３の平衡状態（構造が緩和した状態）が計算された高分子材料モデル１０（図５に示す）が、確実に作成されうる。高分子材料モデル１０は、コンピュータ１に記憶される。

次に、本実施形態のシミュレーション方法では、コンピュータ１が、高分子材料モデル１０を変形させて、高分子材料モデル１０の物理量を計算する（工程Ｓ６）。本実施形態の工程Ｓ６では、高分子材料モデル１０を、予め定められた方向に引っ張る単軸引張試験が計算される。高分子材料モデル１０を引っ張る方向は、適宜選択されうる。

高分子材料モデル１０の変形計算は、例えば、特許文献（特開２０１６－０８１２９７号公報）に記載の手順に基づいて、高分子材料モデル１０の一端側の面７ａ、及び、他端側の面７ｂが互いに離間するように、高分子材料モデル１０の伸長が計算される。高分子材料モデル１０の伸長計算は、分子動力学計算に基づいて計算される。

高分子材料モデル１０の物理量は、適宜求められる。物理量には、高分子材料モデル１０の応力や、セル６に形成された空孔の大きさ等が含まれてもよい。物理量は、コンピュータ１に記憶される。

次に、本実施形態のシミュレーション方法では、コンピュータ１が、高分子材料モデル１０の物理量が良好か否かを判断する（工程Ｓ７）。物理量が良好か否かの判断は、適宜行われうる。本実施形態の工程Ｓ７では、工程Ｓ６で計算された物理量に対して設定された閾値に基づいて、その物理量が、良好か否かが判断される。

工程Ｓ７において、物理量が良好であると判断された場合（工程Ｓ７で、「Ｙ」）、高分子材料モデル１０（図５に示す）に基づいて、高分子材料（図示省略）が製造される（工程Ｓ８）。一方、工程Ｓ７において、物理量が良好でないと判断された場合（工程Ｓ７で、「Ｎ」）、高分子材料モデル１０（図３に示した分子鎖モデル３等）の諸条件が変更されて（工程Ｓ９）、工程Ｓ１～工程Ｓ７が再度実施される。これにより、本実施形態のシミュレーション方法では、良好な物性を有する高分子材料を、効率よく設計及び製造することができる。

これまでの実施形態の緩和工程Ｓ５では、図６に示されるように、第１計算工程Ｓ５１と、第２計算工程Ｓ５２とを含む態様が例示されたが、このような態様に限定されない。緩和工程Ｓ５では、第１計算工程Ｓ５１に先立ち、複数の分子鎖モデル３（図４に示す）を対象に、分子力学計算（Molecular Mechanics ：ＭＭ）を行う工程Ｓ５０がさらに含まれてもよい。図７は、本発明の他の実施形態の緩和工程Ｓ５の処理手順の一例を示すフローチャートである。この実施形態において、これまでの実施形態と同一の構成については、同一の符号を付し、説明を省略することがある。

この実施形態の工程Ｓ５０では、初期配置された複数の分子鎖モデル３（図４に示す）に対する分子力学計算により、粒子モデル４、４間に定義された力場Ｆ１（相互作用のポテンシャルエネルギー）に基づいて、分子鎖モデル３の分子構造が最適化される。これにより、工程Ｓ５０では、分子鎖モデル３の初期配置において、隣接する粒子モデル４、４の距離Ｄ１が必要以上に近くなったとしても、分子鎖モデル３の運動を進展させることなく、それらの粒子モデル４、４を離間させることができる。さらに、工程Ｓ５０では、分子力学計算により、初期配置で不自然な形状になった分子鎖モデル３を、自然な形状に変化させることができる。このような分子力学計算は、上記のシミュレーションソフトウェアが用いられることで、容易に計算することができる。

このように、この実施形態の緩和工程Ｓ５では、第１計算工程Ｓ５１の前に実施される工程Ｓ５０において、分子鎖モデル３の分子構造が最適化される。これにより、緩和工程Ｓ５では、工程Ｓ５０の後に実施される第１計算工程Ｓ５１及び第２計算工程Ｓ５２において、分子鎖モデル３の運動を、より安定して進展させることが可能となる。

これまでの実施形態の緩和工程Ｓ５では、図６及び図７に示されるように、第１計算工程Ｓ５１及び第２計算工程Ｓ５２が含まれることにより、時間刻みが２段階で大きくされたが、例えば、時間刻みが、３段階以上で大きくされてもよい。図８は、本発明のさらに他の実施形態の緩和工程Ｓ５の処理手順の一例を示すフローチャートである。この実施形態において、これまでの実施形態と同一の構成については、同一の符号を付し、説明を省略することがある。

この実施形態の緩和工程Ｓ５では、第２計算工程Ｓ５２の後に、第２の時間刻みよりも大きな第３の時間刻みで、分子鎖モデル３（図４に示す）の運動を進展させる第３計算工程Ｓ５３がさらに含まれる。したがって、この実施形態の緩和工程Ｓ５では、第１計算工程Ｓ５１～第３計算工程Ｓ５３のそれぞれにおいて、分子鎖モデル３の運動の進展が、時間刻み（第１の時間刻み、第２の時間刻み、又は、第３の時間刻み）ごとに繰り返し計算される。

このように、この実施形態の緩和工程Ｓ５では、これまでの実施形態の緩和工程Ｓ５に比べて、より多くの段階（本例では、３段階）で時間刻みを大きくしながら、分子鎖モデル３の運動を徐々に進展させることができる。これにより、緩和工程Ｓ５では、これまでの実施形態の緩和工程Ｓ５に比べて、分子鎖モデル３の構造を、より少しずつ緩和することができるため、上記のエネルギーの発散を効果的に防ぐことができる。

この実施形態において、第１の時間刻みは、これまでの第１の時間刻みと同一範囲に設定されるのが望ましい。

一方、第２の時間刻みは、第１の時間刻みよりも大であり、かつ、第３の時間刻みよりも小であれば、適宜設定されうる。この実施形態の第２の時間刻みは、これまでの第２の時間刻み（０．３～１．０ｆｓ）よりも小さく設定されるのが望ましい。これにより、第２計算工程Ｓ５２では、これまでの実施形態の第２計算工程Ｓ５２に比べて、分子鎖モデル３の構造を少しずつ緩和することが可能となる。

上記のような作用を効果的に発揮させるために、第２の時間刻みは、好ましくは、０．２５～０．７０ｆｓに設定されるのが望ましい。第２の時間刻みが０．７０ｆｓ以下に設定されることにより、分子鎖モデル３の構造を少しずつ緩和することが可能となる。一方、第２の時間刻みが０．２５ｆｓ以上に設定されることにより、第１計算工程Ｓ５１に比べて、分子鎖モデル３の運動を早く進展させることができる。このような観点より、第２の時間刻みは、０．３５ｆｓ以上が望ましく、また、０．６０ｆｓ以下に設定されるのが望ましい。

第３の時間刻みは、第２の時間刻みよりも大きければ、適宜設定されうる。上記のような作用を効果的に発揮させるためには、第３の時間刻みは、０．８～１．２ｆｓに設定されるのが望ましい。第３時間刻みが１．２ｆｓ以下に設定されることにより、従来の構造緩和に用いられている一般的な時間刻みよりも大きくなるのを防ぐことができる。一方、第３の時間刻みが０．８ｆｓ以上に設定されることにより、第１計算工程Ｓ５１及び第２計算工程Ｓ５２に比べて、分子鎖モデル３の運動を早く進展させることができ、構造緩和の計算時間（計算コスト）が必要以上に大きくなるのを防ぐことができる。このような観点より、第３の時間刻みは、１．１ｆｓ以下が望ましく、また、０．９ｆｓ以上が望ましい。

第１計算工程Ｓ５１での繰り返し計算数、第２計算工程Ｓ５２での繰り返し計算数、及び、第３計算工程Ｓ５３での繰り返し計算数は、分子鎖モデル３の運動を進展させて、構造緩和を計算することができれば、適宜設定されうる。これらの繰り返し計算数は、これまでの実施形態と同一範囲に設定されるのが望ましい。

この実施形態の第２計算工程Ｓ５２では、例えば、全ての分子鎖モデル３について、粒子モデル４、４間の距離Ｄ１が、予め定められた閾値（例えば、１～２Å）以上となるまで、分子鎖モデル３の運動の進展が計算されてもよい。これにより、第２計算工程Ｓ５２の後に実施される第３計算工程Ｓ５３では、隣接する粒子モデル４、４間において、エネルギーが発散するのを防ぐことができ、計算落ちの発生をより確実に防ぐことが可能となる。

この実施形態の緩和工程Ｓ５には、第１計算工程Ｓ５１、第２計算工程Ｓ５２、及び、第３計算工程が含まれたが、このような態様に限定されない。緩和工程Ｓ５には、例えば、第３計算工程Ｓ５３の後に、第３の時間刻みよりも大きな第４の時間刻みで、分子鎖モデル３の運動を進展させる第４計算工程（図示省略）が含まれてもよい。さらに、第４計算工程の後に、第４の時間刻みよりも大きな第５の時間刻みで、分子鎖モデル３の運動を進展させる第５計算工程（図示省略）など、複数の計算工程が含まれてもよい。これにより、緩和工程Ｓ５では、より多くの段階で時間刻みを大きくしながら、分子鎖モデル３の運動を徐々に進展させることができるため、構造緩和をさらに安定して計算することができる。

これまでの実施形態の緩和工程Ｓ５において、分子鎖モデル３のみを対象に、構造緩和が計算されたが、このような態様に限定されない。セル６の内部に、例えば、フィラーをモデリングしたフィラーモデル（図示省略）や、カップリング剤をモデリングしたカップリングモデル（図示省略）が含まれる場合には、これらのモデルとともに、構造緩和が計算されてもよい。

以上、本発明の特に好ましい実施形態について詳述したが、本発明は図示の実施形態に限定されることなく、種々の態様に変形して実施しうる。

図２に示した処理手順に基づいて、セルの内部に配置された複数の前記分子鎖モデルを対象に、分子動力学計算に基づく構造緩和を計算する緩和工程が実施された（実施例１～４）。

実施例１～４の緩和工程では、図６に示した処理手順に基づいて、第１の時間刻みで、分子鎖モデルの運動を進展させる第１計算工程と、第１の時間刻みよりも大きな第２の時間刻みで、分子鎖モデルの運動を進展させる第２計算工程とが実施された。緩和工程では、下記の圧力に設定することで、下記の初期密度から緩和計算後の下記の密度に近づけられた。

実施例３の緩和工程では、図７に示した処理手順に基づいて、第１計算工程に先立ち、複数の分子鎖モデルを対象に、分子力学計算を行う工程がさらに実施された。一方、実施例４の緩和工程では、図８に示した処理手順に基づいて、第２計算工程の後に、第３の時間刻みで、分子鎖モデルの運動を進展させる第３計算工程がさらに実施された。

比較のために、一つの時間刻み（１．０ｆｓ）で、分子鎖モデルの運動を進展させる緩和工程が実施された（従来例）。そして、実施例１～４及び従来例において、１０回の緩和計算を実施したときの正常終了した回数がカウントされた。なお、正常終了した回数が、５回以上であれば、良好である。

さらに、正常終了したときの計算コストが、従来例を１００とする指数で示されている。数値が小さいほど、計算コストが小さいことを示している。なお、計算コストの指数が１０５以下であれば、良好である。共通仕様等は、次のとおりである。
分子鎖：スチレンブタジエンゴム
高分子材料モデル：
分子鎖モデルの個数：１０個
初期密度：０．１ｇ／cm³
緩和計算後の密度：０．９４ｇ／cm³
各分子鎖モデル（全原子モデル）：
炭素粒子モデル：１００個
水素粒子モデル：２０２個
構造緩和計算：
温度：３００Ｋ
圧力：１atm
繰り返し計算数：
従来例：１００００００回
実施例１～３：
第１計算工程：１０００回
第２計算工程：１００００００回
実施例４：
第１計算工程：１０００回
第２計算工程：１０００回
第３計算工程：１００００００回
テスト結果が、表１に示される。

テストの結果、実施例１～４は、従来例に比べて、緩和計算を正常に終了させることができた。また、分子力学計算が実施された実施例３は、実施例１～２に比べて、緩和計算をより多く正常終了させることができた。さらに、第３計算工程が実施された実施例４は、実施例１～３に比べて、緩和計算をより多く正常終了させることができた。

Ｓ５緩和工程
Ｓ５１第１計算工程
Ｓ５２第２計算工程

Claims

高分子材料のシミュレーション方法であって、
前記高分子材料の分子鎖に基づいて、数値計算用の分子鎖モデルをコンピュータに入力する工程と、
前記高分子材料の一部に対応する仮想空間であるセルを、前記コンピュータに定義する工程とを含み、
前記コンピュータが、
前記セルの内部に、複数の前記分子鎖モデルを配置する工程と、
複数の前記分子鎖モデルを対象に、分子動力学計算に基づく構造緩和を計算する緩和工程とを含み、
前記緩和工程は、
予め定められた第１の時間刻みで、前記分子鎖モデルの運動を進展させる第１計算工程と、
前記第１計算工程の後に、前記第１の時間刻みよりも大きな第２の時間刻みで、前記分子鎖モデルの運動を進展させる第２計算工程とを含む、
高分子材料のシミュレーション方法。
前記第１の時間刻みは、前記第２の時間刻みの０．１～０．５倍である、請求項１記載の高分子材料のシミュレーション方法。
前記第１の時間刻みは、０．０５～０．２ｆｓであり、
前記第２の時間刻みは、０．３～１．２ｆｓである、請求項１又は２記載の高分子材料のシミュレーション方法。
前記第１計算工程での繰り返し計算数、及び、前記第２計算工程での繰り返し計算数は、１００～１００００００回である、請求項１ないし３のいずれかに記載の高分子材料のシミュレーション方法。
前記緩和工程は、前記第１計算工程に先立ち、複数の前記分子鎖モデルを対象に、分子力学計算を行う工程をさらに含む、請求項１ないし４のいずれかに記載の高分子材料のシミュレーション方法。
前記緩和工程は、前記第２計算工程の後に、前記第２の時間刻みよりも大きな第３の時間刻みで、前記分子鎖モデルの運動を進展させる第３計算工程をさらに含む、請求項１ないし５のいずれかに記載の高分子材料のシミュレーション方法。