(57)【要約】      （修正有）【課題】分子動力学シミュレーションにおいて、高分子中の低分子の運動性を向上させ、高分子中の低分子の配置および分子間相互作用エネルギーの平均値や分布を実用的な計算時間で、かつ統計力学的に信頼性の高い精度で計算する。【解決手段】（１）原子の座標、結合情報、ポテンシャルパラメータタイプ、電荷、質量、スケール値増減条件、出力条件、スケール値の初期値をメモリに格納し、（２）分子内相互作用エネルギー、分子間相互作用エネルギー、分子内力、分子間力を計算し、（３）スケール値増減や出力の有無を切り替え、（４）ｉ、ｊが高分子中の原子と低分子中の原子の組み合わせである場合に、スケール値λｉｊを１未満の値に減少、（５）または１まで増加させ、（６）原子の座標と分子間相互作用エネルギーを出力し、（７）原子の座標と速度を更新するステップとを少なくとも有し、ステップ（２）から（７）を規定回数に到達するまで繰り返し行う。【選択図】図４
2011123874.tif 000002
【特許請求の範囲】
【請求項１】
電子計算機を利用した分子動力学シミュレーションにおける高分子中の低分子の配置および分子間相互作用エネルギーの計算方法であって、（１）原子の座標、結合情報、ポテンシャルパラメータタイプ、電荷、質量、スケール値増減条件、出力条件、スケール値の初期値を電子計算機に入力し、メモリに格納するステップ（ステップ１）と、（２）メモリに格納された原子の座標、結合情報、ポテンシャルパラメータタイプ、電荷、スケール値、および数式（１）で示したポテンシャル関数を用いて、分子内相互作用エネルギー、分子間相互作用エネルギー、分子内力、分子間力を計算し、メモリに格納するステップ（ステップ２）と、
【数１】
2011123874.tif 000009
（ただし、λｉｊ、δはそれぞれスケール値、エネルギーと力の距離依存性を調節するパラメータを表わし、添字ｉ、ｊは原子の番号を表わす。また、δ≧０ｍ２とし、λｉｊの初期値を1とする。）
（３）メモリに格納されたスケール値増減条件、出力条件、および繰り返し回数に基づきスケール値増減や出力の有無を切り替えるステップ（ステップ３）と、（４）ｉ、ｊが高分子中の原子と低分子中の原子の組み合わせである場合にスケール値λｉｊを１未満の値に減少させ、メモリに格納するステップ（ステップ４）と、（５）ｉ、ｊが高分子中の原子と低分子中の原子の組み合わせである場合にスケール値λｉｊを１まで増加させ、メモリに格納するステップ（ステップ５）と、（６）メモリに格納された原子の座標と分子間相互作用エネルギーを出力するステップ（ステップ６）と、（７）メモリに格納された原子の質量、分子内力、分子間力を用いて、原子の座標と速度を更新し、メモリに格納するステップ（ステップ７）とを少なくとも有し、ステップ２から７を規定回数に到達するまで繰り返し行うことを特徴とする計算方法。
【請求項２】
ステップ４およびステップ５におけるスケール値の増減を１ピコ秒あたり１．０以下の速度で変化させる請求項１に記載の計算方法。
【請求項３】
数式（１）において５．０×１０−２１ｍ２≦δ≦５．０×１０−１９ｍ２とする請求項１または２に記載の計算方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、電子計算機を利用した分子動力学シミュレーションにおける、高分子中の低分子の配置および分子間相互作用エネルギーの計算方法に関する。
【背景技術】
【０００２】
分子動力学シミュレーションは、分子や分子集合体の高次構造や物性を詳細に研究するための手法として広く用いられている。近年の電子計算機の急速な進歩に伴って、その適用分野はますます広がっており、合成高分子や生体高分子の極微視的な構造や性質を解析できる手法として、研究に必須の方法論となりつつある。
【０００３】
分子動力学シミュレーションの重要な計算対象の一つとして、イオン交換膜や分離膜といった高分子機能膜が挙げられる。分子動力学シミュレーションを用いて、透過物質と高分子膜との親和性や、透過物質の高分子膜中での運動性、移動経路、存在位置について詳細な知見を得ることで、高分子膜のさらなる高性能化・高機能化が期待できる。
【０００４】
しかし、高分子は分子量が非常に大きいため、分子運動が遅く緩和時間も極めて長大である。そのため、高分子材料中に含まれる低分子は一般に運動性が低く、シミュレーション開始時の位置からほとんど移動しない。従って、一般に高分子シミュレーション研究で実行される数日程度のシミュレーションでは、低分子の配置や分子間相互作用エネルギーについて極めて限定された情報しか得られず、配置を広げて統計力学的に有意なサンプリングを行うためには、数ヶ月以上という非常に長時間のシミュレーションを実行する必要がある。
【０００５】
このように、従来の分子動力学シミュレーションには、実用的なシミュレーション時間で、統計力学的に信頼性の高い精度で高分子中の低分子の配置および分子間相互作用エネルギーの平均値や分布を得ることが困難であるという問題点があった。
【先行技術文献】
【非特許文献】
【０００６】
【非特許文献１】増渕雄一, 高分子系のナノシミュレーション，日本化学会情報化学部会誌，24, 152 (2006).
【発明の概要】
【発明が解決しようとする課題】
【０００７】
高分子中の低分子の運動性は非常に低いため、高分子中の低分子の配置および分子間相互作用エネルギーの平均値や分布を、実用的な計算時間で、かつ統計力学的に信頼性の高い精度で計算する方法は未だ提唱されていない。そこで本発明では、上記問題点を解決した高分子中の低分子の配置および分子間相互作用エネルギーの計算方法を提供することを目的とする。
【課題を解決するための手段】
【０００８】
上記課題を解決するため、本発明は、電子計算機を利用した分子動力学シミュレーションにおける高分子中の低分子の配置および分子間相互作用エネルギーの計算方法であって、（１）原子の座標、結合情報、ポテンシャルパラメータタイプ、電荷、質量、スケール値増減条件、出力条件、スケール値の初期値を電子計算機に入力し、メモリに格納するステップ（ステップ１）と、（２）メモリに格納された原子の座標、結合情報、ポテンシャルパラメータタイプ、電荷、スケール値、および数式（１）で示したポテンシャル関数を用いて、分子内相互作用エネルギー、分子間相互作用エネルギー、分子内力、分子間力を計算し、メモリに格納するステップ（ステップ２）と、
【０００９】
【数１】
2011123874.tif 000003
【００１０】
（ただし、λｉｊ、δはそれぞれスケール値、エネルギーと力の距離依存性を調節するパラメータを表わし、添字ｉ、ｊは原子の番号を表わす。また、δ≧０ｍ２とし、λｉｊの初期値を1とする。）
（３）メモリに格納されたスケール値増減条件、出力条件、および繰り返し回数に基づきスケール値増減や出力の有無を切り替えるステップ（ステップ３）と、（４）ｉ、ｊが高分子中の原子と低分子中の原子の組み合わせである場合にスケール値λｉｊを１未満の値に減少させ、メモリに格納するステップ（ステップ４）と、（５）ｉ、ｊが高分子中の原子と低分子中の原子の組み合わせである場合にスケール値λｉｊを１まで増加させ、メモリに格納するステップ（ステップ５）と、（６）メモリに格納された原子の座標と分子間相互作用エネルギーを出力するステップ（ステップ６）と、（７）メモリに格納された原子の質量、分子内力、分子間力を用いて、原子の座標と速度を更新し、メモリに格納するステップ（ステップ７）とを少なくとも有し、ステップ２から７を規定回数に到達するまで繰り返し行うことを特徴とする計算方法であることを本旨とし、また、種々の改良された態様を提案する。
【発明の効果】
【００１１】
本発明の方法を用いることによって、高分子中で低分子の運動する範囲を飛躍的に拡大させることができ、高分子中の低分子の配置および分子間相互作用エネルギーの平均値や分布を実用的な計算時間で、かつ統計力学的に信頼性の高い精度で計算することが可能になる。
【図面の簡単な説明】
【００１２】
【図１】本発明の方法のフローを説明するための図である。
【図２】スケール値の増減方法を説明するための図である。
【図３】実施例１〜６、比較例１の水分子のｘｙｚ座標の経時変化を示す図である。
【図４】比較例１（ａ）および実施例４（ｂ）で得られた水分子の軌跡を示す図である。
【図５】実施例７〜１６、比較例２の二酸化炭素分子の自己拡散係数を示す図である。
【発明を実施するための形態】
【００１３】
以下、本発明の好ましい実施形態を説明する。
【００１４】
図１に本発明のフローを示す。
【００１５】
初期座標、計算条件を入力するステップ１においては、電子計算機に読み取り可能な形態で原子の座標、結合情報、ポテンシャルパラメータタイプ、電荷、質量、スケール値増減条件、出力条件、スケール値の初期値を入力する。入力には、一般に入手可能な分子モデル作成ソフトウェアやテキストエディタなどが利用できる。次に、当該情報が記述されたファイルを分子動力学計算プログラムを用いて読み取らせ、メモリに格納する。
【００１６】
分子内相互作用エネルギー、分子間相互作用エネルギー、分子内力、分子間力を計算するステップ２においては、メモリに格納された原子の座標、結合情報、ポテンシャルパラメータタイプ、電荷、スケール値、および上記数式（１）で示したポテンシャル関数を用いて、エネルギーと力を計算し、メモリに格納する。
【００１７】
ここで、数式（１）の右辺第１項は結合長、第２項は結合角、第３項は二面角、第４項は反転、第５項はファンデルワールス（ｖｄＷ）力、第６項はクーロン力のそれぞれポテンシャル関数を表わす。なお、数式（１）の第１項から第４項を総称して分子内ポテンシャル関数と呼び、第５項、第６項をまとめて分子間ポテンシャル関数と呼ぶ。また、対応するエネルギーをそれぞれ分子内相互作用エネルギー、分子間相互作用エネルギーと呼び、分子内ポテンシャル関数、分子間ポテンシャル関数を原子の位置で微分した式を用いて計算した力をそれぞれ分子内力、分子間力と呼ぶ。
【００１８】
数式（１）中のパラメータには、ドライディング（ＤＲＥＩＤＩＮＧ）  〔 エス・エル・メイオー （Ｓ. Ｌ. Ｍａｙｏ） 、ビー・ディー・オラフソン （Ｂ. Ｄ. Ｏｌａｆｓｏｎ） 、ダブリュー・エー・ゴダード３世 （Ｗ. Ａ. Ｇｏｄｄａｒｄ III） 、ジャーナル  オブ  フィジカル  ケミストリー （Ｊ. Ｐｈｙｓ. Ｃｈｅｍ.） 、１９９０、９４、ｐ.８８９７.〕 、アンバー （ＡＭＢＥＲ） 〔ダブリュー・ディー・コーネル （Ｗ. Ｄ. Ｃｏｒｎｅｌｌ） 、ピー・シープラック （Ｐ. Ｃｉｅｐｌａｋ） 、シー・アイ・ベイリー （Ｃ. Ｉ. Ｂａｙｌｙ） 、アイ・アール・グルード （Ｉ. Ｒ. Ｇｏｕｌｄ） 、ケー・エム・メルツ・ジュニア （Ｋ. Ｍ. Ｍｅｒｚ ｊｒ） 、ディー・エム・ファーガソン （Ｄ. Ｍ. Ｆｅｒｇｕｓｏｎ） 、ディー・シー・スペルマイアー （Ｄ. Ｃ. Ｓｐｅｌｌｍｅｙｅｒ） 、ティー・フォックス （Ｔ. Ｆｏｘ） 、ジェー・ダブリュー・コールドウェル （Ｊ. Ｗ. Ｃａｌｄｗｅｌｌ）、ピー・エー・コールマン （Ｐ. Ａ.  Ｋｏｌｌｍａｎ） 、ジャーナル  オブ  アメリカン  ケミカル  ソサイエティー （Ｊ. Ａｍ. Ｃｈｅｍ. Ｓｏｃ.）、１９９５、１１７、ｐ.５１７９〕、桑島等のパラメータ 〔桑島聖、野間祐之、逢坂俊郎、“非晶高分子のＭＤ計算用力場パラメータについてｉｉ：ポリオレフィン”、第４回計算化学シンポジウム&ＣＣＳ展  講演要旨集、日本化学プログラム交換機構（１９９４）〕等の公知のパラメータを用いることができる。
【００１９】
計算精度について考慮した場合、以下に述べる方法によってポテンシャルパラメータを決定することがさらに好ましい。
【００２０】
数式（１）の第１項から第４項で表わされる分子内ポテンシャル関数のパラメータは量子化学計算によって決定することが好ましいが、量子化学計算としては、ＨＦ法〔Ａ．ザボ，Ｎ．Ｓ．オストランド，「新しい量子化学」，東京大学出版会（１９８７）．〕を用いることが好ましく、Ｂ３ＬＹＰ法〔エー・ディー・ベッケ（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.〕が好適に用いられる。
【００２２】
数式（１）の第５項で表わされるｖｄＷポテンシャルのパラメータについては、密度と蒸発熱の実験値を再現するように決定することが好ましい。
【００２３】
数式（１）の第６項のクーロンポテンシャルのパラメータ、すなわち電荷については、ＱＥｑ法〔エー・ケー・ラッペ（Ａ．Ｋ．Ｒａｐｐｅ），ダブリュー・エー・ゴダード３世（Ｗ．Ａ．Ｇｏｄｄａｒｄ  ＩＩＩ），ジャーナル  オブ  フィジカル  ケミストリー（Ｊ．Ｐｈｙｓ．Ｃｈｅｍ．），９５，３３５８−３３６３（１９９１）．〕等の学術的に信頼性の高い手法を用いて決定することが好ましく、量子化学計算から静電ポテンシャル（ＥＳＰ）を求め、ＥＳＰフィッティングして電荷を決定することがさらに好ましい。
【００２４】
静電ポテンシャルを求めるための量子化学計算としては、ＨＦ法やＢ３ＬＹＰ法を好適に用いることができ、基底関数としては６−３１Ｇ（ｄ）基底関数、６−３１Ｇ（ｄ，ｐ）基底関数、６−３１＋Ｇ（ｄ，ｐ）基底関数が好適に用いられ、さらに好ましくは、ａｕｇ−ｃｃ−ｐＶＤＺ基底関数、ａｕｇ−ｃｃ−ｐＶＴＺ基底関数、ａｕｇ−ｃｃ−ｐＶＱＺ基底関数〔ティー・エッチ・ダニング・ジュニア（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.〕といった方法を好適に用いることができる。
【００２６】
以上のような量子化学計算を実行するために、有償または無償で入手可能な種々の量子化学計算プログラムパッケージを利用することができる。例えば、ガウシアン（登録商標）９８〔ガウシアン98（A11.3改訂版）（Gaussian 98 (Revision A.11.3)）, エム・ジェー・フリッシュ（M. J. Frisch）他, ガウシアン社（Gaussian, Inc.）, ペンシルベニア州ピッツバーグ（Pittsburgh PA）, 2002.〕、ガウシアン（登録商標）０３〔ガウシアン０３（Ｅ．０１改訂版）（Ｇａｕｓｓｉａｎ  ０３（Ｒｅｖｉｓｉｏｎ  Ｅ．０１）），エム・ジェー・フリッシュ（Ｍ．Ｊ．Ｆｒｉｓｃｈ）他，ガウシアン社（Ｇａｕｓｓｉａｎ，Ｉｎｃ．），コネチカット州ウォーリングフォード（Ｗａｌｌｉｎｇｆｏｒｄ  ＣＴ），２００４．〕、ゲイムス（ＧＡＭＥＳＳ）〔エム・ダブリュー・シュミット（M. W. Schmidt）他, ジャーナル  オブ  コンピューテーショナル  ケミストリー（J. Comput. Chem.）, 1993, 14, p.1347.〕といったプログラムを好適に用いることができる。
【００２７】
数式（１）の第５項、第６項中のλｉｊ、δはそれぞれスケール値、エネルギーと力の距離依存性を調節するパラメータである。従来の一般的な分子動力学計算においては、分子間ポテンシャル関数として、数式（２）に示した関数を用いる。
【００２８】
【数２】
2011123874.tif 000004
【００２９】
すなわち、数式（２）は数式（１）の第５項、第６項において、λｉｊ＝１、δ＝０ｍ２とした式に相当する。数式（２）を用いて、高分子中の低分子の座標と速度を更新した場合、実用的な計算時間では低分子の座標はほとんど変わらない。これは、高分子中の低分子の運動性が低いためである。
【００３０】
本発明者らは、鋭意検討の結果、高分子中の低分子の運動性を向上させるために、スケール値λｉｊを増減するステップ４、５を加え、数式（１）で示したポテンシャル関数を用いて分子内力とスケールされた分子間力を計算し、これらの力を用いて高分子中の低分子の座標と速度を更新することで、低分子の運動性を向上させるとともに分子間相互作用エネルギーの計算精度を確保することに成功し、本発明に到達したものである。
【００３１】
本発明の方法では、スケール値λｉｊを適切に制御することによって、実用的な計算時間で高分子中の低分子に十分な運動性を与えることができる。ここで、λの添字ｉ、ｊは原子の番号を表わし、全てのｉ、ｊについてλｉｊの初期値を１とする。低分子の運動性を向上させるためには、ｉ、ｊが高分子中の原子と低分子中の原子の組み合わせである場合には、スケール値λｉｊは分子動力学計算中に少なくとも１回は１未満の値をとる必要がある。これに対して、ｉ、ｊが高分子中の原子と高分子中の原子の組み合わせである場合にはλｉｊ＝１に固定する。λｉｊ（高分子−低分子の組み合わせ）の最小値としては、高分子中の低分子の種類が水分子のように１個以下の重元素（水素以外の元素）を含む場合には、１．０×１０−３以下とすることが好ましく、高分子中の低分子の種類が二酸化炭素のように２個以上の重元素を含む場合には、１．０×１０−５以下とすることが好ましい。
【００３２】
また、統計力学的に正しい高分子中の低分子の配置および分子間相互作用エネルギーを得るためには、一定時間スケール値λｉｊ（高分子−低分子の組み合わせ）を１に保った後、各原子の座標および高分子中の低分子の分子間相互作用エネルギーを出力することが好ましい。
【００３３】
メモリに格納されたスケール値増減条件、出力条件、および繰り返し回数に基づいて、スケール値の増減や出力の有無を切り替えるステップ３においては、例えば、図２に示す制御手順が好適に用いられる。図中のＡ〜Ｅはスケール値増減条件、出力条件に相当するが、これは分子動力学計算の繰り返し回数に対応する正の整数であり、１≦Ａ≦Ｂ＜Ｃ＜Ｄ＜Ｅ≦分子動力学計算の全繰り返し回数を満たす。図２の制御手順を以下に詳述する。
【００３４】
手順１．｛（繰り返し回数）−１｝をＥで割った余りに１を加えた変数Ｍを計算し、メモリに格納する。
【００３５】
手順２．Ｍ≦Ａならば、ステップ４を実行し、スケール値λｉｊ（高分子−低分子の組み合わせ）を１未満の値に減少させ、メモリに格納する。
【００３６】
手順３．Ｍ＞Ｂ、かつ、Ｍ≦Ｃならば、ステップ５を実行し、スケール値λｉｊ（高分子−低分子の組み合わせ）を１未満の値から１まで増加させ、メモリに格納する。
【００３７】
手順４．Ｍ＞Ｄならば、ステップ６を実行し、各原子の座標と分子間相互作用エネルギーを出力する。
【００３８】
ここで、手順２，３におけるスケール値の増減の速度は、数値積分の安定性を考慮すると、１ピコ秒あたり１．０以下の速度とすることが好ましい。これ以上の速度でスケール値を増減させると、原子同士が衝突するような配置が発生する可能性があるため好ましくない。また、Ａ〜Ｅについては、数値積分の安定性、出力した座標とエネルギーの力学的な正当性、およびシミュレーションに要する時間を考慮すると、高分子中の低分子の種類が水分子のように１個以下の重元素を含む場合には、例えば、以下の値が好適に用いられる。カッコ内の値は分子動力学計算で用いる時間刻みを０．１２５フェムト秒とした場合の時間である。
【００３９】
Ａ：  １６０００〜  ４００００（  ２〜  ５ピコ秒）
Ｂ：  ６４０００〜１６００００（  ８〜２０ピコ秒）
Ｃ：  ８００００〜２０００００（１０〜２５ピコ秒）
Ｄ：１４４０００〜３６００００（１８〜４０ピコ秒）
Ｅ：１６００００〜４０００００（２０〜５０ピコ秒）
一方、高分子中の低分子の種類が二酸化炭素のように２個以上の重元素を含む場合には、以下の値が好適に用いられる。
【００４０】
Ａ：  ４００００〜  １６００００（  ５〜  ２０ピコ秒）
Ｂ：１６００００〜  ６４００００（２０〜  ８０ピコ秒）
Ｃ：２０００００〜  ８０００００（２５〜１００ピコ秒）
Ｄ：３６００００〜１４４００００（４５〜１８０ピコ秒）
Ｅ：４０００００〜１６０００００（５０〜２００ピコ秒）
これらの値が小さすぎるか、または大きすぎると、低分子が高分子中を十分広い範囲に亘って運動することができないため好ましくない。また、これらの値が小さすぎる場合には出力した座標とエネルギーの力学的な正当性が保証されず、大きすぎる場合にはシミュレーションに要する時間が実用的な範囲を超えるため好ましくない。
【００４１】
数式（１）中のδについては、δ＝０ｍ２とした場合、数式（１）の第５項、第６項は従来の一般的な分子動力学計算で用いる数式（２）にスケール値を乗じた式となり、比較的容易に計算プログラムを作成することができる。しかしながら、δ＝０ｍ２とし、かつスケール値が非常に小さい値をとる場合に、数式（１）の第５項で表わされるファンデルワールス（ｖｄＷ）力が急激に変化し、数値積分の安定性が低下する。この傾向は原子間の距離が近距離となった場合に特に顕著になる。この問題を解決するために、スケール値の最小値を非常に小さな値に設定する場合には、δ＞０ｍ２とすることが好ましい。さらに好ましくは、δ≧５．０×１０−２１ｍ２である。一方、δの値が大きすぎると、スケール値が変化した際、ｖｄＷ力の変化が大きくなりすぎるため、δ≦５．０×１０−１９ｍ２とすることが好ましい。このδの値については、発明者らが種々検討した結果、δ＝１．０×１０−２０が好適に用いられることを見出した。
【００４２】
原子の座標と速度を更新するステップ７においては、メモリに格納された原子の質量、分子内力、分子間力を用いて、微小時間後の原子の座標と速度を数値積分によって計算する。
【００４３】
数値積分には、差分近似法や予測子−修正子法〔岡崎進，「コンピューターシミュレーションの基礎」，化学同人（２０００）．〕などの一般的な解法を好適に用いることができる。異なる時定数の力ごとに大きさの異なる時間刻みを設定することのできる多重時間刻み数値積分（ＲＥＳＰＡ法）法〔ジー・ジェー・マルチナ （Ｇ. Ｊ. Ｍａｒｔｙｎａ）、モルキュラー  フィジックス （Ｍｏｌ. Ｐｈｙｓ.）、１９９６、８７、ｐ.１１１７.〕を用いると、シミュレーション時間を短縮できるためさらに好ましい。
【００４４】
以後、ステップ２〜７を規定の全繰り返し回数に到達するまで繰り返し実行し、本発明の手順は終了する。
【実施例】
【００４５】
本発明の具体的実施態様を以下に実施例をもって述べるが、本発明はこれに限定されるものではない。
【００４６】
実施例１
以下のステップ１〜７の手順により、ポリエチレンテレフタレート（ＰＥＴ）４分子／水１分子の系について、分子動力学計算を実施した。なお、分子動力学計算は東レ株式会社で独自に開発したＦＯＲＴＲＡＮプログラムを用いて行い、分子軌道計算は汎用プログラムガウシアン（登録商標）９８を用いて行った。また、計算にはクロック周波数３．０ＧＨｚのＰｅｎｔｉｕｍ（登録商標）４を搭載した電子計算機、およびクロック周波数３．０ＧＨｚのＸｅｏｎ（登録商標）５４７２を搭載した電子計算機を用いた。
【００４７】
ステップ１：初期座標、計算条件の入力
アクセルリス社製の分子設計システムＣｅｒｉｕｓ２（登録商標）上で３Ｄ−Ｓｋｅｔｃｈｅｒを用いて、構造式（１）に示すポリマーを作成した。
【００４８】
【化１】
2011123874.tif 000005
【００４９】
構造式（１）において、ｎ＝１６（分子量＝３０７７）とした。次いで、Ａｍｏｒｐｈｏｕｓ  Ｂｕｉｌｄｅｒを用いて、表１に示した組成で、ＰＥＴ分子と水分子を３次元周期境界条件下でランダムに配置したモデル系を構成し、原子の座標、結合情報、ポテンシャルパラメータタイプ、電荷、質量に関するデータが記載された分子構造ファイルを作成した。次に、テキストエディタを用いて、スケール値増減条件、出力条件を記載した電子ファイルを作成した。作成した分子構造ファイル、スケール値増減条件、出力条件、スケール値の初期値が記載されたファイルを分子動力学計算プログラムを用いて読み取らせ、データをメモリに格納した。ここで、スケール値の初期値は１とした。
【００５０】
【表１】
2011123874.tif 000006
【００５１】
ステップ２：エネルギーと力を計算
メモリに格納された原子の座標、結合情報、ポテンシャルパラメータタイプ、電荷、スケール値、および数式（１）で示したポテンシャル関数を用いて分子内相互作用エネルギー、分子間相互作用エネルギーを計算し、数式（１）を原子の位置で微分した式を用いて分子内力、分子間力を計算した。さらに、これらの分子内力、分子間力と原子の位置から分子内圧力、分子間圧力を計算した。ここで、数式（１）中のδ＝０ｍ２とした。計算した分子内相互作用エネルギー、分子間相互作用エネルギー、分子内力、分子間力、分子内圧力、分子間圧力をメモリに格納した。
【００５２】
ｖｄＷポテンシャルおよび実空間のクーロンポテンシャルの計算はカットオフ半径ｒｃ＝１０オングストロームとし、逆空間のクーロンポテンシャルはα＝０．２１オングストローム−１、｜ｎ｜２ｍａｘ＝５０としてＥｗａｌｄ法を用いて計算した。
【００５３】
ポテンシャルパラメータの設定は以下の通りである。まず、ＰＥＴ分子については、分子内ポテンシャルパラメータとして、ＨＦ／６−３１Ｇ（ｄ，ｐ）レベルの分子軌道計算によって決定したパラメータ、ＡＭＢＥＲパラメータ、および桑島ポテンシャルパラメータを用い、ｖｄＷパラメータにはＡＭＢＥＲパラメータ、および桑島ポテンシャルパラメータを用い、電荷はＢ３ＬＹＰ／６−３１Ｇ（ｄ，ｐ）レベルの分子軌道計算によって決定した値を用いた。水分子については、結合長・結合角の平衡位置にはＴＩＰ３Ｐパラメータ〔ダブリュー・エル・ヨルゲンセン （Ｗ. Ｌ. Ｊｏｒｇｅｎｓｅｎ） 、アール・ダブリュー・インペイ （Ｒ. Ｗ. Ｉｍｐｅｙ） 、ジェー・チャンドラセカール （Ｊ. Ｃｈａｎｄｒａｓｅｋｈａｒ） 、ジェー・ディー・マドゥラ （Ｊ. Ｄ. Ｍａｄｕｒａ） 、エム・エル・クライン （Ｍ. Ｌ. Ｋｌｅｉｎ） 、ジャーナル  オブ  ケミカル  フィジックス （Ｊ. Ｃｈｅｍ. Ｐｈｙｓ.）、１９８３、７９、ｐ.９２６.〕を用い、結合長・結合角のバネ定数はＩＲの実測値を再現するように決定した値を用い、分子間ポテンシャルパラメータにはＴＩＰ３Ｐパラメータを用いた。
【００５４】
ステップ３：スケール値増減や出力の有無を決定
図２に示した手順に従い、表２の条件１を用いて、スケール値の増減や、座標と分子間相互作用エネルギーの出力の有無を決定した。
【００５５】
【表２】
2011123874.tif 000007
【００５６】
ステップ４：スケール値を減少
スケール値を１から０．００１まで等速度で減少させ、メモリに格納した。
【００５７】
ステップ５：スケール値を増加
スケール値を０．００１から１まで等速度で増加させ、メモリに格納した。
【００５８】
ステップ６：データ出力
各原子の座標と分子間相互作用エネルギーを出力した。
【００５９】
ステップ７：原子の座標と速度を更新
メモリに格納された原子の質量、分子内力、分子間力、分子内圧力、分子間圧力を用いて、微小時間後の原子の座標と速度を数値積分によって計算し、メモリに格納した。数値積分には多重時間刻み数値積分法を用い、時間刻みは表３に示した値を用いた。
【００６０】
【表３】
2011123874.tif 000008
【００６１】
また、温度をＮｏｓｅ−Ｈｏｏｖｅｒ法〔エム・タッカーマン（Ｍ．Ｔｕｃｋｅｒｍａｎ），ビー・ジェー・バーン（Ｂ．Ｊ．Ｂｅｒｎｅ），ジー・ジェー・マルチナ（Ｇ．Ｊ．Ｍａｒｔｙｎａ），ジャーナル  オブ  ケミカル  フィジックス（Ｊ．Ｃｈｅｍ．Ｐｈｙｓ．），９７，１９９０（１９９２）．〕を用いて２５℃に制御し、圧力を斜交系セルを用いるＰａｒｒｉｎｅｌｌｏ−Ｒａｈｍａｎ法〔エム・パリネロ （Ｍ. Ｐａｒｒｉｎｅｌｌｏ）、エー・ラーマン （Ａ. Ｒａｈｍａｎ）、ジャーナル  オブ  アプライド  フィジックス  （Ｊ. Ａｐｐｌ. Ｐｈｙｓ.）、１９８１、５２、ｐ.７１８２〕を用いて１ａｔｍに制御してＮＰＴアンサンブルを構成した。
【００６２】
以後、全繰り返し回数を３２００００００回としてステップ２〜７を実行した。この全繰り返し回数は４０００ピコ秒に相当するが、分子動力学計算中には、本発明の効果を調べるために、１ピコ秒に１回の頻度で水分子中の酸素原子のｘｙｚ座標を出力した。
【００６３】
なお、本計算では構造緩和のための予備計算として、上述の計算の前に２５℃のＮＶＴアンサンブルの分子動力学計算を１０ピコ秒、１ａｔｍ／２５℃のＮＰＴアンサンブルの分子動力学計算を９０ピコ秒、１ａｔｍ／２００℃のＮＰＴアンサンブルの分子動力学計算を４００ピコ秒、１ａｔｍ／２５℃のＮＰＴアンサンブルの分子動力学計算を５００ピコ秒実施した。
【００６４】
実施例２
表２の条件２を用いてスケール値の増減の有無を決定し、それ以外は実施例１と同じ条件で分子動力学計算を実施し、水分子中の酸素原子のｘｙｚ座標を出力した。
【００６５】
実施例３
表２の条件３を用いてスケール値の増減の有無を決定し、それ以外は実施例１と同じ条件で分子動力学計算を実施し、水分子中の酸素原子のｘｙｚ座標を出力した。
【００６６】
実施例４
表２の条件４を用いてスケール値の増減の有無を決定し、それ以外は実施例１と同じ条件で分子動力学計算を実施し、水分子中の酸素原子のｘｙｚ座標を出力した。
【００６７】
実施例５
表２の条件５を用いてスケール値の増減の有無を決定し、それ以外は実施例１と同じ条件で分子動力学計算を実施し、水分子中の酸素原子のｘｙｚ座標を出力した。
【００６８】
実施例６
表２の条件６を用いてスケール値の増減の有無を決定し、それ以外は実施例１と同じ条件で分子動力学計算を実施し、水分子中の酸素原子のｘｙｚ座標を出力した。
比較例１
１未満のスケール値を導入したことによる効果を確認するために、実施例１で記載したステップ１〜７の手順からスケール値を制御する手順（ステップ３〜５）を除いて、それ以外は実施例１と同じ条件で分子動力学計算を実施し、水分子中の酸素原子のｘｙｚ座標を出力した。
【００６９】
実施例１〜６、比較例１で出力した水分子中の酸素原子のｘｙｚ座標を図３に示した。この図をみると、実施例１〜６は比較例１よりも水分子の運動性が飛躍的に向上していることがわかる。特に実施例４、５は水分子の運動性が高く、実施例４は最も運動性が高かった。
【００７０】
実施例４で得られた水分子の軌跡を図４（ｂ）に示した。図４（ａ）に示した比較例１の結果と比較すると、実施例４における水分子はＰＥＴ分子鎖内を広く拡散していることがわかる。
【００７１】
実施例７
ポリエチレンテレフタレート（ＰＥＴ）４分子／二酸化炭素１分子の系について、表２の条件２を用いてスケール値の増減の有無を決定し、それ以外は実施例１と同じ条件で分子動力学計算を実施し、二酸化炭素中の炭素原子のｘｙｚ座標を出力した。出力したｘｙｚ座標を用いて、平均２乗変位を計算し、平均２乗変位の２ｐｓ〜１０ｐｓの傾きを用いて見かけの自己拡散係数を計算し、これを出力した。
【００７２】
二酸化炭素分子の力場パラメータについては、結合長・結合角の平衡位置にはＪｅｄｌｏｖｓｚｋｙらの文献〔Ｐ．Ｊｅｄｌｏｖｓｚｋｙ、Ｍ．Ｍｅｚｅｉ、ジャーナル  オブ  アメリカン  ケミカル  ソサイエティ（Ｊ．Ａｍ．Ｃｈｅｍ．Ｓｏｃ．），１２２，５１２５（２０００）．〕記載のパラメータを用い、結合長・結合角の力の定数にはドライディング（ＤＲＥＩＤＩＮＧ）パラメータを用い、ｖｄＷパラメータおよび電荷にはＪｅｄｌｏｖｓｚｋｙらの文献パラメータを用いた。
【００７３】
実施例８
スケール値を増減させるステップ４、５において、スケール値の最小値を１．０×１０−４とし、それ以外は実施例７と同じ条件で二酸化炭素中の炭素原子の自己拡散係数を計算し、これを出力した。
【００７４】
実施例９
スケール値を増減させるステップ４、５において、スケール値の最小値を１．０×１０−５とし、それ以外は実施例７と同じ条件で二酸化炭素中の炭素原子の自己拡散係数を計算し、これを出力した。
【００７５】
実施例１０
スケール値を増減させるステップ４、５において、スケール値の最小値を１．０×１０−６とし、それ以外は実施例７と同じ条件で二酸化炭素中の炭素原子の自己拡散係数を計算し、これを出力した。
【００７６】
実施例１１
スケール値を増減させるステップ４、５において、スケール値の最小値を０．０とし、それ以外は実施例７と同じ条件で二酸化炭素中の炭素原子の自己拡散係数を計算し、これを出力した。
【００７７】
実施例１２
エネルギーと力を計算するステップ２において、数式（１）中のδ＝１．０×１０−２０ｍ２とし、それ以外は実施例７と同じ条件で二酸化炭素中の炭素原子の自己拡散係数を計算し、これを出力した。
【００７８】
実施例１３
エネルギーと力を計算するステップ２において、数式（１）中のδ＝１．０×１０−２０ｍ２とし、それ以外は実施例８と同じ条件で二酸化炭素中の炭素原子の自己拡散係数を計算し、これを出力した。
【００７９】
実施例１４
エネルギーと力を計算するステップ２において、数式（１）中のδ＝１．０×１０−２０ｍ２とし、それ以外は実施例９と同じ条件で二酸化炭素中の炭素原子の自己拡散係数を計算し、これを出力した。
【００８０】
実施例１５
エネルギーと力を計算するステップ２において、数式（１）中のδ＝１．０×１０−２０ｍ２とし、それ以外は実施例１０と同じ条件で二酸化炭素中の炭素原子の自己拡散係数を計算し、これを出力した。
【００８１】
実施例１６
エネルギーと力を計算するステップ２において、数式（１）中のδ＝１．０×１０−２０ｍ２とし、それ以外は実施例１１と同じ条件で二酸化炭素中の炭素原子の自己拡散係数を計算し、これを出力した。
【００８２】
比較例２
スケール値を制御する手順（ステップ３〜５）を除いて、それ以外は実施例７と同じ条件で二酸化炭素中の炭素原子の自己拡散係数を計算し、これを出力した。
実施例７〜１６、比較例２で出力した二酸化炭素中の炭素原子の自己拡散係数を図５に示した。拡散係数の値が大きい程、二酸化炭素の運動性が高いことを表わすが、図５をみると、実施例７〜１０、１２〜１６は比較例２よりも二酸化炭素の運動性が飛躍的に向上していることがわかる。特に実施例１６は水分子の運動性が高かった。一方、スケール値の最小値を０．０、δ＝０ｍ２とした実施例１１は、二酸化炭素分子がＰＥＴに衝突することによって、分子動力学計算プログラムが異常終了した。これに対して、スケール値の最小値を０．０、δ＝１．０×１０−２０ｍ２とした実施例１６は安定に数値積分が実行された。以上の結果からスケール値の最小値が非常に小さい場合には、δ＞０ｍ２とすることが好ましいことがわかった。
【符号の説明】
【００８３】
１：水分子の軌跡
２：ポリエチレンテレフタレート分子のスナップショット
【図１】
2011123874.tif 000010
【図２】
2011123874.tif 000011
【図３】
2011123874.tif 000012
【図４】
2011123874.tif 000013
【図５】
2011123874.tif 000014