(57)【要約】【課題】  高分子と低分子の混合系における、スケールされた非結合性相互を用いた分子動力学シミュレーションにおいて、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞を求め、正確な圧力の計算方法を提供する。【解決手段】  コンピュータに、原子の座標 、スケール値など初期条件を入力し、（１）分子内力・分子間力を計算し、（２）カットオフ半径内のファンデルワールス相互作用から生じる圧力ＰvdW0→rcを計算し、（３）カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞を計算し、（４）ＰvdW0→rcとＰvdWrc→∞との和からＰvdW0→∞計算し、（５）微小時間Δt後の原子の座標、速度を数値積分によって更新し、さらに、（５）の更新値を初期値として（１）〜（５）を再度実行させ、これを規定回数繰り返すことを含む圧力の計算方法。【選択図】図１
2016181257.tif 000002
【特許請求の範囲】
【請求項１】
コンピュータを利用した高分子と低分子の混合系の分子動力学シミュレーションにおいて、下記数式（２）のように、
（i）低分子の分子内ｖｄＷ相互作用Vvdwlow-intra(r, λlow-intra)と、
（ii）高分子の分子内ｖｄＷ相互作用Vvdwpoly-intra(r, λpoly-intra)と、
（iii）低分子‐低分子間ｖｄＷ相互作用Vvdwlow-low(r, λlow-low)と、
（iv）高分子‐高分子間ｖｄＷ相互作用Vvdwpoly-poly(r, λpoly-poly)と、
（v）低分子‐高分子間ｖｄＷ相互作用Vvdwlow-poly(r, λlow-poly)
との和（0≦λ≦１、ただし（i）から（v）のλのうち少なくとも１つは0≦λ<１）で表される、非結合性相互作用のスケール値λを含むスケールファンデルワールス（vdW）相互作用ポテンシャルVvdW(r, λ)を用いた場合の前記混合系のファンデルワールス（vdW）相互作用から生じる圧力の計算方法であって、
【数１】
2016181257.tif 000021
コンピュータに、
・前記高分子と前記低分子に含まれる原子の座標r 
および、必要に応じ
・原子の速度ｖ
の初期値、
ならびに
・結合情報
・原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ
・原子の質量ｍ
・非結合性相互作用のスケール値λ
・相互作用計算のカットオフ半径ｒｃ
を入力した後、コンピュータに、
（ステップ１）原子の座標r、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、非結合性相互作用のスケール値λ、ならびに下記数式（３）で示したポテンシャル関数Vtotalを原子の座標rで微分した式を用いて、前記高分子および前記低分子の分子内力(数式（３）第１項から第４項の微分の和)および分子間力（数式（３）第５項と第６項の微分の和）を計算するステップ；
【数２】
2016181257.tif 000022
（ただし、R、θ、φ、ωはそれぞれ注目する結合の距離、変角の角度、二面角の角度、反転の角度を表す。）
（ステップ２）ステップ１で計算した分子間力の一つであるファンデルワールス（vdW）力（数式（３）第５項の微分）と、原子の座標rからビリアル定理を用いて圧力ＰvdW0→ｒｃを計算するステップ；
（ステップ３）原子の座標r、原子のポテンシャルパラメータタイプ、ポテンシャルパラメータ、非結合性相互作用のスケール値λおよび相互作用計算のカットオフ半径ｒｃから、下記数式(４Ａ)および（４Ｂ）を用いて圧力ＰvdWrc→∞を計算するステップ；
【数３】
2016181257.tif 000023
（ただし、Ｍはポテンシャルパラメータタイプの総数、Ｎα、Ｎβはそれぞれポテンシャルパラメータタイプα、βを持つ原子の数、diag[・・・]は対角行列を表す。）
（ステップ４）下記数式（５）によりファンデルワールス（vdW）相互作用から生じる圧力ＰvdW0→∞を計算するステップ；
【数４】
2016181257.tif 000024
（ステップ５）原子の質量ｍ、分子内力、分子間力を用いて、微小時間Δt後の
・原子の座標r
・原子の速度ｖ
を数値積分によって更新するステップ：
を実行させ、
さらに、ステップ５で更新した微小時間Δt後の原子の座標r、原子の速度ｖを初期値としてステップ１〜ステップ５を再度実行させ、これを繰り返すことを含むファンデルワールス（vdW）相互作用から生じる圧力の計算方法。
【請求項２】
ステップ１において数式（３）におけるスケールされたvdW相互作用ポテンシャルVvdW(r, λ)が数式（８）に示したDreiding型レナードジョーンズポテンシャルであり、数式（４Ａ）、数式（１１Ａ）および（１１Ｂ）を用いてＰvdWrc→∞を計算する、請求項１に記載の圧力の計算方法。
【数５】
2016181257.tif 000025
（ここで注目する二つの原子をi,jとし、δはエネルギーと力の距離依存性を調節するパラメータ、D０、R０はvdWポテンシャルの深さと斥力が生じる位置を決めるパラメータ、r2ijはij間の距離とし、0≦λij≦1、δ≧0ｍ２とする。）
【数６】
2016181257.tif 000026
【請求項３】
非結合性相互作用のスケール値をλlow-intra＝λpoly-intra＝λlow-low＝λpoly-poly＝１、0≦λlow-poly＜1とすることを特徴とする請求項１または２に記載の圧力の計算方法。
【請求項４】
請求項１〜３のいずれかに記載の圧力の計算方法を用いて圧力を計算することを含む高分子-低分子間の相互作用エネルギーまたは高分子中の低分子の配置のサンプリング方法。
【請求項５】
請求項１または２に記載の圧力の計算方法を用いて圧力を計算することを含む高分子-低分子間の相互作用エネルギーまたは高分子中の低分子の配置のサンプリング方法であって、λlow-intra＝λpoly-intra＝λlow-low＝λpoly-poly＝１とし、所定回数のステップ１〜ステップ５の繰り返しに１回、さらに
（ステップ６）λlow-polyを0≦λlow-poly≦1の範囲で増減するステップ；
を実行させ、
さらに、ステップ５で更新した微小時間Δt後の原子の座標rおよび原子の速度ｖ、ならびにステップ６で増減させたλlow-polyを初期値としてステップ１〜ステップ５を再度実行させ、これを繰り返すことを含む、高分子-低分子間の相互作用エネルギーまたは高分子中の低分子の配置のサンプリング方法。
【請求項６】
ステップ１〜５、および所定回数に１回のステップ１〜ステップ６の繰り返しにおいて、ステップ６の実行時に１回以上λlow-poly＝０およびλlow-poly＝１となるよう設定する、請求項５に記載の高分子-低分子間の相互作用エネルギーまたは高分子中の低分子の配置のサンプリング方法。
【請求項７】
ステップ６におけるλlow-polyの増減間隔Δλが０．０１以上０．３以下である、請求項５または６に記載の高分子-低分子間の相互作用エネルギーまたは高分子中の低分子の配置のサンプリング方法。
【請求項８】
請求項１〜３のいずれかに記載の圧力の計算方法または、請求項４〜７のいずれかに記載の高分子-低分子間の相互作用エネルギーまたは高分子中の低分子の配置のサンプリング方法をコンピュータに実行させるためのプログラム。
【請求項９】
請求項８に記載のプログラムをコンピュータに読み取り可能に記録した記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、コンピュータを利用した分子動力学シミュレーションにおける、高分子中に低分子が存在する系におけるファンデルワールス相互作用から生じる圧力の計算方法および、高分子-低分子間の相互作用エネルギーまたは高分子中の低分子の配置のサンプリング方法に関する。
【背景技術】
【０００２】
分子動力学シミュレーションは、分子や分子集合体の高次構造や物性を詳細に解析するための手法として広く用いられる。分子動力学シミュレーションは、近年のコンピュータの発達により、合成高分子や生体高分子といったソフトマターの研究にもよく適用されるようになってきた。
【０００３】
合成高分子の分子動力学シミュレーションの適用例として、高分子分離膜における物質透過性予測が挙げられる（非特許文献１）。物質の膜透過性を予測するにあたり重要な物理量としては、透過物質である低分子の高分子膜への溶解性を示す溶解度係数がある。溶解度係数を分子動力学シミュレーションから求めるためには、低分子と高分子の相互作用エネルギーを高精度に計算する必要がある。しかし高分子の遅い分子運動のため、高分子中の低分子は通常の分子動力学シミュレーションの計算時間ではほとんど移動しない。そのため限定された空間の相互作用エネルギーしかサンプリングできず、統計力学的に高精度に溶解度係数を求めることができない。
【０００４】
低分子の移動を促すため、低分子と高分子間の非結合相互作用を時間とともにスケールし非結合性相互作用を増減させるという方法で、低分子の運動性を向上させ統計精度を高めた高精度な相互作用エネルギー計算手法が考案された（特許文献１）。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】特開2011-123874号公報
【非特許文献】
【０００６】
【非特許文献１】Y. Tamai, H. Tanaka, and K. Nakanishi, Macromolecules, 28, 2544-2554(1995)
【発明の概要】
【発明が解決しようとする課題】
【０００７】
分子動力学シミュレーションにおいて、系の全圧力Ptotalは、下記数式（１）のように表される。
【０００８】
【数１】
2016181257.tif 000003
【０００９】
（ただし、Pintraは分子内相互作用から生じる圧力、PvdW0→∞はファンデルワールス相互作用から生じる圧力、Pcoulomb0→∞はクーロン相互作用から生じる圧力、PKは原子の運動エネルギーより生じる圧力である。）
ファンデルワールス相互作用から生じる圧力ＰvdW0→∞は、カットオフ半径内のファンデルワールス相互作用から生じる圧力ＰvdW0→rcとカットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞の和で表される。
【００１０】
ここで、特許文献１のようにスケールした非結合性相互作用を使用した場合、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞を計算できないという問題があった、すなわち、ファンデルワールス相互作用から生じる圧力ＰvdW0→∞を正確に求める手法は未だ提唱されておらず、そのため、系の全圧力Ｐtotalも正確に求めることができなかった。
【００１１】
一方、十分に大きなカットオフ半径を用いれば圧力ＰvdWrc→∞の項は０に漸近するため、正確な圧力に近い値を得ることができる。しかし計算時間の観点から一般的に好適に用いられるカットオフ半径を用いた時に比べ長時間の計算を要するため、大きなカットオフ半径を用いることは非効率的である。
【００１２】
正確な圧力を計算できないことによる問題点として、分子動力学シミュレーションにおける現実系を再現するために行う、圧力を一定に制御した状態のシミュレーションを正しく行えないことが挙げられる。圧力を制御するために、ある状態における系の瞬間の圧力を計算し、その圧力が設定圧力に近づくようにシミュレーションセルを変形させる方法が用いられるが、正確な圧力を計算できない場合、シミュレーションセルの大きさを過大評価あるいは過小評価し、その結果、密度や空隙率が正しく計算できず、正確に相互作用エネルギーを求めることができなかった。
【００１３】
以上述べたように、高分子と低分子の混合系において、スケールされた非結合性相互を用いた分子動力学シミュレーションでは、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞を求め、正確な圧力を計算する方法は未だ提唱されていない。そこで本発明では、上記問題を解決したスケールされた非結合性相互を用いた場合のファンデルワールス相互作用から生じる圧力の計算方法を提供することを目的とする。
【課題を解決するための手段】
【００１４】
本発明は上記課題を解決するために、次の手段を採用するものである。すなわち、本発明は以下のとおりである。
【００１５】
コンピュータを利用した高分子と低分子の混合系の分子動力学シミュレーションにおいて、下記数式（２）のように、
（i）低分子の分子内ｖｄＷ相互作用Vvdwlow-intra(r, λlow-intra)と、
（ii）高分子の分子内ｖｄＷ相互作用Vvdwpoly-intra(r, λpoly-intra)と、
（iii）低分子‐低分子間ｖｄＷ相互作用Vvdwlow-low(r, λlow-low)と、
（iv）高分子‐高分子間ｖｄＷ相互作用Vvdwpoly-poly(r, λpoly-poly)と、
（v）低分子‐高分子間ｖｄＷ相互作用Vvdwlow-poly(r, λlow-poly)
との和（0≦λ≦１、ただし（i）から（v）のλのうち少なくとも１つは0≦λ<１）で表される、非結合性相互作用のスケール値λを含むスケールファンデルワールス（vdW）相互作用ポテンシャルVvdW(r, λ)を用いた場合の前記混合系のファンデルワールス（vdW）相互作用から生じる圧力の計算方法であって、
【００１６】
【数２】
2016181257.tif 000004
【００１７】
コンピュータに、
・前記高分子と前記低分子に含まれる原子の座標r 
および、必要に応じ
・原子の速度ｖ
の初期値、
ならびに
・結合情報
・原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ
・原子の質量ｍ
・非結合性相互作用のスケール値λ
・相互作用計算のカットオフ半径ｒｃ
を入力した後、コンピュータに、
（ステップ１）原子の座標r、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、非結合性相互作用のスケール値λ、ならびに下記数式（３）で示したポテンシャル関数Vtotalを原子の座標rで微分した式を用いて、前記高分子および前記低分子の分子内力(数式（３）第１項から第４項の微分の和)および分子間力（数式（３）第５項と第６項の微分の和）を計算するステップ；
【００１８】
【数３】
2016181257.tif 000005
【００１９】
（ただし、R、θ、φ、ωはそれぞれ注目する結合の距離、変角の角度、二面角の角度、反転の角度を表す。）
（ステップ２）ステップ１で計算した分子間力の一つであるファンデルワールス（vdW）力（数式（３）第５項の微分）と、原子の座標rからビリアル定理を用いて圧力ＰvdW0→ｒｃを計算するステップ；
（ステップ３）原子の座標r、原子のポテンシャルパラメータタイプ、ポテンシャルパラメータ、非結合性相互作用のスケール値λおよび相互作用計算のカットオフ半径ｒｃから、下記数式(４Ａ)および（４Ｂ）を用いて圧力ＰvdWrc→∞を計算するステップ；
【００２０】
【数４】
2016181257.tif 000006
【００２１】
（ただし、Ｍはポテンシャルパラメータタイプの総数、Ｎα、Ｎβはそれぞれポテンシャルパラメータタイプα、βを持つ原子の数、diag[・・・]は対角行列を表す。）
（ステップ４）下記数式（５）によりファンデルワールス（vdW）相互作用から生じる圧力ＰvdW0→∞を計算するステップ；
【００２２】
【数５】
2016181257.tif 000007
【００２３】
（ステップ５）原子の質量ｍ、分子内力、分子間力を用いて、微小時間Δt後の
・原子の座標r
・原子の速度ｖ
を数値積分によって更新するステップ：
を実行させ、
さらに、ステップ５で更新した微小時間Δt後の原子の座標r、原子の速度ｖを初期値としてステップ１〜ステップ５を再度実行させ、これを繰り返すことを含むファンデルワールス（vdW）相互作用から生じる圧力の計算方法。
【００２４】
また、本発明は、上記の圧力の計算方法を用いた高分子-低分子間の相互作用エネルギーまたは高分子中の低分子の配置のサンプリング方法を含む。
【図面の簡単な説明】
【００２５】
【図１】本発明の圧力の計算方法を示すフローチャートである。
【図２】本発明の圧力の計算方法を利用した高分子-低分子間の相互作用エネルギーや高分子中の低分子の配置のサンプリング方法を示すフローチャートである。
【図３】実施例１〜６、比較例１〜６における非結合性相互作用のスケール値と密度の関係を示す図である。
【図４】実施例１〜６、比較例１〜６における非結合性相互作用のスケール値と空隙率の関係を示す図である。
【図５】実施例５、比較例５、比較例７〜９におけると計算実行時間と密度、カットオフ半径の関係を示す図である。
【発明を実施するための形態】
【００２６】
図１に本発明の計算方法のフローを示す。以下、本発明の好ましい実施形態を、図１のフローに従って説明する。
【００２７】
本発明の計算方法においては、まず、コンピュータに、計算に必要な諸条件を入力する。計算に必要な諸条件を入力するに際しては、コンピュータに読み取り可能な形態で、高分子と低分子に含まれる原子の座標rの初期値、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、原子の質量ｍ、非結合性相互作用のスケール値λ、相互作用計算のカットオフ半径ｒｃを入力する。また、必要に応じて原子の速度ｖの初期値を入力することもできるが、入力しなければマクスウェル分布に従うよう計算された原子の速度が初期速度として与えられるよう構成すればよい。
【００２８】
入力方法は特に限定されないが、一般に入手可能な分子モデル作成ソフトウェアやテキストエディタなどが利用でき、当該情報が記述されたファイルをコンピュータに読み取らせ、メモリに格納することが好ましい。
【００２９】
〔ステップ１〕
高分子および低分子の分子内力(数式（３）第１項から第４項の微分の和)および分子間力（数式（３）第５項と第６項の微分の和）を計算するステップ１においては、メモリに格納された原子の座標r、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、非結合性相互作用のスケール値λ、ならびに数式（３）で示したポテンシャル関数Vtotalを原子の座標rで微分した式を用いて分子内力・分子外力を計算し、メモリに格納する。
【００３０】
ここで、数式（３）の右辺第１項Vbond(R)は結合長、第２項Vangle(θ)は結合角、第３項Vtorsion(φ)は二面角、第４項Vinversion(ω)は反転、第５項Vvdw(r,λ)はファンデルワールス（ｖｄＷ）、第６項Vcoulomb(r,λ)はクーロンのそれぞれポテンシャル関数を表わす。なお、数式（３）の右辺第１項から第４項を総称して分子内ポテンシャル関数と呼び、右辺第５項、第６項をまとめて分子間ポテンシャル関数と呼ぶ。また、分子内ポテンシャル関数、分子間ポテンシャル関数を原子の位置で微分した式を用いて計算した力をそれぞれ分子内力、分子間力と呼ぶ。
【００３１】
数式（３）における各ポテンシャル関数の具体的な表式は下記数式（６）で表される。
【００３２】
【数６】
2016181257.tif 000008
【００３３】
数式（６）の右辺第５項、第６項はスケールされた非結合性相互作用を表し、非結合性相互作用のスケール値λ=1のときスケールされていない一般的な分子動力学シミュレーションで使用するポテンシャル関数に等しい。
【００３４】
スケールされたvdW相互作用である第5項Vvdw(r,λ)の具体的な関数形としては、下記数式（７）に示したレナードジョーンズポテンシャルやモースポテンシャル、バッキンガムポテンシャルが好ましく用いられる。
【００３５】
【数７】
2016181257.tif 000009
【００３６】
数式（７）において、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型でも好適に用いることができる。
【００３７】
【数８】
2016181257.tif 000010
【００３８】
上記数式（６）中の結合長、結合角、二面角、反転おけるそれぞれの力の定数Kbond、Kangle、Ktorsion、Kinversionおよびそれぞれの平衡位置R0、θ0、φ0、ω0、数式（７）中の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.〕が好適に用いられる。
【００４２】
上記数式（３）および数式（６）の第５項Vvdw(r,λ)で表わされるvdWポテンシャルのパラメータについては、密度と蒸発熱の実験値を再現するように決定することが好ましい。
【００４３】
上記数式（３）および数式（６）の第６項Vcoulomb(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）低分子の分子内非結合性相互作用Vlow-intra(r, λlow-intra)、
（ii）高分子の分子内非結合性相互作用Vpoly-intra(r, λpoly-intra)、
（iii）低分子−低分子間非結合相互作用Vlow-low(r, λlow-low)、
（iv）高分子‐高分子間非結合性相互作用Vpoly-poly(r, λpoly-poly)、
（v）低分子‐高分子間非結合性相互作用Vlow-poly(r, λlow-poly)の五種類の和で表される。
【００４６】
【数９】
2016181257.tif 000011
【００４７】
非結合性相互作用のスケール値λの範囲として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ｍ2とし、かつ非結合性相互作用のスケール値λが非常に小さい値をとる場合に、数式（３）および数式（６）の第５項Vvdw(r,λ)で表わされるファンデルワールス（ｖｄＷ）相互作用が急激に変化し、数値積分の安定性が低下する。この傾向は原子間の距離が近距離となった場合に特に顕著になる。この問題を解決するために、非結合性相互作用のスケール値λの最小値を非常に小さな値に設定する場合には、δ＞0ｍ2とすることが好ましい。さらに好ましくは、δ≧5.0×10-21ｍ2である。一方、δの値が大きすぎると、スケール値が変化した際、ｖｄＷ力の変化が大きくなりすぎるため、δ≦5.0×10-19ｍ2とすることが好ましい。このδの値としては、δ＝1.0×10-20ｍ2が好適に用いられる。さらにカットオフ半径rcは、計算時間と計算精度の観点から8Å≦rc≦12Åが好ましく、rc＝10Åが特に好ましい。
【００５１】
〔ステップ２〕
カットオフ半径内のファンデルワールス相互作用から生じる圧力ＰvdW0→ｒｃを計算するステップ２では、ステップ１で計算したファンデルワールス力（数式（３）第５項の微分）と、原子の座標rとから、下記数式（１０Ａ）および（１０Ｂ）で表されるビリアル定理を用いて圧力ＰvdW0→ｒｃを計算する。
【００５２】
【数１０】
2016181257.tif 000012
【００５３】
（ここでFvdWは注目する原子i、j間に働くファンデルワールス力、l, m(=1、2、3)はベクトルの成分を表す。）
〔ステップ３〕
カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞を計算するステップ３では、原子の座標r、原子のポテンシャルパラメータタイプ、ポテンシャルパラメータ、非結合性相互作用のスケール値λおよび相互作用計算のカットオフ半径ｒｃから前記数式(４Ａ)および（４Ｂ）を用いて圧力ＰvdWrc→∞計算する。
【００５４】
数式（４Ａ）の通り圧力ＰvdWrc→∞は対角項のみで、非対角項は０となる。なお数式（４Ｂ）におけるスケールされたファンデルワールス相互作用として前述の様々な関数形を用いることができる。
【００５５】
前述のように、これまで、非結合性相互作用のスケール値λを含むポテンシャルでの圧力ＰvdWrc→∞の計算方法が確立されていなかったが、本発明者らの検討の結果、数式（８）のDreiding型レナードジョーンズポテンシャルを数式（４Ｂ）に代入して得られた解析解である下記数式（１１Ａ）および（１１Ｂ）と数式（４Ａ）から圧力Ｐrc→∞を計算する方法が好ましいことが見出された。
【００５６】
【数１１】
2016181257.tif 000013
【００５７】
〔ステップ４〕
ファンデルワールス相互作用から生じる圧力ＰvdW0→∞を計算するステップ４では、数式（５）の通り、ステップ２およびステップ３で計算したカットオフ半径内のファンデルワールス相互作用から生じる圧力ＰvdW0→ｒｃとカットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力Ｐrc→∞の和から圧力ＰvdW0→∞を計算する。
【００５８】
また数式（１）に、数式（４Ａ）および数式（４Ｂ）、数式（５）、数式（１０Ａ）および数式（１０Ｂ）より求めたＰvdW0→∞を代入することにより全圧力Ｐtotalを求めることができる。なおファンデルワールス相互作用から生じる圧力ＰvdW0→∞以外の圧力の詳細な表式もビリアル定理から導くことができ、参考文献［Glenn J. Martyna, Mark E. Tuckerman, Douglas J. Tobias, Michael L. Klein, Mol. Phys., 1996, 87, p.1117.］に記述されている。この文献の全圧力の表式はクーロン相互作用から生じる圧力Ｐcoulomb0→∞を、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が０．２５フェムト秒のとき４ナノ秒）以上がより好ましい。  
〔結果の出力〕
本発明の計算方法においては、ステップ１〜ステップ５を所定回数繰り返した後に、ファンデルワールス相互作用から生じる圧力ＰvdW0→∞を出力する。ステップ１〜ステップ５の所定回数の繰り返し時点における瞬間値、あるいは一定回数繰り返して得られた平均値のどちらか一方または両方を出力することができる。
【００６５】
また、本発明においては、必要に応じ、上記のように出力した圧力ＰvdW0→∞の値を用いて諸物理量や原子の座標を計算し、出力することができる。出力する諸物理量は特に限定されず、全圧力Ｐtotal、系の密度などが挙げられる。出力する諸物理量は、所定回数の繰り返し時点における瞬間値、あるいは一定回数繰り返して得られた平均値のどちらか一方または両方とすることができる。原子の座標については、分子動力学シミュレーション最大の特徴である原子の時間経過による運動をとらえることができるという性質から、ある繰り返し回数における瞬間値を出力することが好ましい。
【００６６】
以上の手順を実行すれば、高分子中に低分子が存在する系においてスケールされた非結合性相互作用を用いた場合に、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞を求められず正確な圧力を計算できないという課題を解決することができる。
【００６７】
〔相互作用エネルギーまたは配置のサンプリング〕
特に、本発明の計算方法は、高分子-低分子間の相互作用エネルギーや、高分子中の低分子の配置のサンプリングを行うために好適に用いることができる。様々なλlow-polyにおいて、図１のフローに従ったシミュレーションを並列に実行することにより、様々なλlow-polyにおける低分子‐高分子間非結合性相互作用Vlow-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＋（ｍ±１）Δλの状態間のエネルギー差が大きく、正確に自由エネルギーを求めることができない。一方Δλが小さすぎると、自由エネルギーの計算精度は上がるものの、それに伴い計算コストも上昇する。これらの観点から、増減間隔Δλは０．０１以上０．３以下が好ましく、０．０３以上０．２以下がさらに好ましく、０．０５以上０．１５以下が特に好ましい。
【００７５】
ステップ６を実施する間隔Δｎは１以上の値であれば、どのような値でもよい。しかし、Δｎが小さすぎるとサンプリングできる空間は増加するものの局所的な構造のサンプリング効率が低下する、一方Δｎが大きすぎるとサンプリングできる空間が狭くなる。そこで、Δｎは１以上１００００以下が好ましく、１００以上７０００以下がより好ましく、５００以上５０００以下が特に好ましい。
【００７６】
所定回数に１回ステップ６を実施する回を含めた全体の繰り返し回数ｎとステップ６を実施する間隔Δｎを用いると、ステップ６を実施するのはｎをΔｎで割った余りが０に等しいときである。
【００７７】
本サンプリング方法においては、ステップ１〜ステップ５、および所定回数に１回のステップ１〜ステップ６を繰り返した後に、高分子−低分子間の相互作用エネルギーVlow-polyおよび原子の座標ｒを出力する。ここでは、ステップ１〜ステップ５（ステップ６を実施する回を含む）の所定回数の繰り返し時点における瞬間値、あるいは一定回数ごとに得られた平均値のどちらか一方または両方を出力することができる。ただし原子の座標ｒに関しては平均値に物理的な意味が乏しいため瞬間値を出力することが好ましい。また各λlow-polyでのサンプリングという観点から、各λlow-poly毎に高分子−低分子間の相互作用エネルギーVlow-polyおよび原子の座標ｒを出力することが好ましい。
【００７８】
〔プログラムおよび記録媒体〕
本発明のプログラムは、コンピュータの構成要素に上記の圧力の計算方法または相互作用エネルギーと配置のサンプリング方法の各ステップにおける計算手段として動作させるものである。このプログラムを用いれば、上述の計算方法を実行することができる。
【００７９】
また、本発明の記録媒体は、当該プログラムをコンピュータ読み取り可能に記録したものである。この記録媒体のプログラムをコンピュータにて読み取って実行すれば、上述の計算方法を実行することができる。
【実施例】
【００８０】
以下に実施例を挙げて本発明を具体的に説明するが、これらは例示的なものであって限定的なものではない。
【００８１】
［実施例１］
ポリジメチルシロキサン（PDMS）４分子／水１分子の系について、分子動力学シミュレーションを実施した。分子動力学シミュレーションは東レ株式会社で独自に開発したFORTRANプログラムを用いて行い、分子軌道計算は汎用プログラムGaussian（登録商標）09を用いて行った。また、計算にはクロック周波数2.93GHzのCore i7（登録商標）870を搭載したコンピュータ、およびクロック周波数2.7HzのXeon（登録商標）E5-2680を搭載したコンピュータを用いた。
【００８２】
アクセルリス社(現：バイオビア社)製の分子設計システムMaterials Studio（登録商標）7.0上で下記構造式（１）に示すポリマーを作成した。
【００８３】
【化１】
2016181257.tif 000014
【００８４】
構造式(1)において、n=40（分子量＝2996）とした。次いで、Materials Studio（登録商標）7.0 Amorphous Cell Moduleを用いて、表１に示した組成で、PDMS分子と水分子を三次元周期境界条件下でランダムに配置し、系を作成した。Materials Studio（登録商標）7.0の外部出力機能を使用し、作成した系の原子の座標r 、結合情報、原子のポテンシャルパラメータタイプ、原子の質量ｍ、基本セルベクトルhが記載された分子構造ファイルを作成した。
【００８５】
次に、テキストエディタを用いて分子構造ファイルに記載されたポテンシャルパラメータタイプに対応したポテンシャルパラメータを記述したパラメータファイルを作成した。
【００８６】
さらにテキストエディタを用いて非結合性相互作用のスケール値λ、相互作用計算のカットオフ半径ｒｃ、図1計算フローにおけるループの繰り返しを規定する回数n、出力する物理量とその出力間隔、設定温度T0、設定圧力P0、熱浴の慣性量Q、barostatの慣性量Wを記述した計算条件ファイルを作成した。
【００８７】
作成した分子構造ファイル、パラメータファイル、計算条件ファイルを分子動力学シミュレーションプログラムで読み取らせ、データをメモリに格納した。
【００８８】
【表１】
2016181257.tif 000015
【００８９】
計算条件ファイルにおいて、非結合性相互作用のスケール値λを相互作用の種類ごとに設定した。低分子‐高分子間非結合性相互作用のスケール値λlow-poly＝0.1、高分子の分子内非結合性相互作用と、高分子‐高分子間非結合性相互作用のそれぞれのスケール値λpoly-intra＝λpoly-poly＝１とした。本実施例における低分子は水１分子であるため、低分子の分子内非結合性相互作用のスケール値λlow-intra および低分子‐低分子間非結合性相互作用のスケール値λlow-lowは設定しない。
【００９０】
また、計算条件ファイルにおける他のパラメータとして、カットオフ半径ｒｃ＝１０オングストローム、図1計算フローにおけるループの繰り返しを規定する回数を１６００００００回、出力する物理量を原子の座標r、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞および密度とし、その出力間隔を下記ステップ１〜ステップ５のループを４０００回繰り返すごとに１回とした。
【００９１】
温度制御に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〕を用いて、原子数・温度・圧力一定（ＮＰＴ）アンサンブルを構成することにより、系を設定温度T0、設定圧力P0に制御した状態で下記ステップ１〜ステップ５を実施した。なお、設定温度T0=20℃、設定圧力P0=1atm、熱浴の慣性量Q=120.95(kJ・ps2)/mol、barostatの慣性量W=483.79(kJ・ps2)/molとした。
【００９２】
ＮＰＴアンサンブルにおける分子動力学シミュレーションを実施するためには、下記ステップ５において原子の座標r、原子の速度v、熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルh、基本セルベクトルの速度hvを更新することが必要である。そのため分子構造ファイルおよび計算条件ファイルで入力していない、原子の速度v、熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルの速度hvに初期値を与える必要がある。そこで原子の速度vの初期値としてT0=20℃におけるマクスウェル分布に従うよう計算された速度を、熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルの速度hvの初期値として全ての成分に０を与えメモリに格納した。
【００９３】
〔ステップ１〕
メモリに格納された原子の座標r、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、非結合性相互作用のスケール値λ、ならびに数式（８）で示したDreiding型のｖｄＷポテンシャルを含む、数式（６）で示したポテンシャル関数Vtotalを原子の座標rで微分した式を用いて分子内力・分子外力を計算しメモリに格納した。
【００９４】
また、数式（８）中のδ＝1.0×10-20ｍ２とし、数式（６）第６項のクーロンポテンシャルおよびその力の計算にはＥｗａｌｄ法を用い、Ｅｗａｌｄ法のパラメータとして逆空間のクーロンポテンシャルにα＝0.21オングストローム-1、｜ｎ｜２ｍａｘ＝50を用いて計算した。
【００９５】
ポテンシャルパラメータの設定は以下の通りである。まず、PDMS分子については、分子内ポテンシャルパラメータとして、結合長・結合角の平衡位置R0、θ0にはHF/6-31G(d,p)レベルの分子軌道計算によって決定した値を、力の定数Kbond、Kangle、Ktorsion、Kinversionおよび平衡位置φ0、ω0にはＡＭＢＥＲパラメータ、hybrid/UAパラメータおよび桑島ポテンシャルパラメータを用いた。vdWポテンシャルの深さと斥力が生じる位置を決めるパラメータD０、R０にはＡＭＢＥＲパラメータ、および桑島ポテンシャルパラメータを用い、電荷はB3LYP/6-31G(d,p)レベルの分子軌道計算によって決定した値を用いた。
【００９６】
水分子については、結合長・結合角の平衡位置R0、θ0にはＴＩＰ３Ｐパラメータ[W. L. Jorgensen 、R. W. Impey 、J. Chandrasekhar 、 J. D. Madura、M. L. Klein 、J. Chem. Phys.、1983、79、p.926.〕を用い、結合長・結合角の力の定数Kbond、KangleはＩＲの実測値を再現するように決定した値を用い、vdWポテンシャルの深さと斥力が生じる位置を決めるパラメータD０、R０と電荷にはＴＩＰ３Ｐパラメータを用いた。
【００９７】
〔ステップ２〕
ステップ１で計算した分子間力の一つであるvdW力（数式（８）の微分）と、原子の座標rからビリアル定理（数式（１０Ａ）および（１０Ｂ））を用いて圧力ＰvdW0→ｒｃを計算しメモリに格納した。
【００９８】
〔ステップ３〕
原子の座標r、原子のポテンシャルパラメータタイプ、ポテンシャルパラメータ、非結合性相互作用のスケール値λおよび相互作用計算のカットオフ半径ｒｃと数式（４Ａ）、（１１Ａ）、（１１Ｂ）から圧力ＰvdWrc→∞を計算しメモリに格納した。
【００９９】
〔ステップ４〕
ステップ２および３で計算したカットオフ半径内のファンデルワールス相互作用から生じる圧力ＰvdW0→ｒｃとカットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞の和から圧力ＰvdW0→∞を計算しメモリに格納した。
【０１００】
〔ステップ５〕
メモリに格納された原子の質量ｍ、原子の座標r、原子の速度v、熱浴の慣性量Q、熱浴の位置ζ、熱浴の速度ζv、barostatの慣性量W、基本セルベクトルh、基本セルベクトルの速度hv、分子内力、分子間力を用いて、微小時間Δt後の原子の座標r、原子の速度v、熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルh、基本セルベクトルの速度hvを数値積分によって計算し、更新した。数値積分には多重時間刻み数値積分法（ＲＥＳＰＡ法）を用い、時間刻みΔtは表２に示した値を用いた。
【０１０１】
【表２】
2016181257.tif 000016
【０１０２】
ステップ５で更新した微小時間Δt後の原子の座標r、原子の速度v、熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルh、基本セルベクトルの速度hvを初期値としてステップ１〜ステップ５を再度実行させることをn＝１６００００００回繰り返した。ＲＥＳＰＡ法における最小の時間刻みが０．２５フェムト秒であるため、この繰り返し回数は４ナノ秒に相当する。
【０１０３】
なお本実施例ではNose−Hoover法により温度、Parrinello−Rahman法により圧力を制御しているためスッテプ５において熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルh、基本セルベクトルの速度hvを更新した。
【０１０４】
ステップ１〜ステップ５を４０００回（１ピコ秒に相当）繰り返すごとに、原子の座標r、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞および密度を出力した。また計算終了時（ステップ１〜ステップ５を１６００００００回繰り返し終了後）にカットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞および密度の平均値を計算した。
【０１０５】
なお、本実施例の計算では構造緩和のための予備計算として、上述の計算の前に２０℃の原子数・温度・体積一定（ＮＶＴ）アンサンブルの分子動力学シミュレーションを１０ピコ秒、１ａｔｍ／２０℃のＮＰＴアンサンブルの分子動力学シミュレーションを９０ピコ秒、１ａｔｍの分子動力学シミュレーションにおいて温度を２０℃〜５００℃まで２００ピコ秒周期で上昇下降を繰り返し緩和させるシミュレーテッドアニーリングを１４サイクル（２８００ピコ秒に相当）、１ａｔｍ／２０℃のＮＰＴアンサンブルの分子動力学シミュレーションを３００ピコ秒実施した。
【０１０６】
分子動力学シミュレーション終了後、非結合性相互作用のスケール値λの違いによるPDMS構造の変化を明らかにするため、分子動力学シミュレーションで得られた原子の座標rの時系列データとMaterials Studio（登録商標)7.0 Atom Volumes And Surfaces toolを用いてConnolly Surface法[M. L. Connolly, 1983, Science, 221, p.709.]により自由体積を求め、数式（１２）で示される系の空隙率を計算した。
【０１０７】
【数１２】
2016181257.tif 000017
【０１０８】
[実施例２]
実施例１の計算条件ファイルにおいて低分子‐高分子間非結合性相互作用のスケール値λlow-poly＝0.3としてメモリに格納し、それ以外は実施例１と同条件で分子動力学シミュレーションを実施して、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞の平均値、密度の平均値、空隙率を計算した。
【０１０９】
[実施例３]
実施例１の計算条件ファイルにおいて低分子‐高分子間非結合性相互作用のスケール値λlow-poly＝0.4としてメモリに格納し、それ以外は実施例１と同条件で分子動力学シミュレーションを実施して、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞の平均値、密度の平均値、空隙率を計算した。
【０１１０】
[実施例４]
実施例１の計算条件ファイルにおいて低分子‐高分子間非結合性相互作用のスケール値λlow-poly＝0.5としてメモリに格納し、それ以外は実施例１と同条件で分子動力学シミュレーションを実施して、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞の平均値、密度の平均値、空隙率を計算した。
【０１１１】
[実施例５]
実施例１の計算条件ファイルにおいて低分子‐高分子間非結合性相互作用のスケール値λlow-poly＝0.6としてメモリに格納し、それ以外は実施例１と同条件で分子動力学シミュレーションを実施して、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞の平均値、密度の平均値、空隙率を計算した。
【０１１２】
[実施例６]
実施例１の計算条件ファイルにおいて低分子‐高分子間非結合性相互作用のスケール値λlow-poly＝0.8としてメモリに格納し、それ以外は実施例１と同条件で分子動力学シミュレーションを実施して、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞の平均値、密度の平均値、空隙率を計算した。
【０１１３】
[比較例１]
実施例１の計算条件ファイルにおいて低分子‐高分子間非結合性相互作用のスケール値λlow-poly＝0.1としてメモリに格納し、ステップ１およびステップ２を実施した後、ステップ３を実施せず、ステップ４およびステップ５を実施し、以後ステップ１〜２、ステップ４〜５を１６００００００回繰り返し実施し、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞の平均値、密度の平均値、空隙率を計算した。
【０１１４】
ステップ３を実施しないという以外は実施例１と同条件であり、この操作はカットオフ半径から無限遠までのファンデルワールス相互作用から生じる圧力を無視することに相当する。
【０１１５】
[比較例２]
比較例１の計算条件ファイルにおいて低分子‐高分子間非結合性相互作用のスケール値λlow-poly＝0.3としてメモリに格納し、それ以外は比較例１と同条件で分子動力学シミュレーションを実施し、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞の平均値、密度の平均値、空隙率を計算した。
【０１１６】
[比較例３]
比較例１の計算条件ファイルにおいて低分子‐高分子間非結合性相互作用のスケール値λlow-poly＝0.4としてメモリに格納し、それ以外は比較例１と同条件で分子動力学シミュレーションを実施し、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞の平均値、密度の平均値、空隙率を計算した。
【０１１７】
[比較例４]
比較例１の計算条件ファイルにおいて低分子‐高分子間非結合性相互作用のスケール値λlow-poly＝0.5としてメモリに格納し、それ以外は比較例１と同条件で分子動力学シミュレーションを実施し、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞の平均値、密度の平均値、空隙率を計算した。
【０１１８】
[比較例５]
比較例１の計算条件ファイルにおいて低分子‐高分子間非結合性相互作用のスケール値λlow-poly＝0.6としてメモリに格納し、それ以外は比較例１と同条件で分子動力学シミュレーションを実施し、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞の平均値、密度の平均値、空隙率を計算した。
【０１１９】
[比較例６]
比較例１の計算条件ファイルにおいて低分子‐高分子間非結合性相互作用のスケール値λlow-poly＝0.8としてメモリに格納し、それ以外は比較例１と同条件で分子動力学シミュレーションを実施し、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞の平均値、密度の平均値、空隙率を計算した。
【０１２０】
実施例１〜６、比較例１〜６で得られたカットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞の平均値、密度の平均値、空隙率、および同じ低分子‐高分子間非結合性相互作用のスケール値λlow-polyにおける実施例と比較例の空隙率の比を表３にまとめた。また低分子‐高分子間非結合性相互作用のスケール値λlow-polyに対する密度および空隙率を図３、４に示した。
【０１２１】
PDMSの密度の実験値は0.971g/cm3[I. Vinckier, P. Moldenaers, J. Mewis, Journal of Rheology, 1996, 40, p.613.]である。図３より実施例は実験値と同程度の密度を持つが、比較例においては密度が小さくなることがわかる。さらに図４より比較例は実施例より３割程度大きな空隙を持つことがわかる。よって圧力一定のアンサンブルでの分子動力学シミュレーションでは、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞を考慮しなければ、系の密度を過小に評価し、系の空隙率を過大に評価することになる。
【０１２２】
【表３】
2016181257.tif 000018
【０１２３】
上記実施例および比較例では計算時間およびvdW相互作用エネルギーの計算精度の関係から一般的に好ましく用いられるカットオフ半径ｒｃ＝10オングストロームを用いた。
【０１２４】
しかしカットオフ半径を大きくすれば計算時間が増加するが、カットオフ半径内のファンデルワールス相互作用から生じる圧力ＰvdW0→ｒｃの寄与が増加するため、ＰvdW0→ｒｃがＰvdW0→ｒｃ+ＰvdWrc→∞に漸近する。
【０１２５】
[比較例７]
計算条件ファイルにおいてカットオフ半径ｒｃ＝12オングストロームとする以外は比較例５と同条件で分子動力学シミュレーションを実施し、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞の平均値、密度の平均値、空隙率を計算した。
【０１２６】
[比較例８]
計算条件ファイルにおいてカットオフ半径ｒｃ＝15オングストロームとする以外は比較例５と同条件で分子動力学シミュレーションを実施し、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞の平均値、密度の平均値、空隙率を計算した。
【０１２７】
[比較例９]
計算条件ファイルにおいてカットオフ半径ｒｃ＝17オングストロームとする以外は比較例５と同条件で分子動力学シミュレーションを実施し、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞の平均値、密度の平均値、空隙率を計算した。
【０１２８】
実施例５、比較例５、比較例７〜９において分子動力学シミュレーションプログラムを実行するにあたりその計算実行時間を、今回用いたコンピュータの標準シェルであるcshの”time”コマンドを用いて測定し出力した。
【０１２９】
得られた分子動力学シミュレーションプログラムの計算実行時間および密度を図５に示した。カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞を考慮しない比較例５、比較例７〜９においても、カットオフ半径を17オングストローム程度まで大きくすれば、カットオフを行ったことによる影響が小さくなり実験値の密度に近づくが、計算時間も２倍弱に増加した。
【０１３０】
一方、実施例５ではカットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞を考慮したため、短い計算時間で現実系に近い密度再現する高精度な計算を行うことができた。
【０１３１】
[実施例７]
ポリジメチルシロキサン（PDMS）４分子／エチレングリコール（EG）１分子の系を作成し、低分子の分子内非結合性相互作用のスケール値λlow-intra＝1.0と設定する以外実施例１と同条件で分子動力学シミュレーションを実施した。得られた原子の座標rの時系列データと結合情報からエチレングリコールのO-C-C-O二面角の平均値を計算した。
【０１３２】
エチレングリコールの力場パラメータについては分子内ポテンシャルパラメータの力の定数Kbond、Kangle、Ktorsion、Kinversionおよび平衡位置R0、θ0、φ0、ω0として、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とし、低分子内非結合性相互作用をスケールしないことが好ましいといえる。
【０１３７】
【表４】
2016181257.tif 000019
【０１３８】
[実施例９]
実施例１の計算条件ファイルにおいて低分子‐高分子間非結合性相互作用のスケール値λlow-poly＝0.8としてメモリに格納し、それ以外は実施例１と同条件で分子動力学シミュレーションを実施して、ステップ１〜ステップ５のループを４０００回繰り返すごとにＰＤＭＳ−水間の非結合相互作用Vlow-poly、原子の座標を出力した。同様の手順の計算を独立な７種類の初期配置について実施し、初めに計算を行った初期配置を含め８種類の独立な配置から得られたＰＤＭＳ−水間の非結合相互作用Vlow-polyの平均値を計算した。またcshの”time”コマンドを用い、８種類の独立な配置の計算が全て終了するまでの時間を測定した。
【０１３９】
[実施例１０]
ポリジメチルシロキサン（PDMS）４分子／水１分子の系について、図２のフローに従い分子動力学シミュレーションを実施し、様々なλlow-polyにおけるＰＤＭＳ−水間の非結合相互作用Vlow-polyをサンプリングした。
【０１４０】
実施例１と同様に、表１に示した組成で、ＰＤＭＳ分子と水分子を三次元周期境界条件下でランダムに配置し、系を作成した。Materials Studio（登録商標）7.0の外部出力機能を使用し、作成した系の原子の座標r 、結合情報、原子のポテンシャルパラメータタイプ、原子の質量ｍ、基本セルベクトルhが記載された分子構造ファイルを作成した。
【０１４１】
次に、テキストエディタを用いて分子構造ファイルに記載されたポテンシャルパラメータタイプに対応したポテンシャルパラメータを記述したパラメータファイルを作成した。
【０１４２】
さらにテキストエディタを用いて非結合性相互作用のスケール値λ、λlow-polyの増減間隔Δλ、Δλを加減する間隔Δｎ、相互作用計算のカットオフ半径ｒｃ、図２計算フローにおけるループの繰り返しを規定する回数n、出力する物理量とその出力間隔、設定温度T0、設定圧力P0、熱浴の慣性量Q、barostatの慣性量Wを記述した計算条件ファイルを作成した。
【０１４３】
作成した分子構造ファイル、パラメータファイル、計算条件ファイルを分子動力学シミュレーションプログラムで読み取らせ、データをメモリに格納した。
【０１４４】
計算条件ファイルにおいて、非結合性相互作用のスケール値λを相互作用の種類ごとに設定した。低分子‐高分子間非結合性相互作用のスケール値の初期値をλlow-poly＝１、高分子の分子内非結合性相互作用と、高分子‐高分子間非結合性相互作用のそれぞれのスケール値λpoly-intra＝λpoly-poly＝１とした。実施例１同様に低分子の分子内非結合性相互作用のスケール値λlow-intra および低分子‐低分子間非結合性相互作用のスケール値λlow-lowは設定しない。λlow-polyの増減間隔Δλ＝0.1とし、Δλを加減する間隔Δｎ＝４０００回とした。さらに、カットオフ半径ｒｃ＝１０オングストローム、図２計算フローにおけるループの繰り返しを規定する回数ｎ＝１６００００００回とした。出力する物理量をＰＤＭＳ−水間の非結合相互作用Vlow-poly、原子の座標rとしその出力間隔を図２のステップ１〜ステップ５（Mod(n,Δn)=0のときはステップ６も含む）のループを１０００回繰り返すごとに１回とした。
【０１４５】
実施例１同様、温度制御にNose−Hoover法、圧力制御にParrinello−Rahman法を用いて、原子数・温度・圧力一定（ＮＰＴ）アンサンブルを構成することにより、系を設定温度T0、設定圧力P0に制御した状態で図２ステップ１〜ステップ６を実施した。なお、設定温度T0=20℃、設定圧力P0=1atm、熱浴の慣性量Q=120.95(kJ・ps2)/mol、barostatの慣性量W=483.79(kJ・ps2)/molとした。
【０１４６】
また原子の速度vの初期値としてT0=20℃におけるマクスウェル分布に従うよう計算された速度を、熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルの速度hvの初期値として全ての成分に０を与えメモリに格納した。ステップ１〜ステップ５は実施例１と同様に実施した。ステップ１〜ステップ５を４０００回繰り返す毎に、下記ステップ６を１回実行した。
【０１４７】
〔ステップ６〕
メモリに格納されたλlow-polyに増減間隔Δλ（０．１）を次の規則に従い加減した。
【０１４８】
〔Δλ加減規則〕
１．λlow-poly＝０となるまで、ステップ６に到達する毎にΔλを引く。
２．λlow-poly＝０に到達したら、λlow-poly＝１となるまでステップ６に到達する毎にΔλを加える。
３．λlow-poly＝1に到達したら、λlow-poly＝０となるまでステップ６に到達する毎にΔλを引く。
４．上記１〜３の操作を、繰り返しの規定回数ｎに到達するまで繰り返す。
【０１４９】
このようにしてステップ６で更新されたλlow-polyを用いて再びステップ１〜ステップ５を実施した。繰り返しの規定回数に到達後、λlow-poly＝0.8におけるＰＤＭＳ−水間の非結合相互作用Vlow-polyの平均値を計算した。またcshの”time”コマンドを用い、計算が終了するまでの時間を測定した。
【０１５０】
実施例９、実施例１０で計算したＰＤＭＳ−水間の非結合相互作用Vlow-polyおよび計算実行時間を表５にまとめた。表５より実施例９、実施例１０どちらのＰＤＭＳ−水間の非結合相互作用Vlow-polyもほぼ同じ精度で計算できていることがわかる。しかし実施例１０は実施例９の約１／８の計算時間で終了しており、効率が非常に高いといえる。
【０１５１】
【表５】
2016181257.tif 000020
【図１】
2016181257.tif 000027
【図２】
2016181257.tif 000028
【図３】
2016181257.tif 000029
【図４】
2016181257.tif 000030
【図５】
2016181257.tif 000031