JP2007149075A - 粗視化分子シミュレーション用点電荷決定方法 - Google Patents
粗視化分子シミュレーション用点電荷決定方法 Download PDFInfo
- Publication number
- JP2007149075A JP2007149075A JP2006292293A JP2006292293A JP2007149075A JP 2007149075 A JP2007149075 A JP 2007149075A JP 2006292293 A JP2006292293 A JP 2006292293A JP 2006292293 A JP2006292293 A JP 2006292293A JP 2007149075 A JP2007149075 A JP 2007149075A
- Authority
- JP
- Japan
- Prior art keywords
- charge
- molecular
- coarse
- atom
- electrostatic potential
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
【課題】粗視化モデルを用いて分子シミュレーションを高速かつ高精度に実行するための粗視化分子シミュレーション用点電荷決定法を提供する。
【解決手段】分子のファン・デル・ワールス表面における静電ポテンシャルの空間分布を量子化学計算によって求める静電ポテンシャル決定ステップ２と，分子あるいは分子集合体を構成する原子を２つ以上まとめて粗視化粒子とすることにより分子あるいは分子集合体を構成する単位を減少させる粗視化ステップ３と，粗視化粒子の座標と点電荷から静電ポテンシャルの空間分布を求める近似静電ポテンシャル計算ステップ５と，量子化学による静電ポテンシャルと点電荷による近似静電ポテンシャルとを比較して残差を計算する残差計算ステップ６と，残差の収束を判断し，残差が少なくなるよう点電荷の値を調整する電荷最適化ステップ７とを有する，粗視化分子シミュレーション用点電荷決定方法を用いる。
【選択図】図１
【解決手段】分子のファン・デル・ワールス表面における静電ポテンシャルの空間分布を量子化学計算によって求める静電ポテンシャル決定ステップ２と，分子あるいは分子集合体を構成する原子を２つ以上まとめて粗視化粒子とすることにより分子あるいは分子集合体を構成する単位を減少させる粗視化ステップ３と，粗視化粒子の座標と点電荷から静電ポテンシャルの空間分布を求める近似静電ポテンシャル計算ステップ５と，量子化学による静電ポテンシャルと点電荷による近似静電ポテンシャルとを比較して残差を計算する残差計算ステップ６と，残差の収束を判断し，残差が少なくなるよう点電荷の値を調整する電荷最適化ステップ７とを有する，粗視化分子シミュレーション用点電荷決定方法を用いる。
【選択図】図１
Description
本発明は，電子計算機を利用した分子シミュレーションを高速かつ高精度に実施するための点電荷を算出する方法に関する。
分子力学法や分子動力学法といった分子シミュレーション手法は，分子や分子集合体の高次構造や物性を詳細に研究するための手法として広く用いられており，近年の電子計算機の急速な進歩に伴ってますますその適用分野を広げつつある。とりわけ，ポリオレフィンやポリエステルをはじめとする合成高分子，タンパク質やＤＮＡなどの生体高分子といった巨大分子への適用が盛んに試みられている。
ところが，一般に分子シミュレーションの計算負荷は系を構成する粒子数の２乗に比例して増大することが知られており，近年の高速コンピュータをもってしても巨大分子への分子シミュレーション適用は容易ではない。そのため，計算精度を落とすことなく分子シミュレーションを高速化するための手法が求められていた。
このような課題に対して，分子を構成する原子のうち任意の複数をまとめて１つの粗視化粒子とすることにより分子を構成する粒子数を減少させる，粗視化の考え方が提唱されている。粗視化モデルに対し，全ての原子を顕わに扱う手法を全原子モデルという。
粗視化の手法は，分子鎖をひもとして近似するような大幅な粗視化を進めたモデルから，分子の構造式を想起できる程度の軽度の粗視化まで，目的や系の大きさに応じて様々であるが，後者の例として，融合原子モデル（united atom model）と呼ばれるものがある。
この方法は，水素原子を結合する重原子に融合させて１つの粗視化粒子（融合原子）として扱うものであり，例えばメチル基やメチレン基といった原子団は１つの融合原子として計算することにより，系を構成する粒子の数をおよそ半分程度まで減少させることが可能である。そのため，同程度の計算負荷で２倍の大きさの系についてシミュレーションすることが可能になる。また，同じ大きさの系であれば，粒子数が半分になることにより計算負荷はおよそ４分の１まで軽減されるため，同程度の負荷でシミュレーション時間を４倍長く取ることができ，より長時間の物理現象を研究することが可能となる。さらに，運動の周期の短い水素−重元素結合伸縮を考慮する必要がなく，全原子モデルに比べてシミュレーションの時間刻みを粗く取ることができるため，さらに長時間のシミュレーションを実行することができる。
融合原子モデルは，高分子の分子シミュレーションを中心に広く用いられており，例えばポリオレフィンについてリグビーらによる研究例がある（非特許文献１）。
分子シミュレーションにおいては，原子間にはたらく相互作用を結合（分子内）と非結合（分子間）との２つに分類し，さらに非結合相互作用についてはレナード・ジョーンズ型の分散力ポテンシャルと静電相互作用との２成分で表現する近似が広く用いられている。
後者を計算するには，各原子の持つ点電荷間のクーロン相互作用を計算しなければならない。クーロン力は，分散力に比べて遠距離まで作用するため分子シミュレーション全体の計算精度に対して非常に大きな影響を持っており，原子電荷を精度良く算出することは高精度な分子シミュレーションの必須条件である。
粗視化モデルにおいては，複数の原子を粗視化粒子に融合させ顕わには扱わないことから，全原子モデルにおいて個々の原子が持つ分散力ポテンシャルや点電荷は，全原子モデルから精度を落とすことなく粗視化粒子に統合されなければならない。
融合原子モデルを例に取ると，ゴダードらは様々な原子および融合原子について結合および非結合の相互作用パラメータを決定し，ドライディング（Dreiding）という名称で公開している（非特許文献２）。
しかし上記非特許文献２においては，非結合相互作用パラメータのうち分散力パラメータが与えられているだけであり，原子および融合原子の点電荷についてはラッペらの電荷平衡法（非特許文献３）を用いて計算することが推奨されている。
この手法は，電気陰性度やハードネスなどのパラメータを用いて電荷を計算する方法であり，計算式が比較的単純で高速に点電荷を算出できるメリットがある。
しかし，ラッペらのパラメータは１つの元素について１組のパラメータしか用意されておらず，原子や融合原子，粗視化粒子の置かれた化学的環境については全く考慮されていない。例えば酸素については，エーテル酸素やカルボニル酸素といった様々な環境の酸素原子に対して全て同一のパラメータが用いられるため，エステルのような化合物については適切な電荷を得られなかった。また，フッ素原子や塩を含む化合物では分極を過大評価する傾向が知られており，例えばフッ素系アイオノマーのような分子について適切な電荷を算出できないパラメータであるという問題点があった。さらに，全原子モデル・融合原子モデル・粗視化モデルの区別もなされていないため，いずれのモデルに対しても適切な点電荷を計算することができなかった。
一方，全原子モデルに対しては，量子化学計算から得た静電ポテンシャルに基づいて原子電荷を決定する方法が過去に提唱されてきた（非特許文献４〜９）。量子化学計算に基づく原子電荷は，分子周囲の静電ポテンシャルを精度良く再現する他，双極子モーメントの計算精度にも優れており，分子間相互作用の計算に適した電荷である。また，電荷平衡法に見られるような，特定の系について計算精度が異常に低下するという問題点がない。
このように，量子化学計算に基づく方法は全原子モデルに対して非常に優れた原子電荷決定方法であるが，融合原子モデルをはじめとする粗視化モデルについて量子化学に基づき点電荷を決定する方法は未だ提唱されていない。
極めて単純には，最初に全原子モデルについて原子電荷を決定した後，粗視化粒子を構成する原子の電荷を足し合わせて対応する粗視化粒子の点電荷とする方法が考えられる。しかし上記の方法では，分子の全電荷は保存されるものの，双極子モーメントの値は変わってしまう。さらに，双極子モーメントの精度低下に伴って分子周囲の静電ポテンシャルの精度も低下するため，分子間相互作用を高精度に計算できる点電荷決定方法とは言い難い。
ディー・リグビー（D. Rigby）, アール・ジェイ・ロイ（R. J. Roe）, ジャーナル オブ ケミカル フィジックス（J. Chem. Phys.）, 1987, 87, p.7285. エス・エル・メイオー（S. L. Mayo）, ビー・ディー・オラフソン（B. D. Olafson）, ダブリュー・エー・ゴダード３世（W. A. Goddard III）, ジャーナル オブ フィジカル ケミストリー（J. Phys. Chem.）, 1990, 94, p.8897. エー・ケー・ラッペ（A. K. Rappe）, ダブリュー・エー・ゴダード３世（W. A. Goddard III）, ジャーナル オブ フィジカル ケミストリー（J. Phys. Chem.）, 1991, 95, p.3358. ユー・シー・シン（U. C. Singh）, ピー・エー・コールマン（P. A. Kollman）, ジャーナル オブ コンピューテーショナル ケミストリー（J. Comput. Chem.）, 1984, 5, p.129. エス・ジェー・ウェイナー（S. J. Weiner）, ピー・エー・コールマン（P. A. Kollman）, ディー・ティー・グエン（D. T. Nguyen）, ディー・エー・ケース（D. A. Case）, ジャーナル オブ コンピューテーショナル ケミストリー（J. Comput. Chem.）, 1986, 7, p.230. シー・アイ・ベイリー（C. I. Bayly）, ピー・シープラック（P. Cieplak）, ダブリュー・ディー・コーネル（W. D. Cornell）, ピー・エー・コールマン（P. A. Kollman）, ジャーナル オブ フィジカル ケミストリー（J. Phys. Chem.）, 1993, 97, p.10269. ダブリュー・ディー・コーネル（W. D. Cornell）, ピー・シープラック（P. Cieplak）, シー・アイ・ベイリー（C. I. Bayly）, ピー・エー・コールマン（P. A. Kollman）, アメリカ化学会誌（J. Am. Chem. Soc.）, 1993, 115, p.9620. ティー・フォックス（T. Fox）, ピー・エー・コールマン（P. A. Kollman）, ジャーナル オブ フィジカル ケミストリー（J. Phys. Chem.）, 1998, 102, p.8070. ピー・シープラック（P. Cieplak）, ダブリュー・ディー・コーネル（W. D. Cornell）, シー・アイ・ベイリー（C. I. Bayly）, ピー・エー・コールマン（P. A. Kollman）, ジャーナル オブ コンピューテーショナル ケミストリー（J. Comput. Chem.）, 1995, 16, p.1357.
ディー・リグビー（D. Rigby）, アール・ジェイ・ロイ（R. J. Roe）, ジャーナル オブ ケミカル フィジックス（J. Chem. Phys.）, 1987, 87, p.7285. エス・エル・メイオー（S. L. Mayo）, ビー・ディー・オラフソン（B. D. Olafson）, ダブリュー・エー・ゴダード３世（W. A. Goddard III）, ジャーナル オブ フィジカル ケミストリー（J. Phys. Chem.）, 1990, 94, p.8897. エー・ケー・ラッペ（A. K. Rappe）, ダブリュー・エー・ゴダード３世（W. A. Goddard III）, ジャーナル オブ フィジカル ケミストリー（J. Phys. Chem.）, 1991, 95, p.3358. ユー・シー・シン（U. C. Singh）, ピー・エー・コールマン（P. A. Kollman）, ジャーナル オブ コンピューテーショナル ケミストリー（J. Comput. Chem.）, 1984, 5, p.129. エス・ジェー・ウェイナー（S. J. Weiner）, ピー・エー・コールマン（P. A. Kollman）, ディー・ティー・グエン（D. T. Nguyen）, ディー・エー・ケース（D. A. Case）, ジャーナル オブ コンピューテーショナル ケミストリー（J. Comput. Chem.）, 1986, 7, p.230. シー・アイ・ベイリー（C. I. Bayly）, ピー・シープラック（P. Cieplak）, ダブリュー・ディー・コーネル（W. D. Cornell）, ピー・エー・コールマン（P. A. Kollman）, ジャーナル オブ フィジカル ケミストリー（J. Phys. Chem.）, 1993, 97, p.10269. ダブリュー・ディー・コーネル（W. D. Cornell）, ピー・シープラック（P. Cieplak）, シー・アイ・ベイリー（C. I. Bayly）, ピー・エー・コールマン（P. A. Kollman）, アメリカ化学会誌（J. Am. Chem. Soc.）, 1993, 115, p.9620. ティー・フォックス（T. Fox）, ピー・エー・コールマン（P. A. Kollman）, ジャーナル オブ フィジカル ケミストリー（J. Phys. Chem.）, 1998, 102, p.8070. ピー・シープラック（P. Cieplak）, ダブリュー・ディー・コーネル（W. D. Cornell）, シー・アイ・ベイリー（C. I. Bayly）, ピー・エー・コールマン（P. A. Kollman）, ジャーナル オブ コンピューテーショナル ケミストリー（J. Comput. Chem.）, 1995, 16, p.1357.
このように，粗視化モデルを利用して粒子数を削減しながら分子シミュレーションを実行する際に非結合静電相互作用を精度良く記述できる点電荷を決定する方法は，未だ提唱されていない。そこで本発明は，上記問題点を解決し，分子シミュレーションを精度良く実行するための粗視化分子シミュレーション用点電荷決定法を提供することを目的とする。
上記課題を解決するため，本発明は，分子あるいは分子集合体の高次構造あるいは物性の少なくとも一方を求めるための分子シミュレーションにおいて，（１）分子あるいは分子集合体の構造式を電子計算機に入力するステップ（ステップ１）と，（２）分子のファン・デル・ワールス表面における静電ポテンシャルの空間分布を量子化学計算によって求める静電ポテンシャル決定ステップ（ステップ２）と，（３）分子あるいは分子集合体を構成する原子を２つ以上まとめて粗視化粒子とすることにより分子あるいは分子集合体を構成する単位を減少させる粗視化ステップ（ステップ３）と，（４）粗視化粒子に対して点電荷の初期値を与える初期値設定ステップ（ステップ４）と，（５）粗視化粒子の座標と点電荷から静電ポテンシャルの空間分布を求める近似静電ポテンシャル計算ステップ（ステップ５）と，（６）量子化学による静電ポテンシャルと点電荷による近似静電ポテンシャルとを比較して残差を計算する残差計算ステップ（ステップ６）と，（７）残差の収束を判断し，残差が少なくなるよう点電荷の値を調整する電荷最適化ステップ（ステップ７）とを有し，ステップ１を最初に行い，ステップ３〜ステップ７をこの順に行い，ステップ２をステップ６の前に行うことを特徴とする粗視化分子シミュレーション用点電荷決定方法であることを本旨とし，また，種々の改良された態様を提案する。
本発明の粗視化分子シミュレーション用点電荷決定方法により，従来シミュレーションが困難であった巨大な分子または分子集合体の高精度な粗視化分子シミュレーションが可能になり，新規な機能を持った合成高分子の研究開発や生体高分子の機能解明といった分野においてスピードアップ，ブレイクスルーが期待できる。
量子化学計算の結果に基づいて粗視化粒子の点電荷を計算するにあたって，全原子モデルについて決定した原子電荷を経由して粗視化の操作を行うと静電ポテンシャルや双極子モーメントの記述精度が低下するが，本発明者らは，量子化学計算の結果に基づき直接粗視化モデルについて点電荷を決定することにより，静電ポテンシャルや双極子モーメントを精度良く記述できる粗視化粒子の点電荷を決定できることに気づき，さらに計算に必要となる有効半径パラメータの最適化を行って本発明に到達したものである。
図１に本発明のフローを示す。本発明は分子あるいは分子集合体に適用されるものであるが，以下分子の場合について言及する。
分子の構造式を入力するステップ１においては，電子計算機に読み取り可能な形態で分子の構造式を入力する。入力には，一般に入手可能な分子構造式作成ソフトウェアやテキストエディタなどが利用できる。
量子化学計算により静電ポテンシャルを決定するステップ２においては，ステップ１において入力された分子構造式に対して量子化学計算による構造最適化を実行した後，最適化された分子構造について静電ポテンシャルを量子化学的に計算しその空間分布を出力する。
量子化学計算の手法としては，ヒュッケル分子軌道法，半経験的分子軌道法，Ｘ−α法，密度汎関数法，非経験的分子軌道法といった一般的手法を用いることができるが，分子構造および静電ポテンシャルの計算精度の観点から，密度汎関数法または非経験的分子軌道法が好ましく用いられる。さらに好ましくは，密度汎関数法と非経験的分子軌道法とのハイブリッド法を用いるのが良く，例えば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.〕は本発明の目的のため極めて好ましく用いられる。
密度汎関数法，非経験的分子軌道法あるいはハイブリッド法を用いる場合，１電子軌道を展開する基底関数を指定する必要がある。分子構造および静電ポテンシャルについて十分な計算精度を得るためには，二重価電子軌道（valence double zeta）と呼ばれる基底関数に分極関数を加えたものが好ましく用いられる。例えば，ポープルらによる二重価電子軌道基底関数〔ダブリュー・ジェー・ヘーレ（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.〕との組み合わせである6−31G(d)基底関数や6−31G(d,p)基底関数は，本発明の目的のため好適に用いられる。
静電ポテンシャルを量子化学的に計算しその空間分布を出力する手法としては，MK法〔ビー・エッチ・ベスラー（B. H. Besler）, ケイ・エム・メルツ・ジュニア（K. M. Merz Jr）, ピー・エー・コールマン（P. A. Kollman）, ジャーナル オブ コンピューテーショナル ケミストリー（J. Comput. Chem.）, 1990 11, p.431.〕やCHelpG法〔シー・エム・ブレネマン（C. M. Breneman）, ケイ・ビー・ヴィーベリ（K. B. Wiberg）, ジャーナル オブ コンピューテーショナル ケミストリー（J. Comput. Chem.）, 1990, 11, p.361.〕といった方法を好適に用いることができる。これらの方法を用いることにより，分子のファン・デル・ワールス表面上に密なグリッドを設定し，各グリッド点における量子化学的静電ポテンシャルの値と座標との組を出力として得ることができる。
ここで分子のファン・デル・ワールス表面とは，分子を構成する原子の位置に各原子に対応するファン・デル・ワールス半径の球を置いたとき，これらの球が連結して構成する立体の表面である。そのため，MK法やCHelpG法といった静電ポテンシャル計算手法においては，グリッドの設定にあたってファン・デル・ワールス半径のパラメータを必要とする。この値は，論文やソフトウェアが与える初期値のままでも良いが，粗視化の手法や粗視化粒子の種類によって調整することがより好ましい。例えば，融合原子モデルにおけるメチレン基の点電荷をMK法により決定する場合，ファン・デル・ワールス半径の初期値として水素原子に1.2オングストローム，炭素原子に1.5オングストロームの値が設定されるが，炭素原子に水素原子を融合させたものがメチレン融合原子であることを考慮すれば，メチレン融合原子に対応するファン・デル・ワールス表面上の静電ポテンシャルをより正確に算出するには，水素原子のファン・デル・ワールス半径を初期値より小さく，炭素原子のファン・デル・ワールス半径を初期値より大きくするのが好ましい。
以上のような計算を実行するために，有償または無償で入手可能な種々の量子化学計算プログラムパッケージを利用することができる。例えば，ガウシアン98（Gaussian98）〔ガウシアン98（A11.3改訂版）（Gaussian 98 (Revision A.11.3)）, エム・ジェー・フリッシュ（M. J. Frisch）他, ガウシアン社（Gaussian, Inc.）, ペンシルベニア州ピッツバーグ（Pittsburgh PA）, 2002.〕やGAMESS〔エム・ダブリュー・シュミット（M. W. Schmidt）他, ジャーナル オブ コンピューテーショナル ケミストリー（J. Comput. Chem.）, 1993, 14, p.1347.〕といったプログラムは，本発明の実施に必要な機能を全て備えており好適に用いることができる。
分子モデルを粗視化するステップ３においては，分子を構成する原子を２つ以上まとめて粗視化粒子とし，粗視化粒子の座標や有効粒子径など粗視化粒子の特性を表すパラメータを入力する。例えば融合原子モデルを利用する場合であれば，分子モデルから水素原子を削除するとともに，当該水素原子が結合していた重原子に対して融合させた水素原子の数を標識する操作を行う。この操作には，一般に入手可能な分子設計ソフトウェアやテキストエディタなどが利用できる。また，酸素原子，窒素原子，硫黄原子のいずれかに結合した水素原子は他の酸素，窒素，硫黄原子との間に水素結合を形成し強い分子間相互作用の主体となるため，このような水素原子は粗視化の対象としないことが望ましい。
粗視化粒子の点電荷初期値を設定するステップ４においては，全ての粒子に対してゼロ電荷を設定する方法や，全原子モデルについて原子電荷を決定した後，粗視化粒子を構成する原子の電荷を足し合わせて粗視化粒子の点電荷初期値とする方法などが考えられる。
点電荷から近似静電ポテンシャルを計算するステップ５においては，ステップ３において量子化学的静電ポテンシャルの値が計算されたグリッド点のそれぞれについて，数式（１）により近似静電ポテンシャルViの値を計算する。ここで，iはグリッド点の番号，jは粗視化粒子の番号，qｊは粗視化粒子jの持つ点電荷，rｉｊはグリッド点iと粗視化粒子jとの距離である。
量子化学計算結果と比較し残差を計算するステップ６においては，ステップ３において量子化学的に計算された静電ポテンシャルVi quantumとステップ５において近似的に計算された静電ポテンシャルViとの残差２乗和χESP 2を数式（２）により計算する。
収束判定を行うステップ７においては，残差２乗和が閾値未満に収束したか否かを判断する。収束していれば電荷の値を出力するステップ９に進み，本発明の手順は終了する。
収束していなければ，残差が少なくなるよう点電荷の値を調整するステップ８に進む。電荷の調整にあたっては，最急降下法などの一般的な最適化アルゴリズムを用いることができる。電荷の調整が終了すると，再びステップ５から７の手順を繰り返し，これらの手順を残差２乗和が閾値未満に収束するまで続行する。
本発明の具体的実施態様を以下に実施例をもって述べるが，本発明はこれに限定されるものではない。
実施例１
以下のステップ１〜９の手順により，Ｎ−メチルアセトアミドの融合原子モデルについて点電荷を決定した。
以下のステップ１〜９の手順により，Ｎ−メチルアセトアミドの融合原子モデルについて点電荷を決定した。
ステップ１：分子の構造式の入力
Windows（登録商標）標準のテキストエディタを用いて，化学式（１）に示すＮ−メチルアセトアミドの構造を電子計算機に入力した。
Windows（登録商標）標準のテキストエディタを用いて，化学式（１）に示すＮ−メチルアセトアミドの構造を電子計算機に入力した。
ステップ２：量子化学計算による静電ポテンシャルの決定
ステップ１で作成した分子モデルについて量子化学計算ソフトウェア・ガウシアン98による構造最適化および静電ポテンシャル計算を実行した。近似レベルのキーワードにはB3LYP/6−31G(d,p)，構造最適化のキーワードにはOptを用いた。静電ポテンシャル計算のキーワードにはPop=MKおよびiop(6/33=2)を指定した。計算に必要な閾値やパラメータについては，全てソフトウェアの初期設定値を用いた。以上の量子化学計算は，クロック周波数3.0GHzのPentium（登録商標）4を搭載した電子計算機上で実行し，計算に要した時間は約1分30秒であった。計算の結果，構造最適化された分子モデルの原子座標が得られ，さらに662点のグリッド点について静電ポテンシャル値と座標の組が出力された。また，波動関数から計算した厳密な双極子モーメントの値は3.6227デバイであった。双極子モーメントの計算結果を表１に示す。
ステップ１で作成した分子モデルについて量子化学計算ソフトウェア・ガウシアン98による構造最適化および静電ポテンシャル計算を実行した。近似レベルのキーワードにはB3LYP/6−31G(d,p)，構造最適化のキーワードにはOptを用いた。静電ポテンシャル計算のキーワードにはPop=MKおよびiop(6/33=2)を指定した。計算に必要な閾値やパラメータについては，全てソフトウェアの初期設定値を用いた。以上の量子化学計算は，クロック周波数3.0GHzのPentium（登録商標）4を搭載した電子計算機上で実行し，計算に要した時間は約1分30秒であった。計算の結果，構造最適化された分子モデルの原子座標が得られ，さらに662点のグリッド点について静電ポテンシャル値と座標の組が出力された。また，波動関数から計算した厳密な双極子モーメントの値は3.6227デバイであった。双極子モーメントの計算結果を表１に示す。
ステップ３：分子モデルの粗視化
ステップ２で構造最適化した分子モデルから炭素原子に結合した水素原子（以下非極性水素原子という）を削除し，２つのメチル融合原子および水素原子，炭素原子，窒素原子，酸素原子各１つずつを持つ融合原子モデルを作成した。
ステップ２で構造最適化した分子モデルから炭素原子に結合した水素原子（以下非極性水素原子という）を削除し，２つのメチル融合原子および水素原子，炭素原子，窒素原子，酸素原子各１つずつを持つ融合原子モデルを作成した。
ステップ４：粗視化粒子の点電荷初期値を設定
分子モデルを構成する融合原子および各原子について，電荷ゼロの初期値を設定した。
分子モデルを構成する融合原子および各原子について，電荷ゼロの初期値を設定した。
ステップ５〜８：点電荷の最適化
ベイリーらのRESP法（非特許文献６）が実装され，ステップ５，６，７，８の機能を備えた無償のプログラムであるRESPを用い，ステップ３で作成した融合原子モデルの融合原子および原子の座標について点電荷の値を最適化した。データ入力ファイルには，２つの融合原子および水素原子，炭素原子，窒素原子，酸素原子各１つずつの座標を与え，662点のグリッド点について静電ポテンシャル値と座標の組を入力した。また，最適化計算に必要なパラメータは全てソフトウェアの初期値を用いた。
ベイリーらのRESP法（非特許文献６）が実装され，ステップ５，６，７，８の機能を備えた無償のプログラムであるRESPを用い，ステップ３で作成した融合原子モデルの融合原子および原子の座標について点電荷の値を最適化した。データ入力ファイルには，２つの融合原子および水素原子，炭素原子，窒素原子，酸素原子各１つずつの座標を与え，662点のグリッド点について静電ポテンシャル値と座標の組を入力した。また，最適化計算に必要なパラメータは全てソフトウェアの初期値を用いた。
ステップ９：電荷の値の出力
本発明は分子シミュレーションにおいて正確な分子間相互作用を記述する点電荷決定方法を提供するものである。しかし，点電荷は量子力学における観測可能量（オブザーバブル）ではなく，分子の波動関数から量子力学的に厳密な点電荷の値を定義ないし計算することはできない。そこで，点電荷の値そのものではなく，分子間相互作用の記述の正確さをもって点電荷の良否を評価することとした。評価の指標としては，後述のとおり，点電荷と座標から計算される双極子モーメント，分子シミュレーションによる構造最適化を実行して得られる原子間距離や分子間相互作用エネルギーを用いた。
本発明は分子シミュレーションにおいて正確な分子間相互作用を記述する点電荷決定方法を提供するものである。しかし，点電荷は量子力学における観測可能量（オブザーバブル）ではなく，分子の波動関数から量子力学的に厳密な点電荷の値を定義ないし計算することはできない。そこで，点電荷の値そのものではなく，分子間相互作用の記述の正確さをもって点電荷の良否を評価することとした。評価の指標としては，後述のとおり，点電荷と座標から計算される双極子モーメント，分子シミュレーションによる構造最適化を実行して得られる原子間距離や分子間相互作用エネルギーを用いた。
出力された点電荷の値と座標から双極子モーメントを計算したところ3.6430デバイであり，波動関数から計算した値と良好な一致を示した。双極子モーメントの計算結果を表１に示す。
さらに以下のステップＡ〜Ｄの手順により，Ｎ−メチルアセトアミド二量体の分子シミュレーションを実行し，安定な水素結合配置を求めた。
ステップＡ：原子タイプおよび点電荷の設定
ステップ１で作成した分子モデルについて，非極性水素原子を削除してメチル基を融合原子としたモデルを作成し，アクセルリス社製の分子設計ソフトウェアCerius2（登録商標）を用いてドライディング（非特許文献２）型の原子タイプを設定した。この粗視化により，粒子数を12から6に削減できた。さらに，ステップ９の出力から点電荷の数値を取得し，対応する融合原子および原子の点電荷として設定した。
ステップ１で作成した分子モデルについて，非極性水素原子を削除してメチル基を融合原子としたモデルを作成し，アクセルリス社製の分子設計ソフトウェアCerius2（登録商標）を用いてドライディング（非特許文献２）型の原子タイプを設定した。この粗視化により，粒子数を12から6に削減できた。さらに，ステップ９の出力から点電荷の数値を取得し，対応する融合原子および原子の点電荷として設定した。
ステップＢ：一分子の構造最適化
ステップＡで原子タイプおよび点電荷を設定した分子モデルについて，Cerius2（登録商標）を用いて分子力学法による構造最適化計算を実施し，安定な分子構造および一分子のエネルギーを求めた。上記の構造最適化計算は，IBM社製のRS/6000ワークステーション上で実行した。
ステップＡで原子タイプおよび点電荷を設定した分子モデルについて，Cerius2（登録商標）を用いて分子力学法による構造最適化計算を実施し，安定な分子構造および一分子のエネルギーを求めた。上記の構造最適化計算は，IBM社製のRS/6000ワークステーション上で実行した。
ステップＣ：水素結合二量体モデルの作成
ステップＢで構造最適化した分子モデルを２個準備し，片方の分子のアミド基の水素原子がもう片方のカルボニル基の酸素原子と約２オングストロームの距離で水素結合している配置を作成した。
ステップＢで構造最適化した分子モデルを２個準備し，片方の分子のアミド基の水素原子がもう片方のカルボニル基の酸素原子と約２オングストロームの距離で水素結合している配置を作成した。
ステップＤ：水素結合二量体の構造最適化
ステップＣで作成した水素結合二量体モデルについてステップＢと同様の方法で構造最適化計算を実施し，安定構造および二量体のエネルギーを求めた。構造最適化計算に要した時間は数秒以内であった。安定構造において，水素結合している酸素原子と水素原子との距離は1.978オングストロームであり，酸素原子と水素原子が結合している窒素原子との距離は2.943オングストロームであった。さらに，数式（３）を用いて分子間相互作用エネルギーを計算したところ，−8.584キロカロリー毎モルであった。ここで，相互作用エネルギーの符号が負であれば２つの分子の間に引力がはたらいていることを，エネルギーの絶対値が大きければ強い相互作用がはたらいていることを示している。相互作用エネルギーおよび原子間距離の計算結果を表２に示す。
ステップＣで作成した水素結合二量体モデルについてステップＢと同様の方法で構造最適化計算を実施し，安定構造および二量体のエネルギーを求めた。構造最適化計算に要した時間は数秒以内であった。安定構造において，水素結合している酸素原子と水素原子との距離は1.978オングストロームであり，酸素原子と水素原子が結合している窒素原子との距離は2.943オングストロームであった。さらに，数式（３）を用いて分子間相互作用エネルギーを計算したところ，−8.584キロカロリー毎モルであった。ここで，相互作用エネルギーの符号が負であれば２つの分子の間に引力がはたらいていることを，エネルギーの絶対値が大きければ強い相互作用がはたらいていることを示している。相互作用エネルギーおよび原子間距離の計算結果を表２に示す。
（相互作用エネルギー）＝（二量体のエネルギー）−２×（一分子のエネルギー） （３）
比較例１（量子化学計算による精密な構造最適化）
実施例１の結果を評価するにあたり,比較のため量子化学計算を用いて精密な構造最適化計算を実施した。実施例１のステップ１〜２と同様の方法で，Ｎ−メチルアセトアミドのガウシアン98による構造最適化を実施し，安定構造および一分子のエネルギーを求めた。得られた安定構造を用いてステップＣと同様の方法で水素結合二量体モデルを作成し，ガウシアン98による構造最適化を実施して安定構造および二量体のエネルギーを求めた。二量体の構造最適化計算では，近似レベルのキーワードとしてB3LYP/6−31G(d,p)，構造最適化のキーワードとしてOptを用いた。構造最適化計算に要した時間はおよそ１時間２０分であり，実施例１に比べて遙かに長大な時間を要した。安定構造において，酸素−水素原子間距離は1.980オングストロームであり，酸素−窒素原子間距離は2.988オングストロームであった。さらに，分子間相互作用エネルギーは−7.978キロカロリー毎モルであった。
実施例１の結果を評価するにあたり,比較のため量子化学計算を用いて精密な構造最適化計算を実施した。実施例１のステップ１〜２と同様の方法で，Ｎ−メチルアセトアミドのガウシアン98による構造最適化を実施し，安定構造および一分子のエネルギーを求めた。得られた安定構造を用いてステップＣと同様の方法で水素結合二量体モデルを作成し，ガウシアン98による構造最適化を実施して安定構造および二量体のエネルギーを求めた。二量体の構造最適化計算では，近似レベルのキーワードとしてB3LYP/6−31G(d,p)，構造最適化のキーワードとしてOptを用いた。構造最適化計算に要した時間はおよそ１時間２０分であり，実施例１に比べて遙かに長大な時間を要した。安定構造において，酸素−水素原子間距離は1.980オングストロームであり，酸素−窒素原子間距離は2.988オングストロームであった。さらに，分子間相互作用エネルギーは−7.978キロカロリー毎モルであった。
比較例２（全原子モデルによる，粗視化を行わない計算）
実施例１のステップ１および２と同様の方法でＮ−メチルアセトアミド分子モデルの静電ポテンシャルおよびMK法による全原子モデルの原子電荷を計算した。原子電荷から計算した双極子モーメントは3.5878デバイであり，波動関数から計算した値と良好に一致した。上記の手順で電荷を設定し，さらにドライディング型の原子タイプを設定した分子モデルについて，実施例１のステップＢ〜Ｄと同様の手順で安定構造および相互作用エネルギーを求めた。安定構造において，酸素−水素原子間距離は1.978オングストロームであり，酸素−窒素原子間距離は2.947オングストロームであった。さらに，分子間相互作用エネルギーは−9.139キロカロリー毎モルであった。粗視化を行わなかったため，実施例１に比べ２倍の粒子数を持つモデルであった結果，構造最適化計算には５秒程度の時間を必要とした。
実施例１のステップ１および２と同様の方法でＮ−メチルアセトアミド分子モデルの静電ポテンシャルおよびMK法による全原子モデルの原子電荷を計算した。原子電荷から計算した双極子モーメントは3.5878デバイであり，波動関数から計算した値と良好に一致した。上記の手順で電荷を設定し，さらにドライディング型の原子タイプを設定した分子モデルについて，実施例１のステップＢ〜Ｄと同様の手順で安定構造および相互作用エネルギーを求めた。安定構造において，酸素−水素原子間距離は1.978オングストロームであり，酸素−窒素原子間距離は2.947オングストロームであった。さらに，分子間相互作用エネルギーは−9.139キロカロリー毎モルであった。粗視化を行わなかったため，実施例１に比べ２倍の粒子数を持つモデルであった結果，構造最適化計算には５秒程度の時間を必要とした。
比較例３（点電荷を単に足し合わせて計算）
比較例２と同様の方法でＮ−メチルアセトアミド分子モデルの原子電荷を計算した後，非極性水素原子の電荷を結合する炭素原子に足し合わせてメチル融合原子の電荷とした。点電荷から計算した双極子モーメントは3.9118デバイであり，波動関数から計算した厳密な値に比べて過大評価されていた。以上で決定した点電荷およびドライディング型の原子タイプを分子モデルに対して設定した後，ステップＢ〜Ｄと同様の手順で安定構造および相互作用エネルギーを求めた。安定構造において，酸素−水素原子間距離は1.958オングストロームであり，酸素−窒素原子間距離は2.923オングストロームであった。さらに，分子間相互作用エネルギーは−9.829キロカロリー毎モルであった。構造最適化計算に要した時間は数秒以内であり，実施例１とほぼ同等であった。
比較例２と同様の方法でＮ−メチルアセトアミド分子モデルの原子電荷を計算した後，非極性水素原子の電荷を結合する炭素原子に足し合わせてメチル融合原子の電荷とした。点電荷から計算した双極子モーメントは3.9118デバイであり，波動関数から計算した厳密な値に比べて過大評価されていた。以上で決定した点電荷およびドライディング型の原子タイプを分子モデルに対して設定した後，ステップＢ〜Ｄと同様の手順で安定構造および相互作用エネルギーを求めた。安定構造において，酸素−水素原子間距離は1.958オングストロームであり，酸素−窒素原子間距離は2.923オングストロームであった。さらに，分子間相互作用エネルギーは−9.829キロカロリー毎モルであった。構造最適化計算に要した時間は数秒以内であり，実施例１とほぼ同等であった。
比較例４（全原子モデルの電荷を電荷平衡法を用いて決定）
実施例１のステップ１と同様の方法でＮ−メチルアセトアミド分子モデルを作成し，Cerius2（登録商標）に組み込まれた電荷平衡法（非特許文献３）を用いて原子電荷を設定した。点電荷から求めた双極子モーメントは2.7618デバイであり，分極が過小評価されていた。さらにCerius2（登録商標）を用いてドライディング型の原子タイプを設定した後，実施例１のステップＢ〜Ｄと同様の手順で分子力学計算を実行した。安定構造における酸素−水素原子間距離は2.056オングストロームであり，酸素−窒素原子間距離は3.022オングストロームであった。さらに，分子間相互作用エネルギーは−6.560キロカロリー毎モルであった。構造最適化計算には５秒程度の時間を要した。
実施例１のステップ１と同様の方法でＮ−メチルアセトアミド分子モデルを作成し，Cerius2（登録商標）に組み込まれた電荷平衡法（非特許文献３）を用いて原子電荷を設定した。点電荷から求めた双極子モーメントは2.7618デバイであり，分極が過小評価されていた。さらにCerius2（登録商標）を用いてドライディング型の原子タイプを設定した後，実施例１のステップＢ〜Ｄと同様の手順で分子力学計算を実行した。安定構造における酸素−水素原子間距離は2.056オングストロームであり，酸素−窒素原子間距離は3.022オングストロームであった。さらに，分子間相互作用エネルギーは−6.560キロカロリー毎モルであった。構造最適化計算には５秒程度の時間を要した。
比較例５（融合原子モデルの電荷を電荷平衡法を用いて決定）
比較例４と同様の方法でＮ−メチルアセトアミド分子モデルを作成した後，非極性水素原子を削除してメチル基を融合原子とした。その後，Cerius2（登録商標）に組み込まれた電荷平衡法を用いて点電荷を設定した。点電荷から求めた双極子モーメントは2.9328デバイであり，分極は過小評価されていた。さらにドライディング型の原子タイプを設定し，比較例４と同様の手順で分子力学計算を実行した。安定構造における酸素−水素原子間距離は2.008オングストロームであり，酸素−窒素原子間距離は2.973オングストロームであった。さらに，分子間相互作用エネルギーは−6.849キロカロリー毎モルであった。構造最適化計算に要した時間は数秒以内であり，実施例１とほぼ同等であった。
比較例４と同様の方法でＮ−メチルアセトアミド分子モデルを作成した後，非極性水素原子を削除してメチル基を融合原子とした。その後，Cerius2（登録商標）に組み込まれた電荷平衡法を用いて点電荷を設定した。点電荷から求めた双極子モーメントは2.9328デバイであり，分極は過小評価されていた。さらにドライディング型の原子タイプを設定し，比較例４と同様の手順で分子力学計算を実行した。安定構造における酸素−水素原子間距離は2.008オングストロームであり，酸素−窒素原子間距離は2.973オングストロームであった。さらに，分子間相互作用エネルギーは−6.849キロカロリー毎モルであった。構造最適化計算に要した時間は数秒以内であり，実施例１とほぼ同等であった。
比較例６（全ての水素原子を粗視化した融合原子モデルを用いて計算）
粗視化にあたって全ての水素原子を削除した他は実施例１と同様の方法でＮ−メチルアセトアミド分子モデルを作成した。点電荷から求めた双極子モーメントは3.5228デバイであり，厳密な値に近かった。その後，実施例１と同様の手順で二量体モデルの安定構造および相互作用エネルギーを求めた。安定構造における酸素−窒素原子間距離は3.359オングストロームであり，分子間相互作用エネルギーは−5.576キロカロリー毎モルであった。モデルの粒子数は実施例１より１つ少ない５粒子であったが，構造最適化計算の過程で分子構造が大きく変化したため，計算には５秒以上の時間を要した。
粗視化にあたって全ての水素原子を削除した他は実施例１と同様の方法でＮ−メチルアセトアミド分子モデルを作成した。点電荷から求めた双極子モーメントは3.5228デバイであり，厳密な値に近かった。その後，実施例１と同様の手順で二量体モデルの安定構造および相互作用エネルギーを求めた。安定構造における酸素−窒素原子間距離は3.359オングストロームであり，分子間相互作用エネルギーは−5.576キロカロリー毎モルであった。モデルの粒子数は実施例１より１つ少ない５粒子であったが，構造最適化計算の過程で分子構造が大きく変化したため，計算には５秒以上の時間を要した。
表１に示した実施例１および比較例２〜６の結果について以下に考察する。本発明の方法により計算した融合原子モデルの双極子モーメント（実施例１）は，全原子モデルについて量子化学計算に基づき原子電荷を決定した場合（比較例２）と同等の値であり，波動関数から計算した厳密な値をほぼ再現する結果となった。これに対し，全原子モデルについて決定した点電荷に基づき，非極性水素原子の電荷を結合する炭素原子に足し合わせて融合原子モデルの点電荷とした場合（比較例３）は，双極子モーメントを過大評価していた。また，電荷平衡法を用いた場合（比較例４，５）は双極子モーメントが過小評価されており，比較例３〜５の方法では分子間相互作用を正確に評価できないと考えられる。
次に，表２に示した実施例１および比較例１〜６の結果について以下に考察する。比較例１は量子化学計算のみを用いて計算した厳密な計算結果であるが，この結果を得るためには１時間２０分もの計算時間を必要とした。これと比較して，本発明の方法により作成した融合原子モデルによる計算結果（実施例１）および全原子モデルによる結果（比較例２）は分子間相互作用エネルギー，安定構造における原子間距離の両者ともに厳密な結果に近い値を与え，特に本発明の方法による融合原子モデルでは数秒以下のごく短い計算時間で良好な一致が得られた。一方，比較例３の方法は双極子モーメントの過大評価を反映して，分子間相互作用を過大に算出している。また，比較例４，５の方法は相互作用が過小に計算されているが，これは双極子モーメントの過小評価を反映した結果である。さらに，窒素原子に結合した水素原子まで粗視化した融合原子モデル（比較例６）については，双極子モーメントの計算結果こそ厳密な値との差は小さいものの，二量体の計算結果は量子化学計算と大きな差違が生じている。
以上の結果は，本発明の点電荷決定方法を用いれば，計算時間のごく短い分子シミュレーションを実行することにより，双極子モーメントや分子間相互作用を正確に計算できることを示している。
実施例２
以下のステップ１〜９の手順により，ポリフッ化ビニリデン（以下PvDFと略す）の融合原子モデルについて点電荷を決定した。
以下のステップ１〜９の手順により，ポリフッ化ビニリデン（以下PvDFと略す）の融合原子モデルについて点電荷を決定した。
ステップ１：分子の構造式の入力
Cerius2（登録商標）を用い，PvDFの結晶構造〔長谷川他, ポリマージャーナル（Polymer Journal）, 1972, 3, p.600.〕を電子計算機上に入力した。単位セルの大きさはａ＝8.58，ｂ＝4.91，ｃ＝2.56オングストローム（繊維軸）であり，空間群Cm2mの対称性を持っていた。単位セルの分子構造を図２の左側に，対応する構造式を同図の右側に示す。
Cerius2（登録商標）を用い，PvDFの結晶構造〔長谷川他, ポリマージャーナル（Polymer Journal）, 1972, 3, p.600.〕を電子計算機上に入力した。単位セルの大きさはａ＝8.58，ｂ＝4.91，ｃ＝2.56オングストローム（繊維軸）であり，空間群Cm2mの対称性を持っていた。単位セルの分子構造を図２の左側に，対応する構造式を同図の右側に示す。
次に単位セルを繊維軸方向について８回繰り返したモデルを作成し，その中から１本の分子鎖を取り出して両端をメチル基とした分子モデルを作成した。８回繰り返しモデルおよび両端メチル基の分子モデルを図３に示す。
ステップ２：量子化学計算による静電ポテンシャルの決定
ステップ１で作成した両端メチル基の分子モデルについてガウシアン98による構造最適化および静電ポテンシャル計算を実行した。近似レベルのキーワードにはB3LYP/6−31G(d,p)，構造最適化のキーワードにはOpt=ModRedunを用い，炭素原子の座標を固定して構造最適化を実行した。静電ポテンシャル計算のキーワードにはPop=(MK,ReadRadii)およびiop(6/33=2)を指定した。計算に必要な閾値やパラメータについては，全てソフトウェアの初期設定値を用いた。用いたファン・デル・ワールス半径の値を表３に示す。以上の量子化学計算は，クロック周波数3.0GHzのPentium（登録商標）4を搭載した電子計算機上で実行した。計算の結果，構造最適化された分子モデルの原子座標が得られ，さらに2145点のグリッド点について静電ポテンシャル値と座標の組が出力された。また，波動関数から計算した双極子モーメントの値は13.6829デバイであった。
ステップ１で作成した両端メチル基の分子モデルについてガウシアン98による構造最適化および静電ポテンシャル計算を実行した。近似レベルのキーワードにはB3LYP/6−31G(d,p)，構造最適化のキーワードにはOpt=ModRedunを用い，炭素原子の座標を固定して構造最適化を実行した。静電ポテンシャル計算のキーワードにはPop=(MK,ReadRadii)およびiop(6/33=2)を指定した。計算に必要な閾値やパラメータについては，全てソフトウェアの初期設定値を用いた。用いたファン・デル・ワールス半径の値を表３に示す。以上の量子化学計算は，クロック周波数3.0GHzのPentium（登録商標）4を搭載した電子計算機上で実行した。計算の結果，構造最適化された分子モデルの原子座標が得られ，さらに2145点のグリッド点について静電ポテンシャル値と座標の組が出力された。また，波動関数から計算した双極子モーメントの値は13.6829デバイであった。
ステップ３：分子モデルの粗視化
ステップ２で構造最適化した分子モデルから水素原子を削除し，９つの融合原子，８つの炭素原子および１６個のフッ素原子を持つ融合原子モデルを作成した。作成した融合原子モデルの構造を図４に示す。
ステップ２で構造最適化した分子モデルから水素原子を削除し，９つの融合原子，８つの炭素原子および１６個のフッ素原子を持つ融合原子モデルを作成した。作成した融合原子モデルの構造を図４に示す。
ステップ４：粗視化粒子の点電荷初期値を設定
分子モデルを構成する融合原子およびフッ素原子について，電荷ゼロの初期値を設定した。
分子モデルを構成する融合原子およびフッ素原子について，電荷ゼロの初期値を設定した。
ステップ５〜８：点電荷の最適化
実施例１と同様の方法で，RESPプログラムを用いてステップ３で作成した融合原子モデルの点電荷の値を最適化した。データ入力ファイルには，９つの融合原子，８つの炭素原子および１６個のフッ素原子の座標を与え，2145点のグリッド点について静電ポテンシャル値と座標の組を入力した。その他のパラメータは全てソフトウェアの初期値を用いた。
実施例１と同様の方法で，RESPプログラムを用いてステップ３で作成した融合原子モデルの点電荷の値を最適化した。データ入力ファイルには，９つの融合原子，８つの炭素原子および１６個のフッ素原子の座標を与え，2145点のグリッド点について静電ポテンシャル値と座標の組を入力した。その他のパラメータは全てソフトウェアの初期値を用いた。
ステップ９：電荷の値の出力
出力された点電荷の値と座標から双極子モーメントを計算したところ13.5962デバイであり，波動関数から計算した値と良好に一致した。双極子モーメントの計算結果を表４に示す。
出力された点電荷の値と座標から双極子モーメントを計算したところ13.5962デバイであり，波動関数から計算した値と良好に一致した。双極子モーメントの計算結果を表４に示す。
さらに以下のステップＡ〜Ｃの手順により，PvDF結晶の分子シミュレーションを実行した。
ステップＡ：原子タイプおよび点電荷の設定
ステップ１で作成したPvDF結晶の単位セルから水素原子を削除して融合原子モデルを作成し，Cerius2（登録商標）を用いてドライディング型の原子タイプを設定した。次に，図４において破線で囲った部分の融合原子および原子に相当する点電荷をステップ９の出力から取得し，図２の対応する融合原子および原子の点電荷として設定した。さらに，系の全電荷がゼロになるように平均し，点電荷を決定した。決定した点電荷の値を表５に示す。
ステップ１で作成したPvDF結晶の単位セルから水素原子を削除して融合原子モデルを作成し，Cerius2（登録商標）を用いてドライディング型の原子タイプを設定した。次に，図４において破線で囲った部分の融合原子および原子に相当する点電荷をステップ９の出力から取得し，図２の対応する融合原子および原子の点電荷として設定した。さらに，系の全電荷がゼロになるように平均し，点電荷を決定した。決定した点電荷の値を表５に示す。
ステップＢ：シミュレーションモデル系の作成
ステップＡで点電荷を設定した単位セルについて，ａ，ｂ，ｃ軸方向にそれぞれ３，４，１０倍に拡張したモデル系を作成した。セルパラメータは（ａ，ｂ，ｃ）＝（25.74，19.64，25.60）となり，系に含まれる粒子数は960個となった。このモデルについてCerius2（登録商標）を用いた分子力学計算を50ステップ実行し，系の初期歪みを緩和した。
ステップＡで点電荷を設定した単位セルについて，ａ，ｂ，ｃ軸方向にそれぞれ３，４，１０倍に拡張したモデル系を作成した。セルパラメータは（ａ，ｂ，ｃ）＝（25.74，19.64，25.60）となり，系に含まれる粒子数は960個となった。このモデルについてCerius2（登録商標）を用いた分子力学計算を50ステップ実行し，系の初期歪みを緩和した。
ステップＣ：分子シミュレーションの実行
ステップＢで作成したモデル系について周期境界条件を課し，定温定圧の分子動力学シミュレーションを実行した。ここで，結合相互作用については全てドライディングポテンシャルを用いて計算した。非結合相互作用のうち，レナード・ジョーンズ型ポテンシャルについてはドライディングパラメータを用いカットオフ半径rｃ＝10オングストロームとして計算した。静電相互作用については，実空間部分はカットオフ半径rｃ＝10オングストローム，逆空間部分はα＝0.21オングストローム−１，|n|２ ｍａｘ=30とし，エワルド法を用いて計算した。系の温度は能勢・フーバー（Hoover）法〔能勢修一（S. Nose）, エム・エル・クライン（M. L. Klein）, モレキュラー フィジックス（Mol. Phys.）, 1983, 50, p.1055／能勢修一（S. Nose）, モレキュラー フィジックス（Mol. Phys.）, 1984, 52, p.255／能勢修一（S. Nose）, ジャーナル オブ ケミカル フィジックス（J. Chem. Phys.）, 1984, 81, p.511.〕を用いて25℃に制御し，圧力は斜方セルを用いるパリネロ・ラーマン（Parrinello−Rahman）法〔エム・パリネロ（M. Parrinello）, エー・ラーマン（A. Rahman）, ジャーナル オブ アプライド フィジックス（J. Appl. Phys.）, 1981, 52, p.7182.〕を用いて１気圧に制御した。数値積分に用いる時間刻みは0.5fsとし，シミュレーション時間は200ps（400000ステップ）とした。シミュレーションは，クロック周波数3.0GHzのPentium（登録商標）4を搭載した電子計算機上で自作のプログラムを用いて実行した。
ステップＢで作成したモデル系について周期境界条件を課し，定温定圧の分子動力学シミュレーションを実行した。ここで，結合相互作用については全てドライディングポテンシャルを用いて計算した。非結合相互作用のうち，レナード・ジョーンズ型ポテンシャルについてはドライディングパラメータを用いカットオフ半径rｃ＝10オングストロームとして計算した。静電相互作用については，実空間部分はカットオフ半径rｃ＝10オングストローム，逆空間部分はα＝0.21オングストローム−１，|n|２ ｍａｘ=30とし，エワルド法を用いて計算した。系の温度は能勢・フーバー（Hoover）法〔能勢修一（S. Nose）, エム・エル・クライン（M. L. Klein）, モレキュラー フィジックス（Mol. Phys.）, 1983, 50, p.1055／能勢修一（S. Nose）, モレキュラー フィジックス（Mol. Phys.）, 1984, 52, p.255／能勢修一（S. Nose）, ジャーナル オブ ケミカル フィジックス（J. Chem. Phys.）, 1984, 81, p.511.〕を用いて25℃に制御し，圧力は斜方セルを用いるパリネロ・ラーマン（Parrinello−Rahman）法〔エム・パリネロ（M. Parrinello）, エー・ラーマン（A. Rahman）, ジャーナル オブ アプライド フィジックス（J. Appl. Phys.）, 1981, 52, p.7182.〕を用いて１気圧に制御した。数値積分に用いる時間刻みは0.5fsとし，シミュレーション時間は200ps（400000ステップ）とした。シミュレーションは，クロック周波数3.0GHzのPentium（登録商標）4を搭載した電子計算機上で自作のプログラムを用いて実行した。
ステップＤ：シミュレーション結果の解析
200psのシミュレーションを実行するのに要した時間は4時間41分10秒であった。また，全シミュレーション時間を通して系のハミルトニアンは良好に保存されており，運動方程式の解が安定して得られていることを示していた。シミュレーション終了後，後半100psについてセルパラメータの値を平均した。結果を表６に示す。ａ軸方向に膨らみ，ｃ軸方向に押しつぶされた形となった。また，セルベクトルのなす角度が90°から若干外れており，結晶セルがやや崩れる結果となった。
200psのシミュレーションを実行するのに要した時間は4時間41分10秒であった。また，全シミュレーション時間を通して系のハミルトニアンは良好に保存されており，運動方程式の解が安定して得られていることを示していた。シミュレーション終了後，後半100psについてセルパラメータの値を平均した。結果を表６に示す。ａ軸方向に膨らみ，ｃ軸方向に押しつぶされた形となった。また，セルベクトルのなす角度が90°から若干外れており，結晶セルがやや崩れる結果となった。
実施例３
水素原子のファン・デル・ワールス半径として0.7オングストローム，炭素原子のファン・デル・ワールス半径として2.0オングストロームを指定した他は，実施例２のステップ１および２と同様の方法でPvDF分子モデルの静電ポテンシャル計算を実行し，2464個のグリッド点について静電ポテンシャル値と座標の組を得た。この出力を用いて実施例２のステップ４〜９と同様の方法で点電荷を最適化し，双極子モーメントの値として13.8109デバイを得た。この値は波動関数から計算した値と良好な一致を示している。
水素原子のファン・デル・ワールス半径として0.7オングストローム，炭素原子のファン・デル・ワールス半径として2.0オングストロームを指定した他は，実施例２のステップ１および２と同様の方法でPvDF分子モデルの静電ポテンシャル計算を実行し，2464個のグリッド点について静電ポテンシャル値と座標の組を得た。この出力を用いて実施例２のステップ４〜９と同様の方法で点電荷を最適化し，双極子モーメントの値として13.8109デバイを得た。この値は波動関数から計算した値と良好な一致を示している。
上記の手順で得た点電荷の値を用い，実施例２のステップＡ〜Ｄと同様の手順を用いて分子動力学シミュレーションを実行した。シミュレーションに用いた点電荷の値を表５に，セルパラメータの平均値を表６に示す。200psのシミュレーションを実行するのに要した時間は4時間46分22秒であり，全シミュレーション時間を通して系のハミルトニアンは良好に保存されていた。セルはａ，ｂ軸方向にやや膨らみ，ｃ軸方向にやや縮んでいるが，セルベクトルのなす角はほぼ90°を保っており，結晶セルはかなり良好に保存されていた。
実施例４
水素原子のファン・デル・ワールス半径として0.3オングストローム，炭素原子のファン・デル・ワールス半径として2.4オングストロームを指定した他は，実施例３と同様の方法でPvDF分子モデルの静電ポテンシャル計算を実行し，3027個のグリッド点について静電ポテンシャル値と座標の組を得た。この出力を用いて実施例３と同様の方法で点電荷を最適化し，双極子モーメントの値として，波動関数から計算した値と良好に一致する13.7528デバイを得た。
水素原子のファン・デル・ワールス半径として0.3オングストローム，炭素原子のファン・デル・ワールス半径として2.4オングストロームを指定した他は，実施例３と同様の方法でPvDF分子モデルの静電ポテンシャル計算を実行し，3027個のグリッド点について静電ポテンシャル値と座標の組を得た。この出力を用いて実施例３と同様の方法で点電荷を最適化し，双極子モーメントの値として，波動関数から計算した値と良好に一致する13.7528デバイを得た。
上記の手順で得た点電荷の値を用い，実施例３と同様の手順を用いて分子動力学シミュレーションを実行した。シミュレーションに用いた点電荷の値を表５に，セルパラメータの平均値を表６に示す。シミュレーション実行に要した時間は4時間47分35秒であり，ハミルトニアンの保存は良好であった。セルパラメータは実施例３と同様の傾向を示しており，結晶セルはかなり良好に保存されていた。
比較例７（全原子モデルを用いた計算）
実施例２のステップ１および２と同様の方法でPvDF分子モデルの静電ポテンシャルおよびMK法による原子電荷を計算した。原子電荷から計算した双極子モーメントは13.8386デバイであり，波動関数から計算した値と良好に一致した。ガウシアン98の出力から図３右側のモデルにおいて実線で囲った部分の原子に相当する点電荷を取得し，図２の対応する原子の点電荷として設定してPvDFの全原子モデルを作成した。さらに，系の全電荷がゼロになるように平均し，点電荷を決定した。決定した点電荷の値を表５に示す。さらにCerius2（登録商標）を用いてドライディング型の原子タイプを設定した。
実施例２のステップ１および２と同様の方法でPvDF分子モデルの静電ポテンシャルおよびMK法による原子電荷を計算した。原子電荷から計算した双極子モーメントは13.8386デバイであり，波動関数から計算した値と良好に一致した。ガウシアン98の出力から図３右側のモデルにおいて実線で囲った部分の原子に相当する点電荷を取得し，図２の対応する原子の点電荷として設定してPvDFの全原子モデルを作成した。さらに，系の全電荷がゼロになるように平均し，点電荷を決定した。決定した点電荷の値を表５に示す。さらにCerius2（登録商標）を用いてドライディング型の原子タイプを設定した。
上記の手順で作成したPvDF結晶単位セルについて，実施例２のステップＢ〜Ｄと同様の手順で分子動力学シミュレーションを実行した。系に含まれる原子数は1440個であり，シミュレーション実行に要した時間は8時間50分56秒であって実施例２〜３のほぼ２倍の時間を要した。また，全シミュレーション時間を通してハミルトニアンの保存は良好であった。セルパラメータの平均値は表６に示すとおりであり，ａ，ｂ，ｃ軸すべての方向についてやや膨らむ傾向にあるものの，セルの形状はほぼ完全な直方体を保っており，結晶セルは良好に保存されていた。
比較例８（全原子モデルの原子電荷を電荷平衡法を用いて決定）
実施例２のステップ１と同様の方法でPvDF結晶単位セルを作成し，Cerius2（登録商標）に組み込まれた電荷平衡法を用いて原子電荷を設定した。さらにドライディング型の原子タイプを設定した後，この結晶セルについて実施例２のステップＢ〜Ｄと同様の手順で分子動力学シミュレーションを実行した。系に含まれる原子数は1440個であり，シミュレーション実行に要した時間は8時間40分3秒であって比較例７とほぼ同等の時間を要したが，ハミルトニアンの保存は良好であった。セルパラメータの平均値は表６に示すとおりであり，ａ，ｂ軸方向について膨張，ｃ軸方向について圧縮の傾向が見られた。また，セルベクトルのなす角はαが90°からやや外れており，結晶セルは若干歪んでいた。
実施例２のステップ１と同様の方法でPvDF結晶単位セルを作成し，Cerius2（登録商標）に組み込まれた電荷平衡法を用いて原子電荷を設定した。さらにドライディング型の原子タイプを設定した後，この結晶セルについて実施例２のステップＢ〜Ｄと同様の手順で分子動力学シミュレーションを実行した。系に含まれる原子数は1440個であり，シミュレーション実行に要した時間は8時間40分3秒であって比較例７とほぼ同等の時間を要したが，ハミルトニアンの保存は良好であった。セルパラメータの平均値は表６に示すとおりであり，ａ，ｂ軸方向について膨張，ｃ軸方向について圧縮の傾向が見られた。また，セルベクトルのなす角はαが90°からやや外れており，結晶セルは若干歪んでいた。
比較例９（全原子モデルの原子電荷を全てゼロとして計算）
比較例８と同様の方法でPvDF結晶単位セルを作成し，全ての原子について原子電荷をゼロに設定した。続いて静電相互作用の計算を省略した他は比較例８と同様の手順で分子シミュレーションを実行した。系に含まれる原子数は1440個であるが，シミュレーション実行に要した時間は5時間19分15秒であって，比較例１２に比べて大幅に短かった。ハミルトニアンの保存は良好であった。セルパラメータの平均値は表６に示すとおりであり，ａ，ｂ軸方向について膨張しているが，ｃ軸方向については大きな変化が見られなかった。しかし，セルベクトルのなす角はβが90°からやや外れており，結晶セルは若干歪んでいた。さらに，分子鎖がｃ軸方向に滑っているのが観察され，分子鎖の相対位置が安定しない結果となった。
比較例８と同様の方法でPvDF結晶単位セルを作成し，全ての原子について原子電荷をゼロに設定した。続いて静電相互作用の計算を省略した他は比較例８と同様の手順で分子シミュレーションを実行した。系に含まれる原子数は1440個であるが，シミュレーション実行に要した時間は5時間19分15秒であって，比較例１２に比べて大幅に短かった。ハミルトニアンの保存は良好であった。セルパラメータの平均値は表６に示すとおりであり，ａ，ｂ軸方向について膨張しているが，ｃ軸方向については大きな変化が見られなかった。しかし，セルベクトルのなす角はβが90°からやや外れており，結晶セルは若干歪んでいた。さらに，分子鎖がｃ軸方向に滑っているのが観察され，分子鎖の相対位置が安定しない結果となった。
比較例１０（点電荷を単に足し合わせて計算）
比較例７と同様の方法でPvDF分子モデルの静電ポテンシャルおよびMK法による原子電荷を計算し，非極性水素原子の電荷を結合する炭素原子に足し合わせて融合原子の電荷とした。点電荷から計算した双極子モーメントは7.4089デバイであり，波動関数から計算した値に比べておよそ半分の大きさであった。以上で決定した点電荷を，実施例２のステップＡと同様の方法でPvDF結晶単位セルに対して設定し，ステップＢ〜Ｄと同様の手順で分子シミュレーションを実行した。系に含まれる粒子の数は960個であって，シミュレーション実行に要した時間は4時間43分56秒であって実施例２〜３とほぼ同等であり，ハミルトニアンの保存は良好であった。表６のセルパラメータの表に見るとおり，ａ，ｂ軸方向に膨張，ｃ軸方向に収縮していた。
比較例７と同様の方法でPvDF分子モデルの静電ポテンシャルおよびMK法による原子電荷を計算し，非極性水素原子の電荷を結合する炭素原子に足し合わせて融合原子の電荷とした。点電荷から計算した双極子モーメントは7.4089デバイであり，波動関数から計算した値に比べておよそ半分の大きさであった。以上で決定した点電荷を，実施例２のステップＡと同様の方法でPvDF結晶単位セルに対して設定し，ステップＢ〜Ｄと同様の手順で分子シミュレーションを実行した。系に含まれる粒子の数は960個であって，シミュレーション実行に要した時間は4時間43分56秒であって実施例２〜３とほぼ同等であり，ハミルトニアンの保存は良好であった。表６のセルパラメータの表に見るとおり，ａ，ｂ軸方向に膨張，ｃ軸方向に収縮していた。
比較例１１（融合原子モデルの原子電荷を電荷平衡法を用いて決定）
実施例２のステップ１と同様の方法でPvDF結晶単位セルを作成した後水素原子を削除して融合原子モデルとし，Cerius2（登録商標）に組み込まれた電荷平衡法を用いて原子電荷を設定した後，Cerius2（登録商標）を用いてドライディング型の原子タイプを設定した。上記手順により作成した結晶単位セルについて，実施例２のステップＢ〜Ｄと同様の手順で分子動力学シミュレーションを実行した。系に含まれる原子数は960個であり，シミュレーション実行に要した時間は5時間5分56秒であって実施例２〜３に比べてやや長い時間を要したが，ハミルトニアンの保存は良好であった。セルパラメータの平均値は表６に示すとおりであり，セルベクトルのなす角はほぼ90°を保っているが，ａ軸方向について膨張，ｂ，ｃ軸方向について収縮の傾向が見られ，特にｃ軸方向に強く収縮していた。
実施例２のステップ１と同様の方法でPvDF結晶単位セルを作成した後水素原子を削除して融合原子モデルとし，Cerius2（登録商標）に組み込まれた電荷平衡法を用いて原子電荷を設定した後，Cerius2（登録商標）を用いてドライディング型の原子タイプを設定した。上記手順により作成した結晶単位セルについて，実施例２のステップＢ〜Ｄと同様の手順で分子動力学シミュレーションを実行した。系に含まれる原子数は960個であり，シミュレーション実行に要した時間は5時間5分56秒であって実施例２〜３に比べてやや長い時間を要したが，ハミルトニアンの保存は良好であった。セルパラメータの平均値は表６に示すとおりであり，セルベクトルのなす角はほぼ90°を保っているが，ａ軸方向について膨張，ｂ，ｃ軸方向について収縮の傾向が見られ，特にｃ軸方向に強く収縮していた。
比較例１２（融合原子モデルの原子電荷を全てゼロとして計算）
実施例２のステップ１と同様の方法でPvDF結晶単位セルを作成した後水素原子を削除して融合原子モデルとし，全ての原子および融合原子について電荷ゼロを設定した後，Cerius2（登録商標）を用いてドライディング型の原子タイプを設定した。上記手順により作成した結晶単位セルについて，静電相互作用の計算を省略した他は実施例２のステップＢ〜Ｄと同様の手順で分子動力学シミュレーションを実行した。系に含まれる原子数は960個であり，シミュレーション実行に要した時間は2時間40分53秒であって，実施例２〜３のおよそ６割の時間で計算が終了した。ハミルトニアンの保存は良好であった。セルパラメータの平均値は表６に示すとおりであり，角α，βが90°から大きく外れており，結晶セルは強く歪んでいた。さらに，分子鎖がｃ軸方向に滑っているのが観察され，分子鎖の相対位置が安定しない結果となった。
実施例２のステップ１と同様の方法でPvDF結晶単位セルを作成した後水素原子を削除して融合原子モデルとし，全ての原子および融合原子について電荷ゼロを設定した後，Cerius2（登録商標）を用いてドライディング型の原子タイプを設定した。上記手順により作成した結晶単位セルについて，静電相互作用の計算を省略した他は実施例２のステップＢ〜Ｄと同様の手順で分子動力学シミュレーションを実行した。系に含まれる原子数は960個であり，シミュレーション実行に要した時間は2時間40分53秒であって，実施例２〜３のおよそ６割の時間で計算が終了した。ハミルトニアンの保存は良好であった。セルパラメータの平均値は表６に示すとおりであり，角α，βが90°から大きく外れており，結晶セルは強く歪んでいた。さらに，分子鎖がｃ軸方向に滑っているのが観察され，分子鎖の相対位置が安定しない結果となった。
表４に示した実施例２〜４および比較例７，１０の結果について以下に考察する。本発明の方法により計算した融合原子モデルの双極子モーメント（実施例２〜４）は，全原子モデルについて量子化学計算に基づき原子電荷を決定した場合（比較例７）と同等の値であり，波動関数から計算した厳密な値をほぼ再現する結果となった。これに対し，全原子モデルについて決定した点電荷に基づき，非極性水素原子の電荷を結合する炭素原子に足し合わせて融合原子モデルの点電荷とした場合（比較例１０）は，厳密な結果の半分程度の値にしかならず，双極子モーメントの再現性は極めて不十分であった。
次に，表６に示した実施例２〜４および比較例７〜１２の結果について以下に考察する。融合原子モデルを用いたシミュレーション（実施例２〜４および比較例１０〜１２）は，全原子モデルによるシミュレーション（比較例７〜９）に比べて，計算速度の点でおよそ２倍有利であった。しかし，点電荷の決定方法として，非極性水素原子の電荷を融合原子に足し合わせる方法（比較例１０）や電荷平衡法（比較例１１）を用いた場合，さらに電荷をゼロとした場合（比較例１２）では，繊維軸方向への収縮やセルの歪みが見られ，精度の良い計算は実行できなかった。全原子モデルを用いたシミュレーションでは，量子化学計算により原子電荷を決定した場合（比較例７）は結晶セルが良好に保存されていたが，電荷平衡法（比較例８）やゼロ電荷（比較例９）を用いた場合は結晶セルの膨張・収縮や歪みが見られ，電荷を精密に決定する方法の必要性を示している。これらの結果に比較して，実施例２は結晶セル保存の点でやや劣る結果を与えるものの，実施例３および４は比較例７の２倍の計算速度でありながら同等の計算精度であり，本発明の電荷決定法を用いた分子シミュレーションが精度を落とすことなく計算速度を向上できることを示している。
実施例５
数値積分に用いる時間刻みを1.0fsとした他は，実施例３と同様の手順で分子シミュレーションを実行した。全シミュレーション時間を通して系のハミルトニアンは良好に保存されており，運動方程式の解が安定して得られていることを示していた。シミュレーション実行に要した時間は2時間22分41秒で，実施例３のおよそ半分の時間であった。
数値積分に用いる時間刻みを1.0fsとした他は，実施例３と同様の手順で分子シミュレーションを実行した。全シミュレーション時間を通して系のハミルトニアンは良好に保存されており，運動方程式の解が安定して得られていることを示していた。シミュレーション実行に要した時間は2時間22分41秒で，実施例３のおよそ半分の時間であった。
実施例６
数値積分に用いる時間刻みを2.0fsとした他は，実施例３と同様の手順で分子シミュレーションを実行した。ハミルトニアンは良好に保存されており，シミュレーション実行に要した時間は1時間12分9秒で，実施例３のおよそ４分の１の時間であった。
数値積分に用いる時間刻みを2.0fsとした他は，実施例３と同様の手順で分子シミュレーションを実行した。ハミルトニアンは良好に保存されており，シミュレーション実行に要した時間は1時間12分9秒で，実施例３のおよそ４分の１の時間であった。
実施例７
数値積分に用いる時間刻みを4.0fsとした他は，実施例３と同様の手順で分子シミュレーションを実行した。シミュレーション時間の経過とともにハミルトニアンの値が若干増加しており，やや不安定な解が得られていることを示していた。シミュレーション実行に要した時間は36分32秒で，実施例３のおよそ８分の１の時間であったが，シミュレーション結果の信頼性は実施例３および５，６に比べて劣ると考えられる。
数値積分に用いる時間刻みを4.0fsとした他は，実施例３と同様の手順で分子シミュレーションを実行した。シミュレーション時間の経過とともにハミルトニアンの値が若干増加しており，やや不安定な解が得られていることを示していた。シミュレーション実行に要した時間は36分32秒で，実施例３のおよそ８分の１の時間であったが，シミュレーション結果の信頼性は実施例３および５，６に比べて劣ると考えられる。
実施例８
数値積分に用いる時間刻みを1.0fsとした他は，実施例４と同様の手順で分子シミュレーションを実行した。ハミルトニアンは良好に保存されており，シミュレーション実行に要した時間は2時間24分10秒で，実施例４のおよそ半分の時間であった。
数値積分に用いる時間刻みを1.0fsとした他は，実施例４と同様の手順で分子シミュレーションを実行した。ハミルトニアンは良好に保存されており，シミュレーション実行に要した時間は2時間24分10秒で，実施例４のおよそ半分の時間であった。
実施例９
数値積分に用いる時間刻みを2.0fsとした他は，実施例４と同様の手順で分子シミュレーションを実行した。ハミルトニアンは良好に保存されており，シミュレーション実行に要した時間は1時間11分55秒で，実施例４のおよそ４分の１の時間であった。
数値積分に用いる時間刻みを2.0fsとした他は，実施例４と同様の手順で分子シミュレーションを実行した。ハミルトニアンは良好に保存されており，シミュレーション実行に要した時間は1時間11分55秒で，実施例４のおよそ４分の１の時間であった。
実施例１０
数値積分に用いる時間刻みを4.0fsとした他は，実施例４と同様の手順で分子シミュレーションを実行した。シミュレーション時間の経過とともにハミルトニアンの値が若干増加しており，やや不安定な解が得られていることを示していた。シミュレーション実行に要した時間は36分37秒で，実施例４のおよそ８分の１の時間であったが，シミュレーション結果の信頼性は実施例４および８，９より劣ると考えられる。
数値積分に用いる時間刻みを4.0fsとした他は，実施例４と同様の手順で分子シミュレーションを実行した。シミュレーション時間の経過とともにハミルトニアンの値が若干増加しており，やや不安定な解が得られていることを示していた。シミュレーション実行に要した時間は36分37秒で，実施例４のおよそ８分の１の時間であったが，シミュレーション結果の信頼性は実施例４および８，９より劣ると考えられる。
比較例１３（全原子モデルで数値積分の時間刻みを1.0fsとした分子シミュレーション）
数値積分に用いる時間刻みを1.0fsとした他は，比較例７と同様の手順で分子シミュレーションを実行した。ハミルトニアンは良好に保存されており，シミュレーション実行に要した時間は4時間26分11秒で，比較例７のおよそ半分の時間であった。
数値積分に用いる時間刻みを1.0fsとした他は，比較例７と同様の手順で分子シミュレーションを実行した。ハミルトニアンは良好に保存されており，シミュレーション実行に要した時間は4時間26分11秒で，比較例７のおよそ半分の時間であった。
比較例１４（全原子モデルで数値積分の時間刻みを2.0fsとした分子シミュレーション）
数値積分に用いる時間刻みを2.0fsとした他は，比較例７と同様の手順で分子シミュレーションを実行した。シミュレーション時間の経過とともにハミルトニアンの値が若干減少しており，やや不安定な解が得られていることを示していた。シミュレーション実行に要した時間は2時間13分31秒で，比較例７のおよそ４分の１の時間であったが，シミュレーション結果の信頼性は比較例７および１３より劣ると考えられる。
数値積分に用いる時間刻みを2.0fsとした他は，比較例７と同様の手順で分子シミュレーションを実行した。シミュレーション時間の経過とともにハミルトニアンの値が若干減少しており，やや不安定な解が得られていることを示していた。シミュレーション実行に要した時間は2時間13分31秒で，比較例７のおよそ４分の１の時間であったが，シミュレーション結果の信頼性は比較例７および１３より劣ると考えられる。
比較例１５（全原子モデルで数値積分の時間刻みを4.0fsとした分子シミュレーション）
数値積分に用いる時間刻みを4.0fsとした他は，比較例７と同様の手順で分子シミュレーションを実行した。シミュレーション開始直後に運動エネルギーの値が異常に大きくなり，正常にシミュレーションを続行することができなかった。
数値積分に用いる時間刻みを4.0fsとした他は，比較例７と同様の手順で分子シミュレーションを実行した。シミュレーション開始直後に運動エネルギーの値が異常に大きくなり，正常にシミュレーションを続行することができなかった。
実施例３および５〜７，実施例４および８〜１０，比較例７および１３〜１５のシミュレーション結果を表７に示す。融合原子を用いた分子シミュレーションにおいては，時間刻みを全原子モデルの２倍の長さに設定した場合でも運動方程式を安定に解くことが可能であり，粒子数削減の効果と併せておよそ４倍の計算速度を実現できることが示された。
１：比較例７において分子シミュレーション用の電荷を採取した部分
２：実施例２から１０において分子シミュレーション用の電荷を採取した部分
２：実施例２から１０において分子シミュレーション用の電荷を採取した部分
Claims (5)
- 分子あるいは分子集合体の高次構造あるいは物性の少なくとも一方を求めるための分子シミュレーションにおいて，（１）分子あるいは分子集合体の構造式を電子計算機に入力するステップ（ステップ１）と，（２）分子のファン・デル・ワールス表面における静電ポテンシャルの空間分布を量子化学計算によって求める静電ポテンシャル決定ステップ（ステップ２）と，（３）分子あるいは分子集合体を構成する原子を２つ以上まとめて粗視化粒子とすることにより分子あるいは分子集合体を構成する単位を減少させる粗視化ステップ（ステップ３）と，（４）粗視化粒子に対して点電荷の初期値を与える初期値設定ステップ（ステップ４）と，（５）粗視化粒子の座標と点電荷から静電ポテンシャルの空間分布を求める近似静電ポテンシャル計算ステップ（ステップ５）と，（６）量子化学による静電ポテンシャルと点電荷による近似静電ポテンシャルとを比較して残差を計算する残差計算ステップ（ステップ６）と，（７）残差の収束を判断するステップ（ステップ７），（８）残差が少なくなるよう点電荷の値を調整する電荷最適化ステップ（ステップ８），（９）電荷の値を出力するステップ（ステップ９）とを有し，ステップ１を最初に行い，ステップ３〜ステップ９をこの順に行い，ステップ２をステップ６の前に行うことを特徴とする粗視化分子シミュレーション用点電荷決定方法。
- ステップ２をステップ３の後に行い，ステップ２において，ファン・デル・ワールス半径を，ステップ３においてその原子がまとめられる粗視化粒子の種類に応じて与える請求項１に記載の粗視化分子シミュレーション用点電荷決定方法。
- ステップ２において，量子化学計算の手法が非経験的分子軌道法または密度汎関数法またはこれらのハイブリッド法である請求項１または２に記載の点電荷決定方法。
- ステップ３における粗視化手段が融合原子モデルである請求項１から３のいずれかに記載の点電荷決定方法。
- 融合原子モデルにおいて酸素原子，窒素原子，硫黄原子のいずれかに結合した水素原子は粗視化しないことを特徴とする請求項４に記載の点電荷決定方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2006292293A JP2007149075A (ja) | 2005-10-27 | 2006-10-27 | 粗視化分子シミュレーション用点電荷決定方法 |
Applications Claiming Priority (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2005312471 | 2005-10-27 | ||
| JP2006292293A JP2007149075A (ja) | 2005-10-27 | 2006-10-27 | 粗視化分子シミュレーション用点電荷決定方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| JP2007149075A true JP2007149075A (ja) | 2007-06-14 |
Family
ID=38210368
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2006292293A Pending JP2007149075A (ja) | 2005-10-27 | 2006-10-27 | 粗視化分子シミュレーション用点電荷決定方法 |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP2007149075A (ja) |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2010108183A (ja) * | 2008-10-29 | 2010-05-13 | Ihi Corp | シミュレーション方法及びシミュレーション装置 |
| JP2010107314A (ja) * | 2008-10-29 | 2010-05-13 | Ihi Corp | シミュレーション方法及びシミュレーション装置 |
| WO2011148639A1 (ja) * | 2010-05-25 | 2011-12-01 | 株式会社ブリヂストン | 分子間力のシミュレーション手法 |
| CN103539057A (zh) * | 2012-07-17 | 2014-01-29 | 常政 | 一种高危高纯化学液体分配系统 |
| JP2016066284A (ja) * | 2014-09-25 | 2016-04-28 | 日本電気株式会社 | 情報入力手段、情報入力手段形成方法、表計算装置、表計算方法及びプログラム |
| WO2023108622A1 (zh) * | 2021-12-17 | 2023-06-22 | 深圳晶泰科技有限公司 | 获得电荷参数的方法、分子力学模拟结果的方法及装置 |
-
2006
- 2006-10-27 JP JP2006292293A patent/JP2007149075A/ja active Pending
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2010108183A (ja) * | 2008-10-29 | 2010-05-13 | Ihi Corp | シミュレーション方法及びシミュレーション装置 |
| JP2010107314A (ja) * | 2008-10-29 | 2010-05-13 | Ihi Corp | シミュレーション方法及びシミュレーション装置 |
| WO2011148639A1 (ja) * | 2010-05-25 | 2011-12-01 | 株式会社ブリヂストン | 分子間力のシミュレーション手法 |
| CN103539057A (zh) * | 2012-07-17 | 2014-01-29 | 常政 | 一种高危高纯化学液体分配系统 |
| JP2016066284A (ja) * | 2014-09-25 | 2016-04-28 | 日本電気株式会社 | 情報入力手段、情報入力手段形成方法、表計算装置、表計算方法及びプログラム |
| WO2023108622A1 (zh) * | 2021-12-17 | 2023-06-22 | 深圳晶泰科技有限公司 | 获得电荷参数的方法、分子力学模拟结果的方法及装置 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Rowe et al. | Development of a machine learning potential for graphene | |
| Du et al. | Recent progress in robust and quality Delaunay mesh generation | |
| JP2007149075A (ja) | 粗視化分子シミュレーション用点電荷決定方法 | |
| Vogiatzis et al. | Local segmental dynamics and stresses in polystyrene–C60 mixtures | |
| Ghoniem et al. | Multiscale modelling of nanomechanics and micromechanics: an overview | |
| Jirásek et al. | Consistent tangent stiffness for nonlocal damage models | |
| Chen et al. | A novel Volume-Compensated Particle method for 2D elasticity and plasticity analysis | |
| Al-Raoush et al. | Simulation of random packing of polydisperse particles | |
| Nguyen et al. | A Burger model for the effective behavior of a microcracked viscoelastic solid | |
| Duan et al. | Machine-learning assisted coarse-grained model for epoxies over wide ranges of temperatures and cross-linking degrees | |
| Ullah et al. | An adaptive finite element/meshless coupled method based on local maximum entropy shape functions for linear and nonlinear problems | |
| Tam et al. | Moisture effect on interfacial integrity of epoxy-bonded system: a hierarchical approach | |
| WO2011148639A1 (ja) | 分子間力のシミュレーション手法 | |
| Hetenyi et al. | Multiple “time step” Monte Carlo | |
| Martyn et al. | Efficient fully-coherent quantum signal processing algorithms for real-time dynamics simulation | |
| Ohkuma et al. | Equilibrating high-molecular-weight symmetric and miscible polymer blends with hierarchical back-mapping | |
| Bar-Yam | About complex systems | |
| Lu et al. | Molecular dynamics simulations of plastic damage in metals | |
| Schollwöck | Time-dependent density-matrix renormalization-group methods | |
| Sala et al. | Fitting properties from density functional theory based molecular dynamics simulations to parameterize a rigid water force field | |
| Savrasov et al. | Ab initio calculations of lattice dynamics of crystals | |
| JP5312299B2 (ja) | 輸送係数の算出方法、輸送係数の算出装置、および輸送係数の算出プログラム | |
| Guo et al. | Small-data-based machine learning interatomic potentials for graphene grain boundaries enabled by structural unit model | |
| Nguyen et al. | An improved meshless method for finite deformation problem in compressible hyperelastic media | |
| Kay | IVR formulation of Miller's correspondence relations |