JP2009037334A - 分子シミュレーション方法、分子シミュレーション装置、分子シミュレーションプログラム、及び該プログラムを記録した記録媒体 - Google Patents

Abstract

【課題】安定かつ適切に物性値又は物理量を求める。
【解決手段】予め設定された形状の物性値又は物理量をシミュレーションから求める分子シミュレーション方法において、前記形状に原子を配置する原子配置ステップと、前記原子配置ステップにより配置された原子の位置に基づいて原子間ポテンシャルを取得する原子間ポテンシャル取得ステップと、前記原子間ポテンシャル取得ステップにより得られる原子間ポテンシャルに基づいて分子動力学計算を行い、物性値又は物理量を取得する分子動力学計算ステップとを有することにより、上記課題を解決する。
【選択図】図２

Description

本発明は、分子シミュレーション方法、分子シミュレーション装置、分子シミュレーションプログラム、及び該プログラムを記録した記録媒体に係り、特に安定かつ適切に物性値又は物理量を求めるための分子シミュレーション方法、分子シミュレーション装置、分子シミュレーションプログラム、及び該プログラムを記録した記録媒体に関する。

従来、古典力学や量子力学等を基に計算機を用いて物質科学全般の現象を探るための方法論として、分子動力学法に基づく分子シミュレーションが知られている。分子シミュレーションは、分子のポテンシャルエネルギーや最安定構造等、物質の特性を分子レベルで解明することができるものである。

ところで、分子シミュレーションは、主に材料の物性予測に用いられていたため、複雑な形状や一般的な機械構造物等に原子を配置する方法は存在しなかった。なお、物体の複雑な形状を有限の領域に分割してモデル化し、計算する手法としては、例えば、有限要素法がある。また、２次元の計算の場合には、三角形あるいは四角形の領域で物体を満遍なく区切るためにメッシュの生成プログラムを使用することがある（例えば、特許文献１参照。）。
特開２００１−６７４９５号公報

ところで、上述した特許文献１に示されるメッシュ生成手法を、例えば、分子動力学の原子配列に応用すると、三角形要素のみを使用する場合、物体領域を三角形で満たすことを目的としているため，三角形の大きさは一定にはならず、異なってしまう場合が多い。

また、原子を配置する形状によっては、三角形が非常に歪んだ状態となることがある。したがって、単純にメッシュの格子点に原子を配置すると、初期原子配置の粗密の具合によっては極めて大きな応力が発生し、最悪の場合、物理量の計算が破綻することもある。

本発明は、上述した問題点に鑑みなされたものであり、安定かつ適切に物性値又は物理量を求めるための分子シミュレーション方法、分子シミュレーション装置、分子シミュレーションプログラム、及び該プログラムを記録した記録媒体を提供することを目的とする。

上述の目的を達成するために、本発明は、予め設定された形状の物性値又は物理量をシミュレーションから求める分子シミュレーション方法において、前記形状に原子を配置する原子配置ステップと、前記原子配置ステップにより配置された原子の位置に基づいて原子間ポテンシャルを取得する原子間ポテンシャル取得ステップと、前記原子間ポテンシャル取得ステップにより得られる原子間ポテンシャルに基づいて分子動力学計算を行い、物性値又は物理量を取得する分子動力学計算ステップとを有することを特徴とする。これにより、安定かつ適切に物性値又は物理量を求めることができる。したがって、原子配置の情報によって適宜ポテンシャルを修正し、物性値又は物理量が再現されるように工夫することにより、複雑形状を表現しても非物理的な残留応力等が生じず安定した計算を実現することができる。

更に、前記原子配置ステップは、前記形状にメッシュを生成し、生成したメッシュの交点に原子を配置することが好ましい。これにより、メッシュを生成して原子を容易かつ効率的に配置することができる。

更に、前記原子間ポテンシャルは、２個の原子間の距離のみを独立変数とする２体間ポテンシャルであることが好ましい。これにより、最も簡単なポテンシャル形式を用いて効率的に２体間ポテンシャルを計算することができる。

更に、前記原子間ポテンシャル取得ステップは、前記原子配置ステップにより配置された各原子点を用いたボロノイ多面体解析により得られる各ボロノイ点に基づいて、２つの原子間の力の働く断面積及び原子間の結合エネルギーを含むポテンシャルパラメータを取得し、前記ポテンシャルパラメータを用いて原子間ポテンシャルを決定することが好ましい。これにより、適切な物性値又は物理量を算出でき、適宜原子間ポテンシャルを修正することができる。

また本発明は、予め設定された形状の物性値又は物理量をシミュレーションから求める分子シミュレーション装置において、前記形状に原子を配置する原子配置手段と、前記原子配置手段により配置された原子の位置に基づいて原子間ポテンシャルを取得する原子間ポテンシャル取得手段と、前記原子間ポテンシャル取得手段により得られる原子間ポテンシャルに基づいて分子動力学計算を行い、物性値又は物理量を取得する分子動力学計算手段とを有することを特徴とする。これにより、安定かつ適切に物性値又は物理量を求めることができる。したがって、原子配置の情報によって適宜ポテンシャルを修正し、物性値又は物理量が再現されるように工夫することにより、複雑形状を表現しても非物理的な残留応力等が生じず安定した計算を実現することができる。

更に、前記原子配置手段は、前記形状にメッシュを生成し、生成したメッシュの交点に原子を配置することが好ましい。これにより、メッシュを生成して原子を容易かつ効率的に配置することができる。

更に、前記原子間ポテンシャルは、２個の原子間の距離のみを独立変数とする２体間ポテンシャルであることが好ましい。これにより、最も簡単なポテンシャル形式を用いて効率的に２体間ポテンシャルを計算することができる。

更に、前記原子間ポテンシャル取得手段は、前記原子配置手段により配置された各原子点を用いたボロノイ多面体解析により得られる各ボロノイ点に基づいて、２つの原子間の力の働く断面積及び原子間の結合エネルギーを含むポテンシャルパラメータを取得し、前記ポテンシャルパラメータを用いて原子間ポテンシャルを決定することが好ましい。これにより、適切な物性値又は物理量を算出でき、適宜原子間ポテンシャルを修正することができる。

また本発明は、予め設定された形状の物性値又は物理量をシミュレーションから求める分子シミュレーションプログラムにおいて、コンピュータを、前記形状に原子を配置する原子配置手段、前記原子配置手段により配置された原子の位置に基づいて原子間ポテンシャルを取得する原子間ポテンシャル取得手段、及び、前記原子間ポテンシャル取得手段により得られる原子間ポテンシャルに基づいて分子動力学計算を行い、物性値又は物理量を取得する分子動力学計算手段として機能させる。これにより、安定かつ適切に物性値又は物理量を求めることができる。したがって、原子配置の情報によって適宜ポテンシャルを修正し、物性値又は物理量が再現されるように工夫することにより、複雑形状を表現しても非物理的な残留応力等が生じず安定した計算を実現することができる。また、実行プログラムをコンピュータにインストールすることにより、容易に分子シミュレーション処理を実現することができる。

更に、前記原子配置手段は、前記形状にメッシュを生成し、生成したメッシュの交点に原子を配置することが好ましい。これにより、メッシュを生成して原子を容易かつ効率的に配置することができる。

更に、前記原子間ポテンシャルは、２個の原子間の距離のみを独立変数とする２体間ポテンシャルであることが好ましい。これにより、最も簡単なポテンシャル形式を用いて効率的に２体間ポテンシャルを計算することができる。

更に、前記原子間ポテンシャル取得手段は、前記原子配置手段により配置された各原子点を用いたボロノイ多面体解析により得られるボロノイ点に基づいて、２つの原子間の力の働く断面積及び原子間の結合エネルギーを含むポテンシャルパラメータを取得し、前記ポテンシャルパラメータを用いて原子間ポテンシャルを決定することが好ましい。これにより、適切な物性値又は物理量を算出でき、適宜原子間ポテンシャルを修正することができる。

また、本発明は前記請求項９乃至請求項１２に記載された分子シミュレーションプログラムが記録されたコンピュータ読み取り可能な記録媒体を提供する。これにより、記録媒体により他の複数のコンピュータに容易に分子シミュレーションプログラムをインストールすることができる。

本発明によれば、安定かつ適切に物性値又は物理量を求めることができる。これにより、原子配置の情報によって適宜ポテンシャルを修正し、物性値が再現されるように工夫することにより、複雑形状を表現しても非物理的な残留応力等が生じず安定した計算を実現することができる。

<本発明の概要>
本発明では、例えば、汎用のメッシュ生成ソフト等を用いて複雑形状に原子を配置した場合に、原子配置はそのままとし、その周辺原子の情報によってポテンシャルを適宜修正する方法を提供する。なお、このときの計算方法については、例えば、分子動力学法を用いることができるが、本発明においてはこれに限定されるものではなく、例えばモンテカルロ法等、ポテンシャルの安定点を探す方法であればよい。

なお、分子動力学法（Ｍｏｌｅｃｕｌａｒ Ｄｙｎａｍｉｃｓ ｍｅｔｈｏｄ、ＭＤ法）とは、２体（あるいはそれ以上）の原子間ポテンシャルの下に、古典力学におけるニュートン方程式を解いて、系の静的、動的安定構造や、動的過程（ダイナミクス）を解析する手法である。分子動力学により、定温、定圧、定エネルギー、定積、定ケミカルポテンシャル等の様々なアンサンブル（統計集団）の計算が可能である。また、結合長や位置の固定等の様々な拘束条件を付加することも可能である。

以下に、上述したような特徴を有する本発明における分子シミュレーション方法、分子シミュレーション装置、分子シミュレーションプログラム、及び該プログラムを記録した記録媒体を好適に実施した形態について、図面を用いて詳細に説明する。

<本実施例における装置構成>
まず、本実施例における分子シミュレーション装置について図を用いて説明する。図１は、本実施例における分子シミュレーション装置のハードウェア構成の一例を示す図である。

図１における分子シミュレーション装置１０は、入力装置１１と、出力装置１２と、ドライブ装置１３と、補助記憶装置１４と、メモリ装置１５と、各種制御を行うＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）１６と、ネットワーク接続装置１７とを有するよう構成されており、これらはシステムバスＢで相互に接続されている。

入力装置１１は、ユーザ等が操作するキーボードやマウス等のポインティングデバイスを有しており、ユーザ等からのプログラムの実行や処理結果の表示指示、各処理に必要なデータの入力等、各種操作信号を入力する。また、出力装置１２は、本発明における処理を行うためのコンピュータ本体を操作するのに必要な各種ウィンドウやデータ等を表示するディスプレイを有し、ＣＰＵ１６が有する制御プログラムによりプログラムの実行経過や結果等を表示することができる。

ここで、本発明において、コンピュータ本体にインストールされる実行プログラムは、例えばＣＤ−ＲＯＭ等の記録媒体１８等により提供される。プログラムを記録した記録媒体１８は、ドライブ装置１３にセット可能であり、記録媒体１８に含まれる実行プログラムが、記録媒体１８からドライブ装置１３を介して補助記憶装置１４にインストールされる。

補助記憶装置１４は、ハードディスク等のストレージ手段であり、本発明における実行プログラムや、コンピュータに設けられた制御プログラム、本発明における分子シミュレーションに必要な入力データや処理結果等を蓄積し、更に必要に応じて入出力を行うことができる。

メモリ装置１５は、ＣＰＵ１６により補助記憶装置１４から読み出された実行プログラム等を格納する。なお、メモリ装置１５は、ＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）やＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）等からなる。

ＣＰＵ１６は、ＯＳ（Ｏｐｅｒａｔｉｎｇ Ｓｙｓｔｅｍ）等の制御プログラム、メモリ装置１５により読み出され格納されている実行プログラム等に基づいて、各種演算や各ハードウェア構成部とのデータの入出力等、コンピュータ全体の処理を制御して各処理を実現することができる。また、ＣＰＵ１６は、プログラムの実行中に必要な各種情報を補助記憶装置１４から取得することができ、また処理結果等を補助記憶装置１４に格納することができる。

ネットワーク接続装置１７は、通信ネットワーク等と接続することにより、実行プログラムを通信ネットワークに接続されている他の端末等から取得したり、プログラムを実行することで得られた実行結果又は本発明における実行プログラム自体を他の端末等に提供することができる。

<分子シミュレーション装置：機能構成例>
次に、分子シミュレーション装置１０の機能構成例について図を用いて説明する。図２は、分子シミュレーション装置の機能構成の一例を示す図である。分子シミュレーション装置１０は、入力手段２１と、出力手段２２と、蓄積手段２３と、原子配置手段２４と、原子間ポテンシャル取得手段２５と、分子動力学計算手段２６と、画面生成手段２７と、送受信手段２８と、制御手段２９とを有するよう構成されている。

入力手段２１は、ユーザ等から分子シミュレーションの実施に伴う各指示情報等の入力を受け付ける。なお、入力手段２１は、例えばキーボードや、マウス等のポインティングデバイス、マイク等の音声入力インタフェイス等からなる。

また、出力手段２２は、入力手段２１により入力された指示内容や、指示内容に基づいて生成された分子シミュレーションの編集内容、編集結果、実行プログラムの実施指示画面や、原子間ポテンシャルパラメータの途中経過等の内容を表示したり、音声出力を行う。なお、出力手段２２は、ディスプレイやスピーカ等からなる。

また、蓄積手段２３は、分子シミュレーションした原子の配置結果や、原子間ポテンシャル配置結果等の結果、入力データ等を蓄積する。なお、蓄積手段２３は入力データ等を、送受信手段２８を介して通信ネットワーク等により接続された外部装置（図示せず）等から取得して蓄積することもできる。

原子配置手段２４は、例えば、予め設定された形状（単純形状や複雑形状）に対して原子を配置する。なお、適切な位置に原子を配置する場合には、例えば、予め設定された形状に多数のメッシュを生成し、生成されたメッシュの各格子点に原子を配置するか、あるいは、予め設定された形状の枠に原子を配置した後に形状内部に原子を発生させて充填することで原子を配置する手法等を用いることができる。

したがって、例えば、上述のメッシュ生成手法を用いる場合、原子配置手段２４は、具体的には、メッシュ生成により得られるメッシュ交点に原子を配置する。

原子間ポテンシャル取得手段２５は、原子配置手段２４により配置された原子を、原子配置はそのままとし、その周辺原子の情報によってポテンシャルを適宜修正する。具体的には、メッシュ生成により配置した原子に対してボロノイ（Ｖｏｒｏｎｏｉ）多面体解析を行い、ポテンシャルパラメータを決定する。

分子動力学計算手段２６は、原子間ポテンシャル取得手段２５から得られるポテンシャルパラメータに基づいて、分子動力学による物性値又は物理量を計算する。

画面生成手段２７は、原子配置手段２４により予め設定された形状に配置された原子の様子や、原子間ポテンシャル取得手段２５における初期設定や実行結果の表示、分子動力学計算手段２６による計算結果の物性値やその値に対応する表やグラフ、図等の処理結果の表示等を行う。

送受信手段２８は、例えばインターネット等の通信ネットワーク等を介して接続される外部装置等から必要なデータを取得したり、本発明における分子シミュレーション結果等を送信するための通信インタフェイスである。

制御手段２９は、分子シミュレーション装置１０における各機能構成全体の制御を行う。具体的には、制御手段２９は、入力手段２１により入力されたユーザからの入力情報に基づいて原子配置手段２４により原子を配置させたり、原子間ポテンシャル取得手段２５によりポテンシャルパラメータを取得させたり、分子動力学計算手段２６により原子の物性値を取得させる等といった処理を制御する。

これにより、分子動力学を用いて複雑な形状に対しても適切に原子を配置することができる。また、安定かつ適切に物性値又は物理量を求めることができる。したがって、対象とする形状が複雑な場合の分子動力学法の計算において、適当な原子配置であっても安定に計算を行い、かつ適切な物性値又は物理量を算出できるよう原子間ポテンシャルを修正することができる。

<分子シミュレーション処理手順例>
次に、本発明における実行プログラム（分子シミュレーションプログラム）による分子シミュレーション処理手順について説明する。図３は、本実施形態における分子シミュレーション処理手順の一例を示すフローチャートである。

図３に示すように、本実施形態における分子シミュレーションでは、例えば、複雑形状等の予め設定された対象形状に原子を配置する場合、まず、ＣＡＤ（Ｃｏｍｐｕｔｅｒ Ａｉｄｅｄ Ｄｅｓｉｇｎ）データ等からの対象形状を読み込み（Ｓ０１）、メッシュの生成を行う（Ｓ０２）。なお、メッシュの生成は、例えば、上述した特許文献１に示すようなメッシュ生成手法を行ってもよい。

また、メッシュ生成手法については、例えば、デラウニー法（ボロノイ分割に基づく三角形分割）、Ｔｒａｎｓｆｉｎｉｔｅマッピング法（マッピングを４角形以外の領域に拡張）や、バウンダリーフィット法（楕円形偏微分方程式を解く）、アドバンシングフロント法（ある境界から３角形（４面体）を生成）、４分木法、８分木法、ペービング法（ある境界から順に４辺形メッシュを配置していき、矛盾が出た時点で修正）、６面体メッシュ生成法（Ｗｈｉｓｋｅｒ Ｗｅａｖｉｎｇ法（表面の４角形分割データから内部に６面体メッシュを生成していく手法）、Ｇｒｉｄ Ｂａｓｅｄ法（立体の内部に立方体のグリッドを配置して、その表面を処理）、Ｍｅｄｉａｌ Ｓｕｒｆａｃｅを用いた空間分割法（マッピングによる６面体生成が可能な基本領域に分割する。なお、全ての頂点において、稜線が３本である必要がある。））、アダプティブ法（解析結果に基づき離散化誤差の分布を推定し、誤差を減らすようにメッシュを改善する手法等を用いることができる。

また、メッシュ生成においては、例えば、汎用のＣＡＥ（Ｃｏｍｐｕｔｅｒ Ａｉｄｅｄ Ｅｎｇｉｎｅｅｒｉｎｇ）ソフトウェア（例えば、「ＰＡＴＲＡＮ」、「ＦＥＭＡＰ」、「ＡＴＬＡＳ」、「ＭＥＮＴＡＴ」、「ＴＲＵＥＧＲＩＤ」、「ＳＯＬＩＤＷＯＲＫＳ」等）を用いることもできる。

次に、Ｓ０２の処理にて生成したメッシュの各交点に対して原子を配置し（Ｓ０３）、ボロノイ（Ｖｏｒｏｎｏｉ）多面体解析を行う（Ｓ０４）。具体的には、原子の配置に対して平面でない点をどの点に最も近いかによって分割してできるボロノイ図（領域）を求め、その図（ボロノイ領域）の境目の線によりボロノイ境界を求め、その境界の各交点によりボロノイ点を求める。

また、Ｓ０４の処理により、原子間ポテンシャルパラメータを決定する（Ｓ０５）。なお、Ｓ０５の処理では、予めユーザが設定する所望のヤング率等に基づいて、後述する最安定となる原子間距離「ｒ_０」や結合エネルギー「ε」等のパラメータを決定する。

また、Ｓ０５の処理に得られる原子間ポテンシャルパラメータに基づいて、分子動力学における物性値又は物理量を計算する（Ｓ０６）。また、Ｓ０６の処理により得られた処理結果等をディスプレイ等の出力装置１２（出力手段２２）に表示する画面を生成し（Ｓ０７）、生成した画面を表示する（Ｓ０８）。

<Ｓ０４〜Ｓ０６の処理について>
次に、上述したＳ０４（ボロノイ多面体解析）〜Ｓ０６（分子動力学計算）の処理について、具体的に説明する。

まず、ボロノイ多面体について説明する。例えば、平面Ｒ^２上にｎ個の点ｐ_１，ｐ_２，・・・，ｐ_ｎが与えられたときのｐ_ｉの勢力図Ｖ_ｉを以下に示す（１）式で定義する。

このとき、上述の（１）式におけるｘの取り得る領域をボロノイ領域といい、ｐ_ｉはボロノイ領域の母点（本実施形態では原子位置）である。また、このボロノイ領域は、凸多面体となり、それをボロノイ多面体という。Ｓ０４では、このボロノイ多面体に用いて解析を行う。なお、具体的な解析手法の例については後述する。

次に、分子動力学法の概念について説明する。分子動力学法で実際に解くべき運動方程式は、以下に示す（２）式のニュートン（Ｎｅｗｔｏｎ）の運動方程式であり、直接数値積分して解くものである。

なお、上述した（２）式において、ｍ_ｉは原子の質量であり、ｒ_ｉは原子の位置ベクトル、Ｆ_ｉは原子ｉに働く力のベクトル（原子間相互作用の合力）である。また、速度ｖ_ｉは、以下に示す（３）式のように表すことができる。

したがって、時間刻みをΔｔとすると、以下に示す（４）式から速度ｖ_ｉを求めることができる。

また、このときの位置は、以下に示す（５）式から求めることができる。

ここで、原子の初期位置及び初期速度を決めておけば、上述した（４）式，（５）式を用いて原子の位置・速度を順次更新することができる。なお、この数値積分法を蛙飛び法という（詳細は、例えば、岡田勲著「分子シミュレーション入門」、海文堂、１９８９年、Ｐ．４７に記載されている）。

また、原子に働く力Ｆ_ｉは、以下に示す（６）式から求めることができる。

このとき、φは原子間ポテンシャル、Ｎは全原子数である。ここで、上述した（４）式，（５）式を用いて速度・位置を求めるためには、原子間力のみを求めればよく、原子間力を求めるためには原子間ポテンシャルが決定できればよい。

したがって、本実施形態では、原子間ポテンシャルだけ決まれば、分子動力学の計算を行うことが可能となる。しかしながら、このポテンシャルの決定は、経験的要素が強く、容易でない場合が多い。したがって、本発明では、周辺原子の配置に合わせてポテンシャルを決定し、物性を再現する。

なお、本実施形態において、ポテンシャルが２個の原子間の距離のみを独立変数として定義されるものを２体間ポテンシャルといい、最も簡単なポテンシャル形式である。また、この２体間ポテンシャルを式で表すと、一般的には通常（７）式のような形で表すことができる。

なお、「ε」と「σ」は、それぞれ結合エネルギーと原子の大きさを表すポテンシャルパラメータであり、ｒ_０は最安定となる原子間距離であり、ｒ_ｉｊは実際の原子間距離である。

ここで、このポテンシャルを用いて（系のエネルギーが最も低くなる結晶構造である）面心立方格子を保ったまま複雑形状を表現しようとすると、以下に示す問題によりポテンシャルを修正する必要が生じる。図４は、複雑形状を表現する際のポテンシャルの修正理由を説明するための一例を示す図である。

上述したポテンシャルを用いて面心立方体格子を保ったまま、例えば、複雑形状３１に原子３２を配置すると、図４（ａ）に示すように、どうしても表面が凹凸を有する形状になってしまう。また、図４（ｂ）に示すように、最初に複雑形状３１に沿って原子３２を配置し、その後、内部に原子を充填させていくと、確かに表面は滑らかに表現できるが、このままでは残留応力や空孔が発生してしまう。そのため、エネルギーが安定した状態とはならない。したがって、周辺原子の配置に合わせてポテンシャルを修正する必要がある。

そこで、本実施形態では、ヤング率を基にしたパラメータの導出法を行う。以下に示す例では、代表的な２体のポテンシャルであるＬｅｎｎａｒｄ−Ｊｏｎｅｓポテンシャルの修正法を適用した例について説明する。このポテンシャルは、以下に示す（８）式のように表すことができる。

次に、上述した（８）式を２次の項までｒ＝ｒ_０の周りでテイラー（Ｔａｙｌｏｒ）展開すると、弾性領域を記述するポテンシャルが以下に示す（９）式のように表すことができる。

なお、上述の（９）式において、φ^（１），φ^（２）は、それぞれ１回微分、２回微分したことを示している。また、ｒ_０は最安定原子間距離であり、Ｌｅｎｎａｒｄ−Ｊｏｎｅｓポテンシャルでは、以下に示す（１０）式のように表すことができる。

ここで、図５は、Ｌｅｎｎａｒｄ−Ｊｏｎｅｓポテンシャルの例を示す図である。なお、図５において、横軸は原子間距離（ｒ）を示し、縦軸は原子間ポテンシャル（φ）を示している。図５に示すように、ｒ＝ｒ_０のときに原子間ポテンシャルが極小値−εをとる。

また、φ^（１），φ^（２）は、以下に示す（１１）式，（１２）式で表されることから、上述した（１０）式を用いると、以下に示す（１３）式のようになる。

したがって、上述した（１１）式，（１２）式，（１３）式を上述した（９）式に代入すると、以下に示す（１４）式に示すような調和振動子のポテンシャルの形となる。

ここで、バネ定数をｋとすると、上述の（９）式より、以下に示す（１５）式が導出される。

更に、バネ定数とヤング率との関係は以下に示す（１６）式となる。

なお、Ｓは力の働く断面積である。ここで、上述した（１５）式，（１６）式より、結合エネルギーεは、以下に示す（１７）式により求められる。

ここで、ｒ_０は、初期配置時の原子間距離であり本実施形態では最安定原子間距離であるが、Ｓは別途決める必要がある。具体的には、上述したようにボロノイ多面体解析（Ｓ０４）等を行い、原子間距離を表す線分ｒ_０を横切る面積をＳとすればよい。

<ボロノイ多面体解析について>
ここで、ボロノイ多面体解析について、具体的に図を用いて説明する。図６は、本実施形態におけるボロノイ多面体解析を説明するための一例を示す図である。なお、図６（ａ）に示す例では、原子４１−１〜４１−７を有し、原子４１−１を母点として原子４１−１の周りにある全ての粒子（原子）４１−２〜４１−７をそれぞれ直線で結ぶ。また、上述の処理にて得られた線分４２−１〜４２−６に垂直で原子点間距離が等しくなる２等分線（３次元の場合には、２等分面）を生成する。更に、上述の処理にて得られた２等分線４３−１〜４３−５により囲まれる多角形（３次元の場合には、多面体）が、ボロノイ多面体となる。したがって、上述したように定義された多面体から面積Ｓを求めることができる。

ここで、上述した内容を原子配置手段２４等により配置された全ての原子について行うことができればよいが、例えば、対象形状の表面側にある原子については、原子数が足りず上述したボロノイ多面体を定義することができない。

そこで、本実施形態では、図６（ｂ）に示すような処理を行う。なお、図６（ｂ）に示す原子５１−１〜５１−３のうち、原子５１−１，５１−２は、対象形状の表面側に存在する原子（表面原子）とし、また、原子５１−３は、表面原子に対して内部にある原子（内部原子）とする。

このとき、各原子間を結ぶ線分５２−１，５２−２，５２−３の垂直２等分線の交点を点ｐとする。この点ｐは、ボロノイ多面体を形成する頂点の一つとなるため、ボロノイ点とする。

ここで、このボロノイ点ｐから線分５２−１に下ろした垂線５３と、線分５２−１との交点を点ｍとする。また、点ｍに対して垂線５３上に存在し、点ｐ，ｍ間の距離と等しい距離にある点を点ｑとする。このようにして定義した点ｐ，ｑに基づいて、原子５１−１と原子５１−２との相互作用には線分ｐｑで定義された面積Ｓ、具体的には線分ｐｑに奥行きと長さを乗算したものを定義することができる。

なお、３次元の場合には、線分の代わりに面を扱う必要があるが手順は同じである。したがって、シミュレーション処理実行時にユーザ等により予めヤング率を与えることにより煩雑なポテンシャルパラメータの決定に悩む必要がなく、原子間ポテンシャルパラメータを決定することができ、そのポテンシャルパラメータに基づいて、上述した分子動力学の計算を行い、その計算結果に基づいて物性値又は物理量を高精度に決定することができる。したがって、上述にて定めたポテンシャルと、一般的な繰り込み理論を分子動力学法に適用することによりマクロな系の計算が可能になる。

<本実施形態のシミュレーションによる精度結果>
次に、本実施形態のシミュレーションによる精度結果について図を用いて説明する。なお、以下に示す例では、片持ち梁における曲げと振動との関係を例に取って、本実施形態のシミュレーション手法の精度確認を行う。なお、梁の材質は、例えばＳＵＳ等を用いることができる。

図７は、２次元系片持ち梁の曲げと振動のシミュレーション結果の一例を示す図である。また、図８は、３次元系片持ち梁の曲げと振動のシミュレーション結果の一例を示す図である。なお、以下に示す図７，図８は、上述した画面生成手段２７により生成され、出力手段２２により出力される図の一例を示すものでもある。

ここで、図７（ａ）は、２次元における横からの荷重をかけた梁の様子を示し、図７（ｂ）は、梁の曲げ（Ｌｅｎｇｔｈ）とそり（Ｄｅｆｌｅｃｔｉｏｎ）との関係のグラフの一例を示し、図７（ｃ）は、梁の振動（周期）（Ｔｉｍｅ）と振幅（Ａｍｐｌｉｔｕｄｅ）との関係のグラフの一例を示している。

なお、図７（ａ）では、横からの荷重１００［Ｎ］、ヤング率：２０８［ＧＰａ］、密度：７８００［ｋｇ／ｍ^３］、横幅×高さ×奥行き＝１０×１００×１［ｍｍ］の梁の様子を示している。このとき、図７（ｂ）に示すように、繰り込み分子動力学計算により得られた値（ＲＭＤ：図７（ｂ）における点部）が、理論値の値（厳密解）とよく一致しているのがわかる。また、図７（ｃ）においては、理論周期とＲＭＤの周期（図７（ｃ）における実線部）とがよく一致しているのがわかる。

また、図８において、図８（ａ）は３次元における荷重をかけた梁の様子を示し、図８（ｂ）は梁の曲げ（Ｌｅｎｇｔｈ）そり（Ｄｅｆｌｅｃｔｉｏｎ）との関係のグラフの一例を示し、図８（ｃ）は梁の振動（周期）（Ｔｉｍｅ）と振幅（Ａｍｐｌｉｔｕｄｅ）との関係のグラフの一例を示している。

なお、図８（ａ）では、横からの荷重４．２［Ｎ］、ヤング率：１００［ＧＰａ］、密度：７８００［ｋｇ／ｍ^３］、横幅×高さ×奥行き＝１×１０×１［ｍｍ］の梁の様子を示している。このとき、図８（ｂ）に示すように、繰り込み分子動力学計算により得られた値（ＲＭＤ：図８（ｂ）における点部）が、理論値の値（厳密解）とよく一致しているのがわかる。また、図８（ｃ）においては、理論周期と、ＲＭＤの周期（図８（ｃ）における実線部）とがよく一致しているのがわかる。

上述したように本発明によれば、安定かつ適切に物性値又は物理量を求めることができる。これにより、原子配置の情報によって適宜ポテンシャルを修正し、物性値が再現されるように工夫することにより、複雑形状を表現しても非物理的な残留応力等が生じず安定した計算を実現することができる。

つまり、本発明によれば、対象とする形状が複雑な場合の分子動力学法の計算において、適当な原子配置であっても安定に計算を行うことができ、かつ適切な物理量を算出できるよう原子間ポテンシャルを修正することができる。

また、本発明を適用することにより、機構・弾性の動解析における摩擦による発熱、破壊の計算や、原理的には、量子的な現象を除く全ての古典現象を計算することができる。更に、構造解析や流体解析、運動解析、最適設計、電磁気、熱物質の分野等においても本発明を広く適用することができる。

以上、本発明の好ましい実施例について詳述したが、本発明は係る特定の実施形態に限定されるものではなく、特許請求の範囲に記載された本発明の要旨の範囲内において、種々の変形、変更が可能である。

本実施例における分子シミュレーション装置のハードウェア構成の一例を示す図である。 分子シミュレーション装置の機能構成の一例を示す図である。 本実施形態における分子シミュレーション処理手順の一例を示すフローチャートである。 複雑形状を表現する際のポテンシャルの修正理由を説明するための一例を示す図である。 Ｌｅｎｎａｒｄ−Ｊｏｎｅｓポテンシャルの例を示す図である。 本実施形態におけるボロノイ多面体解析を説明するための一例を示す図である。 ２次元系片持ち梁の曲げと振動のシミュレーション結果の一例を示す図である。 ３次元系片持ち梁の曲げと振動のシミュレーション結果の一例を示す図である。

符号の説明

１０ 分子シミュレーション装置
１１ 入力装置
１２ 出力装置
１３ ドライブ装置
１４ 補助記憶装置
１５ メモリ装置
１６ ＣＰＵ
１７ ネットワーク接続装置
１８ 記録媒体
２１ 入力手段
２２ 出力手段
２３ 蓄積手段
２４ 原子配置手段
２５ 原子間ポテンシャル取得手段
２６ 分子動力学計算
２７ 画面生成手段
２８ 送受信手段
２９ 制御手段
３１ 複雑形状
３２，４１，５１ 原子
４２，５２ 線分
４３ ２等分線
５３ 垂線

Claims (13)

1. 予め設定された形状の物性値又は物理量をシミュレーションから求める分子シミュレーション方法において、
   前記形状に原子を配置する原子配置ステップと、
   前記原子配置ステップにより配置された原子の位置に基づいて原子間ポテンシャルを取得する原子間ポテンシャル取得ステップと、
   前記原子間ポテンシャル取得ステップにより得られる原子間ポテンシャルに基づいて分子動力学計算を行い、物性値又は物理量を取得する分子動力学計算ステップとを有することを特徴とする分子シミュレーション方法。
2. 前記原子配置ステップは、
   前記形状にメッシュを生成し、生成したメッシュの交点に原子を配置することを特徴とする請求項１に記載の分子シミュレーション方法。
3. 前記原子間ポテンシャルは、２個の原子間の距離のみを独立変数とする２体間ポテンシャルであることを特徴とする請求項１又は２に記載の分子シミュレーション方法。
4. 前記原子間ポテンシャル取得ステップは、
   前記原子配置ステップにより配置された各原子点を用いたボロノイ多面体解析により得られる各ボロノイ点に基づいて、２つの原子間の力の働く断面積及び原子間の結合エネルギーを含むポテンシャルパラメータを取得し、前記ポテンシャルパラメータを用いて原子間ポテンシャルを決定することを特徴とする請求項１乃至３の何れか１項に記載の分子シミュレーション方法。
5. 予め設定された形状の物性値又は物理量をシミュレーションから求める分子シミュレーション装置において、
   前記形状に原子を配置する原子配置手段と、
   前記原子配置手段により配置された原子の位置に基づいて原子間ポテンシャルを取得する原子間ポテンシャル取得手段と、
   前記原子間ポテンシャル取得手段により得られる原子間ポテンシャルに基づいて分子動力学計算を行い、物性値又は物理量を取得する分子動力学計算手段とを有することを特徴とする分子シミュレーション装置。
6. 前記原子配置手段は、
   前記形状にメッシュを生成し、生成したメッシュの交点に原子を配置することを特徴とする請求項５に記載の分子シミュレーション装置。
7. 前記原子間ポテンシャルは、２個の原子間の距離のみを独立変数とする２体間ポテンシャルであることを特徴とする請求項５又は６に記載の分子シミュレーション装置。
8. 前記原子間ポテンシャル取得手段は、
   前記原子配置手段により配置された各原子点を用いたボロノイ多面体解析により得られる各ボロノイ点に基づいて、２つの原子間の力の働く断面積及び原子間の結合エネルギーを含むポテンシャルパラメータを取得し、前記ポテンシャルパラメータを用いて原子間ポテンシャルを決定することを特徴とする請求項５乃至７の何れか１項に記載の分子シミュレーション装置。
9. 予め設定された形状の物性値又は物理量をシミュレーションから求める分子シミュレーションプログラムにおいて、
   コンピュータを、
   前記形状に原子を配置する原子配置手段、
   前記原子配置手段により配置された原子の位置に基づいて原子間ポテンシャルを取得する原子間ポテンシャル取得手段、及び、
   前記原子間ポテンシャル取得手段により得られる原子間ポテンシャルに基づいて分子動力学計算を行い、物性値又は物理量を取得する分子動力学計算手段として機能させるための分子シミュレーションプログラム。
10. 前記原子配置手段は、
    前記形状にメッシュを生成し、生成したメッシュの交点に原子を配置することを特徴とする請求項９に記載の分子シミュレーションプログラム。
11. 前記原子間ポテンシャルは、２個の原子間の距離のみを独立変数とする２体間ポテンシャルであることを特徴とする請求項９又は１０に記載の分子シミュレーションプログラム。
12. 前記原子間ポテンシャル取得手段は、
    前記原子配置手段により配置された各原子点を用いたボロノイ多面体解析により得られるボロノイ点に基づいて、２つの原子間の力の働く断面積及び原子間の結合エネルギーを含むポテンシャルパラメータを取得し、前記ポテンシャルパラメータを用いて原子間ポテンシャルを決定することを特徴とする請求項９乃至１１の何れか１項に記載の分子シミュレーションプログラム。
13. 前記請求項９乃至請求項１２に記載された分子シミュレーションプログラムが記録されたコンピュータ読み取り可能な記録媒体。