JP3892167B2 - 粒子集団の配置を生成する生成装置および方法 - Google Patents
粒子集団の配置を生成する生成装置および方法 Download PDFInfo
- Publication number
- JP3892167B2 JP3892167B2 JP7498799A JP7498799A JP3892167B2 JP 3892167 B2 JP3892167 B2 JP 3892167B2 JP 7498799 A JP7498799 A JP 7498799A JP 7498799 A JP7498799 A JP 7498799A JP 3892167 B2 JP3892167 B2 JP 3892167B2
- Authority
- JP
- Japan
- Prior art keywords
- potential function
- displacement
- position coordinates
- value
- arrangement
- 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.)
- Expired - Fee Related
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
- G16C10/00—Computational theoretical chemistry, i.e. ICT specially adapted for theoretical aspects of quantum chemistry, molecular mechanics, molecular dynamics or the like
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- Computing Systems (AREA)
- Mathematical Physics (AREA)
- General Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- General Engineering & Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Pure & Applied Mathematics (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Mathematical Analysis (AREA)
- Computational Mathematics (AREA)
- Mathematical Optimization (AREA)
- Algebra (AREA)
- General Health & Medical Sciences (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Description
【０００１】
【発明の属する技術分野】
本発明は、粒子系のシミュレーションにおいて、粒子集団の配置を表す情報を生成する生成装置およびその方法に関する。
【０００２】
【従来の技術】
無機物または有機物で構成される材料の性質の解析・予測を行う方法の１つとして、分子動力学（molecular dynamics，ＭＤ）計算に基づくシミュレーションがある。この計算を行うためには、材料を構成する原子・分子の初期配置を生成する必要がある。ＭＤ計算においては、与えられた初期配置を初期値として、ニュートン（Newton）方程式で表されるような微分方程式を解き、計算結果を用いて材料を表す系の性質を調べ、性質の予測にまでつなげることが試みられている。
【０００３】
しかしながら、１分子の安定状態と考えられる分子構造や、常温・常圧における質量密度の実験データは得られていても、分子集団の配置（各分子を構成するすべての原子の位置座標）の実験データは得られていないことが多い。したがって、分子系のＭＤ計算を行う際、実験的に得られた質量密度を実現する分子集団の初期配置をどのようにして生成したらよいか不明である。
【０００４】
粒子集団の初期配置を生成する簡単な方法として、質量密度の実験データから一定のセルサイズと１セルあたりの粒子数を求め、セル内に粒子をランダムに配置する方法が考えられる。しかし、粒子を全くランダムに配置すると、しばしば、接近した粒子対が生成されてしまう。これらの接近した粒子間には、通常、非常に強い力が働くことが知られている。
【０００５】
このような接近した粒子対が１つでも存在すると、相互作用のためにそれらの粒子の速度が速くなり、その近傍で系が局所的に非常に高温になる。このため、微分方程式を解く際に莫大な数値が発生し、数値積分アルゴリズムが、しばしば、破綻を引き起こす。
【０００６】
このようなシミュレーションの破綻を回避する従来の方法として、次の３つの方法が挙げられる。通常は、これらの方法を組み合わせて用いることが多い。
（１）熱排出法
シミュレーション中に発生した系の熱を適当なタイミングで排出する方法である。具体的には、各粒子の速度を低くするか、または完全に零にしてしまう。
（２）ポテンシャル緩和法
粒子間距離の小さいところで、ポテンシャル関数を変形させて相互作用を弱めることにより、接近した粒子間に非常に強い力が働くのを避ける方法である。この方法では、シミュレーション中にポテンシャル関数の関数形を時々刻々と変化させ、徐々に元の関数形に戻していく。このとき、所定値Ｒｉを用いて各粒子間距離を判定し、すべての粒子間距離がＲｉ以上になれば、ポテンシャル関数を元の関数形に戻す。Ｒｉの値は、あらかじめパラメタとして入力しておく。
（３）セルサイズ変更法
質量密度の実験データから計算されるサイズより大きなセルを設定してから、シミュレーションを開始する方法である。このようなセルを用いれば、各粒子をランダムに配置しても、接近した粒子対が生成される確率は小さくなる。このとき、所定値Ｒｄを用いて各粒子間距離を判定し、Ｒｄ未満の距離の粒子対があれば、ランダム配置をやりなおす。このような試行を繰り返して、すべての粒子間距離がＲｄ以上になれば、通常のシミュレーションを開始する。
【０００７】
【発明が解決しようとする課題】
しかしながら、上述した従来のシミュレーション方法には、次のような問題がある。
（１）熱排出法
与えられた初期配置が、たまたま、数値積分アルゴリズムが破綻を引き起こさない程度の配置であれば機能するが、そのような保証はなく、他の方法と組み合わせる必要がある。
【０００８】
また、熱の排出のタイミングが適切でないと、うまく機能しない。タイミングが遅いと、数値積分アルゴリズムが破綻してしまうか、そうでなくとも、数値誤差が累積し、破綻の可能性が生ずる。逆に、熱の排出を頻繁に行いすぎると、各粒子の速度が遅くなるので、系の時間変化が緩慢になり、計算コストがかかってしまう。
（２）ポテンシャル緩和法
ポテンシャル関数を変形させる際のパラメタ（変形のタイミング、変形後の関数形を定めるパラメタ）を設定するのが簡単ではなく、ノウハウを要する。経験的に定められた設定で間に合うこともあるが、対象とする系によってはそれでは機能せず、パラメタ設定のトライ＆エラーが必要になることが多い。
【０００９】
また、それらのパラメタを効果的に設定するには、原子の種類のペア毎に設定する必要があるが、複雑な系になると、この操作は煩雑になる。
また、粒子間距離の判定条件の設定も適切であることが望まれる。充分に平衡な状態になる前に、元の関数形に戻すｇｏ判定を下してしまうと、余計な計算コストが生ずる場合がある。
【００１０】
例えば、アルカン分子の集合に対して、熱排出法とポテンシャル緩和法を組み合わせて適用した結果、図１２に示すようなシミュレーション結果が得られた。図１２において、横軸は系の物理的な変化を記述する時間（計算時間とは異なる）を表し、Ｕはポテンシャル関数（内部エネルギー）の値（ｇＡ2 ｆｓ-2）を表し、Ｖはセルの体積（Ａ3 ）を表し、Ｔは温度（Ｋ）を表す。
【００１１】
この例では、０．６×１０-1ｐｓ付近でｇｏ判定が下されているが、そのタイミングが適切でなかったため、ｇｏ判定の後に粒子間に大きな力が働き、セルの体積Ｖが膨張している。このため、体積Ｖが実験データに基づく値に近づくまでに、かなりの計算コストを要した。
【００１２】
さらに、相互作用が２体力でない場合は、ポテンシャル関数の変形操作をフォーミュレートしづらいという問題もある。
（３）セルサイズ変更法
セルサイズの大きさが充分でないと、初期配置生成のためのランダム配置の試行の繰り返し回数が多くなり、煩雑である。しかし、セルサイズを大きくとりすぎると、シミュレーションによりそれが元のサイズに戻るのに時間がかかる。高分子の集合のように、相互作用が複雑で粒子数が多い系では、そのコストは大きい。セルサイズから計算される質量密度が実験データに近づくまでに、スーパコンピュータで数時間の計算を要する場合もある。
【００１３】
例えば、２０個の液晶分子の集合に対してセルサイズ変更法を適用した結果、図１３に示すようなシミュレーション結果が得られた。図１３において、横軸の時間とＵ、Ｖ、Ｔについては、図１２と同様である。この例では、体積Ｖが収束するまでに長時間の計算を要した。
【００１４】
このように、（１）〜（３）の方法は、使い易さ（汎用性）または計算コストの点で、必ずしも満足できる方法であるとは言い難い。そこで、莫大な数値によるアルゴリズムの破綻の心配がなく、使い易く、汎用性に富み、また、計算コストが低いシミュレーション方法が望まれる。
【００１５】
本発明の課題は、粒子系のシミュレーション処理を改良するような粒子集団の初期配置を生成する生成装置およびその方法を提供することである。
【００１６】
【課題を解決するための手段】
図１は、本発明の生成装置の原理図である。図１の生成装置は、計算手段１と出力手段２を備え、粒子系のシミュレーションにおける粒子集団の配置を生成する。計算手段１は、任意次元の空間内の粒子集団の平衡配置を計算し、出力手段２は、得られた平衡配置をシミュレーションにおける初期配置として出力する。
【００１７】
粒子集団の平衡配置は、各粒子に働く相互作用が充分に小さいと考えられるような配置に対応し、例えば、粒子系の相互作用を記述するポテンシャル関数または力関数から求めることができる。計算手段１は、与えられた次元の空間内における粒子集団の平衡配置を計算し、出力手段２は、その平衡配置を初期配置として、シミュレーションプログラム等に出力する。
【００１８】
このような生成装置においては、平衡配置を生成する過程で、例えば、一定量で配置を変位させるため、従来のように粒子の加速度・速度をもとにした変位操作とは異なる。したがって、粒子の加速度・速度が原因となって発生する莫大な数値によるアルゴリズムの破綻の心配はない。
【００１９】
したがって、シミュレーションにおいて、従来の熱排出法、ポテンシャル緩和法、およびセルサイズ変更法を用いる必要がなく、これらの方法に必要であった計算コストを削減することができる。具体的には、力学と熱の計算に必要な速度計算、熱排出法におけるタイミング判定と排出処理計算、ポテンシャル緩和法におけるパラメタ調整のための計算、およびセルサイズ変更法におけるセルサイズの調整のための計算が不要になる。
【００２０】
例えば、図１の計算手段１と出力手段２は、後述する図１０のＣＰＵ（中央処理装置）２１とメモリ２２に対応する。
また、このような生成装置を含むシミュレーション装置は、計算手段１、出力手段２のほかに、得られた平衡配置を初期配置として用いて粒子系のシミュレーションを行うシミュレーション手段を備える。この場合、出力手段２は、図１０の出力装置２４に対応し、シミュレーションの結果を出力する。また、シミュレーション手段は、図１０のＣＰＵ２１とメモリ２２に対応する。
【００２１】
このように、本発明のポイントは、粒子集団の平衡配置をシミュレーションの初期配置として用いることである。
本発明の別の生成装置は、計算手段１、出力手段２、および記憶手段を備え、分子系のシミュレーションにおける分子集団の配置を生成する。記憶手段は、任意次元の空間内における分子集団を構成する複数の分子の初期座標と、微小変位幅を記憶する。計算手段１は、各分子の初期座標を位置座標として、その位置座標を正方向および負方向に微小変位幅だけ変位させて、複数の分子の位置座標により記述され、分子系の相互作用を表すポテンシャル関数の値を計算し、正方向の変位後のポテンシャル関数値、負方向の変位後のポテンシャル関数値、および変位前のポテンシャル関数値を比較し、最も小さい値に対応する位置座標を新規の位置座標として記憶手段に格納し、その新規の位置座標を用いて微小変位幅の変位操作とポテンシャル関数の値の計算を繰り返すことで、複数の分子の位置座標からなる分子集団の平衡配置を計算する。出力手段２は、得られた平衡配置をシミュレーションにおける初期配置として出力する。
【００２２】
【発明の実施の形態】
以下、図面を参照しながら、本発明の実施の形態を詳細に説明する。
本実施形態においては、与えられた粒子をランダムに配置するのではなく、粒子間の相互作用を考慮して粒子集団の平衡配置を求め、得られた平衡配置をシミュレーションにおける初期配置として採用する。
【００２３】
まず、粒子間の相互作用がポテンシャル関数により与えられている場合を考える。この場合、Ｎ個の粒子からなる粒子集団の平衡配置は、系のポテンシャル関数Ｕが極小値をとる状態として表される。ここで、Ｕは、Ｒn の部分集合Ｄを定義域とし実数値をとる連続関数として、次式により定義される。
Ｕ：Ｒn ⊃Ｄ→Ｒ，
ｑ＝（ｑ1 ，．．．，ｑn ）→Ｕ（ｑ） ［１］
ここで、ｑはｎ次元ベクトルであり、実数ｑ1 ，．．．，ｑn はｎ個の自由度を表す。また、粒子が配置される空間をｍ次元空間とすると、ｎ＝Ｎｍである。このとき、各粒子の位置ベクトルｒ1 ，．．．，ｒN は、次式により与えられる。
ｒ1 ≡（ｑ1 ，．．．，ｑm ），
ｒ2 ≡（ｑm+1 ，．．．，ｑ2m），
．．．．．．．．．．．．．．．．
ｒN ≡（ｑ(N-1)m+1，．．．，ｑNm） ［２］
また、任意のｉ（ｉ＝１，．．．，ｎ）について、Ｕは次式のように分解することができる。
【００２４】
【数１】
【００２５】
ただし、ｑ(i) は次式で定義されるｎ次元ベクトルである。
ｑ(i) ≡（ｑ1 ，．．．，ｑi ＋Δｑi ，．．．，ｑn ） ［５］
［４］式は、ｑi をΔｑi だけ変化させたときに 外１ （以下では、Ｕ(i)
【００２６】
【外１】
【００２７】
バーと記す）が変化しないことを表している。言い換えれば、［３］式のＵ(i) バーは、Ｕのｑi に依存しない部分に対応し、Ｕ(i) は、Ｕのｑi に依存する部分に対応する。
【００２８】
Ｕが複数の項の和で表され、ｑi に依存しない項を含んでいる場合は、その項をＵ(i) バーとおいて、Ｕ(i) ＝Ｕ−Ｕ(i) バーと定義すればよい。それ以外の場合は、Ｕ(i) ≡Ｕ、Ｕ(i) バー≡０とおけばよい。例えば、よく用いられる２体力のポテンシャル関数は、次式に示すように２つの粒子間のポテンシャルφの総和で表される。
【００２９】
【数２】
【００３０】
この場合、Ｕ(i) は次式で与えられる。
【００３１】
【数３】
【００３２】
ところで、一般に、ＵとＵ(i) の間には次式のような関係が成り立つ。
Ｕ（ｑ(i) ）＜Ｕ（ｑ）⇔Ｕ(i) （ｑ(i) ）＜Ｕ(i) （ｑ） ［８］
したがって、Ｕを小さくするためには、Ｕ(i) を小さくするようにｑを変化させればよい。このような操作を繰り返すことにより、Ｕの極小値に対応する各粒子の位置を容易に求めることができる。
【００３３】
次に、図２から図４までを参照しながら、このような平衡配置生成処理のアルゴリズムを説明する。
図２は、入力データを示している。この入力データにおいて、自由度１１は、与えられた問題を記述する状態変数ｑi の数ｎを表し、粒子の初期座標１２は、ｍ次元空間におけるＮ個の粒子の初期位置を表す。
【００３４】
また、降下法の条件１３には、ステップ数、微小変位の幅、終了条件、出力指定パラメタ等が含まれる。ステップ数は、降下法の反復回数を表し、微小変位の幅は、変数ｑi の変位量を表し、終了条件は、降下法の終了条件を表し、出力指定パラメタは、出力データの出力間隔等を指定するパラメタである。終了条件としては、例えば、次のようなものが用いられる。
（ａ）計算時間または処理ステップ数があらかじめ決められた値に到達したとき、反復処理を終了する。
（ｂ）ポテンシャルＵ（ｑ）の変化があらかじめ決められた判定値を下回ったとき、反復処理を終了する。
【００３５】
また、境界条件１４は、ポテンシャル関数の定義域Ｄに関する境界条件を表す。例えば、トーラスが境界条件として指定されると、領域Ｄがトーラス状に連続しているものとみなして、計算が行われる。
【００３６】
また、カットオフ長１５は、相互作用を簡略化するためのパラメタである。計算の効率化のため、オプション機能として、カットオフ長１５に対応する距離以上離れた遠方のポテンシャルを無視することができる。
【００３７】
図３および図４は、本実施形態の生成装置による平衡配置生成処理のフローチャートである。生成装置は、まず、図２に示したような入力データを設定する（図３のステップＳ１）。このとき、入力データの設定方法をユーザに問い合わせ（ステップＳ２）、ユーザの選択に応じて、所定の方法でデータを自動生成するか（ステップＳ３）、または所定の外部ファイルからデータを読み込む（ステップＳ４）。
【００３８】
ここでは、自由度１１としてｎが入力され、降下法の条件１３として、反復回数、ｎ個の微小変位サイズΔ１，Δ２，．．．，Δｎ、収束判定パラメタε（＞０）が入力される。Δ１，Δ２，．．．，Δｎとしては、一般に、それぞれ異なる値が用いられるが、同一の値を用いてもよい。
【００３９】
次に、自由度の制御変数ｉを１とおき（ステップＳ５）、ポテンシャル関数の境界条件処理を行い（ステップＳ６）、ステップＳ７、Ｓ８、およびＳ９の処理を行う。ステップＳ７では、［３］式のＵ(i) （ｑ1 ，．．．，ｑn ）の値を計算し、得られた値をＵ（ｉ）０とおく。
【００４０】
ステップＳ８では、ｉ番目の自由度に対応するｑi の値を正の方向に微小変位させてｑｉ＋とする。そして、ｑｉ＋を用いてＵ（ｑ1 ，．．．，ｑn ）の値を計算し、得られた値をＵ（ｉ）＋とおく。ステップＳ９では、ｑi の値を負の方向に微小変位させてｑｉ−とする。そして、ｑｉ−を用いてＵ（ｑ1 ，．．．，ｑn ）の値を計算し、得られた値をＵ（ｉ）−とおく。ｑｉ＋およびｑｉ−は、Δｉを用いた任意の計算により生成される。最も簡単には、ｑｉ＋＝ｑｉ＋Δｉ、ｑｉ−＝ｑｉ−Δｉとすればよい。
【００４１】
次に、Ｕ（ｉ）０、Ｕ（ｉ）＋、およびＵ（ｉ）−の値を比較し（図４のステップＳ１０）、Ｕ（ｉ）０が最も小さければ、ｑi を新規のｑi の値として採択する。この場合、ｑi の値は変化しないことになる。また、Ｕ（ｉ）＋が最も小さければ、ｑｉ＋を新規のｑi の値として採択し、Ｕ（ｉ）−が最も小さければ、ｑｉ−を新規のｑi の値として採択する。
【００４２】
次に、ｉ＝ｉ＋１とおき（ステップＳ１１）、ｉをｎと比較する（ステップＳ１２）。ｉがｎを超えていなければ、次の自由度について、図３のステップＳ６以降の処理を繰り返す。このような処理を、ｉ＝１，２，．．．，ｎについて繰り返すことにより、すべての自由度に関する微小変位が行われる。
【００４３】
ステップＳ１２においてｉがｎを超えれば、次に、所定の終了条件が成立したかどうかを判定する（ステップＳ１３）。終了条件が成立しなければ、図３のステップＳ５以降の降下法の処理を繰り返す。そして、終了条件が成立すれば、降下法の処理を終了し（ステップＳ１４）、そのときのｑ1 ，．．．，ｑn 、Ｕ(1) （ｑ1 ，．．．，ｑn ），．．．，Ｕ(n) （ｑ1 ，．．．，ｑn ）、およびＵ（ｑ1 ，．．．，ｑn ）を出力データとして出力して（ステップＳ１５）、処理を終了する。
【００４４】
ステップＳ１３における終了条件としては、例えば、ポテンシャル関数の変化に関する適当な収束条件が満たされるか、または反復回数が入力された回数に達したとき、降下法を終了するという条件が考えられる。また、収束条件としては、例えば、収束判定パラメタεを用いた次式のような条件が用いられる。
｜Ｕ(i) （ｑ1 ，．．．，ｑn ）−Ｕ（ｉ）０｜／｜Ｕ（ｉ）０｜＜ε
ｆｏｒ ｉ＝１，．．．，ｎ ［９］
また、ステップＳ１５においては、得られたｑ1 ，．．．，ｑn を用いて、Ｎ個の粒子の位置がディスプレイ画面上で可視化され、Ｎ個の粒子の平衡配置が表示される。平衡配置を表すデータｑ1 ，．．．，ｑn は、ＭＤ等のシミュレーションにおける初期配置データとして用いられる。
【００４５】
ここでは、降下法の各繰り返しにおいて同じΔｉを用いているが、Δｉを繰り返し毎に変更してもよい。また、収束条件としては、［９］式以外にも任意の条件を用いることができる。
【００４６】
このようなアルゴリズムによれば、降下法の反復回数や収束条件のようなパラメタは簡単に設定することができ、従来のポテンシャル緩和法のような煩雑なパラメタ設定は必要がない。したがって、特別なノウハウを要することなく、容易に平衡配置を生成することが可能になる。
【００４７】
次に、粒子間の相互作用が、次式のような力関数（連続関数）により与えられている場合を考える。
Ｆ：Ｒn ⊃Ｄ→Ｒn ，
ｑ＝（ｑ1 ，．．．，ｑn ）→Ｆ（ｑ） ［１０］
ここで、Ｆ（ｑ）はｎ次元の力ベクトルを表す。この場合、Ｎ個の粒子の平衡配置は、Ｆ（ｑ）＝０となるｎ次元ベクトルｑにより与えられる。そこで、連続狭義単調増加関数ｇを用いて、関数Ｕを次式により定義する。
【００４８】
【数４】
【００４９】
そして、このＵを［３］式のポテンシャル関数Ｕと同様に扱い、その極小値を実現する状態を求める。得られた状態の各粒子の位置が平衡配置に対応する。一般に、相互作用が力関数により与えられる場合は、相互作用がポテンシャル関数により与えられる場合を含んでおり、より多くの問題に適用することができると考えられる。
【００５０】
このように、本実施形態によれば、［６］式のような２体力以外の相互作用（Ｎ体力の相互作用、力関数で与えられる相互作用等）でも、汎用的に扱うことができる。
【００５１】
次に、相互作用がポテンシャル関数により与えられる場合の具体例について説明する。２次元空間（平面）における１００個のアルゴン原子のシミュレーションにおいて、ポテンシャル関数が［６］式により与えられ、［６］式のφが次式により与えられたとする。
φ（ｒ）＝ε（（σ／ｒ）12−（σ／ｒ）6 ） ［１２］
ここで、εおよびσは正の実数であり、パラメタとして設定される。この場合、ｍ＝２、Ｎ＝１００、ｎ＝２００として、上述のアルゴリズムにより平衡配置が生成される。
【００５２】
図５は、ＸＹ平面上において、１００個のアルゴン原子を正方格子の位置から微小にずらして、疑似ランダムに配置した結果を示している。この状態からスタートして、図３および図４に示したアルゴリズムにより平衡配置を計算すると、図６のような配置が得られた。図６においては、ところどころ格子欠損のようなものが生じているが、実験により得られる面心立方（face centered cubic ，ＦＣＣ）構造に近い配置が得られている。
【００５３】
また、図５の配置よりさらにランダムにアルゴン原子を配置すると、図７のような配置が得られる。この場合、非常に接近した原子対が生成されているので、これを初期配置として動力学のシミュレーションを行うと、すぐに局所的な高温が発生し、オーバーフローが生じてしまう。
【００５４】
そこで、この状態からスタートして、図３および図４に示したアルゴリズムにより平衡配置を計算すると、図８のような配置が得られた。図８においては、図６と同様に、実験により得られる構造に近い配置が得られている。また、このとき、ポテンシャル関数の値は、図９に示すように変化した。図９においては、降下法の計算を繰り返すに従って、ポテンシャル関数の値が単調に減少していることが分かる。
【００５５】
以上説明した実施形態においては、比較的簡単な降下法を用いて平衡配置を求めているが、代わりに、シミュレーティド・アニーリング、遺伝子アルゴリズム等の他の探索法を用いてもよい。また、先願の「最適化問題を解決する処理装置および方法」（特願平１１−１６５００）に開示されているように、最適化問題を表す形状モデルを変形しながらコスト値の小さな形状を探索するアルゴリズムを用いて平衡配置を求めてもよい。この場合、各粒子が形状モデルを構成する変形要素として用いられる。
【００５６】
また、シミュレーション対象の粒子集団は原子や分子に限られず、本発明は、あらかじめ相互作用が定義された任意の粒子集団のシミュレーションに適用することができる。
【００５７】
さらに、粒子を配置する空間は２次元空間に限られず、本発明は、任意次元の空間における粒子集団のシミュレーションに適用することができる。通常は、３次元空間内においてシミュレーションが行われることが多い。
【００５８】
ところで、上述した生成装置は、図１０に示すような情報処理装置（コンピュータ）を用いて構成することができる。図１０の情報処理装置は、ＣＰＵ（中央処理装置）２１、メモリ２２、入力装置２３、出力装置２４、外部記憶装置２５、媒体駆動装置２６、およびネットワーク接続装置２７を備え、それらはバス２８により互いに接続されている。
【００５９】
メモリ２２は、例えば、ＲＯＭ（read only memory）、ＲＡＭ（random access memory）等を含み、処理に用いられるプログラムとデータを格納する。ＣＰＵ２１は、メモリ２２を利用してプログラムを実行することにより、必要な処理を行う。
【００６０】
入力装置２３は、例えば、キーボード、ポインティングデバイス、タッチパネル等であり、ユーザからの指示や情報の入力に用いられる。出力装置２４は、例えば、ディスプレイ、プリンタ、スピーカ等であり、ユーザへの問い合わせや処理結果の出力に用いられる。
【００６１】
外部記憶装置２５は、例えば、磁気ディスク装置、光ディスク装置、光磁気ディスク（magneto-optical disk）装置等である。情報処理装置は、この外部記憶装置２５に、上述のプログラムとデータを保存しておき、必要に応じて、それらをメモリ２２にロードして使用することができる。また、外部記憶装置２５は、入力データのファイルを格納するデータベースとしても用いられる。
【００６２】
媒体駆動装置２６は、可搬記録媒体２９を駆動し、その記録内容にアクセスする。可搬記録媒体２９としては、メモリカード、フロッピーディスク、ＣＤ−ＲＯＭ（compact disk read only memory ）、光ディスク、光磁気ディスク等、任意のコンピュータ読み取り可能な記録媒体が用いられる。ユーザは、この可搬記録媒体２９に上述のプログラムとデータを格納しておき、必要に応じて、それらをメモリ２２にロードして使用することができる。
【００６３】
ネットワーク接続装置２７は、ＬＡＮ（local area network）等の任意のネットワーク（回線）を介して外部の装置と通信し、通信に伴うデータ変換を行う。情報処理装置は、必要に応じて、ネットワーク接続装置２７を介して上述のプログラムとデータを外部の装置から受け取り、それらをメモリ２２にロードして使用することができる。
【００６４】
また、平衡配置を初期配置として粒子系のシミュレーションを行うシミュレーション装置は、上述した生成装置を含み、図１０と同様の情報処理装置を用いて構成することができる。この場合、ＣＰＵ２１は、メモリ２２を利用してプログラムを実行することにより、ＭＤ等のシミュレーション処理を行う。
【００６５】
図１１は、図１０の情報処理装置にプログラムとデータを供給することのできるコンピュータ読み取り可能な記録媒体を示している。可搬記録媒体２９や外部のデータベース３０に保存されたプログラムとデータは、メモリ２２にロードされる。そして、ＣＰＵ２１は、そのデータを用いてそのプログラムを実行し、必要な処理を行う。
【００６６】
【発明の効果】
本発明によれば、従来のように粒子の速度を考える必要がなく、莫大な数値によるアルゴリズムの破綻の心配がない。したがって、従来の熱排出法、ポテンシャル緩和法、およびセルサイズ変更法を用いる必要がなく、これらの方法に必要であった計算コストが削減されるとともに、シミュレーションの操作性が向上する。また、様々な相互作用を汎用的に扱うことが可能になる。
【図面の簡単な説明】
【図１】本発明の生成装置の原理図である。
【図２】入力データを示す図である。
【図３】平衡配置生成処理のフローチャート（その１）である。
【図４】平衡配置生成処理のフローチャート（その２）である。
【図５】疑似ランダム配置を示す図である。
【図６】第１の平衡配置を示す図である。
【図７】ランダム配置を示す図である。
【図８】第２の平衡配置を示す図である。
【図９】ポテンシャルの変化を示す図である。
【図１０】情報処理装置の構成図である。
【図１１】記録媒体を示す図である。
【図１２】従来の第１のシミュレーション結果を示す図である。
【図１３】従来の第２のシミュレーション結果を示す図である。
【符号の説明】
１ 計算手段
２ 出力手段
１１ 自由度
１２ 粒子の初期座標
１３ 降下法の条件
１４ 境界条件
１５ カットオフ長
２１ ＣＰＵ
２２ メモリ
２３ 入力装置
２４ 出力装置
２５ 外部記憶装置
２６ 媒体駆動装置
２７ ネットワーク接続装置
２８ バス
２９ 可搬記録媒体
３０ データベース
【発明の属する技術分野】
本発明は、粒子系のシミュレーションにおいて、粒子集団の配置を表す情報を生成する生成装置およびその方法に関する。
【０００２】
【従来の技術】
無機物または有機物で構成される材料の性質の解析・予測を行う方法の１つとして、分子動力学（molecular dynamics，ＭＤ）計算に基づくシミュレーションがある。この計算を行うためには、材料を構成する原子・分子の初期配置を生成する必要がある。ＭＤ計算においては、与えられた初期配置を初期値として、ニュートン（Newton）方程式で表されるような微分方程式を解き、計算結果を用いて材料を表す系の性質を調べ、性質の予測にまでつなげることが試みられている。
【０００３】
しかしながら、１分子の安定状態と考えられる分子構造や、常温・常圧における質量密度の実験データは得られていても、分子集団の配置（各分子を構成するすべての原子の位置座標）の実験データは得られていないことが多い。したがって、分子系のＭＤ計算を行う際、実験的に得られた質量密度を実現する分子集団の初期配置をどのようにして生成したらよいか不明である。
【０００４】
粒子集団の初期配置を生成する簡単な方法として、質量密度の実験データから一定のセルサイズと１セルあたりの粒子数を求め、セル内に粒子をランダムに配置する方法が考えられる。しかし、粒子を全くランダムに配置すると、しばしば、接近した粒子対が生成されてしまう。これらの接近した粒子間には、通常、非常に強い力が働くことが知られている。
【０００５】
このような接近した粒子対が１つでも存在すると、相互作用のためにそれらの粒子の速度が速くなり、その近傍で系が局所的に非常に高温になる。このため、微分方程式を解く際に莫大な数値が発生し、数値積分アルゴリズムが、しばしば、破綻を引き起こす。
【０００６】
このようなシミュレーションの破綻を回避する従来の方法として、次の３つの方法が挙げられる。通常は、これらの方法を組み合わせて用いることが多い。
（１）熱排出法
シミュレーション中に発生した系の熱を適当なタイミングで排出する方法である。具体的には、各粒子の速度を低くするか、または完全に零にしてしまう。
（２）ポテンシャル緩和法
粒子間距離の小さいところで、ポテンシャル関数を変形させて相互作用を弱めることにより、接近した粒子間に非常に強い力が働くのを避ける方法である。この方法では、シミュレーション中にポテンシャル関数の関数形を時々刻々と変化させ、徐々に元の関数形に戻していく。このとき、所定値Ｒｉを用いて各粒子間距離を判定し、すべての粒子間距離がＲｉ以上になれば、ポテンシャル関数を元の関数形に戻す。Ｒｉの値は、あらかじめパラメタとして入力しておく。
（３）セルサイズ変更法
質量密度の実験データから計算されるサイズより大きなセルを設定してから、シミュレーションを開始する方法である。このようなセルを用いれば、各粒子をランダムに配置しても、接近した粒子対が生成される確率は小さくなる。このとき、所定値Ｒｄを用いて各粒子間距離を判定し、Ｒｄ未満の距離の粒子対があれば、ランダム配置をやりなおす。このような試行を繰り返して、すべての粒子間距離がＲｄ以上になれば、通常のシミュレーションを開始する。
【０００７】
【発明が解決しようとする課題】
しかしながら、上述した従来のシミュレーション方法には、次のような問題がある。
（１）熱排出法
与えられた初期配置が、たまたま、数値積分アルゴリズムが破綻を引き起こさない程度の配置であれば機能するが、そのような保証はなく、他の方法と組み合わせる必要がある。
【０００８】
また、熱の排出のタイミングが適切でないと、うまく機能しない。タイミングが遅いと、数値積分アルゴリズムが破綻してしまうか、そうでなくとも、数値誤差が累積し、破綻の可能性が生ずる。逆に、熱の排出を頻繁に行いすぎると、各粒子の速度が遅くなるので、系の時間変化が緩慢になり、計算コストがかかってしまう。
（２）ポテンシャル緩和法
ポテンシャル関数を変形させる際のパラメタ（変形のタイミング、変形後の関数形を定めるパラメタ）を設定するのが簡単ではなく、ノウハウを要する。経験的に定められた設定で間に合うこともあるが、対象とする系によってはそれでは機能せず、パラメタ設定のトライ＆エラーが必要になることが多い。
【０００９】
また、それらのパラメタを効果的に設定するには、原子の種類のペア毎に設定する必要があるが、複雑な系になると、この操作は煩雑になる。
また、粒子間距離の判定条件の設定も適切であることが望まれる。充分に平衡な状態になる前に、元の関数形に戻すｇｏ判定を下してしまうと、余計な計算コストが生ずる場合がある。
【００１０】
例えば、アルカン分子の集合に対して、熱排出法とポテンシャル緩和法を組み合わせて適用した結果、図１２に示すようなシミュレーション結果が得られた。図１２において、横軸は系の物理的な変化を記述する時間（計算時間とは異なる）を表し、Ｕはポテンシャル関数（内部エネルギー）の値（ｇＡ2 ｆｓ-2）を表し、Ｖはセルの体積（Ａ3 ）を表し、Ｔは温度（Ｋ）を表す。
【００１１】
この例では、０．６×１０-1ｐｓ付近でｇｏ判定が下されているが、そのタイミングが適切でなかったため、ｇｏ判定の後に粒子間に大きな力が働き、セルの体積Ｖが膨張している。このため、体積Ｖが実験データに基づく値に近づくまでに、かなりの計算コストを要した。
【００１２】
さらに、相互作用が２体力でない場合は、ポテンシャル関数の変形操作をフォーミュレートしづらいという問題もある。
（３）セルサイズ変更法
セルサイズの大きさが充分でないと、初期配置生成のためのランダム配置の試行の繰り返し回数が多くなり、煩雑である。しかし、セルサイズを大きくとりすぎると、シミュレーションによりそれが元のサイズに戻るのに時間がかかる。高分子の集合のように、相互作用が複雑で粒子数が多い系では、そのコストは大きい。セルサイズから計算される質量密度が実験データに近づくまでに、スーパコンピュータで数時間の計算を要する場合もある。
【００１３】
例えば、２０個の液晶分子の集合に対してセルサイズ変更法を適用した結果、図１３に示すようなシミュレーション結果が得られた。図１３において、横軸の時間とＵ、Ｖ、Ｔについては、図１２と同様である。この例では、体積Ｖが収束するまでに長時間の計算を要した。
【００１４】
このように、（１）〜（３）の方法は、使い易さ（汎用性）または計算コストの点で、必ずしも満足できる方法であるとは言い難い。そこで、莫大な数値によるアルゴリズムの破綻の心配がなく、使い易く、汎用性に富み、また、計算コストが低いシミュレーション方法が望まれる。
【００１５】
本発明の課題は、粒子系のシミュレーション処理を改良するような粒子集団の初期配置を生成する生成装置およびその方法を提供することである。
【００１６】
【課題を解決するための手段】
図１は、本発明の生成装置の原理図である。図１の生成装置は、計算手段１と出力手段２を備え、粒子系のシミュレーションにおける粒子集団の配置を生成する。計算手段１は、任意次元の空間内の粒子集団の平衡配置を計算し、出力手段２は、得られた平衡配置をシミュレーションにおける初期配置として出力する。
【００１７】
粒子集団の平衡配置は、各粒子に働く相互作用が充分に小さいと考えられるような配置に対応し、例えば、粒子系の相互作用を記述するポテンシャル関数または力関数から求めることができる。計算手段１は、与えられた次元の空間内における粒子集団の平衡配置を計算し、出力手段２は、その平衡配置を初期配置として、シミュレーションプログラム等に出力する。
【００１８】
このような生成装置においては、平衡配置を生成する過程で、例えば、一定量で配置を変位させるため、従来のように粒子の加速度・速度をもとにした変位操作とは異なる。したがって、粒子の加速度・速度が原因となって発生する莫大な数値によるアルゴリズムの破綻の心配はない。
【００１９】
したがって、シミュレーションにおいて、従来の熱排出法、ポテンシャル緩和法、およびセルサイズ変更法を用いる必要がなく、これらの方法に必要であった計算コストを削減することができる。具体的には、力学と熱の計算に必要な速度計算、熱排出法におけるタイミング判定と排出処理計算、ポテンシャル緩和法におけるパラメタ調整のための計算、およびセルサイズ変更法におけるセルサイズの調整のための計算が不要になる。
【００２０】
例えば、図１の計算手段１と出力手段２は、後述する図１０のＣＰＵ（中央処理装置）２１とメモリ２２に対応する。
また、このような生成装置を含むシミュレーション装置は、計算手段１、出力手段２のほかに、得られた平衡配置を初期配置として用いて粒子系のシミュレーションを行うシミュレーション手段を備える。この場合、出力手段２は、図１０の出力装置２４に対応し、シミュレーションの結果を出力する。また、シミュレーション手段は、図１０のＣＰＵ２１とメモリ２２に対応する。
【００２１】
このように、本発明のポイントは、粒子集団の平衡配置をシミュレーションの初期配置として用いることである。
本発明の別の生成装置は、計算手段１、出力手段２、および記憶手段を備え、分子系のシミュレーションにおける分子集団の配置を生成する。記憶手段は、任意次元の空間内における分子集団を構成する複数の分子の初期座標と、微小変位幅を記憶する。計算手段１は、各分子の初期座標を位置座標として、その位置座標を正方向および負方向に微小変位幅だけ変位させて、複数の分子の位置座標により記述され、分子系の相互作用を表すポテンシャル関数の値を計算し、正方向の変位後のポテンシャル関数値、負方向の変位後のポテンシャル関数値、および変位前のポテンシャル関数値を比較し、最も小さい値に対応する位置座標を新規の位置座標として記憶手段に格納し、その新規の位置座標を用いて微小変位幅の変位操作とポテンシャル関数の値の計算を繰り返すことで、複数の分子の位置座標からなる分子集団の平衡配置を計算する。出力手段２は、得られた平衡配置をシミュレーションにおける初期配置として出力する。
【００２２】
【発明の実施の形態】
以下、図面を参照しながら、本発明の実施の形態を詳細に説明する。
本実施形態においては、与えられた粒子をランダムに配置するのではなく、粒子間の相互作用を考慮して粒子集団の平衡配置を求め、得られた平衡配置をシミュレーションにおける初期配置として採用する。
【００２３】
まず、粒子間の相互作用がポテンシャル関数により与えられている場合を考える。この場合、Ｎ個の粒子からなる粒子集団の平衡配置は、系のポテンシャル関数Ｕが極小値をとる状態として表される。ここで、Ｕは、Ｒn の部分集合Ｄを定義域とし実数値をとる連続関数として、次式により定義される。
Ｕ：Ｒn ⊃Ｄ→Ｒ，
ｑ＝（ｑ1 ，．．．，ｑn ）→Ｕ（ｑ） ［１］
ここで、ｑはｎ次元ベクトルであり、実数ｑ1 ，．．．，ｑn はｎ個の自由度を表す。また、粒子が配置される空間をｍ次元空間とすると、ｎ＝Ｎｍである。このとき、各粒子の位置ベクトルｒ1 ，．．．，ｒN は、次式により与えられる。
ｒ1 ≡（ｑ1 ，．．．，ｑm ），
ｒ2 ≡（ｑm+1 ，．．．，ｑ2m），
．．．．．．．．．．．．．．．．
ｒN ≡（ｑ(N-1)m+1，．．．，ｑNm） ［２］
また、任意のｉ（ｉ＝１，．．．，ｎ）について、Ｕは次式のように分解することができる。
【００２４】
【数１】
ただし、ｑ(i) は次式で定義されるｎ次元ベクトルである。
ｑ(i) ≡（ｑ1 ，．．．，ｑi ＋Δｑi ，．．．，ｑn ） ［５］
［４］式は、ｑi をΔｑi だけ変化させたときに 外１ （以下では、Ｕ(i)
【００２６】
【外１】
バーと記す）が変化しないことを表している。言い換えれば、［３］式のＵ(i) バーは、Ｕのｑi に依存しない部分に対応し、Ｕ(i) は、Ｕのｑi に依存する部分に対応する。
【００２８】
Ｕが複数の項の和で表され、ｑi に依存しない項を含んでいる場合は、その項をＵ(i) バーとおいて、Ｕ(i) ＝Ｕ−Ｕ(i) バーと定義すればよい。それ以外の場合は、Ｕ(i) ≡Ｕ、Ｕ(i) バー≡０とおけばよい。例えば、よく用いられる２体力のポテンシャル関数は、次式に示すように２つの粒子間のポテンシャルφの総和で表される。
【００２９】
【数２】
この場合、Ｕ(i) は次式で与えられる。
【００３１】
【数３】
ところで、一般に、ＵとＵ(i) の間には次式のような関係が成り立つ。
Ｕ（ｑ(i) ）＜Ｕ（ｑ）⇔Ｕ(i) （ｑ(i) ）＜Ｕ(i) （ｑ） ［８］
したがって、Ｕを小さくするためには、Ｕ(i) を小さくするようにｑを変化させればよい。このような操作を繰り返すことにより、Ｕの極小値に対応する各粒子の位置を容易に求めることができる。
【００３３】
次に、図２から図４までを参照しながら、このような平衡配置生成処理のアルゴリズムを説明する。
図２は、入力データを示している。この入力データにおいて、自由度１１は、与えられた問題を記述する状態変数ｑi の数ｎを表し、粒子の初期座標１２は、ｍ次元空間におけるＮ個の粒子の初期位置を表す。
【００３４】
また、降下法の条件１３には、ステップ数、微小変位の幅、終了条件、出力指定パラメタ等が含まれる。ステップ数は、降下法の反復回数を表し、微小変位の幅は、変数ｑi の変位量を表し、終了条件は、降下法の終了条件を表し、出力指定パラメタは、出力データの出力間隔等を指定するパラメタである。終了条件としては、例えば、次のようなものが用いられる。
（ａ）計算時間または処理ステップ数があらかじめ決められた値に到達したとき、反復処理を終了する。
（ｂ）ポテンシャルＵ（ｑ）の変化があらかじめ決められた判定値を下回ったとき、反復処理を終了する。
【００３５】
また、境界条件１４は、ポテンシャル関数の定義域Ｄに関する境界条件を表す。例えば、トーラスが境界条件として指定されると、領域Ｄがトーラス状に連続しているものとみなして、計算が行われる。
【００３６】
また、カットオフ長１５は、相互作用を簡略化するためのパラメタである。計算の効率化のため、オプション機能として、カットオフ長１５に対応する距離以上離れた遠方のポテンシャルを無視することができる。
【００３７】
図３および図４は、本実施形態の生成装置による平衡配置生成処理のフローチャートである。生成装置は、まず、図２に示したような入力データを設定する（図３のステップＳ１）。このとき、入力データの設定方法をユーザに問い合わせ（ステップＳ２）、ユーザの選択に応じて、所定の方法でデータを自動生成するか（ステップＳ３）、または所定の外部ファイルからデータを読み込む（ステップＳ４）。
【００３８】
ここでは、自由度１１としてｎが入力され、降下法の条件１３として、反復回数、ｎ個の微小変位サイズΔ１，Δ２，．．．，Δｎ、収束判定パラメタε（＞０）が入力される。Δ１，Δ２，．．．，Δｎとしては、一般に、それぞれ異なる値が用いられるが、同一の値を用いてもよい。
【００３９】
次に、自由度の制御変数ｉを１とおき（ステップＳ５）、ポテンシャル関数の境界条件処理を行い（ステップＳ６）、ステップＳ７、Ｓ８、およびＳ９の処理を行う。ステップＳ７では、［３］式のＵ(i) （ｑ1 ，．．．，ｑn ）の値を計算し、得られた値をＵ（ｉ）０とおく。
【００４０】
ステップＳ８では、ｉ番目の自由度に対応するｑi の値を正の方向に微小変位させてｑｉ＋とする。そして、ｑｉ＋を用いてＵ（ｑ1 ，．．．，ｑn ）の値を計算し、得られた値をＵ（ｉ）＋とおく。ステップＳ９では、ｑi の値を負の方向に微小変位させてｑｉ−とする。そして、ｑｉ−を用いてＵ（ｑ1 ，．．．，ｑn ）の値を計算し、得られた値をＵ（ｉ）−とおく。ｑｉ＋およびｑｉ−は、Δｉを用いた任意の計算により生成される。最も簡単には、ｑｉ＋＝ｑｉ＋Δｉ、ｑｉ−＝ｑｉ−Δｉとすればよい。
【００４１】
次に、Ｕ（ｉ）０、Ｕ（ｉ）＋、およびＵ（ｉ）−の値を比較し（図４のステップＳ１０）、Ｕ（ｉ）０が最も小さければ、ｑi を新規のｑi の値として採択する。この場合、ｑi の値は変化しないことになる。また、Ｕ（ｉ）＋が最も小さければ、ｑｉ＋を新規のｑi の値として採択し、Ｕ（ｉ）−が最も小さければ、ｑｉ−を新規のｑi の値として採択する。
【００４２】
次に、ｉ＝ｉ＋１とおき（ステップＳ１１）、ｉをｎと比較する（ステップＳ１２）。ｉがｎを超えていなければ、次の自由度について、図３のステップＳ６以降の処理を繰り返す。このような処理を、ｉ＝１，２，．．．，ｎについて繰り返すことにより、すべての自由度に関する微小変位が行われる。
【００４３】
ステップＳ１２においてｉがｎを超えれば、次に、所定の終了条件が成立したかどうかを判定する（ステップＳ１３）。終了条件が成立しなければ、図３のステップＳ５以降の降下法の処理を繰り返す。そして、終了条件が成立すれば、降下法の処理を終了し（ステップＳ１４）、そのときのｑ1 ，．．．，ｑn 、Ｕ(1) （ｑ1 ，．．．，ｑn ），．．．，Ｕ(n) （ｑ1 ，．．．，ｑn ）、およびＵ（ｑ1 ，．．．，ｑn ）を出力データとして出力して（ステップＳ１５）、処理を終了する。
【００４４】
ステップＳ１３における終了条件としては、例えば、ポテンシャル関数の変化に関する適当な収束条件が満たされるか、または反復回数が入力された回数に達したとき、降下法を終了するという条件が考えられる。また、収束条件としては、例えば、収束判定パラメタεを用いた次式のような条件が用いられる。
｜Ｕ(i) （ｑ1 ，．．．，ｑn ）−Ｕ（ｉ）０｜／｜Ｕ（ｉ）０｜＜ε
ｆｏｒ ｉ＝１，．．．，ｎ ［９］
また、ステップＳ１５においては、得られたｑ1 ，．．．，ｑn を用いて、Ｎ個の粒子の位置がディスプレイ画面上で可視化され、Ｎ個の粒子の平衡配置が表示される。平衡配置を表すデータｑ1 ，．．．，ｑn は、ＭＤ等のシミュレーションにおける初期配置データとして用いられる。
【００４５】
ここでは、降下法の各繰り返しにおいて同じΔｉを用いているが、Δｉを繰り返し毎に変更してもよい。また、収束条件としては、［９］式以外にも任意の条件を用いることができる。
【００４６】
このようなアルゴリズムによれば、降下法の反復回数や収束条件のようなパラメタは簡単に設定することができ、従来のポテンシャル緩和法のような煩雑なパラメタ設定は必要がない。したがって、特別なノウハウを要することなく、容易に平衡配置を生成することが可能になる。
【００４７】
次に、粒子間の相互作用が、次式のような力関数（連続関数）により与えられている場合を考える。
Ｆ：Ｒn ⊃Ｄ→Ｒn ，
ｑ＝（ｑ1 ，．．．，ｑn ）→Ｆ（ｑ） ［１０］
ここで、Ｆ（ｑ）はｎ次元の力ベクトルを表す。この場合、Ｎ個の粒子の平衡配置は、Ｆ（ｑ）＝０となるｎ次元ベクトルｑにより与えられる。そこで、連続狭義単調増加関数ｇを用いて、関数Ｕを次式により定義する。
【００４８】
【数４】
そして、このＵを［３］式のポテンシャル関数Ｕと同様に扱い、その極小値を実現する状態を求める。得られた状態の各粒子の位置が平衡配置に対応する。一般に、相互作用が力関数により与えられる場合は、相互作用がポテンシャル関数により与えられる場合を含んでおり、より多くの問題に適用することができると考えられる。
【００５０】
このように、本実施形態によれば、［６］式のような２体力以外の相互作用（Ｎ体力の相互作用、力関数で与えられる相互作用等）でも、汎用的に扱うことができる。
【００５１】
次に、相互作用がポテンシャル関数により与えられる場合の具体例について説明する。２次元空間（平面）における１００個のアルゴン原子のシミュレーションにおいて、ポテンシャル関数が［６］式により与えられ、［６］式のφが次式により与えられたとする。
φ（ｒ）＝ε（（σ／ｒ）12−（σ／ｒ）6 ） ［１２］
ここで、εおよびσは正の実数であり、パラメタとして設定される。この場合、ｍ＝２、Ｎ＝１００、ｎ＝２００として、上述のアルゴリズムにより平衡配置が生成される。
【００５２】
図５は、ＸＹ平面上において、１００個のアルゴン原子を正方格子の位置から微小にずらして、疑似ランダムに配置した結果を示している。この状態からスタートして、図３および図４に示したアルゴリズムにより平衡配置を計算すると、図６のような配置が得られた。図６においては、ところどころ格子欠損のようなものが生じているが、実験により得られる面心立方（face centered cubic ，ＦＣＣ）構造に近い配置が得られている。
【００５３】
また、図５の配置よりさらにランダムにアルゴン原子を配置すると、図７のような配置が得られる。この場合、非常に接近した原子対が生成されているので、これを初期配置として動力学のシミュレーションを行うと、すぐに局所的な高温が発生し、オーバーフローが生じてしまう。
【００５４】
そこで、この状態からスタートして、図３および図４に示したアルゴリズムにより平衡配置を計算すると、図８のような配置が得られた。図８においては、図６と同様に、実験により得られる構造に近い配置が得られている。また、このとき、ポテンシャル関数の値は、図９に示すように変化した。図９においては、降下法の計算を繰り返すに従って、ポテンシャル関数の値が単調に減少していることが分かる。
【００５５】
以上説明した実施形態においては、比較的簡単な降下法を用いて平衡配置を求めているが、代わりに、シミュレーティド・アニーリング、遺伝子アルゴリズム等の他の探索法を用いてもよい。また、先願の「最適化問題を解決する処理装置および方法」（特願平１１−１６５００）に開示されているように、最適化問題を表す形状モデルを変形しながらコスト値の小さな形状を探索するアルゴリズムを用いて平衡配置を求めてもよい。この場合、各粒子が形状モデルを構成する変形要素として用いられる。
【００５６】
また、シミュレーション対象の粒子集団は原子や分子に限られず、本発明は、あらかじめ相互作用が定義された任意の粒子集団のシミュレーションに適用することができる。
【００５７】
さらに、粒子を配置する空間は２次元空間に限られず、本発明は、任意次元の空間における粒子集団のシミュレーションに適用することができる。通常は、３次元空間内においてシミュレーションが行われることが多い。
【００５８】
ところで、上述した生成装置は、図１０に示すような情報処理装置（コンピュータ）を用いて構成することができる。図１０の情報処理装置は、ＣＰＵ（中央処理装置）２１、メモリ２２、入力装置２３、出力装置２４、外部記憶装置２５、媒体駆動装置２６、およびネットワーク接続装置２７を備え、それらはバス２８により互いに接続されている。
【００５９】
メモリ２２は、例えば、ＲＯＭ（read only memory）、ＲＡＭ（random access memory）等を含み、処理に用いられるプログラムとデータを格納する。ＣＰＵ２１は、メモリ２２を利用してプログラムを実行することにより、必要な処理を行う。
【００６０】
入力装置２３は、例えば、キーボード、ポインティングデバイス、タッチパネル等であり、ユーザからの指示や情報の入力に用いられる。出力装置２４は、例えば、ディスプレイ、プリンタ、スピーカ等であり、ユーザへの問い合わせや処理結果の出力に用いられる。
【００６１】
外部記憶装置２５は、例えば、磁気ディスク装置、光ディスク装置、光磁気ディスク（magneto-optical disk）装置等である。情報処理装置は、この外部記憶装置２５に、上述のプログラムとデータを保存しておき、必要に応じて、それらをメモリ２２にロードして使用することができる。また、外部記憶装置２５は、入力データのファイルを格納するデータベースとしても用いられる。
【００６２】
媒体駆動装置２６は、可搬記録媒体２９を駆動し、その記録内容にアクセスする。可搬記録媒体２９としては、メモリカード、フロッピーディスク、ＣＤ−ＲＯＭ（compact disk read only memory ）、光ディスク、光磁気ディスク等、任意のコンピュータ読み取り可能な記録媒体が用いられる。ユーザは、この可搬記録媒体２９に上述のプログラムとデータを格納しておき、必要に応じて、それらをメモリ２２にロードして使用することができる。
【００６３】
ネットワーク接続装置２７は、ＬＡＮ（local area network）等の任意のネットワーク（回線）を介して外部の装置と通信し、通信に伴うデータ変換を行う。情報処理装置は、必要に応じて、ネットワーク接続装置２７を介して上述のプログラムとデータを外部の装置から受け取り、それらをメモリ２２にロードして使用することができる。
【００６４】
また、平衡配置を初期配置として粒子系のシミュレーションを行うシミュレーション装置は、上述した生成装置を含み、図１０と同様の情報処理装置を用いて構成することができる。この場合、ＣＰＵ２１は、メモリ２２を利用してプログラムを実行することにより、ＭＤ等のシミュレーション処理を行う。
【００６５】
図１１は、図１０の情報処理装置にプログラムとデータを供給することのできるコンピュータ読み取り可能な記録媒体を示している。可搬記録媒体２９や外部のデータベース３０に保存されたプログラムとデータは、メモリ２２にロードされる。そして、ＣＰＵ２１は、そのデータを用いてそのプログラムを実行し、必要な処理を行う。
【００６６】
【発明の効果】
本発明によれば、従来のように粒子の速度を考える必要がなく、莫大な数値によるアルゴリズムの破綻の心配がない。したがって、従来の熱排出法、ポテンシャル緩和法、およびセルサイズ変更法を用いる必要がなく、これらの方法に必要であった計算コストが削減されるとともに、シミュレーションの操作性が向上する。また、様々な相互作用を汎用的に扱うことが可能になる。
【図面の簡単な説明】
【図１】本発明の生成装置の原理図である。
【図２】入力データを示す図である。
【図３】平衡配置生成処理のフローチャート（その１）である。
【図４】平衡配置生成処理のフローチャート（その２）である。
【図５】疑似ランダム配置を示す図である。
【図６】第１の平衡配置を示す図である。
【図７】ランダム配置を示す図である。
【図８】第２の平衡配置を示す図である。
【図９】ポテンシャルの変化を示す図である。
【図１０】情報処理装置の構成図である。
【図１１】記録媒体を示す図である。
【図１２】従来の第１のシミュレーション結果を示す図である。
【図１３】従来の第２のシミュレーション結果を示す図である。
【符号の説明】
１ 計算手段
２ 出力手段
１１ 自由度
１２ 粒子の初期座標
１３ 降下法の条件
１４ 境界条件
１５ カットオフ長
２１ ＣＰＵ
２２ メモリ
２３ 入力装置
２４ 出力装置
２５ 外部記憶装置
２６ 媒体駆動装置
２７ ネットワーク接続装置
２８ バス
２９ 可搬記録媒体
３０ データベース
Claims (5)
- 分子系のシミュレーションにおける分子集団の配置を生成する生成装置であって、
任意次元の空間内における分子集団を構成する複数の分子の初期座標と、微小変位幅を記憶する記憶手段と、
各分子の初期座標を位置座標として、該位置座標を正方向および負方向に前記微小変位幅だけ変位させて、前記複数の分子の位置座標により記述され、前記分子系の相互作用を表すポテンシャル関数の値を計算し、正方向の変位後のポテンシャル関数値、負方向の変位後のポテンシャル関数値、および変位前のポテンシャル関数値を比較し、最も小さい値に対応する位置座標を新規の位置座標として前記記憶手段に格納し、該新規の位置座標を用いて該微小変位幅の変位操作と該ポテンシャル関数の値の計算を繰り返すことで、該複数の分子の位置座標からなる前記分子集団の平衡配置を計算する計算手段と、
得られた平衡配置を前記シミュレーションにおける初期配置として出力する出力手段と
を備えることを特徴とする生成装置。 - 前記分子系の相互作用が力関数により記述されるとき、前記計算手段は、該力関数を用いて前記ポテンシャル関数を生成することを特徴とする請求項１記載の生成装置。
- 任意次元の空間内における分子集団を構成する複数の分子の初期座標と、微小変位幅を記憶する記憶手段と、
各分子の初期座標を位置座標として、該位置座標を正方向および負方向に前記微小変位幅だけ変位させて、前記複数の分子の位置座標により記述され、分子系の相互作用を表すポテンシャル関数の値を計算し、正方向の変位後のポテンシャル関数値、負方向の変位後のポテンシャル関数値、および変位前のポテンシャル関数値を比較し、最も小さい値に対応する位置座標を新規の位置座標として前記記憶手段に格納し、該新規の位置座標を用いて該微小変位幅の変位操作と該ポテンシャル関数の値の計算を繰り返すことで、該複数の分子の位置座標からなる前記分子集団の平衡配置を計算して該記憶手段に格納する計算手段と、
得られた平衡配置を初期配置として用いて前記分子系のシミュレーションを行うシミュレーション手段と、
前記シミュレーションの結果を出力する出力手段と
を備えることを特徴とするシミュレーション装置。 - 分子系のシミュレーションにおける分子集団の配置を生成するコンピュータのためのプログラムを記録した記録媒体であって、
記憶手段に格納された、任意次元の空間内における分子集団を構成する複数の分子の初期座標を位置座標として、該位置座標を正方向および負方向に、該記憶手段に格納された微小変位幅だけ変位させて、該複数の分子の位置座標により記述され、前記分子系の相互作用を表すポテンシャル関数の値を計算するステップと、
正方向の変位後のポテンシャル関数値、負方向の変位後のポテンシャル関数値、および変位前のポテンシャル関数値を比較し、最も小さい値に対応する位置座標を新規の位置座標として前記記憶手段に格納するステップと、
前記新規の位置座標を用いて前記微小変位幅の変位操作と前記ポテンシャル関数の値の計算を繰り返すことで、前記複数の分子の位置座標からなる前記分子集団の平衡配置を計算するステップと、
得られた平衡配置を前記シミュレーションにおける初期配置として出力するステップと
を含む処理を前記コンピュータに実行させるためのプログラムを記録したコンピュータ読み取り可能な記録媒体。 - 計算手段が、記憶手段に格納された、任意次元の空間内における分子集団を構成する複数の分子の初期座標を位置座標として、該位置座標を正方向および負方向に、該記憶手段に格納された微小変位幅だけ変位させて、該複数の分子の位置座標により記述され、分子系の相互作用を表すポテンシャル関数の値を計算し、
前記計算手段が、正方向の変位後のポテンシャル関数値、負方向の変位後のポテンシャル関数値、および変位前のポテンシャル関数値を比較し、最も小さい値に対応する位置座標を新規の位置座標として前記記憶手段に格納し、
前記計算手段が、前記新規の位置座標を用いて前記微小変位幅の変位操作と前記ポテンシャル関数の値の計算を繰り返すことで、前記複数の分子の位置座標からなる前記分子集団の平衡配置を計算して前記記憶手段に格納し、
シミュレーション手段が、得られた平衡配置を初期配置として用いて前記分子系のシミュレーションを行う
ことを特徴とするシミュレーション方法。
Priority Applications (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP7498799A JP3892167B2 (ja) | 1999-03-19 | 1999-03-19 | 粒子集団の配置を生成する生成装置および方法 |
| US09/435,215 US7065477B2 (en) | 1999-03-19 | 1999-11-05 | Device and method for generating an arrangement of a set of particles |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP7498799A JP3892167B2 (ja) | 1999-03-19 | 1999-03-19 | 粒子集団の配置を生成する生成装置および方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JP2000268064A JP2000268064A (ja) | 2000-09-29 |
| JP3892167B2 true JP3892167B2 (ja) | 2007-03-14 |
Family
ID=13563152
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP7498799A Expired - Fee Related JP3892167B2 (ja) | 1999-03-19 | 1999-03-19 | 粒子集団の配置を生成する生成装置および方法 |
Country Status (2)
| Country | Link |
|---|---|
| US (1) | US7065477B2 (ja) |
| JP (1) | JP3892167B2 (ja) |
Families Citing this family (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP5011689B2 (ja) * | 2005-09-15 | 2012-08-29 | 日本電気株式会社 | 分子シミュレーション方法及び装置 |
| JP5169225B2 (ja) * | 2007-05-25 | 2013-03-27 | 横浜ゴム株式会社 | 複合材料モデルの作成装置、シミュレーション装置、複合材料モデルの作成プログラム及び複合材料モデルの作成装置の動作方法 |
| JP5052985B2 (ja) * | 2007-07-31 | 2012-10-17 | 住友重機械工業株式会社 | 分子シミュレーション方法、分子シミュレーション装置、分子シミュレーションプログラム、及び該プログラムを記録した記録媒体 |
| US8443458B2 (en) * | 2008-01-16 | 2013-05-14 | Telefonaktiebolaget Lm Ericsson (Publ) | Method and device for detecting the opening of a cover enclosing the device |
| JP5517964B2 (ja) * | 2011-02-14 | 2014-06-11 | 住友重機械工業株式会社 | 解析装置 |
| US8855982B2 (en) | 2012-02-06 | 2014-10-07 | Sumitomo Heavy Industries, Ltd. | Analysis device and simulation method |
| JP5854067B2 (ja) * | 2014-02-28 | 2016-02-09 | 横浜ゴム株式会社 | 不均質材料のシミュレーションモデルの作成方法、不均質材料のシミュレーション方法、及びプログラム |
Family Cites Families (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2638359B2 (ja) * | 1991-11-26 | 1997-08-06 | 富士通株式会社 | 分子動力学法の拘束条件作成装置 |
| US5241470A (en) * | 1992-01-21 | 1993-08-31 | The Board Of Trustees Of The Leland Stanford University | Prediction of protein side-chain conformation by packing optimization |
| US5596511A (en) * | 1993-01-27 | 1997-01-21 | Fuji Xerox Co., Ltd. | Computing method and apparatus for a many-body problem |
| US5596703A (en) * | 1993-10-22 | 1997-01-21 | Lucent Technologies Inc. | Graphical display of relationships |
| US5553004A (en) * | 1993-11-12 | 1996-09-03 | The Board Of Trustees Of The Leland Stanford Jr. University | Constrained langevin dynamics method for simulating molecular conformations |
| US5600571A (en) * | 1994-01-18 | 1997-02-04 | The Trustees Of Columbia University In The City Of New York | Method for determining protein tertiary structure |
| US5724252A (en) * | 1994-12-09 | 1998-03-03 | Kirin Brewery | System for prediction of protein side-chain conformation and method using same |
| GB9616105D0 (en) * | 1996-07-31 | 1996-09-11 | Univ Kingston | TrkA binding site of NGF |
| US6188965B1 (en) * | 1997-04-11 | 2001-02-13 | California Institute Of Technology | Apparatus and method for automated protein design |
-
1999
- 1999-03-19 JP JP7498799A patent/JP3892167B2/ja not_active Expired - Fee Related
- 1999-11-05 US US09/435,215 patent/US7065477B2/en not_active Expired - Fee Related
Also Published As
| Publication number | Publication date |
|---|---|
| US7065477B2 (en) | 2006-06-20 |
| JP2000268064A (ja) | 2000-09-29 |
| US20020049572A1 (en) | 2002-04-25 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Wang et al. | Isogeometric analysis for parameterized LSM-based structural topology optimization | |
| Grimberg et al. | Mesh sampling and weighting for the hyperreduction of nonlinear Petrov–Galerkin reduced‐order models with local reduced‐order bases | |
| Kedward et al. | Efficient and exact mesh deformation using multiscale RBF interpolation | |
| Farhat et al. | Computational bottlenecks for PROMs: precomputation and hyperreduction | |
| König et al. | A flexible C++ framework for the partitioned solution of strongly coupled multifield problems | |
| Brüls et al. | The generalized-α scheme as a linear multistep integrator: toward a general mechatronic simulator | |
| Bauchau et al. | Scaling of constraints and augmented Lagrangian formulations in multibody dynamics simulations | |
| JP3892167B2 (ja) | 粒子集団の配置を生成する生成装置および方法 | |
| Butcher et al. | Transition curve analysis of linear fractional periodic time-delayed systems via explicit harmonic balance method | |
| Liu et al. | Optimal design of nonlinear multimaterial structures for crashworthiness using cluster analysis | |
| Meyer et al. | Co-simulation: Error estimation and macro-step size control | |
| Feng et al. | Subspace recycling accelerates the parametric macro‐modeling of MEMS | |
| JP5176895B2 (ja) | Ｓｒａｍ形状パラメータ等の多目的最適化設計支援装置、方法、及びプログラム | |
| Rastegarzadeh et al. | Neural network-assisted design: A study of multiscale topology optimization with smoothly graded cellular structures | |
| Sobota et al. | Implicit dynamic analysis using an isogeometric Reissner–Mindlin shell formulation | |
| Liu et al. | A sequential sampling generation method for multi-fidelity model based on Voronoi region and sample density | |
| Ji et al. | Two-step multi-resolution reconstruction-based compact gas-kinetic scheme on tetrahedral mesh | |
| Banerjee et al. | Graph-theoretic sensitivity analysis of multibody systems | |
| Wu et al. | Unstructured adaptive triangular mesh generation techniques and finite volume schemes for the air bearing problem in hard disk drives | |
| Eichmeir et al. | Time-optimal control of dynamic systems regarding final constraints | |
| Senecal | Efficient coupling algorithms and reduced-order methods for high-fidelity multiphysics simulations of nuclear reactors | |
| Mastrippolito et al. | RBF-based mesh morphing improvement using Schur complement applied to rib shape optimization | |
| Srinivasan et al. | Latency tolerance through parallelization of time in scientific applications | |
| Sun et al. | A novel node-to-segment algorithm in smoothed finite element method for contact problems | |
| Rastiello et al. | Path-following methods for unstable structural responses induced by strain softening: a critical review |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20040324 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20060912 |
|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20061109 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20061205 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20061206 |
|
| R150 | Certificate of patent or registration of utility model |
Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20091215 Year of fee payment: 3 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20101215 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20111215 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20111215 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20121215 Year of fee payment: 6 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20121215 Year of fee payment: 6 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20131215 Year of fee payment: 7 |
|
| LAPS | Cancellation because of no payment of annual fees |