WO2011021279A1 - 物質の応力を算出する情報処理装置 - Google Patents

Abstract

情報処理装置は、物質中の原子位置を含む原子構造を記憶する手段と、原子を含む結晶中の原子位置を含む原子構造を記憶する手段と、物質を複数に分割した分割部分の原子位置と原子を含む結晶の原子位置とを比較し、分割部分と結晶との間で対応する原子間の相対距離の評価値が最も小さくなるように、分割部分の原子位置と結晶の原子位置とを対応付け、結晶の単位格子に相当する、物質の分割部分を特定する分割手段と、分割部分の頂点と平行六面体の頂点との間で相対距離の評価値が最も小さくなるような平行六面体を決定する平行六面体形成手段と、平行六面体を繰り返した仮想結晶構造において平行六面体に加わる平均応力を算出する平均応力算出手段と、部分特定手段、平行六面体形成手段、平均応力算出手段を繰り返し機能させることによって物質のそれぞれの分割部分の応力を特定する制御手段と、を備える。

Description

物質の応力を算出する情報処理装置

　本発明は、物質の応力を算出する情報処理装置に関するものである。

　ナノデバイスの開発を行う際に、量子論に基づくシミュレーションを用いる傾向が年々高まっている。なぜなら、ナノデバイスの微細化が進む現在において、目的にあったデバイス開発には、古典力学ではなく量子力学（量子論）が適用されるナノスケールでの物理現象を理解することが有効と考えられるからである。ナノスケールでの物理現象の一例として、物質中の欠陥・亀裂の発生が挙げられる。ナノスケールの欠陥や亀裂の発生の理解には、量子論が適用される微細な領域の原子構造の安定性をより詳細に調べることが有効である。いずれにしても、ナノスケールの欠陥や亀裂の発生の理解、あるいは、微細な領域の原子構造の安定性調査のためには、局所応力の計算が有用である。

　量子論に基づいて物質のシミュレーションを行う際の計算手法には、タイトバインディング法など経験的パラメータを用いる計算手法と、経験的パラメータを用いない計算手法、例えば第一原理計算が存在する。

　経験的パラメータを用いた量子論に基づく計算手法の１つであるタイトバインディング（Tight-Binding）法とは以下のような計算手法である。

　通常、タイトバインディング法では、全系のエネルギーを以下のように表す。

　右辺の第一項は電子と電子、電子とイオンの相互作用によるエネルギーを示す。また、第二項はイオンとイオンの相互作用によるエネルギーと第一項の補正を表す。また、ε_nは一電子ハミルトニアンの固有値である。すなわち、一電子のシュレーディンガー方程式

の固有値である。タイトバインディング法では、例えば、一電子のシュレーディンガー方程式の固有関数である波動関数を数３のように、原子軌道の和で表す。一電子のシュレーディンガー方程式とは、シュレーディンガー方程式を平均場近似と呼ばれる手法で近似した式である。したがって、一電子のシュレーディンガー方程式は、電子１個の系を取り扱う式という意味ではない。

　ここで、φ_iは原子内の電子軌道、nは固有値のモード、Ｒ（太字）は並進ベクトル、τ（太字）_iは原子の位置、ｋ（太字）は波動関数の波数である。また、数式上太字で示されている、Ｒ（太字）、τ（太字）_i、ｋ（太字）を、明細書文中では、単に、Ｒ、τ_i、ｋのように通常のフォントで表すことにする。

　このように波動関数を原子軌道で展開することによって、一電子のシュレーディンガー方程式（微分方程式）を解く問題を、原子軌道の係数を求める問題（固有値問題）に置き換える。具体的には、シュレーディンガー方程式の両辺に左から、Σ_l0exp(-ikR_l0)φ_j ^*（ｒ－τ_j－Ｒ_l0）を掛け、全空間で積分することによって、一電子のシュレーディンガー方程式が固有値問題に置き換えられる。すなわち、

　ここで、Ｈ_ij（τ_i，τ_j），Ｓ_ij（τ_i，τ_j）を以下のように定義した。

　タイトバインディング法では、距離の関数であるＨ_ij（τ_i，τ_j），Ｓ_ij（τ_i，τ_j）を、固有値などが実験値を再現するように人為的に決定する。このとき、経験的パラメータを用いる。また、タイトバインディング法の特徴として、求める係数ｃ_iの個数が（全系の原子数）×（1原子当たりの価電子数）で済むため計算負荷が軽いという特徴を持つ。例えば、１００原子のシリコン原子の最外殻の電子を考慮してシミュレーションを行うためには、s軌道とp軌道の電子のみ扱う場合には４００次元の固有値問題を解けばよいことになる。

　一方、経験的パラメータを用いない量子論に基づく計算手法の１つとして、密度汎関数理論に基づく第一原理計算がある。その手法は以下の通りである。

　第一原理計算では、原子単位系での全系のエネルギーは

　と表すことができる。第一項目は電子の運動エネルギー、第二項目は電子間のクーロンエネルギー、第三項目は電子の交換相関項、第四項目は電子とイオンの相互作用エネルギー、第五項目はイオン間の相互作用エネルギーである。また、ｖ_ext（ｒ）はイオンが作るポテンシャルである。ここで、上記の|_i（ｒ）とｎ（ｒ）は電子間の相互作用がない仮想系における一電子の波動関数と電子密度である。|_i（ｒ）とｎ（ｒ）は、下記数７の方程式（Kohn Sham方程式）を自己無撞着に解くことによって計算することができる。

　Kohn Sham 方程式を解く際には、波動関数を平面波で展開する方法が広く使われている。波動関数を平面波で展開した形式を数８に例示する。

　このように波動関数を平面波で展開することによって、Kohn Sham方程式（微分方程式）を解く問題を、平面波の係数を求める問題（固有値問題）に置き換える。固有値問題を解く際は、タイトバインディング法と異なり経験的パラメータを用いずに行列要素を求める。また、第一原理計算の特徴として、平面波の係数の個数が非常に多いという特徴を持つ。これは、イオンがいるサイトに局在した波動関数を広がった平面波の重ね合わせで表現しなければならないためである。例えば、２原子を有するシリコン原子のシミュレーションに、数百次元の固有値問題を解く必要がある。したがって、経験的パラメータを用いない場合には、取り扱うことができる原子数がタイトバインディング法などの経験的パラメータを用いた場合と比較して少なくなる。

　現在、周期的構造を持った結晶中の応力の平均値を求める方法として、両手法に共通して適用できる理論が存在する（非特許文献１）。具体的な方法は、以下の通りである。
（１）シュレーディンガー方程式を解き、結晶のエネルギーと原子に働く力を求め、原子構造最適化、もしくは分子動力学計算を行う。
（２）（１）で求めた原子座標とシュレーディンガー方程式の解である波動関数を用いて、応力の平均値を求める。応力の平均値は、以下の数９の式で表される。

　ここで、Tαβは応力、rは原子座標、Ψは波動関数、ｍｉは原子ｉの質量、ｐにハット（＾）を付した要素は運動量の演算子、Vはポテンシャルである。ｉは、原子を識別するラベル、Σは、系に含まれる原子についての加算である。なお、αβは、ｘ，ｙ，ｚのいずれか２つを意味する。したがって、Tαβは、Ｔｘｘ，Ｔｘｙ（＝Ｔｙｘ），Ｔｘｚ（＝Ｔｚｘ），Ｔｙｙ，Ｔｙｚ(=Ｔｚｙ)のいずれかである。例えば、Ｔｘｙは、ｘ軸に垂直な面（ＹＺ平面）に働くｙ方向の応力を意味する。

　また、経験的パラメータを用いない計算手法においては、局所応力の計算手法も開発されている（非特許文献２）。開発されている手法の一例として密度汎関数理論に基づく第一原理計算を用いた局所応力の計算手法を挙げる。まず、基礎方程式であるKohn Sham方程式を解いて一電子波動関数と電子密度を求める。一電子波動関数と電子密度は三次元座標の関数として求めることができるので、密度汎関数理論では、全系のエネルギーが一電子波動関数および電子密度の汎関数で与えられる。したがって、エネルギーの微分量である応力も三次元座標の関数として得ることができ、局所応力を計算することができる。

特開２００３－３４７３０１

O. H. Nielsen and R. M. Martin, Phys. Rev. ，米国，B 32, 3780 (1985) A. Filippetti and V. Fiorentini, Phys. Rev. ，米国，B 61, 8433 (2000)

　すでに説明したように、量子論に基づいた物質のシミュレーションには、経験的パラメータを用いた手法と経験的パラメータを用いない手法がある。現在、応力の平均値の計算については、両者の手法を用いて実行することが可能である。一方、応力分布すなわち局所応力の計算については、経験的パラメータを用いない手法に対して手法開発が行われている。しかし、経験的パラメータを用いない手法は、取り扱うことができる原子数が数十原子～数百原子程度であるという欠点を持つ。したがって、精度は落ちるものの高速計算が可能である経験的パラメータを用いた量子論に基づく手法で、局所応力を計算する方法の開発が課題となる。

　ところで、経験的パラメータを用いない計算手法の場合には、数６に示したように、エネルギーが三次元座標の関数として与えられるので、エネルギーを三次元座標で微分することによって応力も三次元座標の関数として得ることができる。しかし、経験的パラメータを用いた計算手法の場合は、数１に示したように、エネルギーを三次元座標の関数をして表示することができないので、経験的パラメータを用いない手法のように微分操作により局所応力を得ることはできないという課題がある。

　タイトバインディング法に代表される経験的パラメータを用いた量子論に基づく手法で局所応力を算出することができれば、現実的な数の原子を含むナノ構造、例えば、数千個以上の原子を含む物質における局所応力の計算が可能になると期待される。

　そこで、開示の技術の目的は、経験的パラメータを用いた計算手法で得られた、三次元座標に対するエネルギー分布が得られない物質内の原子配置のシミュレーション結果に対して、物質内の応力分布を算出する技術を提供することにある。

　開示の技術の一側面は、情報処理装置として例示できる。本情報処理装置は、物質中の原子位置を含む原子構造を記憶する手段と、原子を含む結晶中の原子位置を含む原子構造を記憶する手段と、物質を複数に分割した分割部分の原子位置と原子を含む結晶の原子位置とを比較し、分割部分と結晶との間で対応する原子間の相対距離の評価値が最も小さくなるように、分割部分の原子位置と結晶の原子位置とを対応付け、結晶の単位格子に相当する、物質の分割部分を特定する分割手段と、分割部分の頂点と平行六面体の頂点との間で相対距離の評価値が最も小さくなるような平行六面体を決定する平行六面体形成手段と、平行六面体を繰り返し配置した仮想結晶構造において平行六面体に加わる平均応力を算出する平均応力算出手段と、分割手段、平行六面体形成手段、平均応力算出手段を繰り返し機能させることによって物質のそれぞれ分割部分の応力を特定する制御手段と、を備える。

　開示の技術によれば、経験的パラメータを用いた計算手法で得られた、三次元座標に対するエネルギー分布が得られない物質内の原子配置のシミュレーション結果に対して、物質内の応力分布を算出することができる。

本情報処理装置による処理概要の一側面を例示する図である。 原子構造を決定する場合のシミュレーションの手順を例示する図である。 物質の原子構造から、微小領域を求めた処理例である。 微小領域から平行六面体を作成する処理例である。 微小領域の境界面の応力から、任意の位置の応力を補間する処理例である。 不適切な分割の例である。 望ましいと考えられる微小領域の作成例である。 原子構造を微小領域に分割する処理を例示するフローチャートである。 選択した原子および選択した近傍の原子の原子種に対してリストアップされる結晶構造を例示する図である。 原子間の距離の二乗平均平方根（RMS）の算出例を示す図である。 分割する微小領域の決定手順を例示する図である。 結晶の基本格子内に２個の原子が存在する場合の処理例である。 結晶の基本格子内に２個の原子が存在する場合の処理例である。 結晶の基本格子内に２個の原子が存在する場合の処理例である。 結晶構造のデータベースのフォーマット例である。 単純立方格子の結晶構造を示す図である。 現実の物質、あるいは、シミュレーションによって求めた仮想系の原子構造を例示する図である。 六方最密充填構造を持つ結晶構造の例である。 シミュレーションで求めた仮想系の原子構造の例である。 微小領域から平行六面体への変形処理を例示するフローチャートである。 図形を用いて微小領域から平行六面体への変形処理を例示する図である。 ３頂点を選択する処理を例示する図である。 平行四辺形作成処理を例示するフローチャートである。 原子にラベル付けを例示する図である。 変化させた平行四辺形の例を示す図である。 平行四辺形の条件式を解いたときに得られる２種類の平行四辺形を例示する図である。 平行六面体作成処理を例示するフローチャートである。 平行六面体形成時に移動前の平行四辺形を例示する図である。 平行六面体形成時の移動前後の平行四辺形を例示する図である。 並進ベクトルに課せられる条件を例示する図である。 情報処理装置のハードウェア構成を例示する図である。 情報処理装置の機能構成を例示する図である。

　以下、図面を参照してを実施形態に係る情報処理装置について説明する。以下の実施形態の構成は例示であり、本情報処理装置は実施形態の構成には限定されない。

　本情報処理装置は、材料設計支援アプリケーションにおける、局所応力の計算手法に関する。特に、本情報処理装置は、ターゲット物質を微小領域に分割し周期構造を持った仮想系を想定することにより高速演算が可能な経験的パラメータを用いた量子論に基づいて局所応力を計算する。

　経験的パラメータを用いた手法では、全系のエネルギーを三次元座標で表される一電子波動関数の汎関数で得ることができないので、経験的パラメータを用いない量子論による手法と同じ手段では、エネルギーの微分量である応力を三次元座標の関数として得ることができない。

　そこで、本実施形態の情報処理装置は、物質の原子座標を求めるシミュレーション、例えば、原子構造最適化、もしくは分子動力学計算を行った後、計算した原子座標を用いて物質を微小領域に分割し、分割された物質の微小領域を平行六面体で近似する。そして、本情報処理装置は、その平行六面体を単位格子とする周期的構造を持った仮想系を想定することによって、微小領域の応力の平均値を計算する。最後に、本情報処理装置は、各微小領域の応力を補間することによって、応力分布を求める。

　＜処理手順＞
　図１に、本情報処理装置による処理概要の一側面を例示する。本情報処理装置は、まず、経験的パラメータを用いたシミュレーションによって、原子の座標を含む原子構造を決定する（Ｆ１）。図２に、原子構造を決定する場合のシミュレーションの手順を例示する。ここでは、本情報処理装置による応力分布算出の前提として、物質内の原子構造、すなわち、原子座標を求めるシミュレーションの概要を説明する。図２は、物質内の原子構造を求める処理フローチャートの例である。情報処理装置のＣＰＵは、メモリ上に実行可能に展開されたコンピュータプログラムを実行し、図２に例示したシミュレーションを実行する。以下、ＣＰＵがコンピュータプログラムを実行するときの処理を単に情報処理装置の処理として説明する。

　図２の処理では、まず、情報処理装置は、ユーザの指定にしたがって、シミュレーション対象の物質に含まれる原子の種類と、原子の位置の初期値を設定する（Ｆ１１）。原子の種類と、原子の位置の初期値は、情報処理装置が読み出すハードディスク等の外部記憶装置上のパラメータファイルに設定しておいてもよい。パラメータは、例えば、（原子の種類、Ｘ，Ｙ，Ｚ）のような形式で設定するようにしてもよい。ここで、Ｘ，Ｙ，Ｚは、原子の中心を設定する３次元座標の値である。

　また、例えば、ディスプレイ上に、既存の物質の結晶構造を表示し、ユーザがその結晶構造をマニュアル操作で変形するようにしてもよい。ここで、マニュアル操作とは、ディスプレイに表示された結晶構造を、マウス等のポインティングデバイスで変形する操作をいう。例えば、マウスカーソルでディスプレイに表示された結晶構造をドラッグできるようにしてもよい。また、既存の結晶構造の一部に、切断、亀裂等を生じさせるメニューを設けて、ユーザのメニュー選択に応じて、結晶構造をディスプレイ上で編集できるようにしてもよい。また、結晶構造中に、その結晶中の原子と同種類の原子、または他の種類の原子を挿入できるメニューを設けてもよい。ユーザが、結晶中に、原子を挿入する位置と、挿入される原子の種類を選択できるようにすればよい。以上のようなＦ１の処理にしたがい、パラメータファイルからのデータの読み出し、あるいは、既存の結晶構造の変形によって、ユーザが所望する原子の位置の初期値が設定される。

　次に、現在の原子の位置と当該原子の原子種を基に、数５の行列要素を経験的パラメータを用いて設定する（F 12）。例えば、τ_iにシリコン原子、τ_jに水素原子が配置している場合には、ｒ＝｜τ_i－τ_j｜だけ離れたシリコン原子と水素原子の電子的性質を再現するように経験的パラメータを用いて行列要素を求める。

　次に、情報処理装置は、Ｆ１２で求めた行列要素、Ｈ_ij（τ_i，τ_j），Ｓ_ij（τ_i，τ_j）を用い、数４の連立方程式を解くことで、ｃ_iを求める。そして、情報処理装置は、求めたｃ_iを用いて、数３と同様の形式でシュレーディンガー方程式の解である波動関数を求める。

　次に、情報処理装置は、Ｆ３で得られた波動関数を基に、原子間に作用する力を求める（Ｆ１４）。原子間に作用する力としては、例えば、原子の最外殻の電子による結合力を例示できる。

　次に、情報処理装置は、原子間に作用する力にしたがって、各原子に作用する力の和を求める。そして、情報処理装置は、それぞれの原子について、力の和が所定の許容値未満か否かを判定する。情報処理装置は、シミュレーションの対象の原子のすべてを確認し、原子に作用する力が、所定の許容値ＥＰＳ未満であるか、否かを判定する（Ｆ１５）。そして、シミュレーションの対象のすべての原子に対して、作用する力が所定の許容値ＥＰＳ未満である場合、情報処理装置は、シミュレーションを終了する。

　一方、シミュレーションの対象のいずれかの原子に作用する力が、所定の許容値ＥＰＳ以上である場合、情報処理装置は、それぞれの原子を力の方向に所定量Δ移動する。なお、原子の移動方向はそれぞれの原子に作用する力の方向とし、一方、原子を移動する所定量Δは、シミュレーション対象の原子に共通の値でもよい。また、原子を移動する所定量Δについても、原子に作用する力に応じた値としてもよい。そして、情報処理装置は、制御をＦ２に戻し、再度、原子の周囲のポテンシャルを求める。

　そして、原子に作用する力が所定の許容値ＥＰＳ未満となるまで、情報処理装置は、Ｆ２からＦ６の処理を繰り返す。このような手順で、情報処理装置は、Ｆ１で設定された原子の位置の初期値に対して、安定な原子構造を求める。図２のような原子の安定な位置を求める手法は原子構造最適化と呼ばれており、一般に、最急降下法や共役勾配法が用いられている。このとき得られる原子構造は、絶対零度における内部エネルギーが最低となる構造である。また、ある特定の温度・圧力下において原子構造を求める処理は、分子動力学計算と呼ばれる手法を適用したコンピュータシステムによっても実行されている。

　次に、再び図１に戻って説明する。図２のような処理によって原子構造が決定された後、決定された原子構造中の原子座標を用いて局所応力を求めるために、Ｆ２以下の３段階の手続きを行う。まず、情報処理装置は、原子構造が決定された物質を微小領域に分割する（Ｆ２）。情報処理装置は、原子が結晶を組んだときの単位格子を基準にして、分割した微小領域の形状を決定する。“基準にして決定する“とは、例えば、原子が結晶を組んだときの既知の単位格子と、Ｆ１で得られた原子構造を比較し、既知の結晶の単位格子から最もずれ量が小さい原子の組み合わせを取り出す処理をいう。図３は、物質の原子構造から、微小領域を求めた処理例である。図３の例では、円で囲まれた部分に含まれる原子の組み合わせを微小領域として取り出している。

　次に、情報処理装置は、分割した微小領域に働く応力を計算する。本情報処理装置は、分割した微小領域を平行六面体に近似し、分割した微小領域が周期的に並んでいる仮想系を設定する。周期的構造を持った結晶中の応力の平均値の計算手法については、経験的パラメータを用いる手法と用いない手法の両者に対して共通した理論が確立されている。そこで、本情報処理装置は、分割した微小領域が周期的に並んでいる仮想系を設定することで、確立された理論に基づいて量子計算を行い、仮想系における応力の平均値を求めることを可能にする。

　そこで、微小領域に働く応力算出のため、本情報処理装置は、Ｆ２で求めた微小領域から平行六面体を作成する（Ｆ３）。図４は、微小領域から平行六面体を作成する処理例である。本情報処理装置は、Ｆ２の処理により原子構造から取り出した微小領域を、底面が平行四辺形で所定の並進ベクトルを有する平行六面体に置き換える。図４は、平行六面体の１つの軸方向から見た平面図を示している。

　そして、情報処理装置は、平行六面体に働く応力を算出する（Ｆ４）。応力の計算では、情報処理装置は、Ｆ３で求めた平行六面体が無限に繰り返される３次元構造を仮定する。そして、情報処理装置は、規則的な平行六面体の繰り返し構造中で、１つの平行六面体に作用する応力の平均値を求める。そして、情報処理装置は、求めた応力の平均値を微小領域の応力として決定する。そして、情報処理装置は、Ｆ３とＦ４の処理をＦ２で求めた微小領域の数だけ繰り返す（Ｆ５）。

　Ｆ３－Ｆ５の処理によって、物質中の微小領域ごとの応力が求められる。微小領域ごとの応力は、微小領域の境界面の応力ということができる（図５参照）。次に情報処理装置は、それぞれの微小領域ごとの応力を基に、微小領域の境界面以外の応力を補間し、応力分布を求める（Ｆ６）。図５に、微小領域の境界面の応力から、任意の位置の応力を補間する処理例を示す。以上の処理によって、情報処理装置は、Ｆ１の処理で得られた物質の応力分布を求める。
＜原子構造の微小領域への分割＞
　ここでは、図１のＦ２の処理の詳細を説明する。ナノスケールの構造を持った物質、すなわち完全には結晶構造を形成していない物質を微小領域に分割する方法を以下に示す。分割を単に碁盤の目状に行うと、現実の系と微小領域を周期的に並べた仮想系との間で、応力計算に影響を与える原子構造が異なる可能性がある。

　図６に不適切な分割の例を示す。図６で現実の系は、実際の物質の原子構造を例示している。また、図６で仮想系は、図２のようなシミュレーションで決定された原子構造である。図６のように、現実の系と、シミュレーションによって決定された仮想系の原子構造とでは、必ずしも一致しない箇所が存在し得る。

　今、仮想系での微小領域の分割の仕方として、ＳＱ１とＳＱ２で示される微小領域に分割した場合を考える。ここで、微小領域ＳＱ１、ＳＱ２は、ともに、仮想系の原子構造を固定のメッシュに切断することで求めた微小領域である。

　微小領域ＳＱ１とＳＱ２とに対して、図１で説明した平行六面体を形成し、応力の平均を求める手順をとった場合、ＳＱ１に含まれる原子構造と、ＳＱ１の周囲の原子構造との関係が、現実の系とは一致しない。また、微小領域ＳＱ２についても、ＳＱ２に含まれる原子構造と、ＳＱ２の周囲の原子構造との関係が、現実の系とは一致しない。したがって、微量領域ＳＱ１、ＳＱ２について、それぞれ単純に平行六面体を形成し、平均応力を計算しても、現実の物質に生じる応力には整合しない可能性が高い。

　そこで、シミュレーションで求めた仮想系の原子構造を微小領域に分割する場合に、極力、現実の物質中の原子構造に合わせて、微小領域を作成することが望ましいと考えられる。例えば、現実の物質が、結晶構造の歪み、格子欠陥、転移、粒界、不純物原子の混入等を含み、原子構造が本来の物質の原子間の位置関係からずれていたとしても、極力、現実の物質の原子構造に合わせて微小領域を作成することが望ましい。現実の物質中の原子構造に合わせた微小領域を作成することで、仮想系の原子構造が現実の系と一致していなくても、仮想系での応力分布の計算結果が現実の系の応力分布に類似する可能性が高まるからである。

　図７に、望ましいと考えられる微小領域の作成例を示す。図７の例では、現実の規則的な原子構造に近いＳＱ３の部分と、ＳＱ３に隣接して、原子間の距離が圧縮されたＳＱ４の部分を含む。このような原子構造に対して、本情報処理装置は、図６のように固定領域に分割する代わりに、ＳＱ３、ＳＱ４のように、本来の原子構造中の基本格子に含まれる原子間の位置関係に相当する関係を優先して、微小領域を作成する。そのため、物質の原子種について知られている結晶構造を現実の物質の原子構造の参照モデルとして利用する。

　すなわち、本情報処理装置は、固定の領域に分割するのではなく、例えば、現実の結晶格子中の単位格子、あるいは、基本格子と、仮想系の原子構造との関係を比較し、現実の結晶格子中の単位格子、あるいは、基本格子に相当する原子の組み合わせを抽出し、微小領域とする。したがって、仮想系での原子構造が歪んだ状態であっても、本情報処理装置は、現実の結晶格子中の単位格子、あるいは、基本格子に相当する原子構造が歪んだ状態のものとなっている箇所を見いだし、微小領域とする。

　ここで、単位格子とは、結晶が有する繰り返し構造のうち、通常、辺の長さが最も短くなる格子をいう。単位格子は、単位セルと呼ばれる。また、結晶固有の対象性をよりよく示すため、複数の単位格子を組み合わせた繰り返し構造を、繰り返しの単位に定義することがある。複数の単位格子を組み合わせたものを含めて、結晶の繰り返し構造は、基本格子、あるいは、基本セルとも呼ばれる。以下、本実施例中では、単位格子および基本格子の双方の意味で、単位格子という用語を使用する。

　本情報処理装置は、図７のように、シミュレーションされた仮想系の原子構造中で、現実の結晶格子中の単位格子、あるいは、基本格子に相当する原子構造を抽出する。このような単位格子に相当する原子構造では、着目する微小領域と、微小領域の周囲の原子構造との関係は、現実の結晶格子の単位格子と周囲の原子構造との関係と類似する。したがって、微小領域を個々に抽出し、平行六面体に置き換え、規則的な平行六面体の並びから、応力の平均値を求めたとしても、微小領域と微小領域の周囲の原子構造との関係は、現実の物質に近いものとなり得る。

　そこで、本情報処理装置は、まず注目した原子の原子種と、注目した原子近傍の原子の原子種とから、注目した原子および近傍の原子が結晶を形成したときの結晶構造を１以上特定する。そして、本情報処理装置は、特定した結晶構造の基本格子を基準にして微小領域の形状を決定し、物質を微小領域に分割する。単位格子を用いる理由は、図１のＦ３の処理において微小領域を周期的に並べるために適しているためである。

　図８から図１４に、仮想系を微小領域に分割する具体的な手順を示す。図８は、仮想系を微小領域に分割する処理手順を例示するフローチャートである。図９から図１４は、本情報処理装置による原子構造の微小領域への分割処理の例である。図８－図１４に示す処理も、本情報処理装置のＣＰＵがメモリ上に実行可能に展開されたコンピュータプログラムを実行することによりなされる処理である。さらに、図１５に、本情報処理装置が有する結晶構造のデータベースのフォーマットを例示する。

　本情報処理装置は、物質の原子構造の微小領域への分割のため、予め結晶構造と原子種に対応する格子定数、単位ベクトル、原子座標のデータを用意し、結晶構造のデータベースに格納している。例えば、図１５の例では、原子種としてシリコン原子が例示されている。そして、シリコン原子、すなわち、当該原子種がとり得る結晶構造の数＝２が定義されている。

　そして、当該原子種がとり得る結晶構造の数＝２に対応して、２つの結晶構造が定義される。それぞれの結晶構造として、格子定数、結晶構造の単位ベクトル、および原子座標が定義される。ここで、格子定数は、格子の長さである。例えば、立方格子の場合、格子定数は１つである。一方、立方格子以外の格子では、一般には、格子定数は、３つ定義される。例えば、第１の軸方向の格子定数ａ、第２の軸方向の格子定数ｂ、第３の軸方向の格子定数ｃ等である。ここで、第１の軸方向、第２の軸方向、第３の軸方向とは、それぞれ、結晶格子の方向をいう。

　なお、物質の結晶中に六面体を仮定した場合、８頂点のすべての位置に、原子が存在するように六面体を定義できない物質もあり得る。しかしながら、例えば、六面体の８頂点のうち、５頂点に原子が存在する六面体を定義できれば、５頂点の原子座標から、原子の存在しない頂点の位置は特定できる。そこで、本情報処理装置は、六面体の頂点の一部に原子が存在しない物質については、原子が存在する六面体の頂点から原子の存在しない頂点位置を推定する。

　また、例えば、六面体の中心に原子が存在するような結晶構造の単位格子を定義したとしても、格子定数ａ、ｂ、ｃを基に六面体の頂点を推定できる。例えば、原子の位置（０，０，０）に対して、ａ／２、ｂ／２、ｃ／２を加算あるいは減算した位置として、頂点を特定できる。したがって、原子の座標と、六面体の頂点座標との関数関係を維持することで、原子の位置から結晶構造の六面体の頂点を決定できる。

　また、単位ベクトルとは、結晶格子の３軸の単位ベクトルとである。図１５では、例えば、（ａ１ｘ，ａ１ｙ，ａ１ｚ）、（ａ２ｘ，ａ２ｙ，ａ２ｚ）、（ａ３ｘ，ａ３ｙ，ａ３ｚ）の３つの単位ベクトルで示される。

　また、原子座標とは、結晶中の原子の位置を示す座標である。なお、原子座標の原点は、結晶中の所定の位置、例えば、結晶の中心、結晶格子の１つの頂点等に設定される。結晶構造のデータベースに、原子座標が含まれるのは、単位ベクトルと格子定数で結晶格子を作成できたとしても、格子の頂点に原子座標があるとは限らないからである。

　次に、図８を参照して、原子構造を微小領域に分割する処理を説明する。まず、本情報処理装置は、シミュレーションにより求めた原子座標のデータから原子座標を１つ選択する（Ｆ２１）。

　そして、本情報処理装置は、選択した原子と、選択した近傍の原子の原子種を調べる（Ｆ２２）。ここで、近傍とは、例えば、選択した原子と、原子構造中で選択した原子に隣接する位置に置かれた原子とすればよい。そして、本情報処理装置は、選択した原子と選択した近傍の原子の原子種とが結晶を形成する際にとり得る結晶構造を結晶構造のデータベースからリストアップする。そして、本情報処理装置は、各結晶構造の格子定数、単位ベクトル、原子座標を結晶構造のデータベースから取り出す（Ｆ２３）。

　例えば、注目した原子とその近傍の原子の原子種がすべてシリコン原子であるならば、シリコン原子の結晶構造はダイヤモンド構造かβスズ構造である。また、シリコン原子の結晶構造として、ダイヤモンド構造とβスズ構造は、それぞれ、結晶構造のデータベースに図１５のフォーマットで定義されている。そこで、本情報処理装置は、例えば、シリコン原子に対する２つの格子定数、単位ベクトル、原子座標を結晶構造のデータベースから取り出す。

　図９に、Ｆ２１で選択した原子および選択した近傍の原子の原子種に対して、リストアップされる結晶構造を例示する。なお、図９では、Ｆ２１で選択した原子は、注目する原子とされている。

　次に、本情報処理装置は、リストアップしたいくつかの結晶構造から、１つの結晶構造を選択する（Ｆ２４）。そして、本情報処理装置は、Ｆ２１で選択した原子を中心とし、格子定数の２倍を半径とする球内の原子座標をシミュレーションにより求めた物質の原子座標データから取り出す（Ｆ２５）。

　そして、シミュレーションから求めた原子座標と選択した結晶構造のデータベースにある原子構造中の原子座標の差の二乗平均平方根（RMS）を計算する（Ｆ２６）。格子定数の２倍を半径とする球内の原子座標の数が、複数ある場合には、原子構造の原子座標の差の二乗平均平方根（RMS）は、最小値を選択する。シミュレーションから求めた原子座標と選択した結晶構造のデータベースにある原子構造中の原子座標の差の二乗平均平方根（RMS）が相対距離の評価値に相当する。例えば、格子定数の２倍を半径とする球内に、Ｆ２１で選択した原子の座標以外に、Ｔ１－Ｔｋの原子座標がある場合を仮定する。

　まず、本情報処理装置は、Ｆ２１で選択した原子の位置を結晶格子中のいずれかの原子の位置に一致させる。そして、本情報処理装置は、結晶中で、Ｆ２１で選択した原子が位置する結晶中の位置以外の原子の位置（例えば、Ｐ１からＰｋまで）を抽出する。そして、本情報処理装置は、シミュレーションから求めた原子座標Ｔ１からＴｋまでの原子座標と、結晶中の原子の位置Ｐ１－Ｐｋとの間で、原子間の距離の二乗平均平方根（RMS）を以下の式にしたがって計算する。

　ＲＭＳ＝（（１／Ｎ）×ΣＤｉ×Ｄｉ）^1/2  ；ここで、Ｄｉは、対応する原子対ｉにおける原子間の距離であり、Σは、ｉ＝１からｋまでの和であり、ｋは、ＲＭＳ計算の対象となる原子の数である。

　すでに説明したように、原子種によって、あるいは、結晶構造によっては、格子の格子点に原子が存在するとは限らない。しかしながら、図１５に示したように、結晶構造のデータベースは、原子種の結晶構造ごとに、原子座標を保持するので、結晶構造がリストアップできれば、本情報処理装置は、当該原子種における当該結晶構造での原子の位置を特定することができる。

　そして、ミュレーションから求めた物質中の原子座標Ｔ１からＴｋまでと、結晶中の原子座標Ｐ１からＰｋまでとの間で、原子間の距離の二乗平均平方根（RMS）が最小となる原子の位置の対応関係を求める。そして、本情報処理装置は、原子間の距離の二乗平均平方根（RMS）の最小値を保持する。

　図１０に、原子間の距離の二乗平均平方根（RMS）の算出例を示す。図１０では、Ｆ２１で選択した、注目する原子の位置が、結晶構造の中央の原子の位置に設定されている。そして、結晶中の中央の原子の周囲の原子との間で、８個の距離Ｒｉ（ｉ＝１～８）が算出されている。

　本情報処理装置は、Ｆ２２－Ｆ２６の一連の手順、すなわち、結晶構造の選択、シミュレーションから求めた原子座標からのデータの取り出し、RMSの計算を、リストアップしたすべての結晶構造について繰り返す（Ｆ２７）。

　そして、本情報処理装置は、各結晶構造で計算したRMSを比較し、最もRMSが小さい結晶構造の単位格子を決定する。そして、本条法処理装置は、最もRMSが小さい結晶構造の単位格子を基準にして、分割する領域の形状を決定する（Ｆ２８）。すなわち、情報処理装置は、最もRMSが小さい結晶構造の単位格子中の原子に相当する、シミュレーションから求めた原子座標を求める。そして、最もRMSが小さい結晶構造の単位格子中の原子に相当する、シミュレーションから求めた原子座標が形成する微小領域を求める。

　図１１に、分割する微小領域の決定手順を例示する。例えば、図１１において、結晶中の原子の位置Ｐ１－Ｐ４と、シミュレーションから求めた物質中の原子座標Ｔ１－Ｔ４との対応付けで、原子間の距離の差の二乗平均平方根（RMS）が最小となる場合に、情報処理装置は、原子座標Ｔ１－Ｔ４によって得られる領域を微小領域とする。

　なお、すでに述べたように、複数の頂点の一部に、原子が存在しない場合には、原子が存在する頂点から、原子が存在市内頂点の位置を推定する。例えば、８頂点中３頂点に原子が存在しない場合、原子が存在する５頂点についてＲＭＳが最も小さくなる六面体を作成すればよい。また、結晶中の単位格子の８頂点以外の位置に原子座標がある場合には、まず、結晶構造中の原子と、シミュレーションされた物質中の原子との間で、座標の差分値のＲＭＳが最も小さくなる対応関係を求める。そして、結晶中での原子座標に対して、格子定数から得られる六面体の位置を特定する座標間の関数関係を求めておく。そして、結晶中の原子座標に対応する物質のシミュレーションされた原子座標に対して、上記座標の関数関係を適用し、物質の微小領域を求めればよい。いずれにしても、本情報処理装置によれば、微小領域は、六面体の形状である。

　なお、最小のRMSが所定の限界値より大きくなる場合には、注目した原子近傍での応力計算を行わない（Ｆ２９）。ここで、所定の限界値は、例えば、所定の整数Ｎを用いて、格子定数のＮ倍のように、情報処理装置のパラメータとして、メモリに設定しておけばよい。本情報処理装置は、Ｆ２２－Ｆ２９の操作を、シミュレーションで得た原子座標に対して行い（Ｆ２Ａ）、領域の分割を行う。

　図１２－図１４に、他の結晶格子を基準にして分割した処理例を示す。図１２－１４では、結晶の基本格子内に２個の原子が存在する場合の処理例である。図１２は、Ｆ２１で選択した原子と、選択した近傍の原子の原子種と、リストアップされる結晶構造の例を示す。

　図１２のように、本情報処理装置は、Ｆ２１で選択した原子と、選択した近傍の原子の原子種を基に、結晶構造のデータベースを検索し、結晶構造をリストアップする。

　そして、図１３のように、結晶中の原子の位置と、シミュレーションから求めた原子座標を対応付けて、原子間の距離の差の二乗平均平方根（RMS）が最小となる場合を求める。そして、図１４のように、各結晶構造のうち、原子間の距離の差の二乗平均平方根（RMS）が最小となる対応関係を基に、結晶の基本格子中の原子に対応するシミュレーションから求めた原子座標を求める。

　なお、結晶の並進対称性を特徴付けるものとして、ブラベー格子がある。ブラベー格子の種類は全部で１４種類存在する。具体的には、１４種類のブラベー格子は、単純三斜、単純単斜、底心単斜、単純斜方、体心斜方、面心斜方、底心斜方、単純六方、単純菱面体、単純正方、体心正方、単純立方、体心立方、面心立方である。

　例えば、注目した原子種が単純立方のブラベー格子を有し、単純立方構造の結晶構造をとるときにRMSが最小になるのであれば、本情報処理装置は、微小領域が立方体に近い形状になるように分割する。図１６Ａに、単純立方格子の結晶構造を示す。図１６Ａは、理想の結晶構造ということができる。また、図１６Ｂに、図１６Ａと同一の原子種で、現実の物質、あるいは、シミュレーションによって求めた仮想系の原子構造を例示する。図１６Ｂのように、現実の物質の原子構造、あるいは、シミュレーションで求めた仮想系の原子構造は、理想の結晶構造からずれた歪んだ構造をとる場合が多い。したがって、現実の物質の原子構造、あるいは、シミュレーションで求めた仮想系の原子構造を、寸法、形状が固定された微小領域に分割しても、原子種特有の結晶構造に対応する微小領域を得ることはできない。一方、図８に示した手順にしたがって、結晶構造中の原子と、現実の物質の原子構造、あるいは、シミュレーションで求めた仮想系の原子構造中の原子との間で、原子間の距離の差の二乗平均平方根（RMS）が最小となる場合を求めることで、原子種本来の結晶構造が歪んだ状態を取得できる。

　図１７Ａは、六方最密充填構造を持つ結晶構造の例である。図１７のように、原子の座標は、六角柱の頂点、端面の中心、および、六角柱の端面で、端面の中心に位置する原子に位置する。さらに、端面の６角形の頂点のうちの矩形を形成する４頂点と６角形の中心を結ぶ線で形成される２つの矩形と、２つの三角形とを形成したときに、２つの矩形を端面とする四角柱、および２つの三角形を端面とする三角柱にが形成される。このうち、四角柱の中心部にも原子が配置される。また、三角柱は、隣接する六角柱に含まれる三角柱と組み合わせことで、四角柱を形成し、形成された四角柱の中心部にも原子が配置される。このような六方最密充填構造を持つ結晶構造においても、図１７Ａの矢印で「結晶の単位格子」の文字で示される位置に、六面体を定義できる。

　また、図１７Ｂは、図１７Ａの結晶と同一の原子種による現実の物質の原子構造推定のため、シミュレーションで求めた仮想系の原子構造の例である。単純六方のブラベー格子を有し、六方最密充填構造の結晶構造の原子と、仮想系の原子構造の原子との間で、ＲＭＳが最小になるのであれば、本情報処理装置は、微小領域として底面が菱形の四角柱に近い形状になるように分割する。分割のサイズは、結晶における単位格子の体積程度とすればよい。

　＜微小領域の平行六面体への変形＞
　図１のＦ３およびＦ４で説明したように、分割した微小領域を周期的に並べるためには、領域を平行六面体に変形する必要がある。その理由は、微小領域の応力を計算するためには微小領域の並進操作で作られる周期的構造を用意する必要があり、並進操作によって３次元空間を埋めることができる多面体は平行六面体であるためである。また、周期的構造は、一般的に知られている量子論に基づく計算手法を適用するために用意した。本情報処理装置は、変形による局所応力の変化を最小にするために、元の原子配置からの変位の二乗平均平方根（RMS）が最小になるような変形を行う。

　図１８Ａは、微小領域から平行六面体への変形処理を例示するフローチャートである。また、図１８Ｂに、図形を用いて微小領域から平行六面体への変形処理を例示する。図１８Ａ、図１８Ｂの処理も、本情報処理装置のＣＰＵがメモリ上に実行可能に展開されたコンピュータプログラムを実行することによりなされる処理である。

　まず、情報処理装置は、図１のＦ２の処理で取り出した微小領域の８つの頂点のうち、１つの頂点を選択する（Ｆ３１）。選択した頂点の位置を図１８Ｂに例示する。

　次に、本情報処理装置は、Ｆ３１で選択した頂点を含む平行四辺形の作成に用いられる残り３頂点を選択する（Ｆ３２）。平行六面体に変形される微小領域としては、六面体が仮定されているので、Ｆ３２の処理は、Ｆ３１で選択した頂点とともに四角形を形成する他の３頂点を選択する処理となる。

　図１９に、３頂点を選択する処理を例示する。７つの頂点から３つの頂点を選ぶときには、図１のＦ２の処理において領域の形を決定する上で基準にした空間格子の単位ベクトル（ａ１->、ａ２->、ａ３->）を用いる。なお、明細書文章中では、ベクトルａ１を“ａ１->”のように表示する。

　例えば、選択した原子を原点としたときは、下記の３通りの選び方がある。
（１）ａ１->、ａ２->、ａ１->＋ａ２->に位置する格子点に相当する３つの原子；
（２）ａ２->、ａ３->、ａ２->＋ａ３->に位置する格子点に相当する３つの原子；
（３）ａ３->、ａ１->、ａ３->＋ａ１->に位置する格子点に相当する３つの原子；
以上の３パターンについて、平行四辺形を作成する。

　そして、本情報処理装置は、選択した頂点とその他の３つの頂点に囲まれた領域が平行四辺形になるようにその他の３つの頂点を変位させる（Ｆ３３）。Ｆ３３の処理では、平行四辺形は任意の辺の長さを有し、平行四辺形がなす平面の角度にも制限はない。Ｆ３３の処理は、図２０－図２３においてさらに詳細に説明する。

　次に、本情報処理装置は、平行六面体を作成する。本情報処理装置は、平行六面体として、「６つの平面によって囲まれ、２つの相対する平面が平行四辺形で合同かつ平行であり、残りの４つの平面が平行四辺形」という定義にしたがった処理を実行する。すなわち、本情報処理装置は、Ｆ３３で作成した平行四辺形を平行移動させた平行六面体の候補を作成する。ここでは、Ｆ３２の処理で選択され、Ｆ３３で平行四辺形となるように移動された４つの頂点以外の微小領域中の残りの４頂点を移動する。すなわち、本情報処理装置は、平行六面体の候補の頂点と、微小領域の頂点が重なるように、残りの４つの頂点を変位させ、平行六面体を作成する（Ｆ３４）。Ｆ３４の処理は、図２４においてさらに詳細に説明する。

　本情報処理装置は、選択した頂点以外の３頂点を変更し、Ｆ３３の処理を繰り返す（Ｆ３５）。さらに、本情報処理装置は、選択する頂点を変更し、領域にある８つの頂点に対して行うＦ３２－Ｆ３５の操作を実行する（Ｆ３６）。

　そして、本情報処理装置は、動かした７つの頂点、すなわち原子の変位のRMSを計算し、RMSが最も小さい平行六面体を仮想系の単位格子として用いる（Ｆ３７）。上記の手順では、多くの平行六面体を作成するが、「動かした原子の変位のRMSが最小になる平行六面体を選ぶ」という条件により、唯一の平行六面体を得ることができる。

　図２０に、図１８ＡのＦ３３の処理、すなわち、平行四辺形作成処理を例示する。以下、Ｆ３１あるいはＦ３６の処理で選択した頂点とその他の３つの頂点に囲まれた領域が平行四辺形になるように、選択した頂点以外の３つの頂点を本情報処理装置が変位させるときの手順を説明する。ここで、図２１のように原子にラベル付けを行う。

　選択した原子（原子１）とその他の３つの原子（原子２，３，４）を頂点とする平行四辺形を作成するためには、４つの原子が同一平面上に存在するという条件が満たされる。そこで、原子１を含み法線ベクトル（a,b,c）を持つ任意の平面を考え、その平面上に原子２，３，４を変位させる。この平面を満たす方程式は、以下の通りである。

　ａ（ｘ－ｘ１）＋ｂ（ｙ－ｙ１）＋ｃ（ｚ－ｚ１）＝０（式１）
ここで、（ｘ１、ｙ１、ｚ１）は原子１の座標である。

　用意した平面を面Ａとすると、原子２，３，４は面Ａ上に変位することになるから、以下の式が成り立つ。

　ａ（ｘ２－ｘ１）＋ｂ（ｙ２－ｙ１）＋ｃ（ｚ２－ｚ１）＝０（式２）
　ａ（ｘ３－ｘ１）＋ｂ（ｙ３－ｙ１）＋ｃ（ｚ３－ｚ１）＝０（式３）
　ａ（ｘ４－ｘ１）＋ｂ（ｙ４－ｙ１）＋ｃ（ｚ４－ｚ１）＝０（式４）
ここで、（ｘｉ、ｙｉ、ｚｉ）は原子iの座標である。

　原子１、２，３，４を頂点とする四角形は、平行四辺形でなければならないので、対辺の長さが等しいという条件から式５および式６が成立する。

　（ｘ２－ｘ１）²＋（ｙ２－ｙ１）²＋（ｚ２－ｚ１）²
 ＝（ｘ３－ｘ４）²＋（ｙ３－ｙ４）²＋（ｚ３－ｚ４）² （式５）
　（ｘ４－ｘ１）²＋（ｙ４－ｙ１）²＋（ｚ４－ｚ１）²
 ＝（ｘ３－ｘ２）²＋（ｙ３－ｙ２）²＋（ｚ３－ｚ２）² （式６）
　式（２）－（６）より、任意の法線ベクトル（a,b,c）を与えたとき、未定数ｘ２，ｘ３，ｘ４，ｙ２，ｙ３，ｙ４，ｚ２，ｚ３，ｚ４のうち４つの値を与えれば、平面Ａ上に原子１、２、３、４の原子座標を頂点とする平行四辺形を作成できる。なお、原子１は変形させないので、座標（ｘ１、ｙ１、ｚ１）はそのまま使うため既知である。

　ただし、例えば４つの変数、ｘ２，ｘ３，ｙ２，ｙ３の値を与えたとき作成される平行四辺形は２種類存在する。つまり式５、式６を解くとの解が２つ存在するので、情報処理装置は、原子４の変位が小さいほうを採用する。図２３に、式５、式６を解いたときに得られる２種類の平行四辺形を例示する。

　上記の操作を法線ベクトルの値と４つの変数ｘ２，ｘ３，ｙ２，ｙ３を変化させて、様々な平行四辺形を用意する。ｘ２，ｘ３，ｙ２，ｙ３を変化させる手順を以下に例示する。まず、ｘ２，ｘ３，ｙ２，ｙ３の変化の範囲は、元の原子位置からユーザが設定した半径を持つ球内とする。変化の刻みは球を等間隔で刻んだメッシュ幅とする。また、法線ベクトル（a,b,c）を変化させる手段として、その変化の範囲は、極座標表示で、θは0°<θ<180°、φは0°＜φ＜360°とする。刻み幅は、極座標表示でのθ、φに対して、ユーザが設定した刻み幅とする。以上の範囲で以上の変化の刻みで、ｘ２，ｘ３，ｙ２，ｙ３、ａ，ｂ，ｃを変化させたときの平行四辺形を用意する。図２２に、変化させた平行四辺形の例を示す。

　図２０は、本情報処理装置による平行四辺形作成処理を例示するフローチャートである。本情報処理装置は、原子１を含む平面Ａの法線ベクトル（ａ，ｂ，ｃ）の初期値を与える（Ｆ３２１）。初期値は、例えば、極座標で（θ，φ）＝（Δθ，Δφ）としてもよい。ここで、Δθ，Δφは、それぞれ、θおよびφを変化させるときの刻みである。また、原子１、原子２および原子３の中心座標が含まれる平面の法線ベクトルとしてもよい。

　次に、本情報処理装置は、４つの変数（ｘ２，ｘ３，ｙ２，ｙ３）の初期値を与える（Ｆ３２２）。４つの変数（ｘ２，ｘ３，ｙ２，ｙ３）の初期値は、例えば、ユーザが設定した半径を持つ球内で、等間隔で刻んだメッシュの端点とすればよい。

　次に、本情報処理装置は、式２－式６を解き、（ｘ４，ｘ４，ｚ２，ｚ３，ｚ４）を求め、平面Ａ上で平行四辺形を作成する。このとき、（Ｘ４，Ｙ４）の解は２つ存在する。そこで、上述のように、原子４からの変位の少なく方を選択する（Ｆ３２３）。

　次に、本情報処理装置は、４つの変数（ｘ２，ｘ３，ｙ２，ｙ３）を所定の変化の刻み幅だけ変化させる（Ｆ３２４）。そして、本情報処理装置は、制御をＦ３２３に戻す。このようにして、情報処理装置は、元の原子位置から、ユーザが設定した半径を持つ球内で、Ｆ３２３の処理を実行する。

　次に、本情報処理装置は、平面の法線ベクトル（a,b,c）を変化させる（Ｆ３２５）。そして、本情報処理装置は、制御をＦ３２２に戻す。このようにして、情報処理装置は、極座標表示で、θは0°<θ<180°、φは0°＜φ＜360°の範囲で法線ベクトルを変化させ、Ｆ３２２、Ｆ３２３の処理を繰り返し実行する。

　図２４は、図１８ＡのＦ３４の処理、すなわち、平行六面体作成処理を例示するフローチャートである。図２５に、平行六面体形成前、すなわち、移動前の平行四辺形を例示する。図２６に、平行六面体形成前、すなわち、移動前後の平行四辺形を例示する。

　まず、本情報補処理装置は、作成した平行四辺形を平行移動させるときの並進ベクトル(p,q,r)を決定する（Ｆ３４０）。並進ベクトル作成のため、本情報処理装置は、平行四辺形を含む平面Ａの法線ベクトルの方向がｚ軸方向になるように新たに座標系を用意する（Ｆ３４１）。

　そして、並進ベクトル(p,q,r)を以下の式７、式８および式９の範囲内で設定する（Ｆ３４２）。今、平行四辺形の４頂点を（ｘｉ、ｙｉ，ｚｉ）（ｉ＝１から４）とする。また、微小領域の残りの４頂点を（ｘｉ、ｙｉ，ｚｉ）（ｉ＝５から８）とする。

　MIN(x5,x6,x7,x8)-MIN(x1,x2,x3,x4) < MAX(x5,x6,x7,x8)-MAX(x1,x2,x3,x4)のとき、
MIN(x5,x6,x7,x8)-MIN(x1,x2,x3,x4)<p<MAX(x5,x6,x7,x8)-MAX(x1,x2,x3,x4);
　MAX(x5,x6,x7,x8)-MAX(x1,x2,x3,x4) < MIN(x5,x6,x7,x8)-MIN(x1,x2,x3,x4)のとき、
MAX(x5,x6,x7,x8)-MAX(x1,x2,x3,x4)<p<MIN(x5,x6,x7,x8)-MIN(x1,x2,x3,x4);(式７)
　図２７に、式７で示される関係を例示する。本情報処理装置は、平行四辺形を平面Ａから移動する際、移動の並進ベクトルのｘ成分ｐとして、移動前の平行四辺形のｘ座標x1,x2,x3,x4と、移動後の平行四辺形のｘ座標x5,x6,x7,x8との間で、座標の最小値側の変化量Δ１と、最大値側の変化量Δ２を算出する。

　Δ１＝MIN(x5,x6,x7,x8)-MIN(x1,x2,x3,x4);
　Δ２＝MAX(x5,x6,x7,x8)-MAX(x1,x2,x3,x4);
そして、本情報処理装置は、Δ１とΔ２で設定される範囲で、ｐを変化させる。並進ベクトルのｙ成分ｑの算出も、ｘ成分と同様である。変化の刻みは、ユーザ指定値とすればよい。

　MIN(y5,y6,y7,y8)-MIN(y1,y2,y3,y4) < MAX(y5,y6,y7,y8)-MAX(y1,y2,y3,y4)のとき、
MIN(y5,y6,y7,y8)-MIN(y1,y2,y3,y4)<q<MAX(y5,y6,y7,y8)-MAX(y1,y2,y3,y4);
　MAX(y5,y6,y7,y8)-MAX(y1,y2,y3,y4) < MIN(y5,y6,y7,y8)-MIN(y1,y2,y3,y4)のとき、
MAX(y5,y6,y7,y8)-MAX(y1,y2,y3,y4)<q<MIN(y5,y6,y7,y8)-MIN(y1,y2,y3,y4);（式８）
また、並進ベクトルのｚ成分は、微小領域の残りの４頂点の最小値と最大値の範囲で変化させる。変化の刻みは、ユーザ指定値とすればよい。

　MIN(z5,z6,z7,z8) <r< MAX(z5,z6,z7,z8)  （式９）
　そして、本情報処理装置は、平行四辺形を並進ベクトル（ｐ、ｑ、ｒ）にしたがって平行移動させ、平行六面体とする（Ｆ３４３）。そして、本情報処理装置は、平行移動させた平行四辺形の頂点と重なるように原子５，６，７，８を変位させる。以上の作業を各並進ベクトルに対して行い、変位量のRMSが最小になる並進ベクトル(p,q,r)を採用する（Ｆ３４４）。そして、本情報処理装置は、採用した並進ベクトル(p,q,r)にしたがって、平行六面体を作成する。

　＜応力分布の計算＞
　本情報処理装置は、図１のＦ３の処理で得られた平行六面体を周期的に配列し、周期的に並んでいる仮想的な結晶に対して、経験的パラメータを用いた量子論に基づく手法で、応力の平均値を求める。平均応力は、数９にしたがって求める。数９にしたがった応力の平均値算出処理の詳細は、以下のとおりである。

　一般に、応力は全エネルギーを歪みの量ε_αβで微分することによって、数１０のように求めることができる。

タイトバインディング法を用いる場合は、全エネルギーが数１で表すことができるので、応力は以下のようになる。

ここで、数５を用いると、

今、数５のＳ_ijが単位行列の場合を示した。
Ｒは歪みがある状況では、下記のように表すことができる。

従って、

一方、

　また、本情報処理装置は、物質の微小領域を平行六面体に置き換え、平行六面体をＸ、ＹおよびＺ方向に無限に繰り返す構造を仮定する。平行六面体をＸ、ＹおよびＺ方向に無限に繰り返す構造では、物質が均一であることが仮定できるので、平均を意味する＜＞は、単純に１個の平行六面体内での積分記号に変形できる。そこで、数１０は、平行六面体内の応力を波動関数Ψで平均する式となっている。

　数１０で、波動関数Ψで平均される応力中の第１項、すなわち、運動量の演算子の積で示される第１項は、原子の運動エネルギーを１回微分した式である。この第１項は、原子の運動エネルギーから生じる圧力の項である。また、波動関数Ψで平均される応力中の第２項は、ポテンシャル、例えば、クーロン力等により生じる圧力の項である。

　そして、平行六面体の周期構造を採用することで、並進ベクトルRを用いて周期構造を仮定した上での応力を求めることができる。また、波動関数Ψは、数３に示したように、原子軌道の和とすることができ、数３の係数は、数４、数５の形式で解くことができる。また、数５は、経験的パラメータを用いて、現実の物質の値と一致するように設定すればよい。したがって、数９によって、平行六面体の応力の平均値を求めることができる。

　以上のように、微小領域を平行六面体に置き換え、無限に繰り返す構造を仮定することで、数９あるいは数１０の応力の平均の式を、平行六面体内の積分に帰着できる。

　以上のように、応力の平均値は、最終的には数１４、数１５の形で表すことができるので、HijやErepをRで微分し、固有値問題の解であるCi、Cjを用いることによって計算できる。

　そして、数９の式で得られる平行六面体内の平均応力を物質の微小領域の応力とする。微小領域は、微小領域の中心座標、例えば、重心（Ｘｉ，Ｙｉ、Ｚｉ）と対応付けて、（Ｘｉ，Ｙｉ、Ｚｉ、Ｔｘｘｉ，Ｔｘｙｉ＝Ｔｙｚｉ，Ｔｘｚｉ＝Ｔｚｘｉ，Ｔｙｙｉ，Ｔｙｚｉ＝Ｔｚｙｉ，Ｔｚｚｉ）の形式で記憶装置のファイルに格納すればよい。したがって、図１に示したＦ３とＦ４とを微小領域の数だけ繰り返すことで、物質中の各微小領域中の応力を算出し、ファイルに格納できる。

　図１のＦ６の処理において、各領域で求めた応力を線形補間することにより局所応力を求める。各微小領域で計算した応力の値は、平行六面体を形成する８原子の重心座標の値とする。そして、線形補間を行うことによって全空間における応力の関数（応力分布）を求める。

　例えば、微小領域Ｖ１において、重心が（Ｘ１，Ｙ１、Ｚ１）であり、Ｘ軸方向の応力Ｔｙｚ１が得られているとする。また、微小領域Ｖ１に隣接する微小領域Ｖ２の重心が（Ｘ１＋Ｄ，Ｙ１、Ｚ１）であり、Ｘ軸方向の応力Ｔｙｚ２が得られているとする。

　この場合に、（Ｘ１，Ｙ１、Ｚ１）と（Ｘ１＋Ｄ，Ｙ１、Ｚ１）との間のＸ軸方向の応力は、Ｔｘｙ（ｄ）＝（（Ｄ－ｄ）Ｔｘｙ１＋ｄＴｘｙ２）／Ｄ；とすればよい。ここで、ｄは、微小領域Ｖ１の重心からの距離である。

　以上述べたように、本情報処理装置は、経験的パラメータを用いた量子論に基づく計算手法において局所応力の計算および計算された応力の分布の計算を可能にする。したがって、求めることができる物理量の範囲を拡張することができる。

　すなわち、本情報処理装置は、物質に含まれる原子の結晶構造を基準に、物質の原子構造を微小領域に分割する。そして、本情報処理装置は、分割によって得られた微小領域を平行六面体に置き換え、平行六面体の繰り返しによる仮想結晶構造を作成する。仮想結晶構造は、均一な物質と考えてよいので、本情報処理装置は、仮想結晶構造中の１つの平行六面体内の積分計算で、応力の平均値を求める。そして、本情報処理装置は、仮想結晶構造で得られた応力の平均値を物質の微小領域の応力とする。以上のような処理を繰り返すことで、本情報処理装置は、物質の原子構造のそれぞれの微小部分の応力を計算する。したがって、本情報処理装置は、物質中の応力分布、すなわち、局所応力を求めることができる。

　以上のような手順において、本情報処理装置は、物質に含まれる原子の結晶構造を基準にして物質を微小部部分に分割するので、微小部分と微小部分の周囲の原子構造との関係を、極力本来の原子構造に類似させた状態で、仮想結晶構造を作成できる。したがって、仮想結晶構造から得られた平均応力は、物質中の微小部分に作用する応力に近いものと考えてよい。本情報処理装置は、以上のような物質本来の物理的性質を極力維持した状態で、仮想結晶構造による平均応力の計算を可能とする。

　そして、本情報処理装置は、計算負荷が比較的少ない経験的パラメータを用いた量子論を用いるので、ナノデバイスの開発において、物質中の欠陥・亀裂の発生を原子レベルで理解するために重要な局所応力を、数千原子以上からなる実用的な規模の系について、実用的な時間で計算できるようになる。

　さらに、以下の三点の特長を持つ。
（１）上記の計算は微小領域ごとに独立して行えるため、プロセッサの並列計算による高速化が容易である。（２）本発明による局所応力計算では、シミュレーションする物質を微小領域に分解するので、原子座標さえ与えられれば、計算負荷の重い経験的パラメータを用いない量子論に基づく計算手法を用いても、微小領域の応力を計算できるという特長を持つ。すなわち、応力計算の前提となる物質の原子構造あるいは原子の位置を求めるシミュレーションに制限はない。例えば、図２の処理以外の処理による原子構造に対しても、本情報処理装置は応力計算可能である。（３）必要な領域だけの局所応力計算が可能であり、要件に応じて計算コストを削減することができる。

　図２８および図２９を参照して、実施例２に係る情報処理装置を説明する。図２８は、実施例２に係る情報処理装置のハードウェア構成を例示する。本情報処理装置は、ネットワーク３０に接続可能なコンピュータである。本情報処理装置は、ＣＰＵ１１、メモリ１２、ハードディスクドライブ１３等の外部記憶装置、ディスプレイ１４、操作部１５、通信部１６、可搬型記憶媒体入出力装置１７を有する。

　ＣＰＵ１１は、メモリ１２に実行可能に展開されたコンピュータプログラムを実行し、情報処理装置の機能を提供する。ＣＰＵ１１は、１個に限定されず、複数のコアを含むものでもよい。

　メモリ１２は、ＣＰＵ１１が実行するコンピュータプログラム、ＣＰＵ１１が処理するデータ等を記憶する。メモリ１１は、不揮発性のＲＯＭ（Read Only Memory）と、揮発性のＤＲＡＭ（Dynamic Random Access Memory）とを含む。

　ハードディスクドライブ１３に駆動されるハードディスクには、メモリ１２に展開されるコンピュータプログラム、あるいは、ＣＰＵ１１が処理するデータ等が格納される。ハードディスクドライブ１３に代えて、フラッシュメモリ等、ＳＳＤ（Solid State Drive）を用いてもよい。ハードディスクドライブ１３、あるいは、ＳＳＤ等の外部記憶装置は、インターフェース１３Ａを通じて、ＣＰＵ２１に接続される。

　インターフェース１３Ａは、例えば、ＵＳＢ（Universal Serial Bus）、ＩＤＥ（Integrated Drive Electronics）、ＳＣＳＩ（Small Computer System Interface）、ＦＣ（Fibre Channel）等のインターフェースである。

　ディスプレイ１４は、例えば、液晶ディスプレイ、エレクトロルミネッセンスパネル等である。ディスプレイ１４は、インターフェース１４Ａを通じて、ＣＰＵ２１に接続される。インターフェース１４Ａは、例えば、ＶＧＡ（Video Graphics Array）等のグラフィックスモジュール、ＤＶＩ（Digital Visual Interface）等のインターフェースである。

　操作部１５は、キーボード、マウス、タッチパネル、静電パッド等の入力装置である。静電パッドは、平面パッドを指等でなぞるユーザ操作を検知し、ユーザ操作に応じてディスプレイ上のカーソルの位置と移動状態とを制御するために使用される装置である。例えば、平面パッド下の電極の静電容量の変化によって、ユーザの指の動きが検知される。操作部１５は、インターフェース１５Ａを通じて、ＣＰＵ１１に接続される。インターフェース１５Ａは、例えば、ＵＳＢのインターフェースである。

　通信部１６は、ＮＩＣ(Network Interface Card)とも呼ばれる。通信部１６は、例えば、ＬＡＮ（Local Area Network）のインターフェースである。通信部１６は、インターフェース１６Ａを通じて、ＣＰＵ１１に接続される。インターフェース１６Ａは、例えば、ＣＰＵ１１の内部バスに接続する拡張スロットである。

　可搬型記憶装置１７は、例えば、ＣＤ（Compact Disc）、ＤＶＤ（Digital Versatile Disk）、ブルーレイディスク、フラッシュメモリカード等の入出力装置である。可搬型記憶装置１７は、インターフェース１７Ａを通じて、ＣＰＵ１１に接続される。インターフェース１７Ａは、例えば、ＵＳＢ、ＳＣＳＩ等のインターフェースである。

　なお、図２８では、情報処理装置として、単独のコンピュータを例示した。しかし、情報処理装置は、複数のコンピュータが違いに連携し、処理を分担して実行する複数のコンピュータであってもよい。

　図２９に、実施例２に係る情報処理装置の機能構成を例示する。本情報処理装置は、ハードディスク等の外部記憶装置上に、結晶構造のデータベース２４および物質の原子構造ファイル２５を有する。結晶構造のデータベース２４が結晶中の原子位置を含む原子構造を記憶する手段に相当する。物質の原子構造ファイル２５が物質中の原子位置を含む原子構造を記憶する手段に相当する。

　ただし、結晶構造のデータベース２４および物質の原子構造ファイル２５の少なくとも一方は、情報処理装置の外部記憶装置以外の装置、例えば、メモリ１２にあってもよい。また、結晶構造のデータベース２４および物質の原子構造ファイル２５の少なくとも一方は、ネットワーク上の記憶装置、例えば、ＳＡＮ(Storage Area Network)にあってもよい。また、結晶構造のデータベース２４および物質の原子構造ファイル２５の少なくとも一方は、ネットワークを通じて、アクセス可能なネットワーク上の他のコンピュータ、例えばデータベースサーバ上にあってもよい。

　結晶構造のデータベース２４のフォーマットは、実施例１の図１５に示したものと同様である。また、物質の原子構造ファイル２５は、実施例１で説明した経験的パラメータを用いた原子構造最適化処理、例えば、最急降下法や共役勾配法で求められる。ここで、物質の原子構造ファイル２５は、例えば、原子種（ATYPEi）、原子のＸ座標（Ｘｉ）、Ｙ座標（Ｙｉ）、Ｚ座標（Ｚｉ）を含むレコードを多数含む。

　そして、情報処理装置は、ＣＰＵ１１が、メモリ１２に実行可能に展開されたプログラムを実行し、制御手段２１、分割手段２２および物質の分割部分平均応力算出手段２３として機能する。

　すなわち、分割手段２２は、物質の原子構造ファイル２５から物質の原子構造中に含まれる原子の位置を読み出す。そして、分割手段２２は、物質を複数に分割した分割部分中の原子位置と原子を含む結晶中の原子位置とを比較する。そして、分割手段２２は、物質を複数に分割した分割部分の原子位置と、当該物質の原子を含む結晶の原子位置とを比較し、分割部分と結晶との間で対応する原子間の相対距離の評価値が最も小さくなるように、分割部分の原子位置と結晶の原子位置とを対応付ける。このような分割部分の原子と、結晶構造中の原子との対応付けの手順は、実施例１において、図８において例示した処理と同様である。そして、分割手段２２は、結晶の単位格子に相当する、物質の分割部分を特定する。ここで、分割部分は、例えば、実施例１でいう物質の微小部分に相当する。

　平行六面体形成手段２３は、分割手段２２によって分割された物質の分割部分を取得する。分割部分は、例えば、六面体である。平行六面体形成手段２３は、分割部分の頂点と平行六面体の頂点との間で相対距離の評価値が最も小さくなるような平行六面体を決定する。この場合の平行六面体形成手段２３の処理手順は、実施例１に示した図１８Ａ－図２４の処理と同様である。

　平均応力算出手段２４は、平行六面体形成手段２３が作成した平行六面体を繰り返した仮想結晶構造を作成する。そして、平均応力算出手段２４は、実施例１に示した式１０と同様の計算式にしたがって、仮想結晶構造の平均応力を計算する。

　制御手段２１は、分割手段２２、平行六面体形成手段２３、および平均応力算出手段２４を用いて、物質の原子構造ファイル２５を処理する。すなわち、物質の原子構造ファイル２５で付与される原子構造を分割部分に分割し、分割部分を平行六面体に変更し、仮想結晶構造を生成し、仮想結晶構造から得られる平均応力を当該分割部分の応力とする。制御手段２１は、分割手段２２、平行六面体形成手段２３、および平均応力算出手段２４を繰り返し実行することによって、物質の原子構造ファイル２５で付与される原子構造で生じる応力分布を作成する。応力分布は、例えば、分割部分の重心（Ｘｉ，Ｙｉ，Ｚｉ）を用いて、（Ｘｉ，Ｙｉ，Ｚｉ，Ｔｘｙｉ，Ｔｙｚｉ，Ｔｚｘｉ）のように表し、メモリ１２に格納することができる。

　なお、均応力算出手段２４で得られる平均応力は、分割手段２２で分割された分割部分に作用する応力、すなわち、分割部分と他の分割部分との境界面に作用するである。一方、分割部分と他の分割部分との境界面以外の位置の応力については、制御手段２１は、境界面からの距離にしたがった線形補間によって求める。このようにして、制御手段２１は、物質の原子構造ファイル２５で付与される物質中の応力分布を出力する。図２９では、省略されているが、応力分布の出力先は、例えば、外部記憶装置中のファイル、ディスプレイ１４Ａ、メモリ１２中の他のアプリケーションプログラムの変数領域等である。
＜コンピュータが読み取り可能な記録媒体＞
　コンピュータその他の機械、装置（以下、コンピュータ等）に上記いずれかの機能を実現させるプログラムをコンピュータ等が読み取り可能な記録媒体に記録することができる。そして、コンピュータ等に、この記録媒体のプログラムを読み込ませて実行させることにより、その機能を提供させることができる。

　ここで、コンピュータ等が読み取り可能な記録媒体とは、データやプログラム等の情報を電気的、磁気的、光学的、機械的、または化学的作用によって蓄積し、コンピュータ等から読み取ることができる記録媒体をいう。このような記録媒体のうちコンピュータ等から取り外し可能なものとしては、例えばフレキシブルディスク、光磁気ディスク、ＣＤ－ＲＯＭ、ＣＤ－Ｒ／Ｗ、ＤＶＤ、ブルーレイディスク、ＤＡＴ、８ｍｍテープ、フラッシュメモリなどのメモリカード等がある。また、コンピュータ等に固定された記録媒体としてハードディスクやＲＯＭ（リードオンリーメモリ）等がある。

　本実施形態で説明した技術を用いることにより、新材料・デバイス開発において、物質中の欠陥・亀裂の発生を原子レベルで解明するができる。そのため、目的にあったデバイスの原子構造的な安定性を、実際に試作する前にシミュレーションを行うことによってある程度予測することができる。これは、開発時間と経費の削減、開発による環境負荷の低減につながる。本実施形態で説明した技術は、新材料・デバイス開発の技術分野、新材料・デバイス開発を支援する情報処理技術の分野に適用できる。

　１１　　ＣＰＵ
　１２　　メモリ
　２１　　制御手段
　２２　　分割手段
　２３　　平行六面体形成手段
　２４　　物質の分割部分の平均応力算出手段
　２４　　結晶構造のデータベース
　２５　　物質の原子構造のファイル

Claims (3)

1. 　物質中の原子位置を含む原子構造を記憶する手段と、
   　前記原子を含む結晶中の原子位置を含む原子構造を記憶する手段と、
   　前記物質を複数に分割した分割部分の原子位置と前記物質の原子を含む結晶の原子位置とを比較し、前記分割部分と前記結晶との間で対応する原子間の相対距離の評価値が最も小さくなるように、前記分割部分の原子位置と前記結晶の原子位置とを対応付け、前記結晶の単位格子に相当する、前記物質の分割部分を特定する分割手段と、
   　前記分割部分の頂点と平行六面体の頂点との間で相対距離の評価値が最も小さくなるような平行六面体を決定する平行六面体形成手段と、
   　前記平行六面体を繰り返した仮想結晶構造において前記平行六面体に加わる平均応力を算出する平均応力算出手段と、
   　前記分割手段、平行六面体形成手段、平均応力算出手段を繰り返し機能させることによって前記物質のそれぞれの分割部分の応力を特定する制御手段と、を備える情報処理装置。
2. 　コンピュータが、物質中の原子位置を含む原子構造を記憶手段から読み出すステップと、
   　前記原子を含む結晶中の原子位置を含む原子構造を記憶手段から読み出すステップと、
   　前記物質を複数に分割した分割部分の原子位置と前記物質の原子を含む結晶の原子位置とを比較し、前記分割部分と前記結晶との間で対応する原子間の相対距離の評価値が最も小さくなるように、前記分割部分の原子位置と前記結晶の原子位置とを対応付け、前記結晶の単位格子に相当する、前記物質の分割部分を特定する分割ステップと、
   　前記分割部分の頂点と平行六面体の頂点との間で相対距離の評価値が最も小さくなるような平行六面体を決定する平行六面体形成ステップと、
   　前記平行六面体を繰り返した仮想結晶構造において前記平行六面体に加わる平均応力を算出する平均応力算出ステップと、
   　前記分割ステップ、平行六面体形成ステップ、および平均応力算出ステップを繰り返し実行することによって前記物質のそれぞれの分割部分の応力を特定する制御ステップと、を実行する情報処理方法。
3. 　コンピュータに、物質中の原子位置を含む原子構造を記憶手段から読み出すステップと、
   　前記原子を含む結晶中の原子位置を含む原子構造を記憶手段から読み出すステップと、
   　前記物質を複数に分割した分割部分の原子位置と前記物質の原子を含む結晶の原子位置とを比較し、前記分割部分と前記結晶との間で対応する原子間の相対距離の評価値が最も小さくなるように、前記分割部分の原子位置と前記結晶の原子位置とを対応付け、前記結晶の単位格子に相当する、前記物質の分割部分を特定する分割ステップと、
   　前記分割部分の頂点と平行六面体の頂点との間で相対距離の評価値が最も小さくなるような平行六面体を決定する平行六面体形成ステップと、
   　前記平行六面体を繰り返した仮想結晶構造において前記平行六面体に加わる平均応力を算出する平均応力算出ステップと、
   　前記分割ステップ、平行六面体形成ステップ、および平均応力算出ステップを繰り返し実行することによって前記物質のそれぞれの分割部分の応力を特定する制御ステップと、を実行させることを特徴とするプログラム。

Images