JP2016181257A - ファンデルワールス相互作用から生じる圧力の計算方法および、高分子−低分子間の相互作用エネルギーまたは高分子中の低分子の配置のサンプリング方法 - Google Patents

Abstract

【課題】
高分子と低分子の混合系における、スケールされた非結合性相互を用いた分子動力学シミュレーションにおいて、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞を求め、正確な圧力の計算方法を提供する。
【解決手段】
コンピュータに、原子の座標 、スケール値など初期条件を入力し、（１）分子内力・分子間力を計算し、（２）カットオフ半径内のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^0→rcを計算し、（３）カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞を計算し、（４）Ｐ_vdW ^0→rcとＰ_vdW ^rc→∞との和からＰ_vdW ^0→∞計算し、（５）微小時間Δt後の原子の座標、速度を数値積分によって更新し、さらに、（５）の更新値を初期値として（１）〜（５）を再度実行させ、これを規定回数繰り返すことを含む圧力の計算方法。
【選択図】図１

Description

本発明は、コンピュータを利用した分子動力学シミュレーションにおける、高分子中に低分子が存在する系におけるファンデルワールス相互作用から生じる圧力の計算方法および、高分子-低分子間の相互作用エネルギーまたは高分子中の低分子の配置のサンプリング方法に関する。

分子動力学シミュレーションは、分子や分子集合体の高次構造や物性を詳細に解析するための手法として広く用いられる。分子動力学シミュレーションは、近年のコンピュータの発達により、合成高分子や生体高分子といったソフトマターの研究にもよく適用されるようになってきた。

合成高分子の分子動力学シミュレーションの適用例として、高分子分離膜における物質透過性予測が挙げられる（非特許文献１）。物質の膜透過性を予測するにあたり重要な物理量としては、透過物質である低分子の高分子膜への溶解性を示す溶解度係数がある。溶解度係数を分子動力学シミュレーションから求めるためには、低分子と高分子の相互作用エネルギーを高精度に計算する必要がある。しかし高分子の遅い分子運動のため、高分子中の低分子は通常の分子動力学シミュレーションの計算時間ではほとんど移動しない。そのため限定された空間の相互作用エネルギーしかサンプリングできず、統計力学的に高精度に溶解度係数を求めることができない。

低分子の移動を促すため、低分子と高分子間の非結合相互作用を時間とともにスケールし非結合性相互作用を増減させるという方法で、低分子の運動性を向上させ統計精度を高めた高精度な相互作用エネルギー計算手法が考案された（特許文献１）。

特開2011-123874号公報

Y. Tamai, H. Tanaka, and K. Nakanishi, Macromolecules, 28, 2544-2554(1995)

分子動力学シミュレーションにおいて、系の全圧力P_totalは、下記数式（１）のように表される。

（ただし、P_intraは分子内相互作用から生じる圧力、P_vdW ^0→∞はファンデルワールス相互作用から生じる圧力、P_coulomb ^0→∞はクーロン相互作用から生じる圧力、P_Kは原子の運動エネルギーより生じる圧力である。）
ファンデルワールス相互作用から生じる圧力Ｐ_vdW ^0→∞は、カットオフ半径内のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^0→rcとカットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞の和で表される。

ここで、特許文献１のようにスケールした非結合性相互作用を使用した場合、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞を計算できないという問題があった、すなわち、ファンデルワールス相互作用から生じる圧力Ｐ_vdW ^0→∞を正確に求める手法は未だ提唱されておらず、そのため、系の全圧力Ｐ_totalも正確に求めることができなかった。

一方、十分に大きなカットオフ半径を用いれば圧力Ｐ_vdW ^rc→∞の項は０に漸近するため、正確な圧力に近い値を得ることができる。しかし計算時間の観点から一般的に好適に用いられるカットオフ半径を用いた時に比べ長時間の計算を要するため、大きなカットオフ半径を用いることは非効率的である。

正確な圧力を計算できないことによる問題点として、分子動力学シミュレーションにおける現実系を再現するために行う、圧力を一定に制御した状態のシミュレーションを正しく行えないことが挙げられる。圧力を制御するために、ある状態における系の瞬間の圧力を計算し、その圧力が設定圧力に近づくようにシミュレーションセルを変形させる方法が用いられるが、正確な圧力を計算できない場合、シミュレーションセルの大きさを過大評価あるいは過小評価し、その結果、密度や空隙率が正しく計算できず、正確に相互作用エネルギーを求めることができなかった。

以上述べたように、高分子と低分子の混合系において、スケールされた非結合性相互を用いた分子動力学シミュレーションでは、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞を求め、正確な圧力を計算する方法は未だ提唱されていない。そこで本発明では、上記問題を解決したスケールされた非結合性相互を用いた場合のファンデルワールス相互作用から生じる圧力の計算方法を提供することを目的とする。

本発明は上記課題を解決するために、次の手段を採用するものである。すなわち、本発明は以下のとおりである。

コンピュータを利用した高分子と低分子の混合系の分子動力学シミュレーションにおいて、下記数式（２）のように、
（i）低分子の分子内ｖｄＷ相互作用V^vdw _low-intra(r, λ_low-intra)と、
（ii）高分子の分子内ｖｄＷ相互作用V^vdw _poly-intra(r, λ_poly-intra)と、
（iii）低分子‐低分子間ｖｄＷ相互作用V^vdw _low-low(r, λ_low-low)と、
（iv）高分子‐高分子間ｖｄＷ相互作用V^vdw _poly-poly(r, λ_poly-poly)と、
（v）低分子‐高分子間ｖｄＷ相互作用V^vdw _low-poly(r, λ_low-poly)
との和（0≦λ≦１、ただし（i）から（v）のλのうち少なくとも１つは0≦λ<１）で表される、非結合性相互作用のスケール値λを含むスケールファンデルワールス（vdW）相互作用ポテンシャルV^vdW(r, λ)を用いた場合の前記混合系のファンデルワールス（vdW）相互作用から生じる圧力の計算方法であって、

コンピュータに、
・前記高分子と前記低分子に含まれる原子の座標r
および、必要に応じ
・原子の速度ｖ
の初期値、
ならびに
・結合情報
・原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ
・原子の質量ｍ
・非結合性相互作用のスケール値λ
・相互作用計算のカットオフ半径ｒ_ｃ
を入力した後、コンピュータに、
（ステップ１）原子の座標r、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、非結合性相互作用のスケール値λ、ならびに下記数式（３）で示したポテンシャル関数V^totalを原子の座標rで微分した式を用いて、前記高分子および前記低分子の分子内力(数式（３）第１項から第４項の微分の和)および分子間力（数式（３）第５項と第６項の微分の和）を計算するステップ；

（ただし、R、θ、φ、ωはそれぞれ注目する結合の距離、変角の角度、二面角の角度、反転の角度を表す。）
（ステップ２）ステップ１で計算した分子間力の一つであるファンデルワールス（vdW）力（数式（３）第５項の微分）と、原子の座標rからビリアル定理を用いて圧力Ｐ_vdW ^0→ｒｃを計算するステップ；
（ステップ３）原子の座標r、原子のポテンシャルパラメータタイプ、ポテンシャルパラメータ、非結合性相互作用のスケール値λおよび相互作用計算のカットオフ半径ｒ_ｃから、下記数式(４Ａ)および（４Ｂ）を用いて圧力Ｐ_vdW ^rc→∞を計算するステップ；

（ただし、Ｍはポテンシャルパラメータタイプの総数、Ｎ_α、Ｎ_βはそれぞれポテンシャルパラメータタイプα、βを持つ原子の数、diag[・・・]は対角行列を表す。）
（ステップ４）下記数式（５）によりファンデルワールス（vdW）相互作用から生じる圧力Ｐ_vdW ^0→∞を計算するステップ；

（ステップ５）原子の質量ｍ、分子内力、分子間力を用いて、微小時間Δt後の
・原子の座標r
・原子の速度ｖ
を数値積分によって更新するステップ：
を実行させ、
さらに、ステップ５で更新した微小時間Δt後の原子の座標r、原子の速度ｖを初期値としてステップ１〜ステップ５を再度実行させ、これを繰り返すことを含むファンデルワールス（vdW）相互作用から生じる圧力の計算方法。

また、本発明は、上記の圧力の計算方法を用いた高分子-低分子間の相互作用エネルギーまたは高分子中の低分子の配置のサンプリング方法を含む。

本発明の圧力の計算方法を示すフローチャートである。 本発明の圧力の計算方法を利用した高分子-低分子間の相互作用エネルギーや高分子中の低分子の配置のサンプリング方法を示すフローチャートである。 実施例１〜６、比較例１〜６における非結合性相互作用のスケール値と密度の関係を示す図である。 実施例１〜６、比較例１〜６における非結合性相互作用のスケール値と空隙率の関係を示す図である。 実施例５、比較例５、比較例７〜９におけると計算実行時間と密度、カットオフ半径の関係を示す図である。

図１に本発明の計算方法のフローを示す。以下、本発明の好ましい実施形態を、図１のフローに従って説明する。

本発明の計算方法においては、まず、コンピュータに、計算に必要な諸条件を入力する。計算に必要な諸条件を入力するに際しては、コンピュータに読み取り可能な形態で、高分子と低分子に含まれる原子の座標rの初期値、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、原子の質量ｍ、非結合性相互作用のスケール値λ、相互作用計算のカットオフ半径ｒ_ｃを入力する。また、必要に応じて原子の速度ｖの初期値を入力することもできるが、入力しなければマクスウェル分布に従うよう計算された原子の速度が初期速度として与えられるよう構成すればよい。

入力方法は特に限定されないが、一般に入手可能な分子モデル作成ソフトウェアやテキストエディタなどが利用でき、当該情報が記述されたファイルをコンピュータに読み取らせ、メモリに格納することが好ましい。

〔ステップ１〕
高分子および低分子の分子内力(数式（３）第１項から第４項の微分の和)および分子間力（数式（３）第５項と第６項の微分の和）を計算するステップ１においては、メモリに格納された原子の座標r、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、非結合性相互作用のスケール値λ、ならびに数式（３）で示したポテンシャル関数V^totalを原子の座標rで微分した式を用いて分子内力・分子外力を計算し、メモリに格納する。

ここで、数式（３）の右辺第１項V^bond(R)は結合長、第２項V^angle(θ)は結合角、第３項V^torsion(φ)は二面角、第４項V^inversion(ω)は反転、第５項V^vdw(r,λ)はファンデルワールス（ｖｄＷ）、第６項V^coulomb(r,λ)はクーロンのそれぞれポテンシャル関数を表わす。なお、数式（３）の右辺第１項から第４項を総称して分子内ポテンシャル関数と呼び、右辺第５項、第６項をまとめて分子間ポテンシャル関数と呼ぶ。また、分子内ポテンシャル関数、分子間ポテンシャル関数を原子の位置で微分した式を用いて計算した力をそれぞれ分子内力、分子間力と呼ぶ。

数式（３）における各ポテンシャル関数の具体的な表式は下記数式（６）で表される。

数式（６）の右辺第５項、第６項はスケールされた非結合性相互作用を表し、非結合性相互作用のスケール値λ=1のときスケールされていない一般的な分子動力学シミュレーションで使用するポテンシャル関数に等しい。

スケールされたvdW相互作用である第5項V^vdw(r,λ)の具体的な関数形としては、下記数式（７）に示したレナードジョーンズポテンシャルやモースポテンシャル、バッキンガムポテンシャルが好ましく用いられる。

数式（７）において、m1、m2は型を決定する乗数であり、LJm1/m2型と表記するとき、LJ12/6型、LJ9/6型、LJ12/4型、LJ6/4型、LJ12/10型が好適に用いられる。さらにLJ12/6型には定数A,BやvdWポテンシャルの深さと斥力が生じる位置を決めるパラメータD^０、R^０の置き方により下記数式（８）で示したDreiding型〔 S. L. Mayo 、B. D. Olafson 、W. A. Goddard III 、J. Phys. Chem.、1990、94、p.8897.〕やAmber型〔W. D. Cornell、P. Cieplak 、C. I. Bayly 、I. R. Gould 、K. M. Merz jr 、D. M. Ferguson 、D. C. Spellmeyer 、T. Fox 、J. W. Caldwell、P. A. Kollman 、J. Am. Chem. Soc.、1995、117、p.5179〕、OPLS型［W. L. Jorgensen、D. S. Maxwell、Julian Tirado-Rives、J. Am. Chem. Soc.、1996、118、p.11225］、CHARMM型［B. R. Brooks、R. E. Bruccoleri、B. D. Olafson、D. J. States、S. Swaminathan、M. Karplus、J. Comput. Chem.、1983、4、p.187］などいくつかの形式が存在し、どのLJ12/6型でも好適に用いることができる。

上記数式（６）中の結合長、結合角、二面角、反転おけるそれぞれの力の定数K^bond、K^angle、K^torsion、K^inversionおよびそれぞれの平衡位置R⁰、θ⁰、φ⁰、ω⁰、数式（７）中のvdWポテンシャルの深さと斥力が生じる位置を決めるパラメータD^０、R^０には、DREIDING、AMBER、hybrid/UA[A. L. Frischknecht、J. G. Curro、Macromolecules、2003、36、p.2122]、桑島らのパラメータ〔桑島聖、野間祐之、逢坂俊郎、“非晶高分子のＭＤ計算用力場パラメータについてｉｉ：ポリオレフィン”、第４回計算化学シンポジウム&ＣＣＳ展 講演要旨集、日本化学プログラム交換機構（１９９４）〕等の公知のパラメータを用いることができる。

計算精度について考慮した場合、以下に述べる方法によってポテンシャルパラメータを
決定することがさらに好ましい。

上記数式（３）および数式（６）の第１項から第４項で表わされる分子内ポテンシャル関数のパラメータは量子化学計算によって決定することが好ましい。量子化学計算としては、HF法〔Ａ．ザボ，Ｎ．Ｓ．オストランド，「新しい量子化学」，東京大学出版会（１９８７）．〕を用いることが好ましく、B3LYP法〔A. D. Becke, J. Chem. Phys., 1993, 98, p.5648／C. Lee, W. Yang, R. G. Parr, Phys. Rev. B, 1988, 37, p.785／B. Miehlich, A. Savin, H. Stoll,H. Preuss, Chem. Phys. Lett., 1989, 157,p.200.〕を用いることがさらに好ましい。

ＨＦ法やＢ３ＬＹＰ法を用いる場合、１電子軌道を展開する基底関数を指定する必要がある。計算時間と計算精度を考慮すると、６−３１Ｇ（ｄ）基底関数、６−３１Ｇ（ｄ，ｐ）基底関数、６−３１＋Ｇ（ｄ，ｐ）基底関数〔W. J.Hehre, R. Ditchfield, J.A. Pople, J. Chem. Phys. 1972, 56, p.2257， P. C. Hariharan, J.A. Pople, Theoret. Chimica Acta, 1973, 28, p.213.〕が好適に用いられる。

上記数式（３）および数式（６）の第５項V^vdw(r,λ)で表わされるvdWポテンシャルのパラメータについては、密度と蒸発熱の実験値を再現するように決定することが好ましい。

上記数式（３）および数式（６）の第６項V^coulomb(r,λ)のクーロンポテンシャルのパラメータ、すなわち電荷については、QEq法〔A.K.Rappe,W.A.Goddard III,J.Phys.Chem.,95,3358-3363(1991).〕等の学術的に信頼性の高い手法を用いて決定することが好ましく、量子化学計算から静電ポテンシャル（ＥＳＰ）を求め、ＥＳＰフィッティングして電荷を決定することがさらに好ましい。静電ポテンシャルを求めるための量子化学計算としては、ＨＦ法やＢ３ＬＹＰ法を好適に用いることができ、基底関数としては６−３１Ｇ（ｄ）基底関数、６−３１Ｇ（ｄ，ｐ）基底関数、６−３１＋Ｇ（ｄ，ｐ）基底関数が好適に用いられ、さらに好ましくは、ａｕｇ−ｃｃ−ｐＶＤＺ基底関数、ａｕｇ−ｃｃ−ｐＶＴＺ基底関数、ａｕｇ−ｃｃ−ｐＶＱＺ基底関数〔T.H. Dunning, Jr.， J. Chem. Phys.，90，1007 (1989).〕が用いられる。ＥＳＰフィッティングして電荷を求める手法としては、ＭＫ法〔B. H. Besler, K. M. Merz Jr, P. A. Kollman, J. Comput. Chem., 1990 11, p.431.〕やＣＨｅｌｐＧ法〔C. M. Breneman, K. B. Wiberg, J. Comput. Chem., 1990, 11, p.361.〕といった方法を好適に用いることができる。

以上のような量子化学計算を実行するために、有償または無償で入手可能な種々の量子化学計算プログラムパッケージを利用することができる。例えば、Gaussian （登録商標）09 Revision C.01, M. J. Frisch他, ガウシアン社（Gaussian, Inc.）, Pittsburgh PA, 2002.〕、ＧＡＭＥＳＳ〔M. W. Schmidt 他, J. Comput. Chem., 1993, 14, p.1347.〕といったプログラムを好適に用いることができる。

高分子分離膜のように高分子中に低分子が存在する系において、上記数式（３）および数式（６）の第５項・第６項における非結合性相互作用は、下記数式（９）のように
（i）低分子の分子内非結合性相互作用V_low-intra(r, λ_low-intra)、
（ii）高分子の分子内非結合性相互作用V_poly-intra(r, λ_poly-intra)、
（iii）低分子−低分子間非結合相互作用V_low-low(r, λ_low-low)、
（iv）高分子‐高分子間非結合性相互作用V_poly-poly(r, λ_poly-poly)、
（v）低分子‐高分子間非結合性相互作用V_low-poly(r, λ_low-poly)の五種類の和で表される。

非結合性相互作用のスケール値λの範囲として0≦λ≦1を用いることができるが、（i）から（v）いずれかのλのうち1つ以上0≦λ<１である必要がある。(i）から（v）それぞれのλの範囲として、λ_poly-intra＝λ_low-low=λ_poly-poly＝１、0≦λ_low-intra≦1、0≦λ_low-poly＜1とすることが好ましい。これは高分子内、高分子間および低分子間に関する非結合性相互作用をスケールしないことを示す。このような取り扱いにより、高分子の平衡構造を維持しながら低分子の運動性を上昇させることができる。なお、低分子を１分子のみとして計算を行う場合には、λ_low-lowは設定不要であり、水やメタンなど分子内非結合相互作用を持たない低分子の計算を行う場合には、λ_low-intraの設定は不要である。一般的に分子動力学シミュレーションでは、分子内非結合相互作用は１−４相互作用以上を扱うため、１−２、１−３相互作用までしか存在しない水やメタンは、分子内非結合相互作用を計算しない。なお１−Ｘとは基準となる原子を１とし、共有結合でつながっている隣の原子を２、さらに２と共有結合でつながっている原子を３と順に番号付けしたときの相互作用ペアを表す。

さらにλ_poly-intra＝λ_low-low=λ_poly-poly＝λ_low-intra＝１、0≦λ_low-poly＜1とすることが特に好ましい。λ_low-intra＝1とすることによって、分子内水素結合を持つ低分子の配座を維持しながら低分子の運動性を上昇させることができる。エチレングリコールのような分子内水素結合を持つ低分子は、分子内非結合性をスケールすると水素結合が弱まるため本来とるべき配座を再現しない。このように分子内非結合性相互作用が低分子のコンフォメーションに大きく寄与する場合、低分子内の非結合性相互作用をスケールしなければ低分子内の配座を維持しながら運動性を上昇させることができる。

λ_low-polyについては、十分に低分子に高分子中を運動させるために、λ_low-poly≦0.3が好ましく、λ_low-poly≦0.2がより好ましく、λ_low-poly≦0.1が特に好ましい。

数式（７）、数式（８）におけるエネルギーと力の距離依存性を調節するパラメータδについては、δ＝0ｍ²とし、かつ非結合性相互作用のスケール値λが非常に小さい値をとる場合に、数式（３）および数式（６）の第５項V^vdw(r,λ)で表わされるファンデルワールス（ｖｄＷ）相互作用が急激に変化し、数値積分の安定性が低下する。この傾向は原子間の距離が近距離となった場合に特に顕著になる。この問題を解決するために、非結合性相互作用のスケール値λの最小値を非常に小さな値に設定する場合には、δ＞0ｍ²とすることが好ましい。さらに好ましくは、δ≧5.0×10^-21ｍ²である。一方、δの値が大きすぎると、スケール値が変化した際、ｖｄＷ力の変化が大きくなりすぎるため、δ≦5.0×10^-19ｍ²とすることが好ましい。このδの値としては、δ＝1.0×10^-20ｍ²が好適に用いられる。さらにカットオフ半径r_cは、計算時間と計算精度の観点から8Å≦r_c≦12Åが好ましく、r_c＝10Åが特に好ましい。

〔ステップ２〕
カットオフ半径内のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^0→ｒｃを計算するステップ２では、ステップ１で計算したファンデルワールス力（数式（３）第５項の微分）と、原子の座標rとから、下記数式（１０Ａ）および（１０Ｂ）で表されるビリアル定理を用いて圧力Ｐ_vdW ^0→ｒｃを計算する。

（ここでF^vdWは注目する原子i、j間に働くファンデルワールス力、l, m(=1、2、3)はベクトルの成分を表す。）
〔ステップ３〕
カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞を計算するステップ３では、原子の座標r、原子のポテンシャルパラメータタイプ、ポテンシャルパラメータ、非結合性相互作用のスケール値λおよび相互作用計算のカットオフ半径ｒ_ｃから前記数式(４Ａ)および（４Ｂ）を用いて圧力Ｐ_vdW ^rc→∞計算する。

数式（４Ａ）の通り圧力Ｐ_vdW ^rc→∞は対角項のみで、非対角項は０となる。なお数式（４Ｂ）におけるスケールされたファンデルワールス相互作用として前述の様々な関数形を用いることができる。

前述のように、これまで、非結合性相互作用のスケール値λを含むポテンシャルでの圧力Ｐ_vdW ^rc→∞の計算方法が確立されていなかったが、本発明者らの検討の結果、数式（８）のDreiding型レナードジョーンズポテンシャルを数式（４Ｂ）に代入して得られた解析解である下記数式（１１Ａ）および（１１Ｂ）と数式（４Ａ）から圧力Ｐ^rc→∞を計算する方法が好ましいことが見出された。

〔ステップ４〕
ファンデルワールス相互作用から生じる圧力Ｐ_vdW ^0→∞を計算するステップ４では、数式（５）の通り、ステップ２およびステップ３で計算したカットオフ半径内のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^0→ｒｃとカットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ^rc→∞の和から圧力Ｐ_vdW ^0→∞を計算する。

また数式（１）に、数式（４Ａ）および数式（４Ｂ）、数式（５）、数式（１０Ａ）および数式（１０Ｂ）より求めたＰ_vdW ^0→∞を代入することにより全圧力Ｐ_totalを求めることができる。なおファンデルワールス相互作用から生じる圧力Ｐ_vdW ^0→∞以外の圧力の詳細な表式もビリアル定理から導くことができ、参考文献［Glenn J. Martyna, Mark E. Tuckerman, Douglas J. Tobias, Michael L. Klein, Mol. Phys., 1996, 87, p.1117.］に記述されている。この文献の全圧力の表式はクーロン相互作用から生じる圧力Ｐ_coulomb ^0→∞を、Ewald法やParticle Mesh Ewald(PME)法といった無限遠のクーロン相互作用の寄与を考慮できる方法を用いた場合の表式である。本発明においてもEwald法やPME法を用いることが好ましく、このような取り扱いをすることで無限遠までの相互作用を考慮した正確な全圧力Ｐ_totalを得ることができる。

〔ステップ５〕
微小時間Δt後の原子の座標r、原子の速度ｖを数値積分によって更新するステップ５においては、原子の質量ｍ、分子内力、分子間力を用いて、微小時間Δt後の原子の座標と速度を数値積分によって計算し、更新する。

数値積分には、差分近似法や予測子−修正子法〔岡崎進，「コンピューターシミュレーションの基礎」，化学同人（2000）．〕などの一般的な解法を好適に用いることができる。異なる時定数の力ごとに大きさの異なる時間刻みΔtを設定することのできる多重時間刻み数値積分（ＲＥＳＰＡ法）法〔G. J. Martyna、Mol. Phys.、1996、87、p.1117.〕を用いると、シミュレーション時間を短縮できるためさらに好ましい。

〔ステップ１〜ステップ５の繰り返し〕
以後、ステップ５で更新した微小時間Δt後の原子の座標r、原子の速度ｖを初期値としてステップ１〜ステップ５を再度実行させることを繰り返し、あらかじめ規定した回数実行する。

ステップ１〜ステップ５の繰り返し回数は、本発明の最初に行われる諸条件の入力に際し初期条件として入力するか、またはあらかじめプログラム内に記述しておくことで規定する。

微小時間Δtは、最も速い時定数の運動をとらえることができるよう十分小さな値でなければならない。一般的に最も速い時定数を持つ運動は水素原子を含む結合の伸縮運動である。この伸縮運動を精度良くとらえることができ、かつ計算時間の観点から微小時間Δtは、０．１フェムト秒以上０．３フェムト秒以下であることが好ましく、０．２５フェムト秒がより好ましい。

また繰り返し回数は、系の平衡構造における構造の揺らぎが十分にシミュレーション時間中にあらわれる程度の回数が必要である。系の構造によって必要な繰り返し回数ｎは異なるが、高分子と低分子の混合系であれば８００００００回（微小時間Δtが０．２５フェムト秒のとき２ナノ秒）以上が好ましく、１６００００００回（微小時間Δtが０．２５フェムト秒のとき４ナノ秒）以上がより好ましい。
〔結果の出力〕
本発明の計算方法においては、ステップ１〜ステップ５を所定回数繰り返した後に、ファンデルワールス相互作用から生じる圧力Ｐ_vdW ^0→∞を出力する。ステップ１〜ステップ５の所定回数の繰り返し時点における瞬間値、あるいは一定回数繰り返して得られた平均値のどちらか一方または両方を出力することができる。

また、本発明においては、必要に応じ、上記のように出力した圧力Ｐ_vdW ^0→∞の値を用いて諸物理量や原子の座標を計算し、出力することができる。出力する諸物理量は特に限定されず、全圧力Ｐ_total、系の密度などが挙げられる。出力する諸物理量は、所定回数の繰り返し時点における瞬間値、あるいは一定回数繰り返して得られた平均値のどちらか一方または両方とすることができる。原子の座標については、分子動力学シミュレーション最大の特徴である原子の時間経過による運動をとらえることができるという性質から、ある繰り返し回数における瞬間値を出力することが好ましい。

以上の手順を実行すれば、高分子中に低分子が存在する系においてスケールされた非結合性相互作用を用いた場合に、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞を求められず正確な圧力を計算できないという課題を解決することができる。

〔相互作用エネルギーまたは配置のサンプリング〕
特に、本発明の計算方法は、高分子-低分子間の相互作用エネルギーや、高分子中の低分子の配置のサンプリングを行うために好適に用いることができる。様々なλ_low-polyにおいて、図１のフローに従ったシミュレーションを並列に実行することにより、様々なλ_low-polyにおける低分子‐高分子間非結合性相互作用V_low-polyや高分子中の低分子の配置をサンプリング可能である。しかしλ_low-polyが大きい場合は高分子中を低分子が移動し難いため、十分なサンプリングを行うには独立な初期配置を複数用意する必要があり、計算時間が長くなる。

そこで、前述の所定回数のステップ１〜ステップ５の繰り返しの後に、λ_low-polyを0≦λ_low-poly≦1の範囲で増減させるステップ６を加えることにより、様々なλ_low-polyにおける低分子‐高分子間非結合性相互作用や高分子中の低分子の配置を、短時間かつ高い統計精度でサンプリングすることができる。

以下、本サンプリング方法の好ましい実施形態を、図２のフローに従って説明する。なお、図１と同様の部分は省略し、本サンプリング方法特有の操作に絞り説明する。

本サンプリング方法においては、非結合性相互作用のスケール値λの範囲として0≦λ≦1を用いることができるが、λ_poly-intra＝λ_low-low=λ_poly-poly＝λ_low-intra＝１、0≦λ_low-poly≦1とすることが好ましい。λ_poly-intra＝λ_low-low=λ_poly-poly＝λ_low-intra＝１とする理由は前述の通りである。λ_low-polyの初期値は0≦λ_low-poly≦1の範囲であればどのような値でもよいが、数値処理の行いやすさの観点から０または１が好ましい。

図２のステップ１〜図２ステップ５までは、図１のステップ１〜図１ステップ５までと同様である。

〔ステップ６〕
ステップ６は、所定回数（後述する「Δｎ」回）のステップ１〜ステップ５の繰り返しに１回、ステップ１〜５にさらに加えて実行される。λ_low-polyの値を増減させるステップ６においては、λ_low-polyの増減間隔Δλをλ_low-polyに加減する。なおλ_low-polyの増減の規則は、λ_low-polyが０から１まで変化するようなものであればどのようなものでも構わないが、様々なλ_low-polyの状態の相互作用をサンプリングするという観点から、ステップ１〜５、および所定回数に１回のステップ１〜ステップ６の繰り返しにおいて、ステップ６の実行時に１回以上λ_low-poly＝０およびλ_low-poly＝１となるよう設定することが好ましい。

さらに、ステップ６の実行時に１回以上λ_low-poly＝０およびλ_low-poly＝１となるよう設定する場合、各λ_low-polyにおける低分子‐高分子間非結合性相互作用と、熱力学的積分法、自由エネルギー摂動法、ベネット受容比法またはマルチステイトベネット受容比法を組み合わせることにより、高分子中に低分子が溶解するときの自由エネルギーを求めることができる。なお、熱力学的積分法の場合は相互作用をλで微分した値を用いる。

増減間隔Δλは、０＜Δλ≦１となる値であればどのような値であってもよい。ステップ６を実行したときに１回以上λ_low-poly＝０およびλ_low-poly＝１なるようにするためには、λ_low-polyにｍ回（ｍは自然数）加えるまたは引いたときに１または０となる増減間隔Δλを設定することが好ましい。また自由エネルギー計算まで考慮に入れると、Δλが大きすぎる場合、隣り合うλ_low-poly＋ｍΔλの状態とλ_low-poly＋（ｍ±１）Δλの状態間のエネルギー差が大きく、正確に自由エネルギーを求めることができない。一方Δλが小さすぎると、自由エネルギーの計算精度は上がるものの、それに伴い計算コストも上昇する。これらの観点から、増減間隔Δλは０．０１以上０．３以下が好ましく、０．０３以上０．２以下がさらに好ましく、０．０５以上０．１５以下が特に好ましい。

ステップ６を実施する間隔Δｎは１以上の値であれば、どのような値でもよい。しかし、Δｎが小さすぎるとサンプリングできる空間は増加するものの局所的な構造のサンプリング効率が低下する、一方Δｎが大きすぎるとサンプリングできる空間が狭くなる。そこで、Δｎは１以上１００００以下が好ましく、１００以上７０００以下がより好ましく、５００以上５０００以下が特に好ましい。

所定回数に１回ステップ６を実施する回を含めた全体の繰り返し回数ｎとステップ６を実施する間隔Δｎを用いると、ステップ６を実施するのはｎをΔｎで割った余りが０に等しいときである。

本サンプリング方法においては、ステップ１〜ステップ５、および所定回数に１回のステップ１〜ステップ６を繰り返した後に、高分子−低分子間の相互作用エネルギーV_low-polyおよび原子の座標ｒを出力する。ここでは、ステップ１〜ステップ５（ステップ６を実施する回を含む）の所定回数の繰り返し時点における瞬間値、あるいは一定回数ごとに得られた平均値のどちらか一方または両方を出力することができる。ただし原子の座標ｒに関しては平均値に物理的な意味が乏しいため瞬間値を出力することが好ましい。また各λ_low-polyでのサンプリングという観点から、各λ_low-poly毎に高分子−低分子間の相互作用エネルギーV_low-polyおよび原子の座標ｒを出力することが好ましい。

〔プログラムおよび記録媒体〕
本発明のプログラムは、コンピュータの構成要素に上記の圧力の計算方法または相互作用エネルギーと配置のサンプリング方法の各ステップにおける計算手段として動作させるものである。このプログラムを用いれば、上述の計算方法を実行することができる。

また、本発明の記録媒体は、当該プログラムをコンピュータ読み取り可能に記録したものである。この記録媒体のプログラムをコンピュータにて読み取って実行すれば、上述の計算方法を実行することができる。

以下に実施例を挙げて本発明を具体的に説明するが、これらは例示的なものであって限定的なものではない。

［実施例１］
ポリジメチルシロキサン（PDMS）４分子／水１分子の系について、分子動力学シミュレーションを実施した。分子動力学シミュレーションは東レ株式会社で独自に開発したFORTRANプログラムを用いて行い、分子軌道計算は汎用プログラムGaussian（登録商標）09を用いて行った。また、計算にはクロック周波数2.93GHzのCore i7（登録商標）870を搭載したコンピュータ、およびクロック周波数2.7HzのXeon（登録商標）E5-2680を搭載したコンピュータを用いた。

アクセルリス社(現：バイオビア社)製の分子設計システムMaterials Studio（登録商標）7.0上で下記構造式（１）に示すポリマーを作成した。

構造式(1)において、n=40（分子量＝2996）とした。次いで、Materials Studio（登録商標）7.0 Amorphous Cell Moduleを用いて、表１に示した組成で、PDMS分子と水分子を三次元周期境界条件下でランダムに配置し、系を作成した。Materials Studio（登録商標）7.0の外部出力機能を使用し、作成した系の原子の座標r 、結合情報、原子のポテンシャルパラメータタイプ、原子の質量ｍ、基本セルベクトルhが記載された分子構造ファイルを作成した。

次に、テキストエディタを用いて分子構造ファイルに記載されたポテンシャルパラメータタイプに対応したポテンシャルパラメータを記述したパラメータファイルを作成した。

さらにテキストエディタを用いて非結合性相互作用のスケール値λ、相互作用計算のカットオフ半径ｒ_ｃ、図1計算フローにおけるループの繰り返しを規定する回数n、出力する物理量とその出力間隔、設定温度T₀、設定圧力P₀、熱浴の慣性量Q、barostatの慣性量Wを記述した計算条件ファイルを作成した。

作成した分子構造ファイル、パラメータファイル、計算条件ファイルを分子動力学シミュレーションプログラムで読み取らせ、データをメモリに格納した。

計算条件ファイルにおいて、非結合性相互作用のスケール値λを相互作用の種類ごとに設定した。低分子‐高分子間非結合性相互作用のスケール値λ_low-poly＝0.1、高分子の分子内非結合性相互作用と、高分子‐高分子間非結合性相互作用のそれぞれのスケール値λ_poly-intra＝λ_poly-poly＝１とした。本実施例における低分子は水１分子であるため、低分子の分子内非結合性相互作用のスケール値λ_low-intra および低分子‐低分子間非結合性相互作用のスケール値λ_low-lowは設定しない。

また、計算条件ファイルにおける他のパラメータとして、カットオフ半径ｒ_ｃ＝１０オングストローム、図1計算フローにおけるループの繰り返しを規定する回数を１６００００００回、出力する物理量を原子の座標r、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞および密度とし、その出力間隔を下記ステップ１〜ステップ５のループを４０００回繰り返すごとに１回とした。

温度制御にNose−Hoover法〔M.Tuckerman,B.J.Berne,G.J.Martyna,J.Chem.Phys.,1992, 97,p.1990.〕、圧力制御にParrinello−Rahman法〔M. Parrinello、A. Rahman、J. Appl. Phys.、1981、52、p.7182〕を用いて、原子数・温度・圧力一定（ＮＰＴ）アンサンブルを構成することにより、系を設定温度T₀、設定圧力P₀に制御した状態で下記ステップ１〜ステップ５を実施した。なお、設定温度T₀=20℃、設定圧力P₀=1atm、熱浴の慣性量Q=120.95(kJ・ps²)/mol、barostatの慣性量W=483.79(kJ・ps²)/molとした。

ＮＰＴアンサンブルにおける分子動力学シミュレーションを実施するためには、下記ステップ５において原子の座標r、原子の速度v、熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルh、基本セルベクトルの速度hvを更新することが必要である。そのため分子構造ファイルおよび計算条件ファイルで入力していない、原子の速度v、熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルの速度hvに初期値を与える必要がある。そこで原子の速度vの初期値としてT₀=20℃におけるマクスウェル分布に従うよう計算された速度を、熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルの速度hvの初期値として全ての成分に０を与えメモリに格納した。

〔ステップ１〕
メモリに格納された原子の座標r、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、非結合性相互作用のスケール値λ、ならびに数式（８）で示したDreiding型のｖｄＷポテンシャルを含む、数式（６）で示したポテンシャル関数V^totalを原子の座標rで微分した式を用いて分子内力・分子外力を計算しメモリに格納した。

また、数式（８）中のδ＝1.0×10^-20ｍ^２とし、数式（６）第６項のクーロンポテンシャルおよびその力の計算にはＥｗａｌｄ法を用い、Ｅｗａｌｄ法のパラメータとして逆空間のクーロンポテンシャルにα＝0.21オングストローム^-1、｜ｎ｜^２ _ｍａｘ＝50を用いて計算した。

ポテンシャルパラメータの設定は以下の通りである。まず、PDMS分子については、分子内ポテンシャルパラメータとして、結合長・結合角の平衡位置R⁰、θ⁰にはHF/6-31G(d,p)レベルの分子軌道計算によって決定した値を、力の定数K^bond、K^angle、K^torsion、K^inversionおよび平衡位置φ⁰、ω⁰にはＡＭＢＥＲパラメータ、hybrid/UAパラメータおよび桑島ポテンシャルパラメータを用いた。vdWポテンシャルの深さと斥力が生じる位置を決めるパラメータD^０、R^０にはＡＭＢＥＲパラメータ、および桑島ポテンシャルパラメータを用い、電荷はB3LYP/6-31G(d,p)レベルの分子軌道計算によって決定した値を用いた。

水分子については、結合長・結合角の平衡位置R⁰、θ⁰にはＴＩＰ３Ｐパラメータ[W. L. Jorgensen 、R. W. Impey 、J. Chandrasekhar 、 J. D. Madura、M. L. Klein 、J. Chem. Phys.、1983、79、p.926.〕を用い、結合長・結合角の力の定数K^bond、K^angleはＩＲの実測値を再現するように決定した値を用い、vdWポテンシャルの深さと斥力が生じる位置を決めるパラメータD^０、R^０と電荷にはＴＩＰ３Ｐパラメータを用いた。

〔ステップ２〕
ステップ１で計算した分子間力の一つであるvdW力（数式（８）の微分）と、原子の座標rからビリアル定理（数式（１０Ａ）および（１０Ｂ））を用いて圧力Ｐ_vdW ^0→ｒｃを計算しメモリに格納した。

〔ステップ３〕
原子の座標r、原子のポテンシャルパラメータタイプ、ポテンシャルパラメータ、非結合性相互作用のスケール値λおよび相互作用計算のカットオフ半径ｒ_ｃと数式（４Ａ）、（１１Ａ）、（１１Ｂ）から圧力Ｐ_vdW ^rc→∞を計算しメモリに格納した。

〔ステップ４〕
ステップ２および３で計算したカットオフ半径内のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^0→ｒｃとカットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞の和から圧力Ｐ_vdW ^0→∞を計算しメモリに格納した。

〔ステップ５〕
メモリに格納された原子の質量ｍ、原子の座標r、原子の速度v、熱浴の慣性量Q、熱浴の位置ζ、熱浴の速度ζv、barostatの慣性量W、基本セルベクトルh、基本セルベクトルの速度hv、分子内力、分子間力を用いて、微小時間Δt後の原子の座標r、原子の速度v、熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルh、基本セルベクトルの速度hvを数値積分によって計算し、更新した。数値積分には多重時間刻み数値積分法（ＲＥＳＰＡ法）を用い、時間刻みΔtは表２に示した値を用いた。

ステップ５で更新した微小時間Δt後の原子の座標r、原子の速度v、熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルh、基本セルベクトルの速度hvを初期値としてステップ１〜ステップ５を再度実行させることをn＝１６００００００回繰り返した。ＲＥＳＰＡ法における最小の時間刻みが０．２５フェムト秒であるため、この繰り返し回数は４ナノ秒に相当する。

なお本実施例ではNose−Hoover法により温度、Parrinello−Rahman法により圧力を制御しているためスッテプ５において熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルh、基本セルベクトルの速度hvを更新した。

ステップ１〜ステップ５を４０００回（１ピコ秒に相当）繰り返すごとに、原子の座標r、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞および密度を出力した。また計算終了時（ステップ１〜ステップ５を１６００００００回繰り返し終了後）にカットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞および密度の平均値を計算した。

なお、本実施例の計算では構造緩和のための予備計算として、上述の計算の前に２０℃の原子数・温度・体積一定（ＮＶＴ）アンサンブルの分子動力学シミュレーションを１０ピコ秒、１ａｔｍ／２０℃のＮＰＴアンサンブルの分子動力学シミュレーションを９０ピコ秒、１ａｔｍの分子動力学シミュレーションにおいて温度を２０℃〜５００℃まで２００ピコ秒周期で上昇下降を繰り返し緩和させるシミュレーテッドアニーリングを１４サイクル（２８００ピコ秒に相当）、１ａｔｍ／２０℃のＮＰＴアンサンブルの分子動力学シミュレーションを３００ピコ秒実施した。

分子動力学シミュレーション終了後、非結合性相互作用のスケール値λの違いによるPDMS構造の変化を明らかにするため、分子動力学シミュレーションで得られた原子の座標rの時系列データとMaterials Studio（登録商標)7.0 Atom Volumes And Surfaces toolを用いてConnolly Surface法[M. L. Connolly, 1983, Science, 221, p.709.]により自由体積を求め、数式（１２）で示される系の空隙率を計算した。

[実施例２]
実施例１の計算条件ファイルにおいて低分子‐高分子間非結合性相互作用のスケール値λ_low-poly＝0.3としてメモリに格納し、それ以外は実施例１と同条件で分子動力学シミュレーションを実施して、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞の平均値、密度の平均値、空隙率を計算した。

[実施例３]
実施例１の計算条件ファイルにおいて低分子‐高分子間非結合性相互作用のスケール値λ_low-poly＝0.4としてメモリに格納し、それ以外は実施例１と同条件で分子動力学シミュレーションを実施して、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞の平均値、密度の平均値、空隙率を計算した。

[実施例４]
実施例１の計算条件ファイルにおいて低分子‐高分子間非結合性相互作用のスケール値λ_low-poly＝0.5としてメモリに格納し、それ以外は実施例１と同条件で分子動力学シミュレーションを実施して、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞の平均値、密度の平均値、空隙率を計算した。

[実施例５]
実施例１の計算条件ファイルにおいて低分子‐高分子間非結合性相互作用のスケール値λ_low-poly＝0.6としてメモリに格納し、それ以外は実施例１と同条件で分子動力学シミュレーションを実施して、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞の平均値、密度の平均値、空隙率を計算した。

[実施例６]
実施例１の計算条件ファイルにおいて低分子‐高分子間非結合性相互作用のスケール値λ_low-poly＝0.8としてメモリに格納し、それ以外は実施例１と同条件で分子動力学シミュレーションを実施して、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞の平均値、密度の平均値、空隙率を計算した。

[比較例１]
実施例１の計算条件ファイルにおいて低分子‐高分子間非結合性相互作用のスケール値λ_low-poly＝0.1としてメモリに格納し、ステップ１およびステップ２を実施した後、ステップ３を実施せず、ステップ４およびステップ５を実施し、以後ステップ１〜２、ステップ４〜５を１６００００００回繰り返し実施し、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞の平均値、密度の平均値、空隙率を計算した。

ステップ３を実施しないという以外は実施例１と同条件であり、この操作はカットオフ半径から無限遠までのファンデルワールス相互作用から生じる圧力を無視することに相当する。

[比較例２]
比較例１の計算条件ファイルにおいて低分子‐高分子間非結合性相互作用のスケール値λ_low-poly＝0.3としてメモリに格納し、それ以外は比較例１と同条件で分子動力学シミュレーションを実施し、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞の平均値、密度の平均値、空隙率を計算した。

[比較例３]
比較例１の計算条件ファイルにおいて低分子‐高分子間非結合性相互作用のスケール値λ_low-poly＝0.4としてメモリに格納し、それ以外は比較例１と同条件で分子動力学シミュレーションを実施し、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞の平均値、密度の平均値、空隙率を計算した。

[比較例４]
比較例１の計算条件ファイルにおいて低分子‐高分子間非結合性相互作用のスケール値λ_low-poly＝0.5としてメモリに格納し、それ以外は比較例１と同条件で分子動力学シミュレーションを実施し、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞の平均値、密度の平均値、空隙率を計算した。

[比較例５]
比較例１の計算条件ファイルにおいて低分子‐高分子間非結合性相互作用のスケール値λ_low-poly＝0.6としてメモリに格納し、それ以外は比較例１と同条件で分子動力学シミュレーションを実施し、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞の平均値、密度の平均値、空隙率を計算した。

[比較例６]
比較例１の計算条件ファイルにおいて低分子‐高分子間非結合性相互作用のスケール値λ_low-poly＝0.8としてメモリに格納し、それ以外は比較例１と同条件で分子動力学シミュレーションを実施し、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞の平均値、密度の平均値、空隙率を計算した。

実施例１〜６、比較例１〜６で得られたカットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞の平均値、密度の平均値、空隙率、および同じ低分子‐高分子間非結合性相互作用のスケール値λ_low-polyにおける実施例と比較例の空隙率の比を表３にまとめた。また低分子‐高分子間非結合性相互作用のスケール値λ_low-polyに対する密度および空隙率を図３、４に示した。

PDMSの密度の実験値は0.971g/cm³[I. Vinckier, P. Moldenaers, J. Mewis, Journal of Rheology, 1996, 40, p.613.]である。図３より実施例は実験値と同程度の密度を持つが、比較例においては密度が小さくなることがわかる。さらに図４より比較例は実施例より３割程度大きな空隙を持つことがわかる。よって圧力一定のアンサンブルでの分子動力学シミュレーションでは、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞を考慮しなければ、系の密度を過小に評価し、系の空隙率を過大に評価することになる。

上記実施例および比較例では計算時間およびvdW相互作用エネルギーの計算精度の関係から一般的に好ましく用いられるカットオフ半径ｒ_ｃ＝10オングストロームを用いた。

しかしカットオフ半径を大きくすれば計算時間が増加するが、カットオフ半径内のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^0→ｒｃの寄与が増加するため、Ｐ_vdW ^0→ｒｃがＰ_vdW ^0→ｒｃ+Ｐ_vdW ^rc→∞に漸近する。

[比較例７]
計算条件ファイルにおいてカットオフ半径ｒ_ｃ＝12オングストロームとする以外は比較例５と同条件で分子動力学シミュレーションを実施し、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞の平均値、密度の平均値、空隙率を計算した。

[比較例８]
計算条件ファイルにおいてカットオフ半径ｒ_ｃ＝15オングストロームとする以外は比較例５と同条件で分子動力学シミュレーションを実施し、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞の平均値、密度の平均値、空隙率を計算した。

[比較例９]
計算条件ファイルにおいてカットオフ半径ｒ_ｃ＝17オングストロームとする以外は比較例５と同条件で分子動力学シミュレーションを実施し、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞の平均値、密度の平均値、空隙率を計算した。

実施例５、比較例５、比較例７〜９において分子動力学シミュレーションプログラムを実行するにあたりその計算実行時間を、今回用いたコンピュータの標準シェルであるcshの”time”コマンドを用いて測定し出力した。

得られた分子動力学シミュレーションプログラムの計算実行時間および密度を図５に示した。カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞を考慮しない比較例５、比較例７〜９においても、カットオフ半径を17オングストローム程度まで大きくすれば、カットオフを行ったことによる影響が小さくなり実験値の密度に近づくが、計算時間も２倍弱に増加した。

一方、実施例５ではカットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐ_vdW ^rc→∞を考慮したため、短い計算時間で現実系に近い密度再現する高精度な計算を行うことができた。

[実施例７]
ポリジメチルシロキサン（PDMS）４分子／エチレングリコール（EG）１分子の系を作成し、低分子の分子内非結合性相互作用のスケール値λ_low-intra＝1.0と設定する以外実施例１と同条件で分子動力学シミュレーションを実施した。得られた原子の座標rの時系列データと結合情報からエチレングリコールのO-C-C-O二面角の平均値を計算した。

エチレングリコールの力場パラメータについては分子内ポテンシャルパラメータの力の定数K^bond、K^angle、K^torsion、K^inversionおよび平衡位置R⁰、θ⁰、φ⁰、ω⁰として、HF/6-31G(d,p)レベルの分子軌道計算によって決定したパラメータ、ＡＭＢＥＲパラメータ、hybrid/UAパラメータおよび桑島ポテンシャルパラメータを用い、ｖｄＷパラメータにはＡＭＢＥＲパラメータ、および桑島ポテンシャルパラメータを用い、電荷はB3LYP/6-31G(d,p)レベルの分子軌道計算によって決定した値を用いた。なおO-C-C-O二面角、H-O-C-C二面角パラメータは、分子軌道計算により仮パラメータを設定した後、NMRにより測定されたバルク溶液中のO-C-C-O二面角の平均値φ＝72°[G. Chidichimo, D. Imbardelli, M. Longeri, A. Saupe, Mol. Phys., 1988, 65, 1143]を再現するよう調整したものを用いた。

[実施例８]
ポリジメチルシロキサン（PDMS）４分子／エチレングリコール（EG）１分子の系について、低分子の分子内非結合性相互作用のスケール値λ_low-intra＝0.1とする以外実施例７と同条件で分子動力学シミュレーションを実施し、エチレングリコールのO-C-C-O二面角の平均値を計算した。

実施例７、実施例８で計算したエチレングリコールのO-C-C-O二面角の平均値およびNMR測定による実験値を表４にまとめた。表４より実施例７はバルク状態の実験値に近い値でありゴーシュ体が支配的である。バルク中とポリマー中のエチレングリコールのコンフォメーションは必ずしも一致する必要はないが、実施例７、実施例８におけるPDMSのようにエチレングリコールに対し特異的に相互作用する官能基を有していない場合は、バルク状態のコンフォメーションとの差異はあまり大きくないと考えられるため、O-C-C-O二面角の平均値はバルク中の実験値と近いはずである。

一方、実施例８では実験値と異なりトランス体が支配的である。これは低分子の分子内非結合性相互作用のスケール値λ_low-intra＝0.1としたことにより、エチレングリコールの分子内水素結合が弱まり、分子内のコンフォメーションが崩れてしまったためである。

よってエチレングリコールなど低分子内の非結合性相互作用がコンフォメーションに大きく寄与している低分子では、低分子の分子内非結合性相互作用のスケール値λ_low-intra＝１.0とし、低分子内非結合性相互作用をスケールしないことが好ましいといえる。

[実施例９]
実施例１の計算条件ファイルにおいて低分子‐高分子間非結合性相互作用のスケール値λ_low-poly＝0.8としてメモリに格納し、それ以外は実施例１と同条件で分子動力学シミュレーションを実施して、ステップ１〜ステップ５のループを４０００回繰り返すごとにＰＤＭＳ−水間の非結合相互作用V_low-poly、原子の座標を出力した。同様の手順の計算を独立な７種類の初期配置について実施し、初めに計算を行った初期配置を含め８種類の独立な配置から得られたＰＤＭＳ−水間の非結合相互作用V_low-polyの平均値を計算した。またcshの”time”コマンドを用い、８種類の独立な配置の計算が全て終了するまでの時間を測定した。

[実施例１０]
ポリジメチルシロキサン（PDMS）４分子／水１分子の系について、図２のフローに従い分子動力学シミュレーションを実施し、様々なλ_low-polyにおけるＰＤＭＳ−水間の非結合相互作用V_low-polyをサンプリングした。

実施例１と同様に、表１に示した組成で、ＰＤＭＳ分子と水分子を三次元周期境界条件下でランダムに配置し、系を作成した。Materials Studio（登録商標）7.0の外部出力機能を使用し、作成した系の原子の座標r 、結合情報、原子のポテンシャルパラメータタイプ、原子の質量ｍ、基本セルベクトルhが記載された分子構造ファイルを作成した。

次に、テキストエディタを用いて分子構造ファイルに記載されたポテンシャルパラメータタイプに対応したポテンシャルパラメータを記述したパラメータファイルを作成した。

さらにテキストエディタを用いて非結合性相互作用のスケール値λ、λ_low-polyの増減間隔Δλ、Δλを加減する間隔Δｎ、相互作用計算のカットオフ半径ｒ_ｃ、図２計算フローにおけるループの繰り返しを規定する回数n、出力する物理量とその出力間隔、設定温度T₀、設定圧力P₀、熱浴の慣性量Q、barostatの慣性量Wを記述した計算条件ファイルを作成した。

作成した分子構造ファイル、パラメータファイル、計算条件ファイルを分子動力学シミュレーションプログラムで読み取らせ、データをメモリに格納した。

計算条件ファイルにおいて、非結合性相互作用のスケール値λを相互作用の種類ごとに設定した。低分子‐高分子間非結合性相互作用のスケール値の初期値をλ_low-poly＝１、高分子の分子内非結合性相互作用と、高分子‐高分子間非結合性相互作用のそれぞれのスケール値λ_poly-intra＝λ_poly-poly＝１とした。実施例１同様に低分子の分子内非結合性相互作用のスケール値λ_low-intra および低分子‐低分子間非結合性相互作用のスケール値λ_low-lowは設定しない。λ_low-polyの増減間隔Δλ＝0.1とし、Δλを加減する間隔Δｎ＝４０００回とした。さらに、カットオフ半径ｒ_ｃ＝１０オングストローム、図２計算フローにおけるループの繰り返しを規定する回数ｎ＝１６００００００回とした。出力する物理量をＰＤＭＳ−水間の非結合相互作用V_low-poly、原子の座標rとしその出力間隔を図２のステップ１〜ステップ５（Mod(n,Δn)=0のときはステップ６も含む）のループを１０００回繰り返すごとに１回とした。

実施例１同様、温度制御にNose−Hoover法、圧力制御にParrinello−Rahman法を用いて、原子数・温度・圧力一定（ＮＰＴ）アンサンブルを構成することにより、系を設定温度T₀、設定圧力P₀に制御した状態で図２ステップ１〜ステップ６を実施した。なお、設定温度T₀=20℃、設定圧力P₀=1atm、熱浴の慣性量Q=120.95(kJ・ps²)/mol、barostatの慣性量W=483.79(kJ・ps²)/molとした。

また原子の速度vの初期値としてT₀=20℃におけるマクスウェル分布に従うよう計算された速度を、熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルの速度hvの初期値として全ての成分に０を与えメモリに格納した。ステップ１〜ステップ５は実施例１と同様に実施した。ステップ１〜ステップ５を４０００回繰り返す毎に、下記ステップ６を１回実行した。

〔ステップ６〕
メモリに格納されたλ_low-polyに増減間隔Δλ（０．１）を次の規則に従い加減した。

〔Δλ加減規則〕
１．λ_low-poly＝０となるまで、ステップ６に到達する毎にΔλを引く。
２．λ_low-poly＝０に到達したら、λ_low-poly＝１となるまでステップ６に到達する毎にΔλを加える。
３．λ_low-poly＝1に到達したら、λ_low-poly＝０となるまでステップ６に到達する毎にΔλを引く。
４．上記１〜３の操作を、繰り返しの規定回数ｎに到達するまで繰り返す。

このようにしてステップ６で更新されたλ_low-polyを用いて再びステップ１〜ステップ５を実施した。繰り返しの規定回数に到達後、λ_low-poly＝0.8におけるＰＤＭＳ−水間の非結合相互作用V_low-polyの平均値を計算した。またcshの”time”コマンドを用い、計算が終了するまでの時間を測定した。

実施例９、実施例１０で計算したＰＤＭＳ−水間の非結合相互作用V_low-polyおよび計算実行時間を表５にまとめた。表５より実施例９、実施例１０どちらのＰＤＭＳ−水間の非結合相互作用V_low-polyもほぼ同じ精度で計算できていることがわかる。しかし実施例１０は実施例９の約１／８の計算時間で終了しており、効率が非常に高いといえる。

Claims (9)

1. コンピュータを利用した高分子と低分子の混合系の分子動力学シミュレーションにおいて、下記数式（２）のように、
   （i）低分子の分子内ｖｄＷ相互作用V^vdw _low-intra(r, λ_low-intra)と、
   （ii）高分子の分子内ｖｄＷ相互作用V^vdw _poly-intra(r, λ_poly-intra)と、
   （iii）低分子‐低分子間ｖｄＷ相互作用V^vdw _low-low(r, λ_low-low)と、
   （iv）高分子‐高分子間ｖｄＷ相互作用V^vdw _poly-poly(r, λ_poly-poly)と、
   （v）低分子‐高分子間ｖｄＷ相互作用V^vdw _low-poly(r, λ_low-poly)
   との和（0≦λ≦１、ただし（i）から（v）のλのうち少なくとも１つは0≦λ<１）で表される、非結合性相互作用のスケール値λを含むスケールファンデルワールス（vdW）相互作用ポテンシャルV^vdW(r, λ)を用いた場合の前記混合系のファンデルワールス（vdW）相互作用から生じる圧力の計算方法であって、
   

   

   コンピュータに、
   ・前記高分子と前記低分子に含まれる原子の座標r
   および、必要に応じ
   ・原子の速度ｖ
   の初期値、
   ならびに
   ・結合情報
   ・原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ
   ・原子の質量ｍ
   ・非結合性相互作用のスケール値λ
   ・相互作用計算のカットオフ半径ｒ_ｃ
   を入力した後、コンピュータに、
   （ステップ１）原子の座標r、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、非結合性相互作用のスケール値λ、ならびに下記数式（３）で示したポテンシャル関数V^totalを原子の座標rで微分した式を用いて、前記高分子および前記低分子の分子内力(数式（３）第１項から第４項の微分の和)および分子間力（数式（３）第５項と第６項の微分の和）を計算するステップ；
   

   

   （ただし、R、θ、φ、ωはそれぞれ注目する結合の距離、変角の角度、二面角の角度、反転の角度を表す。）
   （ステップ２）ステップ１で計算した分子間力の一つであるファンデルワールス（vdW）力（数式（３）第５項の微分）と、原子の座標rからビリアル定理を用いて圧力Ｐ_vdW ^0→ｒｃを計算するステップ；
   （ステップ３）原子の座標r、原子のポテンシャルパラメータタイプ、ポテンシャルパラメータ、非結合性相互作用のスケール値λおよび相互作用計算のカットオフ半径ｒ_ｃから、下記数式(４Ａ)および（４Ｂ）を用いて圧力Ｐ_vdW ^rc→∞を計算するステップ；
   

   

   （ただし、Ｍはポテンシャルパラメータタイプの総数、Ｎ_α、Ｎ_βはそれぞれポテンシャルパラメータタイプα、βを持つ原子の数、diag[・・・]は対角行列を表す。）
   （ステップ４）下記数式（５）によりファンデルワールス（vdW）相互作用から生じる圧力Ｐ_vdW ^0→∞を計算するステップ；
   

   

   （ステップ５）原子の質量ｍ、分子内力、分子間力を用いて、微小時間Δt後の
   ・原子の座標r
   ・原子の速度ｖ
   を数値積分によって更新するステップ：
   を実行させ、
   さらに、ステップ５で更新した微小時間Δt後の原子の座標r、原子の速度ｖを初期値としてステップ１〜ステップ５を再度実行させ、これを繰り返すことを含むファンデルワールス（vdW）相互作用から生じる圧力の計算方法。
2. ステップ１において数式（３）におけるスケールされたvdW相互作用ポテンシャルV^vdW(r, λ)が数式（８）に示したDreiding型レナードジョーンズポテンシャルであり、数式（４Ａ）、数式（１１Ａ）および（１１Ｂ）を用いてＰ_vdW ^rc→∞を計算する、請求項１に記載の圧力の計算方法。
   

   

   （ここで注目する二つの原子をi,jとし、δはエネルギーと力の距離依存性を調節するパラメータ、D^０、R^０はvdWポテンシャルの深さと斥力が生じる位置を決めるパラメータ、r² _ijはij間の距離とし、0≦λ_ij≦1、δ≧0ｍ^２とする。）
   

   
3. 非結合性相互作用のスケール値をλ_low-intra＝λ_poly-intra＝λ_low-low＝λ_poly-poly＝１、0≦λ_low-poly＜1とすることを特徴とする請求項１または２に記載の圧力の計算方法。
4. 請求項１〜３のいずれかに記載の圧力の計算方法を用いて圧力を計算することを含む高分子-低分子間の相互作用エネルギーまたは高分子中の低分子の配置のサンプリング方法。
5. 請求項１または２に記載の圧力の計算方法を用いて圧力を計算することを含む高分子-低分子間の相互作用エネルギーまたは高分子中の低分子の配置のサンプリング方法であって、λ_low-intra＝λ_poly-intra＝λ_low-low＝λ_poly-poly＝１とし、所定回数のステップ１〜ステップ５の繰り返しに１回、さらに
   （ステップ６）λ_low-polyを0≦λ_low-poly≦1の範囲で増減するステップ；
   を実行させ、
   さらに、ステップ５で更新した微小時間Δt後の原子の座標rおよび原子の速度ｖ、ならびにステップ６で増減させたλ_low-polyを初期値としてステップ１〜ステップ５を再度実行させ、これを繰り返すことを含む、高分子-低分子間の相互作用エネルギーまたは高分子中の低分子の配置のサンプリング方法。
6. ステップ１〜５、および所定回数に１回のステップ１〜ステップ６の繰り返しにおいて、ステップ６の実行時に１回以上λ_low-poly＝０およびλ_low-poly＝１となるよう設定する、請求項５に記載の高分子-低分子間の相互作用エネルギーまたは高分子中の低分子の配置のサンプリング方法。
7. ステップ６におけるλ_low-polyの増減間隔Δλが０．０１以上０．３以下である、請求項５または６に記載の高分子-低分子間の相互作用エネルギーまたは高分子中の低分子の配置のサンプリング方法。
8. 請求項１〜３のいずれかに記載の圧力の計算方法または、請求項４〜７のいずれかに記載の高分子-低分子間の相互作用エネルギーまたは高分子中の低分子の配置のサンプリング方法をコンピュータに実行させるためのプログラム。
9. 請求項８に記載のプログラムをコンピュータに読み取り可能に記録した記録媒体。
   

Images