JP3335881B2 - 材料設計支援システム、およびコンピュータシステムに系の安定構造や動的挙動を原子や分子レベルで解析させこれにより材料の設計支援を行うためのコンピュータプログラムを格納した記憶媒体 - Google Patents

Description

【発明の詳細な説明】

【０００１】

【発明の属する技術分野】本発明は計算機を利用して物
質および材料の設計を行うための支援システム（ＣＡＤ
システム）に係り、特にＭＭ（分子力学）計算、ＭＤ
（分子動力学）計算に関するものである。

【０００２】

【従来の技術】計算機を利用して物質および材料の設計
を行う場合、結晶系の構造安定性や動的挙動等を、ＭＭ
（分子力学）計算やＭＤ（分子動力学）計算を用いて解
析することが一般に知られている。

【０００３】ＭＤ計算とは有限温度での原子や分子の運
動の時間変化をニュ−トン方程式に従って追跡する手法
であり、系の徴視的な挙動を調べる手段として物質・材
料設計における必要不可欠な技術の一つとなっている
（文献［１］：「上田顯，コンピュータシミュレーショ
ン，朝倉書店(1990)」）。

【０００４】また、ＭＭ計算とは、系の最安定構造を解
析するための手法であり、絶対零度（Ｔ＝０Ｋ）でのＭ
Ｄ計算に相当する（文献［２］：「岡田勲，大澤映二
編：分子シミュレーション入門，海文堂(1989)」）。

【０００５】ＭＭ計算及びＭＤ計算は無限系（例えば結
晶）と孤立系（例えば分子）の両者に対して適用可能で
あるが、本発明におけるＭＭ／ＭＤ計算は、無限系を対
象とする。

【０００６】ここで「無限系」とは、原子・分子レベル
で見た時に十分な広がりを持つ系、すなわち結晶・アモ
ルファス・液晶等の固体および液体を指す。このような
無限系を計算機シミュレ−ションで扱う場合、そのまま
では取り扱えないため、無限系内に、ユニットセルと呼
ばれる所定の大きさのセル（単位セル）を設定する。そ
して、無限系を、上記ユニットセルが３次元的に周期的
に連続する系として把握するようにする。

【０００７】このようにすれば、ユニットセルについ
て、周期的境界条件を課して計算を行うだけで、無限系
についての計算機シュミレーションを行うことができ
る。

【０００８】このような周期的境界条件を用いた取り扱
いは、もともと周期性を有する結晶系に対しては妥当な
手法と言える。しかしながら、不純物等の存在で周期性
が消滅してしまった系に対しても近似法として適用がな
されているということがある。

【０００９】周期性を消滅させる原因（乱れ）として
は、まず不純物が挙げられる。ここでいう不純物とは、
置換型・侵入型不純物原子及び欠陥の３種類（以下、こ
の３種類を「不純物原子」と称する）のことであり、１
個以上の不純物原子が結晶中で比較的近傍に存在してい
る集合のことを「不純物」と呼ぶ。また、乱れとして
は、不純物の他にも破壊現象等を挙げることができる
が、この明細書では、このような乱れが比較的低濃度で
分散して存在している系（乱れ同士の相互作用がほとん
どない系）を「非周期系」として定義する。

【００１０】このような「非周期系」は、前記乱れが不
純物である場合には、「孤立した不純物を含む結晶系」
ということもできる。この「孤立した不純物を含む結晶
系」とは、結晶中にこのような「不純物」が唯１つしか
存在していない場合、或いは複数個存在している場合で
も、それら不純物同士が十分離れて位置しており、不純
物間の相互作用がほとんどない状態をいう。

【００１１】このような系は、別の表現によれば「不純
物」の含有率が非常に低い系のことであり、例えば（化
合物）半導体中の欠陥を例として挙げることができる。
これは、ＤＸセンターとして知られ、半導体のエネルギ
ーギャップの中央付近に現れる深い不純物準位を与え、
半導体レーザなどの素子特性に大きな影響を及ぼすもの
である。

【００１２】このような「非周期系」の構造安定性や動
的挙動等を前述のＭＭ計算やＭＤ計算によって解析を行
う場合において、前記「周期的境界条件」を用い１つの
ユニットセルに注目してシュミレーションを行う方法
は、「スーパーセル法」として一般に知られている。

【００１３】ス−パ−セル法は、ユニットセルに乱れを
１つだけ配置し、ユニットセルのサイズを十分大きく設
定することで乱れ同士の相互作用がない状態を近似的に
記述する方法である。すなわち、スーパーセル法では、
前述した非周期系（乱れが「不純物」である場合には
「孤立した不純物を含む結晶系」）の状態を実現するた
めに、通常、ユニットセル中に唯１つの乱れを設定す
る。しかし、周期境界条件によってこの乱れはは周囲の
ユニットセル（これを「イメージセル」という）にも自
動的に存在するので、それらとの相互作用を小さくして
正しい計算を行うためには、ユニットセルのサイズをで
きるだけ大きく設定する必要がある。

【００１４】ただし、スーパーセル法では、シミュレー
ションを開始する前に予めユニットセルのサイズを決め
ておく必要があり（シミュレーション中のサイズ変更は
不可能）、そのサイズは計算精度と計算時間との兼ね合
いで（小さめに）決められる場合が多い。

【００１５】このス−パ−セル法を用いることによっ
て、周期系に対する計算手法と同じ枠組み内で非周期系
も取り扱うことができるため、ス−パ−セル法は、従来
より、ＭＤ計算やＭＭ計算に対する便宜的な手法として
ごく一般的に用いられているのである。

【００１６】ところで、このスーパーセル法を用いてＭ
Ｍ計算やＭＤ計算を行うには、系の各構成原子に働く力
を、原子間相互作用に基づいてＭＭ計算における繰り返
し計算の各ステップやＭＤ計算の各時刻（繰り返し計算
の１ステップ）ごとに計算する必要がある。また同時に
系のポテンシャルエネルギ−も原子間相互作用に基づい
て計算される場合が多い。この原子間相互作用が長距離
相互作用である場合、力やエネルギ−の計算の収束が遅
く、特にｒ⁻¹ 関数であるク−ロン相互作用に関して
は、原子のペアについて単純に総和計算を行ったのでは
精度の良い計算が困難であることが知られている。その
ためク−ロン相互作用による力やエネルギ−の計算を上
述した無限系（周期的境界条件が課された系）に対して
行う場合は、エワルド（Ｅｗａｌｄ）法と呼ばれる手法
を用いることが一般的である（前述の文献［１］）。

【００１７】このエワルド法はク−ロン相互作用だけで
なく、ｒ^−ｎ関数で与えられる相互作用（多重極相互作
用やファン・デル・ワ−ルス力）の格子和（異なるユニ
ットセルに属する等価な原子に関する和）の計算に適用
可能であり、格子和の計算を実空間と虚数空間の二つに
分けてそれぞれの収束をうまく速めることでｒ^−ｎ長距
離相互作用の力やエネルギ−の計算を効率化する優れた
手法である。

【００１８】しかしながら、従来の手法による周期系お
よびス−パ−セル法を適用した非周期系に対するＭＭ／
ＭＤ計算は以下のような問題を抱えている。

【００１９】すなわち、ＭＭ／ＭＤ計算の実行時間のほ
とんどは力（およびエネルギ−）の計算によって費やさ
れる。原子間相互作用が長距離相互作用である場合、考
慮しなければならないペアの数が膨大なものとなり、特
にク−ロン相互作用に関してはエワルド法を適用しても
依然計算時間が深刻な間題として残る。

【００２０】さらに、ＭＤ計算に関しては、統計力学的
に許容される範囲内に温度ゆらぎをおさえるにはユニッ
トセル当たりの原子数Ｎが数百個程度以上の系に対して
ＭＤ計算を行う必要があるため、力に関する計算量を少
しでも減らすことが重要な課題となる。

【００２１】特に物質・材料設計を目的とする場合、何
種類もの系に対して温度等のシミュレ−ション条件をい
ろいろ変えてＭＤ計算を繰り返し行うことが要求される
ため、計算時間の短縮、すなわち力の計算の効率化が計
算機による現実的な新材料設計のための大きな鍵となっ
ている。

【００２２】またＭＤ計算における繰り返し計算の１ス
テップは現実の系での１ピコ秒以下（０．０１ＰＳ程
度）にしか相当しないため、物理的に意味のある範囲で
系の時間変化を解析するには何万から何十万ステップ以
上もの繰り返し計算を行う必要がある。そのため計算機
能力の向上にも拘わらず、１回のＭＤ計算に何十日もの
時間を要するという現状が未だ残されている。

【００２３】以上の計算時間に関する問題点は、周期系
およびス−パ−セル法で取り扱う非周期系の両者に対し
て長距離相互作用を考慮する場合に生じるものである
が、非周期系に対するＭＭ／ＭＤ計算ではさらに次のよ
うな問題が生じる。

【００２４】すなわち、非周期系をス−パ−セル法で取
り扱う場合、上述したように異なるユニットセルに含ま
れる乱れ同士の相互作用を無視し得るよう、十分に大き
なサイズのユニットセルを用いて計算を行う必要があ
る。

【００２５】しかし原子間相互作用が長距離相互作用で
ある場合、計算時間との兼ね合いから十分なサイズのユ
ニットセルで計算を実行することは現実問題として困難
であり、その結果セルサイズに依存した計算誤差が生じ
ることになる。すなわち、通常のシミュレ―ションでは
計算時間を削減するために程々の大きさのユニットセル
を用いて計算を実行することが多いが、その場合、ユニ
ットセルの大きさによっては、周囲のイメージセルに含
まれる乱れとの相互作用が無視できないことになり、そ
れに応じた計算誤差が生じることになるのである。これ
は、いいかえれば、前に定義した「非周期系」という状
況を正しく扱っていないことになるまたユニットセルの
サイズ依存性を確認するには、ユニットセルのサイズを
変えてＭＭ／ＭＤ計算を行い結果の比較を行う必要があ
るが、最も一般的な相互作用であるク−ロン相互作用は
ｒ⁻¹ 関数であるため、例えば１０倍の精度でエネルギ
−の収束を確かめるためには、単純には原子数Ｎが１０
^３倍のユニットセルを用いてＭＭ／ＭＤ計算を行わなけ
ればならない。

【００２６】しかし上述したように計算時間がＮの二乗
に比例して増大するためク−ロン相互作用に関してはユ
ニットセルのサイズ依存性の碓認さえ難しい現状になっ
ている。

【００２７】さらに非周期系でク−ロン相互作用（電
荷）を考慮している場合、ス−パ−セル法では次のよう
な問題が生じる。すなわち、局所的な乱れの有する電荷
の総和が０でない場合、そのような乱れ（例えば置換型
・侵入型不純物イオン等）が導入されることでユニット
セル当たりの電荷の総和自体も０からずれることにな
る。周期的境界条件が仮定されているス−パ−セル法で
は、エネルギ−の発散を防ぐためにユニットセル当たり
の電荷の総和を必ず０にしなければならないという、数
値計算の都合でテクニカルに要求される制約が課され
る。

【００２８】この制約の結果、電気的中性を破る乱れを
含む系をス−パ−セル法で取り扱う場合においてユニッ
トセルの電荷の総和が０でない場合には、これを０にす
るために人為的な電荷（通常一様電荷）を各構成原子に
付加しなければならないということがある。これでは、
系の正しい電荷分布での計算が実現できないことにな
る。また、この人為的な電荷の考慮は、特にポテンシャ
ルエネルギ−の値に大きな影響を及ぼすものである。

【００２９】

【発明が解決しようとする課題】以上をまとめると、従
来の周期系及び非周期系をスーパーセル法を用いて解析
する方法では、次のような問題が提起される。

【００３０】（１）まず、ＭＭ／ＭＤ計算においては、
ＭＭ／ＭＤ計算の実行時間のほとんどは力（およびエネ
ルギ−）の計算によって費やされる。原子間相互作用が
長距離相互作用である場合、考慮しなければならないペ
アの数が膨大なものとなり、特にク−ロン相互作用に関
してはエワルド法を適用しても依然計算時間が深刻な間
題として残る。

【００３１】（２）非周期系の計算をスーパセル法に基
づいて行う場合には、異なるユニットセルに含まれる乱
れ同士の相互作用を無視し得るよう、十分に大きなサイ
ズのユニットセルを用いて計算を行う必要がある。

【００３２】しかし原子間相互作用が長距離相互作用で
ある場合、計算時間との兼ね合いから十分なサイズのユ
ニットセルで計算を実行することは現実問題として困難
であり、その結果セルサイズに依存した計算誤差が生じ
ることになる。

【００３３】またユニットセルのサイズ依存性を確認す
るには、ユニットセルのサイズを変えてＭＭ／ＭＤ計算
を行い結果の比較を行う必要があるが、ペアポテンシャ
ルでは、計算時間が原子数Ｎの二乗に比例して増大する
ため、ｒ^−ｎ長距離相互作用に関してはユニットセルの
サイズ依存性の確認さえ難しい現状になっている。

【００３４】（３）電気的中性が破れた系に対してスー
パーセル法を使った解析を行う場合には、ユニットセル
当たりの総電荷を０にするために電荷分布を補正する必
要があるが、この場合、系の全ポテンシャルエネルギー
・平衡格子定数・平衡原子座標等が変化してしまうとい
う不都合が生じ、電気的中性が破れた系の電荷分布を厳
密に取り扱うことができない。

【００３５】この発明は、このような事情に鑑みてなさ
れたものであり、その目的とするところは、無限系の安
定構造や動的挙動を、系が周期系であると非周期系であ
るとを問わず、効率的に解析することができる材料設計
支援システムを提供することにある。

【００３６】

【課題を解決するための手段】

（１）請求項１に記載された発明 （構成）請求項１に記載された発明は、材料の設計支援
に用いられるシステムであって、局所的な乱れを有する
系（非周期系Ｓ² ）の構造安定性や動的挙動を原子や分
子レベルで解析する材料設計支援システムにおいて、乱
れを有さない系（リファレンス系Ｓ⁰ ）に対する物理量
を処理するリファレンス系Ｓ⁰ 処理部と、前記非周期系
Ｓ² の構造をリファレンス系Ｓ⁰ からのずれで記述する
ために、前記非周期系Ｓ² に対する物理量、および前記
リファレンス系Ｓ⁰ と前記非周期系Ｓ² の構成原子を１
対１に対応付けるためのリストを処理する非周期系Ｓ²
処理部と、原子間または分子間相互作用を記述する関数
形を処理する原子間ポテンシャル処理部と、前記リファ
レンス系Ｓ⁰ のエネルギ−および各原子または各分子に
働く力を計算するリファレンス系Ｓ⁰ 計算部と、前記リ
ファレンス系Ｓ⁰ 計算部で計算されたリファレンス系Ｓ
⁰ のエネルギ−および構成原子に働く力を用いて、非周
期系Ｓ² のエネルギ−および構成原子に働く力を計算す
る非周期系Ｓ² 計算部と、非周期系Ｓ² 計算部の計算結
果を用いて、非周期系Ｓ² に対する安定構造あるいは動
的挙動の計算（ＭＭ計算あるいはＭＤ計算）を実行する
ＭＭ／ＭＤ計算部とを有することを特徴とするものであ
る。

【００３７】（作用）上述の如く構成すれば、非周期系
Ｓ² に対する高速なＭＭ／ＭＤ計算を以下のように実現
することができる。

【００３８】この発明では非周期系Ｓ² に対するＭＭ／
ＭＤ計算を実行する際に、リファレンス系Ｓ⁰ （構成原
子が非周期系Ｓ² の構成原子と１対１に対応した系）を
同時に取り扱う。

【００３９】このリファレンス系Ｓ⁰ は、有限温度の下
で温度効果による格子振動を有する系であっても、絶対
零度の下で温度効果による格子振動のない系のいずれで
あっても良い。なお、前者の場合には、リファレンス系
Ｓ⁰ を、請求項２，４の発明におけるようにリファレン
ス系Ｓ⁰ （格子振動を考慮しない系）と周期系Ｓ１（格
子振動を考慮する系）とに分けて扱うこともできる。

【００４０】この発明では、まず、リファレンス系Ｓ⁰
の構造およびポテンシャルパラメ−タ等の処理をリファ
レンス系Ｓ⁰ 処理部において行い、非周期系Ｓ² の構造
およびポテンシャルパラメ−タ等の処理およびリファレ
ンス系Ｓ⁰ の構成原子との１対１対応付けを非周期系Ｓ
² 処理部において行う。

【００４１】次に、リファレンス系Ｓ⁰ に対する力およ
びエネルギ−の計算がリファレンス系Ｓ⁰ 計算部で行わ
れる。リファレンス系Ｓ⁰ は周期性を有する系であるた
め、構成原子に働く力の計算には後述の式(18)を適用さ
れ、さらにｒ^−ｎ長距離相互作用の格子和に対してはエ
ワルド法による計算の効率化も可能である。

【００４２】このリファレンス系Ｓ⁰ が、温度効果によ
る格子振動を考慮しない系である場合には、ユニットセ
ルＵ⁰ 当たりの原子数Ｎ_{Ｔ= ０}が比較的小さいことか
ら、力の計算量少なくて済むという利点を持っている。

【００４３】なお、ＭＭ計算は、絶対零度（Ｔ＝０Ｋ）
での解析であるため、常にリファレンス系Ｓ⁰ に対する
力の計算量は少なくなる。

【００４４】すなわち、この発明は、非周期系Ｓ² の構
成原子に対する力を求める際に、計算量の少ないリファ
レンス系Ｓ⁰ の構成原子に対する力を用いることで計算
の効率化を図るものである。

【００４５】すなわち、この発明では、非周期系Ｓ² の
構造をリファレンス系Ｓ⁰ からのずれとして記述するこ
とで、従来法（ス−パ−セル法）で扱うような人為的な
周期的境界条件を用いることなく、完全に孤立した乱れ
を含む無限系に対するＭＭ／ＭＤ計算を実現する。

【００４６】これを行うため本発明は、非周期系Ｓ² 処
理部を有する。この非周期系Ｓ² 処理部では、局所的な
乱れによって生じる格子緩和を考慮する領域（実施例で
はＲ_relax と記述）を設定し、この領域内の原子に関す
る座標やポテンシャルパラメ−タ等、およびこれらの原
子とリファレンス系Ｓ⁰ の原子とを１対１に対応付ける
対応リストを処理する。そして非周期系Ｓ² 計算部にお
いて非周期系Ｓ² に対する力（エネルギー）を求める際
に、リファレンス系Ｓ⁰ 計算部で計算されたリファレン
ス系Ｓ⁰ に対する力の計算結果を利用することで計算の
効率化を図り、このようにして求めた力を用いてＭＭ／
ＭＤ計算部においてＭＭ／ＭＤ計算を実行することで、
非周期系Ｓ² に対するス−パ−セル法よりも高精度かつ
高速なＭＭ／ＭＤ計算を実現するものである。

【００４７】（２）請求項２に記載された発明 （構成）請求項２に記載された発明は、請求項１記載の
材料設計支援システムにおいて、前記非周期系Ｓ² の動
的挙動の解析（ＭＤ計算）を有限温度（温度Ｔ＞０Ｋ）
で行う場合に、このシステムは、さらに、前記非周期系
Ｓ² の構造を、局所的な乱れを有さないが温度効果によ
る原子の振動を考慮した系（周期系Ｓ¹ ）からのずれで
記述する手段を有し、この手段は、周期系Ｓ¹ の構造を
リファレンス系Ｓ⁰ （温度効果による原子の振動を考慮
しない系）からのずれで記述するために、この周期系Ｓ
¹ に対する物理量及びリファレンス系Ｓ⁰ と周期系Ｓ¹
の構成原子を１対１に対応付けるためのリストを処理す
る周期系Ｓ¹ 処理部と、前記リファレンス系Ｓ⁰ 計算部
で計算されたリファレンス系Ｓ⁰ のエネルギ及び構成原
子に働く力を用いて、周期系Ｓ¹ のエネルギ及び構成原
子に働く力を計算する周期系Ｓ¹ 計算部と、非周期系Ｓ
² の動的挙動を、周期系Ｓ¹ からのずれとして解析する
ために、周期系ＭＤ計算を実行する周期系ＭＤ計算部と
を有し、前記非周期系ＭＭ／ＭＤ計算部は、前記周期系
ＭＤ計算部及び非周期系Ｓ² 計算部の計算結果を用い
て、非周期系Ｓ² に対する動的挙動の計算（ＭＤ計算）
を実行することを特徴とするものである。

【００４８】（作用）この発明は、請求項１に記載の発
明において、非周期系Ｓ^２の動的挙動の解析（ＭＤ計
算）を有限温度（温度Ｔ＞０Ｋ）で行う場合に、リファ
レンス系Ｓ^０のユニットセルが温度効果による格子振動
のために比較的大きなサイズを有する場合にも、非周期
系Ｓ^２に対する効率的なＭＤ計算を実現しうる手段を提
供するものである。

【００４９】すなわち、請求項２記載の発明では、温度
効果による格子振動のために比較的大きなサイズを有す
るリファレンス系Ｓ^０を周期系Ｓ^１と呼び、リファレン
ス系Ｓ^０は温度効果による格子振動を考慮しない系とみ
なす。そして、周期系Ｓ１の構造をリファレンス系の構
成原子に働く力を用いることで計算の効率を図り、ひい
ては非周期系Ｓ^２に対する高精度かつ高速なＭＤ計算を
実現するものである。

【００５０】リファレンス系Ｓ^０及び周期系Ｓ^１の２つ
の系を取り扱うことによる計算の効率化については、後
述する請求項４に記載された発明の説明において説明す
る。なお、請求項１記載の発明においてＭＤ計算部と称
した計算部を請求項２に記載の発明では、周期系ＭＤ計
算部と非周期系ＭＤ計算部とに分けて取り扱っている。

【００５１】この請求項２記載の発明は、非周期系Ｓ^２
に対するＭＤ計算を行う際、リファレンス系Ｓ^０，周期
系Ｓ^１，非周期系Ｓ^２という３つの系を取り扱うことを
意味するものである。

【００５２】（３）請求項３に記載された発明 （構成）請求項３に記載された発明は、請求項１記載の
材料設計支援システムにおいて、非周期系Ｓ^２処理部
は、非周期系Ｓ^２の格子緩和を考慮する領域をシュミレ
ーション中に変更する手段を有することを特徴とするも
のである。

【００５３】（作用）このような構成によれば、格子緩
和（構造変形）が生じる範囲がＭＭ／ＭＤ計算開始前に
設定された領域よりも大きくなるような場合に関しても
非周期系Ｓ^２の安定構造や動的挙動を正確かつ効率的に
計算することが可能になる。

【００５４】（４）請求項４に記載された発明 （構成）請求項４に記載された発明は、材料の設計支援
に用いられるシステムであって、局所的な乱れは有さな
いが、温度効果による原子の振動が存在する系（周期系
Ｓ¹ ）の原子又は分子の動的挙動を解析する材料設計支
援システムにおいて、温度効果による原子の振動を考慮
しない系（リファレンス系Ｓ⁰ ）に対する物理量を処理
するリファレンス系Ｓ⁰ 処理部と、前記周期系Ｓ¹ の構
造をリファレンス系Ｓ^０からのずれで記述するために、
周期系Ｓ¹ に対する物理量およびリファレンス系Ｓ⁰ と
周期系Ｓ¹ の構成原子を１対１に対応付けるためのリス
トを処理する周期系Ｓ¹ 処理部と、リファレンス系Ｓ⁰
のエネルギ−および各原子または各分子に働く力を計算
するリファレンス系Ｓ⁰ 計算部と、前記リファレンス系
Ｓ^０計算部で計算されたリファレンス系Ｓ^０のエネルギ
−および構成原子に働く力を用いて、周期系Ｓ¹ のエネ
ルギ−および構成原子に働く力を計算する周期系Ｓ¹ 計
算部と、前記周期系Ｓ¹ 計算部の計算結果を用いて、周
期系Ｓ¹ に対する動的挙動の計算（ＭＤ計算）を実行す
る周期系ＭＤ計算部とを有することを特徴とするもので
ある。

【００５５】（作用）上述の如く構成すれば、周期系Ｓ
¹ に対する高速なＭD 計算を以下のように実現すること
ができる。

【００５６】この発明では周期系Ｓ¹ に対するＭＤ計算
を実行する際に、リファレンス系Ｓ⁰ （構成原子が周期
系Ｓ¹ の構成原子と１対１に対応した系）を同時に取り
扱う。

【００５７】このリファレンス系Ｓ⁰ は周期系Ｓ¹ の絶
対零度（Ｔ＝０Ｋ）での構造に相当し、温度効果による
格子振動のない系となっている。そのためリファレンス
系Ｓ⁰ のユニットセルＵ⁰ 当たりの原子数Ｎ_{Ｔ= ０}は、
周期系Ｓ¹ のユニットセルＵ¹ 当たりの原子数Ｎ_{Ｔ> ０}
と比較して通常Ｎ_{Ｔ= ０}《Ｎ_{Ｔ> ０}となっている。な
お、アモルファス構造等、もともとの最小単位胞がそれ
ほど小さくない系では、Ｎ_{Ｔ= ０}〜Ｎ_{Ｔ> ０}という状況
が起こり得る。そのような系に対しては、請求項１にお
いて、リファレンス系Ｓ^０を温度効果による格子振動を
考慮した系と見なした計算方法によって高精度かつ高速
なＭＤ計算の実現が可能である。

【００５８】この発明では、まず、リファレンス系Ｓ⁰
の構造およびポテンシャルパラメ−タ等の処理をリファ
レンス系Ｓ⁰ 処理部において行い、周期系Ｓ¹ の構造お
よびポテンシャルパラメ−タ等の処理およびリファレン
ス系Ｓ⁰ の構成原子との１対１対応付けを周期系Ｓ¹ 処
理部において行う。

【００５９】次に、リファレンス系Ｓ⁰ に対する力およ
びエネルギ−の計算がリファレンス系Ｓ⁰ 計算部で行わ
れる。リファレンス系Ｓ⁰ は周期性を有する系であるた
め、構成原子に働く力の計算には後述の式(69)を適用す
ることができ、さらにｒ^−ｎ長距離相互作用の格子和に
対してはエワルド法による計算の効率化も可能である。

【００６０】このリファレンス系Ｓ⁰ はユニットセルＵ
⁰ 当たりの原子数Ｎ_{Ｔ= ０}が比較的小さいことから、力
の計算量少なくて済むという利点を持っている。

【００６１】すなわち、この発明は、周期系Ｓ¹ の構成
原子に対する力を求める際に、計算量の少ないリファレ
ンス系Ｓ⁰ の構成原子に対する力を用いることで計算の
効率化を図るものである。

【００６２】具体的には、周期系Ｓ¹ に対する力（およ
びエネルギ−）の計算を周期系Ｓ¹計算部で行う際に、
以下のように処理する。すなわち、周期系Ｓ¹ の構成原
子（ｋ０）（基本セルξ= ０に含まれるｋ番目の原子）
に働く力を計算する際、原子（ｋ０）を中心とするある
領域Ｒ_cut を想定し、この領域Ｒ_cut 内の原子（ｋ'ξ'
）との相互作用としては周期系Ｓ¹ における正しい値
を計算し、領域Ｒ_cut外の原子（ｋ' ξ' ）との相互作
用としては、温度効果による格子振動を無視した場合、
すなわちリファレンス系Ｓ⁰ での相互作用を近似的に用
いる。

【００６３】原子間の相互作用は近距離にある原子から
最も強く影響を受けるため、遠距離にある原子との相互
作用のみを近似的に扱うこのような取り扱いは妥当であ
るといえる。

【００６４】原子（ｋ０）と相互作用している原子の数
を領域Ｒ_cut の内部と外部で比較した場合、領域Ｒ_cut
外の原子数の方が圧倒的に多く（正確に言えばそうなる
ように領域Ｒ_cut を設定する）相互作用の計算に時間を
要するため、この領域Ｒ_cut外の原子との相互作用の計
算に計算量の少ないリファレンス系Ｓ⁰ の計算結果を利
用することで周期系Ｓ¹ に対する力の計算の効率化が可
能となっている。

【００６５】本発明では、このようにして周期系Ｓ¹ 計
算部で求めた力を用いて周期系ＭＤ計算部においてＭＤ
計算を実行することで、周期系Ｓ¹ に対する従来法より
も高速なＭＤ計算を実現できる。

【００６６】（５）請求項５に記載された発明 （構成）請求項５に記載された発明は、請求項１あるい
は４記載の材料設計支援システムにおいて、リファレン
ス系Ｓ^０処理部は、リファレンス系Ｓ^０にダミー原子又
はダミー分子を設定する手段を有することを特徴とする
ものである。

【００６７】（作用）このような構成によれば、リファ
レンス系Ｓ^０にダミー原子あるいは分子を設置すること
で、格子緩和（構造変形）が任意の大きさで生じるよう
な場合に関しても非周期系Ｓ² あるいは周期系Ｓ^１の安
定構造や動的挙動を精度良く計算することが可能にな
る。

【００６８】（６）請求項６に記載された発明 請求項７に記載された発明は、コンピュータシステム
に、局所的な乱れを有する結晶系（非周期系Ｓ² ）の安
定構造性や動的挙動を原子や分子レベルで解析させこれ
により材料の設計支援を行うためのコンピュータプログ
ラムを格納した記憶媒体において、乱れを有さない系
（リファレンス系Ｓ⁰ ）に対する物理量を処理するリフ
ァレンス系Ｓ⁰ 処理手順と、前記非周期系Ｓ² の構造を
リファレンス系Ｓ⁰ からのずれで記述するために、前記
非周期系Ｓ² に対する物理量、および前記リファレンス
系Ｓ⁰ と前記非周期系Ｓ² の構成原子を１対１に対応付
けるためのリストを処理する非周期系Ｓ² 処理手順と、
原子間または分子間相互作用を記述する関数形を処理す
る原子間ポテンシャル処理手順と、前記リファレンス系
Ｓ⁰ のエネルギ−および各原子または各分子に働く力を
計算するリファレンス系Ｓ⁰ 計算手順と、前記周期系Ｓ
¹ 計算部で計算された周期系Ｓ¹ のエネルギ−および構
成原子に働く力を用いて、非周期系Ｓ² のエネルギ−お
よび構成原子に働く力を計算する非周期系Ｓ² 計算手順
と、非周期系Ｓ² 計算部の計算結果を用いて、非周期系
Ｓ² に対する安定構造あるいは動的挙動の計算（ＭＭ計
算あるいはＭＤ計算）を実行するＭＭ／ＭＤ計算手順と
を実行するための指令を格納することを特徴とするもの
である。

【００６９】この発明の作用は、請求項１の作用に対応
する。

【００７０】（７）請求項７に記載された発明 請求項７に記載された発明は、請求項６記載の記憶媒体
において、前記非周期系Ｓ² の動的挙動の解析（ＭＤ計
算）を有限温度（温度Ｔ＞０Ｋ）で行う場合に、この記
憶媒体には、さらに、前記非周期系Ｓ² の構造を、局所
的な乱れを有さないが温度効果による原子の振動を考慮
した系（周期系Ｓ¹ ）からのずれで記述する処理手順が
格納され、この処理手順は、周期系Ｓ¹ の構造をリファ
レンス系Ｓ⁰ （温度効果による原子の振動を考慮しない
系）からのずれで記述するために、この周期系Ｓ¹ に対
する物理量及びリファレンス系Ｓ⁰ と周期系Ｓ¹ の構成
原子を１対１に対応付けるためのリストを処理する周期
系Ｓ¹ 処理手順と、前記リファレンス系Ｓ⁰ 計算部で計
算されたリファレンス系Ｓ⁰ のエネルギ及び構成原子に
働く力を用いて、周期系Ｓ¹ のエネルギ及び構成原子に
働く力を計算する周期系Ｓ¹ 計算手順と、非周期系Ｓ²
の動的挙動を、周期系Ｓ¹ からのずれとして解析するた
めに、周期系ＭＤ計算を実行する周期系ＭＤ計算手順と
を有し、前記非周期系ＭＭ／ＭＤ計算手順では、前記周
期系ＭＤ計算手順及び非周期系Ｓ² 計算手順の計算結果
を用いて、非周期系Ｓ² に対する動的挙動の計算（ＭＤ
計算）を実行することを特徴とするものである。

【００７１】この発明の作用は、請求項２の作用に対応
する。

【００７２】（８）請求項８に記載された発明 請求項８に記載された発明は、請求項６記載の記憶媒体
において、非周期系Ｓ² 処理手順は、非周期系Ｓ² の格
子緩和を考慮する領域をシュミレーション中に変更する
手順を含むことを特徴とするものである。

【００７３】この発明の作用は、請求項３の作用に対応
する。

【００７４】（９）請求項９に記載された発明 請求項９に記載された発明は、コンピュータシステム
に、局所的な乱れは有さないが、温度効果による原子の
振動が存在する系（周期系Ｓ¹ ）での原子又は分子の動
的挙動を解析させこれにより材料の設計支援を行うため
のコンピュータプログラムを格納した記憶媒体におい
て、温度効果による原子の振動を考慮しない系（リファ
レンス系Ｓ⁰ ）に対する物理量を処理するリファレンス
系Ｓ⁰ 処理手順と、前記周期系Ｓ¹ の構造をリファレン
ス系Ｓ^０からのずれで記述するために、周期系Ｓ¹ に対
する物理量およびリファレンス系Ｓ⁰ と周期系Ｓ¹ の構
成原子を１対１に対応付けるためのリストを処理する周
期系Ｓ¹ 処理手順と、リファレンス系Ｓ⁰ のエネルギ−
および各原子または各分子に働く力を計算するリファレ
ンス系Ｓ⁰ 計算手順と、前記リファレンス系Ｓ^０計算手
順で計算されたリファレンス系Ｓ^０のエネルギ−および
構成原子に働く力を用いて、周期系Ｓ^１のエネルギ−お
よび構成原子に働く力を計算する周期系ＭＤ計算手順
と、前記周期系Ｓ¹ 計算部の計算結果を用いて、周期系
Ｓ¹ に対する動的挙動の計算（ＭＤ計算）を実行する周
期系ＭＤ計算手順とを有することを特徴とするものであ
る。

【００７５】この発明の作用は、請求項４の作用に対応
する。

【００７６】（１０）請求項１０に記載された発明 請求項１０に記載された発明は、請求項６あるいは９記
載の記憶媒体において、リファレンス系Ｓ^０処理手順
は、リファレンス系Ｓ^０にダミー原子又はダミー分子を
設定する手順を含むことを特徴とするものである。

【００７７】この発明の作用は、請求項５の作用に対応
する。

【００７８】

【発明の実施の形態】以下、本発明の第１〜第７の実施
形態を図面を参照して説明する。

【００７９】（第１の実施形態）非周期系Ｓ² に対する
ＭＭ／ＭＤ計算 まず、この発明の第１の実施形態を図面を参照して説明
する。

【００８０】この実施形態は、請求項１記載の発明に対
応する。この実施形態では、「乱れ」が「不純物」であ
る場合を例にとって説明を行う。また、そのような場
合、非周期系Ｓ² を「孤立不純物系」とも呼ぶ。

【００８１】図１は、本発明の第１の実施形態に係わる
材料設計支援システムの構成を示したブロック図であ
る。なお、このシステム構成は、後述する第２、第３及
び第４の実施形態においても基本的に同様である。

【００８２】本実施形態の材料設計支援システムは、リ
ファレンス系Ｓ⁰ 処理部［１−１］、非周期系Ｓ² 処理
部［１−２］、原子間ポテンシャル処理部［１−３］、
リファレンス系Ｓ⁰ 計算部［１−４］、非周期系Ｓ² 計
算部［１−５］、ＭＭ計算部［１−６］、ＭＤ計算部
［１−７］から構成される。本実施形態では、まず図１
の各ブロックにおける処理動作を説明し、その後フロー
チャートを用いて材料設計支援システムにおける処理の
流れを説明する。

【００８３】本実施形態の計算方法では、計算対象とし
ている孤立した不純物を含む系（非周期系Ｓ² ）の原子
座標を、不純物を一切配置していない系（リファレンス
系Ｓ⁰ ）からのずれで表すことを最大の特徴としている
ので、それを図２を用いて説明する。まず、図２におい
て実線で表示されている原子がリファレンス系Ｓ⁰ を構
成する原子であり、ボックスがリファレンス系Ｓ⁰ のユ
ニットセルを表している。

【００８４】図２には２次元的に９個のユニットセルが
表示されているが、いま中央のユニットセルを基本セル
（ζ＝０）と呼び、この基本セル（ζ＝０）のほぼ中心
に位置している原子を他の元素（黒丸で表示）に置換す
ることを考える。これが置換型不純物原子であり、この
ような不純物の存在によってその周囲に格子緩和（構造
変形）が起こる。この格子緩和は、通常、不純物を中心
とする有限な領域のみで生じ、図２ではその領域を点線
で示した。この領域内部の原子は不純物からの相互作用
によって平衡原子座標が変化するため、例えば点線（灰
色丸）で示した原子のような配置をとる。従って、本実
施形態が計算対象としている「孤立した不純物を含む結
晶系」としての非周期系Ｓ² は、リファレンス系Ｓ⁰ の
ある一部の領域のみの原子座標が異なっているような系
に相当する。

【００８５】このような構造を記述する際、リファレン
ス系Ｓ⁰ の結晶構造とリファレンス系Ｓ⁰ からのずれ
（図２では矢印で示されており、以下このずれを緩和量
とも呼ぶ）によって表すことが非常に有効である。何故
なら、まずリファレンス系Ｓ⁰は不純物を一切含まない
ことからユニットセルを比較的小さく設定し得るため、
リファレンス系Ｓ⁰ の結晶構造の記述に必要な情報量
（ユニットセルに含まれる原子の座標やポテンシャルパ
ラメータ等）が少なくて済み、さらにはポテンシャルエ
ネルギーや力の計算量も少なく抑えることができるから
である。

【００８６】次に、リファレンス系Ｓ⁰ からのずれに関
してであるが、これは非周期系Ｓ²の不純物を中心とす
るある有限な範囲内でのみ存在するものであるため、保
有しておかなければならない情報量（各原子のずれを表
すベクトル）はそれほど多くはない。さらに、非周期系
Ｓ² のポテンシャルエネルギーや力などの計算はこのず
れに基づいて計算されるため、この計算量もさほど多く
はならない。

【００８７】以上のようなメリットを持つ本実施形態に
よる計算方法と比較して、従来法であるスーパーセル法
では、例えば図２のような格子緩和を計算するにはリフ
ァレンス系Ｓ⁰ のユニットセルの９倍（３次元的には２
７倍）以上の大きさのユニットセルを用いて計算を実行
しなければならないことになり、これは上述の（発明が
解決しようとする課題）で説明したような不都合を含む
ことになる。

【００８８】このように本実施形態では、リファレンス
系Ｓ⁰ と非周期系Ｓ² という関係付けられた２つの系を
扱う。この２つの系をそれぞれに関する物理量（原子座
標、ポテンシャルパラメータ、系のポテンシャルエネル
ギー、各原子に働く力など）は明確に区別しておく必要
があるため、以下ではリファレンス系Ｓ⁰ の物理量に関
しては添字０を、非周期系Ｓ² の物理量に関しては添字
２を付加することによってそれらを区別する。

【００８９】まず、図１のリファレンス系Ｓ⁰ 処理部
［１−１］で行う処理を説明する。上述したようにリフ
ァレンス系Ｓ⁰ は不純物を一切ドープしていない系であ
るため、比較的小さいサイズのユニットセルを取り得
る。例えば、リファレンス系Ｓ⁰として完全結晶を仮定
した場合は、結晶の最小単位胞をリファレンス系Ｓ⁰ の
ユニットセルと見做すことができる。また、リファレン
ス系Ｓ⁰ としては完全結晶だけでなく、周期性を有する
ものであればアモルファス構造を仮定してもよい。

【００９０】リファレンス系Ｓ⁰ 処理部［１−１］に
は、リファレンス系Ｓ⁰ の格子定数、ユニットセル当た
りの原子数Ｎ⁰ 、Ｎ⁰ 個の原子の原子座標ｒ⁰ （h
0）_α、およびこれらの原子のポテンシャルパラメータ
Ａ⁰ （h ０）が入力され、それらの値が保有される。

【００９１】周期境界条件の下でユニットセルは全て等
価であるため、ここではζ＝０というユニットセルに着
目し、この基本セル（ζ＝０）に対して原子座標やポテ
ンシャルパラメータを指定している。ｒ⁰(h0) _α は、
基本セル（ζ＝０）に含まれるｈ番目の原子のα方向の
座標を表している。任意のユニットセルに含まれる原子
の座標ｒ⁰(h ζ) _α は、格子定数（基本並進ベクトル
ａ_λα）と原子の相対座標（ｘ_λ(h) ）から、下記の
（１）〜（３）式のように計算される。

【００９２】

【数１】

【００９３】ここで、λは３本の基本並進ベクトルを区
別するための添字である。リファレンス系Ｓ⁰ 処理部
［１−１］では、原子座標ｒ⁰(h0) _α の代わりに相対
座標ｘ_λ(h) を入力し、上記の式に従ってｒ⁰(h0) _α
を設定する方式を採用したのでもよい。

【００９４】次に、ポテンシャルパラメータＡ⁰(h0) は
原子間ポテンシャル処理部［１−３］で保持されている
原子間ポテンシャルのパラメータであり、例えばクーロ
ン相互作用の場合は各原子の電荷に相当している。原子
座標ｒ⁰(h0) _α とポテンシャルパラメータＡ⁰(h0) に
関しては、周期境界条件のため任意のζに対して、下記
の（４）（５）が成り立っている。

【００９５】

【数２】

【００９６】以上のリファレンス系Ｓ⁰ の結晶構造を記
述するための取扱いは、通常のスーパーセル法での取扱
いと全く同じである。ポテンシャルパラメータＡ⁰(h0)
の取扱いは図１の原子間ポテンシャル処理部［１−３］
での処理動作とも関係することであるため、原子間ポテ
ンシャル処理部［１−３］の説明において詳述する。

【００９７】次に、図１の非周期系Ｓ² 処理部［１−
２］での処理を説明する。本実施形態では、上述したリ
ファレンス系Ｓ⁰ のいくつかの原子を不純物原子と置き
換えることによって、非周期系Ｓ² の原子を１対１に対
応付けて取り扱う必要がある。また非周期系Ｓ² では、
リファレンス系Ｓ⁰ が有していた周期性が不純物を局所
的に配置したことによって消滅する（図２参照）。その
ため、周期性を表す（hζ）という原子番号の表記方法
はもはや適当ではなく、非周期系Ｓ² の原子に対しては
通し番号(i) を設定する必要がある。

【００９８】以上の処理は対応リスト［Ｓ⁰ ←→Ｓ² ］
を用いて行うことができる。この対応リスト［Ｓ⁰ ←→
Ｓ² ］の一例を図４に示した。この表は上述した図２に
示した結晶構造に対する対応リスト［Ｓ⁰ ←→Ｓ² ］に
なっている。また、図２のリファレンス系Ｓ⁰ と非周期
系Ｓ² の構成原子に原子番号を割り振ったものが図３で
あり、この図３と対応リスト［Ｓ⁰ ←→Ｓ² ］（図４）
の原子番号が対応しているため、以降は結晶構造として
図３を参照することにする。

【００９９】図４の対応リスト［Ｓ⁰ ←→Ｓ² ］にはユ
ニットセル中に４個（Ｎ⁰ ＝４）の構成原子を含むリフ
ァレンス系Ｓ⁰ の原子番号（h ζ）が登録されており、
それに対応する非周期系Ｓ² の構成原子として２１個の
通し番号(i) が登録されている。この２１個の原子は図
３の点線で囲まれた領域に含まれる原子に相当してい
る。図４の対応リスト［Ｓ⁰ ←→Ｓ² ］によって、非周
期系Ｓ² とリファレンス系Ｓ⁰ の構成原子に１対１の対
応付けが行われていることが分かる。この対応付けは例
えばｉ＝ｈ＋ζＮ⁰ という関係式によって一意的に行う
ことも可能である。以下では説明を分かりやすいものと
するため、非周期系Ｓ² の不純物を通し番号(i) の代わ
りに不純物番号（σ）でも表すこととする。

【０１００】以下、非周期系Ｓ² 処理部［１−２］での
処理動作の流れを主にＭＭ／ＭＤ計算の初期設定部分に
関して説明する。非周期系Ｓ² の処理部［１−２］に
は、非周期系Ｓ² に含まれる不純物原子の個数Ｎ^imp 、
Ｎ^imp 個の不純物原子の番号、Ｎ^imp 個の不純物原子の
ポテンシャルパラメータＡ^imp （σ）、及び格子緩和を
考慮する領域Ｒ_relax が入力される。ここで不純物原子
の番号とは、不純物原子との置き換えを行うリファレン
ス系Ｓ⁰ 中の原子の番号である。即ち、不純物原子をＮ
^imp 個発生させる場合は、｛（ｈ₁,ζ₁ ），（ｈ₂ , ζ
₂ ），…，（ｈ_{Ni mp}，ζ_Nimp）｝という組を指定する。

【０１０１】図３に示された例では不純物原子は１個
（Ｎ^imp ＝１）のみになっており｛（ｈ＝２、ζ＝
０）｝が指定される。Ｎ^imp としては任意の個数の不純
物原子を取り扱うことが可能であり、また指定する不純
物原子はリファレンス系Ｓ⁰ の基本セル（ζ＝０）に含
まれるものでもよいし、イメージセル（ζ≠０）に含ま
れるものでもよい。さらに、この不純物原子を指定する
組み｛（ｈ₁,ζ₁ ），（ｈ₂, ζ₂ ），…，（ｈ_Nimp，
ζ_Nimp）｝は、対応リスト［Ｓ⁰ ←→Ｓ² ］によって不
純物番号の組み｛（σ₁ ），（σ₂ ），…，
（σ_Nimp）｝と対応付けられている。図４の対応リスト
［Ｓ⁰ ←→Ｓ² ］では不純物原子は１個であり、非周期
系Ｓ² の原子（ｉ＝２）に不純物番号（σ＝１）が割り
振られている。

【０１０２】次に、格子緩和を考慮する領域Ｒ_relax に
ついて説明する。非周期系Ｓ² では、不純物が導入され
たことによって生じる格子緩和が不純物の周辺で生じ
る。この格子緩和は不純物より遠方では殆ど０であるた
め、不純物を中心とする有限な領域内のみで格子緩和を
計算すればよい。この格子緩和を考慮する領域として、
例えば不純物原子の重心を中心とする半径Ｒ_relax の球
を仮定することができる。この領域は必ずしも球形であ
る必要はなく、本実施形態では格子緩和を考慮する領域
を一般的にＲ_relax で表す。この領域Ｒ_relax は、ユー
ザによって任意に指定されてもよいし、本実施形態の材
料設計支援システム側で設定する方式でもよい。非周期
系Ｓ² 処理部［１−２］では、領域Ｒ_relax に含まれる
原子数Ｎ²を数え上げる（Ｎ² は、Ｎの二乗ではなく、
本発明では、Ｎの二乗はＮ×Ｎと記述する）。当然なが
ら、このＮ² 個の原子にはＮ^imp 個の不純物原子が全て
含まれている。図３に示した例ではＮ² ＝２１となって
いる。次に非周期系Ｓ² 処理部［１−２］では、Ｎ² 個
の原子の座標ｒ²(i)_α を計算し、その値を保有する。
ここでｒ²(i)_α は下記の（６）式によって計算され
る。

【０１０３】

【数３】

【０１０４】上述したように、非周期系Ｓ² の原子(i)
は対応リスト［Ｓ⁰ ←→Ｓ² ］によってリファレンス系
Ｓ⁰ の原子（h ζ）と１対１に対応付けられている。従
って、式（６）のｒ⁰(i)_α に対しては対応リスト［Ｓ
⁰ ←→Ｓ² ］とリファレンス系Ｓ⁰ の周期性から

【数４】

【０１０５】が成り立っており、その値はリファレンス
系Ｓ⁰ 処理部［１−１］において既に保有されているｒ
⁰(h0) _α から簡単に計算される。

【０１０６】次に、式（６）の第２項であるΔｒ(i) _α
の説明を行う。本実施形態では、非周期系Ｓ² の結晶
構造をリファレンス系Ｓ⁰ からのずれによって記述する
ことを既に述べたが、そのずれがΔｒ(i) _α に相当
し、これは物理的には不純物がドープされたことによる
リファレンス系Ｓ⁰ からの格子緩和を表している。例え
ば、本実施形態の計算方法に基づいて非周期系Ｓ² の安
定構造をＭＭ計算によって求める際は、この緩和量Δｒ
(i) _α を計算する。Δｒ(i) _α にはＭＭ／ＭＤ計算
の初期設定段階で初期値が設定される。この初期値はユ
ーザが任意に指定する方式でもよいし、本実施形態の材
料設計支援システムが一括して０に設定する方式でもよ
い。

【０１０７】次に、非周期系Ｓ² の原子のポテンシャル
パラメータの取扱いを説明する。上述したように非周期
系Ｓ² 処理部［１−２］には、Ｎ^imp 個の不純物原子
（σ）のポテンシャルパラメータＡ^imp （σ）が入力さ
れる。これを受けて非周期系Ｓ² 処理部［１−２］で
は、領域Ｒ_relax に含まれるＮ² 個の原子(i) のポテン
シャルパラメータＡ²(i) を

【数５】

【０１０８】と設定する。このポテンシャルパラメータ
Ａ²(i) の取扱いは、リファレンス系Ｓ⁰ のポテンシャ
ルパラメータＡ⁰(h0) と同様に原子間ポテンシャル処理
部［１−３］での処理動作と関連しているため、以下原
子間ポテンシャル処理部［１−３］の説明を行う。

【０１０９】図１の原子間ポテンシャル処理部［１−
３］では、原子間相互作用を表す関数形を保持する。本
実施形態では原子(i) と原子（ｊ）に働く原子間相互作
用を

【数６】

【０１１０】という関数形によって記述する。ここでＡ
(i) は、リファレンス系Ｓ⁰ 処理部［１−１］や非周期
系Ｓ² 処理部［１−２］で述べたポテンシャルパラメー
タＡ⁰(h0) 及びＡ²(i) に相当している。例えば、φ
（ij）がクーロン相互作用を表す場合、Ａ(i) ＝ｑ
_i （ｑ_i は原子(i) の電荷）、ψ（ｒ｜ij）＝１／｛ｒ
（ij）｝（ｒ（ij）は原子(i) と原子（ｊ）の原子間距
離）となる。

【０１１１】式（10）のような関数形で記述される原子
間ポテンシャルとして、クーロン相互作用の他に、ファ
ン・デル・ワ−ルス力、多重極相互作用等のｒ^−ｎで記
述される原子間相互作用がある。これらの相互作用は通
常のシミュレ―ションにおいて基本的に考慮されるもの
であるため、式（10）は原子間ポテンシャルの一般的な
関数形を表しているといえる。本実施形態の材料設計支
援システムでは、例えばクーロン相互作用＋ファン・デ
ル・ワ−ルス力というように、同時に複数の原子間相互
作用を扱うことができる。その場合、原子(i) と原子
（ｊ）に働く原子間相互作用φ（ij）は

【数７】

【０１１２】で記述される。原子間相互作用を１種類の
み考慮する場合も、複数種類考慮する場合も取扱いは同
じであるため、上述したリファレンス系Ｓ⁰ 処理部［１
−１］、非周期系Ｓ² 処理部［１−２］、及び以下の説
明では添字μを省略する。また、本実施形態では原子間
ポテンシャルに関する説明を行うが、本実施形態による
計算方法では分子性結晶に対する取扱いも可能であり、
その場合は分子間ポテンシャルが原子間ポテンシャルと
同様に処理される。

【０１１３】非周期系Ｓ² 処理部［１−２］の説明での
べたように、リファレンス系Ｓ⁰ と非周期系Ｓ² の原子
は対応リスト［Ｓ⁰ ←→Ｓ² ］によって１対１に対応付
けられており、それらのポテンシャルパラメータは式
（８）及び式（９）によって関係付けられている。式
（８）と式（９）によると、リファレンス系Ｓ⁰ と非周
期系Ｓ² のポテンシャルパラメータを比較した場合、非
周期系Ｓ² の不純物原子はリファレンス系Ｓ⁰ の対応す
る原子とは異なるポテンシャルパラメータの値を持つ
が、それ以外の非周期系Ｓ² の原子はリファレンス系Ｓ
⁰ の対応する原子のポテンシャルパラメータと全く同じ
値を持つ。

【０１１４】さらに本実施形態の材料設計支援システム
では、不純物原子として置換型・侵入型及び欠陥の３種
類を取り扱うことが可能であり、それらはリファレンス
系Ｓ⁰ と非周期系Ｓ² のポテンシャルパラメータの関係
によって、図５のように区別される。図５に示した３種
類の不純物原子のうち、取扱いに特に注意を要するのが
侵入型不純物原子である。侵入型不純物原子はそれに対
応する原子がリファレンス系Ｓ⁰ にもともと存在しな
い。よって対応リスト［Ｓ⁰←→ Ｓ² ］で１対１の対応
付けを行うため、リファレンス系Ｓ⁰ に予めダミー原子
を設定しておく必要がある。このダミー原子のポテンシ
ャルパラメータＡ⁰(h0) は図５に示してあるように常に
０になっているため、ダミー原子が設定されることによ
ってリファレンス系Ｓ⁰ のポテンシャルエネルギーや各
原子に働く力が変化することは一切ない。

【０１１５】本実施形態の材料設計支援システムでは、
非周期系Ｓ² の安定構造や動的挙動をＭＭ計算やＭＤ計
算によって求めるために、リファレンス系Ｓ⁰ 及び非周
期系Ｓ² のポテンシャルエネルギーや各原子に働く力を
計算する。以下この方法を説明する。

【０１１６】まず、図１のリファレンス系Ｓ⁰ 計算部
［１−４］での処理動作を説明する。ここでの計算対象
であるリファレンス系Ｓ⁰ は周期性を有しているため、
その取扱いは従来法であるスーパーセル法における計算
方法と全く同じである。リファレンス系Ｓ⁰ 計算部［１
−４］では、リファレンス系Ｓ⁰ のポテンシャルエネル
ギー及び各原子に働く力が計算される。そのために２原
子間の位置ベクトルｒ⁰(ｈζ, ｈ’ζ')_α を本実施形
態では下記の（13）式のように定義する。

【０１１７】

【数８】

【０１１８】以下で説明を行うポテンシャルエネルギー
や力の計算式の符号は、この式（13）の定義に基づいて
いる。まずポテンシャルエネルギーに関しては、式（1
0）で表される原子間ポテンシャルを用いて原子（h
０）に関して

【数９】

【０１１９】クーロン相互作用を計算する場合、式（1
4）と式（15）に関してはエワルド法（参考文献１）を
用いて計算の効率化を図ることも可能である。リファレ
ンス系Ｓ⁰ のユニットセル当たりのポテンシャルエネル
ギーＥ⁰ _cell は式（14）を用いて

【数１０】

【０１２０】のように求められる。従来法であるスーパ
ーセル法では、Ｅ⁰ _cell 或いはＥ⁰ _{cel l} ／Ｎ⁰ （１原子
当たりのポテンシャルエネルギー）が計算される。リフ
ァレンス系Ｓ⁰ に含まれる全てのユニットセル分だけＥ
⁰ _cell を足し合わせたものが、全ポテンシャルエネルギ
ーＥ⁰ _whole に相当する。次に、リファレンス系Ｓ⁰ の
基本セル（ζ＝０）中の原子（h ０）に働く力のα成分
Ｆ⁰(h0) _α は、式（15）を用いて

【数１１】

【０１２１】によって計算される。従来法であるスーパ
ーセル法では、このＦ⁰(h0) _α を用いてＭＭ計算及び
ＭＤ計算が実行される。

【０１２２】次に、図１の非周期系Ｓ² 計算部［１−
５］の動作処理を説明する。非周期系Ｓ² の安定構造等
を求めるため、リファレンス系Ｓ⁰ からの緩和量Δｒ
(i) _αを計算する点が本実施形態の最も大きな特徴であ
る。このΔｒ(i) _α を求めるために本実施形態では、
非周期系Ｓ² に対してΔＥ 及びΔＦ(i) _α を計算す
る。これら２つの量は

【数１２】

【０１２３】によって定義され、非周期系Ｓ² の物理量
（右辺第１項）からリファレンス系Ｓ⁰ の物理量（右辺
第２項）を差し引いた値になっている。具体的には、Δ
Ｅ はリファレンス系Ｓ⁰ に不純物がドープされたこと
によって生じる格子緩和によるポテンシャルエネルギー
の変化分を表している。Ｅ⁰ _whole もＥ² _whole も系全
体に対するポテンシャルエネルギーであるため発散量で
あるが、その差ΔＥ は有限な値となる。これは局所的
に配置された不純物による格子緩和が有限な範囲のみで
生じていることに対応している。

【０１２４】またＦ⁰(i)_α とＦ²(i)_α は、それぞれ
リファレンス系Ｓ⁰ と非周期系Ｓ²の原子に働く力であ
る。ΔＥ 及びΔＦ(i) _α の定義は上の式で与えられ
るが、実施の計算は以下のように行う。ΔＥ 及びΔＦ
(i) _α はそれぞれ２つの項で

【数１３】

【０１２５】と表すことができる。先ずエネルギーに関
しては、

【数１４】

【０１２６】と表される。次に力は、

【数１５】

【０１２７】で与えられる。ここで、式（23）と式（2
8）のｅ⁰(i) 及びｆ⁰(i)_α に関しては、対応リスト
［Ｓ⁰←→ Ｓ²] とリファレンス系Ｓ⁰ の周期性から

【数１６】

【０１２８】上述した一連の式は、数学的にはテーラー
展開を表している。非周期系Ｓ² の物理量（Ｅ² _whole
及びＦ²(i)_α ）をリファレンス系Ｓ⁰ からのずれ（δ
ｒ（ij）或いはΔｒ(ij)_α）によってテーラー展開し、
得られた第０次項からリファレンス系Ｓ⁰ の物理量（Ｅ
⁰ _whole 或いはＦ²(i)_α ）を差し引いたものが式（2
1）と式（22）の第１項に相当する。そして、テーラー
展開の高次の項が式（21）及び式（22）の第２項に相当
する。この高次の項を与えている式（24）と式（29）で
は２次までの項が書き下されており、Ｏ［（δｒ）³ ］
及びＯ［（Δｒ）³ ］は３次以上の任意の次数までの項
を表している。

【０１２９】また、上記の式（２４）と式（２９）の代
わりに、後述する式（１０６）と式（１０８）におい
て、周期系Ｓ^１に対する量をリファレンス系Ｓ^０に対す
る量に置き換えた式を用いても良い。

【０１３０】次に、図１のＭＭ計算部［１−６］での動
作処理を説明する。このＭＭ計算部［１−６］では、本
実施形態による計算方法に基づいて非周期系Ｓ² の安定
構造をＭＭ計算によって求める。ＭＭ計算とは系のポテ
ンシャルエネルギーを最小化することによって安定構造
を求める方法であり、これを実現する手段としていくつ
かの計算手法が存在するが、本実施形態では最急降下法
を例にとって説明する。本実施形態の材料設計支援シス
テムによると、最急降下法による繰り返し計算でＭＭ計
算を実行する場合、非周期系Ｓ² の領域Ｒ_relax に含ま
れるＮ² 個の原子(i) のｎステップ目のずれΔｒ(i|n)
_α が

【数１７】

【０１３１】によって計算される。ここでΔＦ(i) _α
は、非周期系Ｓ² 計算部［１−５］において式（22）に
従って計算されるものである。またεは定数であり、Ｍ
Ｍ計算を効率良く実行するためにユーザが任意に指定し
てもよいし、本実施形態の材料設計支援システム側で設
定したのでもよい。上述の式に従ってΔｒ(i) _α がＭ
Ｍ計算部［１−６］で計算されると、それを受けて非周
期系Ｓ² 処理部［１−２］においてΔｒ(i) _α の値が
更新されると共に、非周期系Ｓ² の原子座標ｒ²(i)_α
が計算される。以上の処理をΔＦ(i|n) _α がある判定
条件の下で殆ど０になるまで繰り返し計算を行うこと
で、計算目的である非周期系Ｓ² の安定構造を求めるこ
とができる。

【０１３２】本実施形態の計算方法によると、式（37）
によって非周期系Ｓ² の安定構造を求めると同時に、リ
ファレンス系Ｓ⁰ の原子座標ｒ⁰(h0) _α に関しても通
常のＭＭ計算を行うことが可能である。但し、リファレ
ンス系Ｓ⁰ として予め安定構造（平衡構造）を入力して
おく、或いは非周期系Ｓ² に対するＭＭ計算を開始する
前にリファレンス系Ｓ⁰ に対するＭＭ計算を別途実行
し、その安定構造を求めておくことが計算の効率アップ
を図るためには望ましい。このＭＭ計算実行中、非周期
系Ｓ² 計算部［１−５］では、任意のステップにおいて
ΔＥ を式（21）に基づいて計算することもできる。こ
のΔＥ は、上述したように物理的にはリファレンス系
Ｓ⁰ に不純物がドープされたことによって生じる格子緩
和によるポテンシャルエネルギーの変化分である。各ス
テップでΔＥ を計算しその値をモニタした場合、非周
期系Ｓ² のポテンシャルエネルギーがＭＭ計算によって
最小化されていく過程を追跡することが可能である。

【０１３３】次に、図１のＭＤ計算部［１−７］での処
理動作を説明する。ＭＤ計算では系の動的な時間変化を
解析するために、各原子の座標と共に速度の計算も行
う。よってＭＤ計算を実行する場合、リファレンス系Ｓ
⁰ 処理部［１−１］及び非周期系Ｓ² 処理部［１−２］
では、各原子の速度を表すｖ⁰(h0) _α ないしｖ²(i)_α
が追加保有される。また、同時に各原子の質量ｍ⁰(h0)
ないしｍ²(i) も保有される。そのためここで、まずリ
ファレンス系Ｓ⁰ 処理部［１−１］及び非周期系Ｓ² 処
理部［１−２］に対する補足説明を行い、その後ＭＤ計
算部［１−７］での処理動作を説明する。リファレンス
系Ｓ⁰ 処理部［１−１］で保有されるｖ⁰(h0) _α とｍ
⁰(h0) に関しては、対応リスト［Ｓ⁰ ←→Ｓ² ］とリフ
ァレンス系Ｓ⁰ の周期性から

【数１８】

【０１３４】が成り立っている。これらの値は、ＭＤ計
算の初期設定段階においてユーザによって指定されるも
のである。もしくは、ある基準に従って本実施形態の材
料設計支援システム側が設定する方式でもよい。ここ
で、原子の質量はポテンシャルパラメータと同様の取扱
いが可能であり、ポテンシャルパラメータの１種として
処理することもできる。次に非周期系Ｓ² 処理部［１−
２］では、まずｖ²(i)_αに関してはｒ²(i)_α に対する
処理と同様の処置が取られる。即ち、式（６）と同様な

【数１９】

【０１３５】という記述を行い、Δｖ(i) _α に初期値
を設定する。この初期値はユーザが任意に指定してもよ
いし、本実施形態による材料設計支援システム側で一括
して０に設定したのでもよい。次に原子の質量ｍ²(i)
は上述したようにポテンシャルパラメータの１種として
みなすことが可能であるため、Ａ²(i) と同様の処理が
行われる。即ち、非周期系Ｓ² 処理部［１−２］にＮ
^imp 個の不純物原子の質量ｍ^imp （σ）が入力され、そ
れを受けて式（８）と式（９）に従ってＮ² 個の原子
(i) の質量ｍ²(i)が設定される。以上がリファレンス系
Ｓ⁰ 処理部［１−１］及び非周期系Ｓ² 処理部［１−
２］に対する補足説明である。

【０１３６】それではＭＤ計算部［１−７］での処理動
作を説明する。本実施形態ではＭＤ計算のアルゴリズム
としてベルレの方法を例にとって説明する。まずリファ
レンス系Ｓ⁰ に対しては、基本セル（ζ＝０）に含まれ
るＮ⁰ 個の原子（h ０）の時刻ｔでの位置と速度が

【数２０】

【０１３７】によって計算され、リファレンス系Ｓ⁰ 処
理部［１−１］のそれらの値が更新される。ここで、Ｆ
⁰(h0) _α はリファレンス系Ｓ⁰ 計算部［１−４］で式
（18）に従って計算される原子（h ０）のα方向の力で
ある。またΔｔはＭＤ計算における時間刻みであり、こ
の値はユーザが任意に指定してもよいし、本実施形態の
材料設計支援システム側で設定してもよい。式（41）と
式（42）によるリファレンス系Ｓ⁰ に対するＭＤ計算
は、従来法であるスーパーセル法と全く同じ取扱いにな
っている。

【０１３８】次に非周期系Ｓ² に対しては、領域Ｒ
_relax に含まれるＮ² 個の原子(i) の時刻ｔでのΔｒ
(i) _α 及びΔｖ(i) _α が

【数２１】

【０１３９】によって計算される。これを受けて非周期
系Ｓ² 処理部［１−２］でそれらの値が更新されると共
に、ｒ²(i)_α 及びｖ²(i)_α が式（６）或いは式（4
0）に従って計算される。ここでΔF'(i) _α は、式（2
8）と式（29）を若干変更した

【数２２】

【０１４０】によって非周期系Ｓ² 計算部［１−５］で
計算されるものである。式（46）及び式（47）において
ｍ⁰(i) ＝１、ｍ²(i) ＝１と置くと、それぞれ式（2
8）及び式（29）と等価となるため、非周期系Ｓ² 計算
部［１−５］には式（22）の代わりに式（45）を用意し
ておき、ＭＭ計算に関してはｍ⁰(i) ＝１、ｍ²(i) ＝
１として計算を実行したのでもよい。

【０１４１】通常のＭＤ計算実行中は、一般に系の全エ
ネルギー（運動エネルギー＋ポテンシャルエネルギー）
をモニタしエネルギーの保存則を確認するが、本実施形
態の計算方法においても、任意の時刻ｔにおいてポテン
シャルエネルギーを計算することが可能である。その
際、非周期系Ｓ² に関してはＥ⁰ _cell を、リファレンス
系Ｓ⁰ に関してはΔＥ を計算する。ＭＤ計算部［１−
７］では、式（41）〜（44）に基づく繰り返し計算を時
刻ｔ_end まで実行し、非周期系Ｓ² の動的な時間変化を
解析する。時刻ｔ_end はユーザが任意に指定できる値で
ある。

【０１４２】以上、本実施形態の材料設計支援システム
の構成を示したブロック図（図１）における処理動作を
各ブロック毎に説明してきた。次に、ブロック間の接続
状態を明確にするため、まずＭＭ計算を行う場合につい
て図６に示したフローチャートを用いて説明する。初め
に、リファレンス系Ｓ⁰ 処理部［１−１］にリファレン
ス系Ｓ⁰ に関する物理量が入力される（ステップ［２−
１］）。この物理量は、リファレンス系Ｓ⁰ の格子定
数、ユニットセル当たりの原子数Ｎ⁰ 、Ｎ⁰ 個の原子の
座標ｒ⁰(h0) _α 及びポテンシャルパラメータＡ⁰(h0)
である。

【０１４３】次に、非周期系Ｓ² 処理部［１−２］に非
周期系Ｓ² に関する物理量が入力される（ステップ［２
−２］）。この物理量は、非周期系Ｓ² に含まれる不純
物原子の個数Ｎ^imp 、Ｎ^imp 個の不純物原子（σ）の番
号、不純物原子（σ）のポテンシャルパラメータＡ^imp
（σ）、及び格子緩和を考慮する領域Ｒ_relax である。

【０１４４】これらの物理量が入力されると、非周期系
Ｓ² 処理部［１−２］において対応リスト［Ｓ⁰ ←→Ｓ
² ］が作成され、非周期系Ｓ² の領域Ｒ_relax に含まれ
るＮ² 個の原子(i) の座標ｒ²(i)_α 及びポテンシャル
パラメータＡ²(i)が設定され、Δｒ(i) _α の初期値が
設定される（ステップ［２−３］）。

【０１４５】次いで、リファレンス系Ｓ⁰ 計算部［１−
４］において、リファレンス系Ｓ⁰に対する部分ポテン
シャルエネルギーｅ⁰(h0) 及び部分力ｆ⁰(h0) _α が計
算される（ステップ［２−４］）。

【０１４６】以上までのステップが初期設定に当たり、
これ以降のステップがＭＭ計算に従う繰り返し計算の処
理動作に当たる。まず、非周期系Ｓ² の計算部［１−
５］において、ΔＥ 及びΔＦ(i) _α が計算され（ス
テップ［２−５］）、このΔＦ(i) _α を用いてＭＭ計
算部［１−６］で緩和量Δｒ(i) _α が計算される（ス
テップ［２−６］）。これを受けて、非周期系Ｓ² 処理
部［１−２］でΔｒ(i)_α の値が更新され、非周期系
Ｓ² の原子座標ｒ²(i)_α が計算される（ステップ［２
−７］）。

【０１４７】以上ステップ［２−５］からステップ［２
−７］に関して、ある判定条件の下でΔＦ(i) _α が十
分小さくなったと判定されるまで繰り返し計算を実行す
る。そして判定条件が満たされる場合、ｒ²(i)_α を非
周期系Ｓ² の安定構造（緩和構造）として出力または表
示する（ステップ［２−８］）。上述したステップ［２
−４］及びステップ［２−５］では、ポテンシャルエネ
ルギーｅ⁰(h0) ないしΔＥ を必ずしも計算する必要は
なく、ユーザの目的に合わせて任意に計算を実行するこ
とが可能である。

【０１４８】次に、ＭＤ計算に関して図７に示したフロ
ーチャートを用いて説明する。まず図７のステップ［３
−１］からステップ［３−３］までは、先に図６で説明
したＭＭ計算に関する処理と全く同じ処理を実行する。
但しこのＭＤ計算に関する処理では、ステップ［３−
１］でリファレンス系Ｓ⁰ の原子の速度ｖ⁰(h0) _α と
質量ｍ⁰(h0) とが追加入力される。また、ステップ［３
−２］では不純物原子の質量ｍ^imp （σ）も入力され、
ステップ［３−３］では上述したステップ［２−３］で
の処理に加えて、ｖ²(i)_α とｍ²(i) が設定される。

【０１４９】ＭＤ計算では、ステップ［３−４］以降が
繰り返し計算処理に相当する。まず、リファレンス系Ｓ
⁰ 計算部［１−４］でｅ⁰(h0) 及びｆ⁰(h0) _α を計算
し（ステップ［３−４］）、非周期系Ｓ² 計算部［１−
５］でΔＥ 及びΔＦ(i) _αを計算する（ステップ［３
−５］）。次いでＭＤ計算部［１−７］において、リフ
ァレンス系Ｓ⁰ に関してはｒ⁰(h0) _α とｖ⁰(h0) _α
を、非周期系Ｓ² に関してはΔｒ(i) _α とΔｖ(i) _α
とを計算する（ステップ［３−６］）。そして、リフ
ァレンス系Ｓ⁰ 処理部［１−１］と非周期系Ｓ² 処理部
［１−２］においてそれぞれｒ⁰(h0) _α とΔｒ(i) _α
の値が更新され、非周期系Ｓ² 処理部［１−２］でｒ
²(i)_α が計算される（ステップ［３−７］）。

【０１５０】以上のステップ［３−４］からステップ
［３−７］までの処理を時刻ｔ_end まで繰り返し実行
し、時刻ｔ_end での非周期系Ｓ² の構造を出力或いは表
示する（ステップ［３−８］）。本実施形態によるＭＤ
計算では、繰り返し計算実行中の任意の時刻において、
リファレンス系Ｓ⁰ に関してｒ⁰(h0) _α 及びｖ⁰(h0)
_αを、非周期系Ｓ² に関してｒ²(i)_α 及びｖ²(i)_α
を出力或いは表示することも可能である。また、ステッ
プ［３−４］とステップ［３−５］ではｅ⁰(h0) ないし
ΔＥ を必ずしも毎回計算する必要はなく、ユーザの目
的に合わせて任意の時刻で計算することが可能である。

【０１５１】以上までが本実施形態による材料設計支援
システムの骨格をなす処理である。

【０１５２】以下、本実施形態における材料設計支援シ
ステムに関する変形例を第２、第３及び第４の実施形態
として説明する。

【０１５３】（第２の実施形態）この第２の実施形態
は、請求項３記載の発明に対応した実施形態である。

【０１５４】まず、非周期系Ｓ² 処理部［１−２］で設
定されている格子緩和を考慮する領域Ｒ_relax に関する
取扱いを説明する。

【０１５５】第１の実施形態の材料設計支援システムで
は、不純物が領域Ｒ_relax のほぼ中心に位置するように
この領域を設定する。そのため、緩和量Δｒ(i) _α は
通常、領域Ｒ_relax の中心部分で大きな値を持ち、領域
Ｒ_relax の端では殆ど０になる。ＭＭ計算及びＭＤ計算
では、繰り返し計算の各ステップないしは各時刻におい
て非周期系Ｓ² 処理部［１−２］で緩和量Δｒ(i) _α
が更新される。緩和量Δｒ(i) _α がシミュレ―ション
実行途中に領域Ｒ_relax の端でもある程度の大きさの値
を持つようになった場合、それは実際に格子緩和が起こ
る範囲が領域Ｒ_relax よりも大きく、そのままシミュレ
―ションを続けたのでは正しい計算結果が得られないこ
とを意味している。

【０１５６】これを避けるため、本実施形態の計算方法
ではシミュレ―ションの途中で領域Ｒ_relax をより大き
なサイズのものに設定変更することが可能である。この
変更は以下のようにして実現される。まずＭＭ計算に関
しては、フローチャートが図６から図８に変更され、ス
テップ［２−９］とステップ［２−１０］が新たに追加
される。図８では、既に説明したようにステップ［２−
７］で緩和量Δｒ(i)_α の値が繰り返し計算の各ステ
ップ毎に更新される。ステップ［２−９］では、領域Ｒ
_relax の端での緩和量Δｒ(i) _α の値をある基準の下
で判定する。そして、領域Ｒ_relax の端での緩和量Δｒ
(i) _α が基準値Δｒ^max _edge よりも大きくなっている
と判定された場合、ステップ［２−１０］に進む。ステ
ップ［２−１０］では、まず非周期系Ｓ² 処理部［１−
２］で領域Ｒ_relax が再設定される。

【０１５７】この再設定は、例えば領域の半径Ｒ_relax
を半径Ｒ´_relax ＝Ｒ_relax ＋Ｒ_{ad d} に置き換えること
によって実行することができる。この新しく設定された
領域Ｒ´_relax には、それまでに取り扱っていたＮ² 個
の原子を含むＮ´² （＞Ｎ²）個の原子が含まれる。そ
のため、続いてステップ［２−１０］では新しく数え上
げられた（Ｎ´² −Ｎ² ）個の原子に関する情報を対応
リスト［Ｓ⁰ ←→Ｓ²］に追加登録すると共に、これら
の原子のΔｒ(i) _α の初期値を設定する。以上の処理
をステップ［２−１０］で行った後、通常の繰り返し計
算に戻る。上述した基準値Δｒ^max _edge 及びＲ_add は、
計算の効率化が計られるようにユーザが任意に指定して
もよいし、本実施形態の材料設計支援システム側で予め
設定しておいてもよい。

【０１５８】以上のような手続きによって領域Ｒ_relax
を適宜変更することで、常に正しい計算結果（緩和構
造）を得ることが可能である。次にＭＤ計算に関して
は、フローチャートが図７から図９に変更される。図９
において追加されているステップ［３−９］とステップ
［３−１０］では、上述したステップ［２−９］及びス
テップ［２−１０］と全く同じ処理が行われる。

【０１５９】（第３の実施形態）この第３の実施形態
は、請求項５記載の発明に対応した実施形態である。

【０１６０】ここでは、緩和量Δｒ(i) _α がシミュレ
―ション実行途中で領域Ｒ_relax の中心部分でも大きな
値を持つようになった場合の処理を第３の実施形態とし
て説明する。

【０１６１】第１の実施形態の材料設計支援システムで
は、非周期系Ｓ² の原子座標をリファレンス系Ｓ⁰ から
のずれΔｒ(i) _α で記述し、非周期系Ｓ² の物理量
（ポテンシャルエネルギー及び原子に働く力）をΔｒ
(i) _α に関するテーラー展開によって近似的に求める
ことを特徴としている。

【０１６２】一般に、テーラー展開を行う場合、展開変
数が小さい場合（展開点の近傍）は近似の精度は保証さ
れているが、展開変数が大きくなった場合（展開点から
離れた点）の計算は精度が悪くなることが知られてい
る。第１の実施形態の材料設計支援システムでは、リフ
ァレンス系Ｓ⁰ が展開点に相当しており、展開点からの
ずれ、即ち緩和量Δｒ(i) _α が大きくなった場合、テ
ーラー展開を用いて計算している物理量（ΔＥ 及びΔ
Ｆ(i) _α ）の計算精度が次第に悪くなる。

【０１６３】これを回避するため本実施形態では、以下
のような手続きによって常に精度のよい計算を実行する
ことができる。まず、ＭＭ計算に関してはフローチャー
トが図６から図１０に変更される。図１０のステップ
［２−１１］では、領域Ｒ_{rela x} に含まれるＮ² 個の原
子に関するΔｒ(i) _α の値が基準値Δｒ^max _whole よ
りも大きくなっているかどうかを判定し、大きい場合は
ステップ［２−１２］に進む。ここで基準値Δｒ^max
_whole はユーザが任意に指定してもよいし、本実施形
態の材料設計支援システム側で予め設定しておいてもよ
い。以下、Ｎ_big （＜Ｎ² ）個の原子の緩和量Δｒ(i)
_α が基準値Δｒ^max _whole よりも大きくなっていると
判定された場合について説明する。また、このＮ_big 個
の原子を（ｉ_{bi g} ）で表すことにする。そして、これら
の非周期系Ｓ² の原子を対応リスト［Ｓ⁰ ←→Ｓ² ］で
対応付けられているリファレンス系Ｓ⁰ の原子を（ｈ
_big ，ζ_{bi g} ）で表す。

【０１６４】以下の処理は、ステップ［２−１１］での
判定の結果、ステップ［２−１２］が初めて実行される
場合に関してのものである。ステップ［２−１２］では
リファレンス系Ｓ⁰ 処理部［１−１］においてリファレ
ンス系Ｓ⁰ に対する追加・変更処理を行う。まず、リフ
ァレンス系Ｓ⁰ の基本セル（ζ＝０）にＮ_big 個のダミ
ー原子（ｈ_dmy ０）（以下、（ｈ_dmy ；ζ＝０）と記述
する）を追加設定する。この処理により、リファレンス
系Ｓ⁰ のユニットセル当たりの原子数はＮ⁰ ＋Ｎ_big と
なる。そして、新たに設定したリファレンス系Ｓ⁰ のダ
ミー原子（ｈ_{dm y} ；ζ＝０）のポテンシャルパラメータ
及び原子位置を

【数２３】

【０１６５】と設定する。

【０１６６】次にステップ［２−１３］では、非周期系
Ｓ² 処理部［１−２］において非周期系Ｓ² に対する追
加・変更処理を行う。ステップ［２−１２］においてリ
ファレンス系Ｓ⁰ にダミー原子が配置されたことによ
り、必然的に非周期系Ｓ² に含まれる原子数も増えるこ
とになる。そのため、まず非周期系Ｓ² の領域Ｒ_relax
に含まれる原子の数え上げを再度行い、対応リスト［Ｓ
⁰ ←→Ｓ² ］にリファレンス系Ｓ⁰ のダミー原子（ｈ
_dmy ，ζ）とそれに対応する非周期系Ｓ² のダミー原
子（ｉ_dmy ）との対応リストを書き加える。対応リスト
［Ｓ⁰ ←→Ｓ² ］に書き加えるペアの数をＮ_add とする
と、Ｎ_add ＞Ｎ_big となり、具体的なＮ_addの値はリフ
ァレンス系Ｓ⁰ のユニットセルと領域Ｒ_relax の大きさ
との兼ね合いから決まる。

【０１６７】また、ステップ［２−１３］以降では非周
期系Ｓ² の領域Ｒ_relax に含まれる原子数はＮ² ＋Ｎ
_add となり、もともとの不純物原子（σ）に加えて原子
（ｉ_{bi g} ）とダミー原子（ｉ_dmy ）もこのステップ以降
は不純物原子として取り扱われる。そのために、続いて
原子（ｉ_big ）とダミー原子（ｉ_dmy ）のポテンシャル
パラメータの書き換えないしは新規設定を行う。まず原
子（ｉ_big ）に関しては、

【数２４】

【０１６８】と変更する。次いで、ダミー原子
（ｉ_dmy ）のポテンシャルパラメータの新規設定を行う
が、それに関しては以下のことを注意しなければならな
い。

【０１６９】上述のステップ［２−１２］でリファレン
ス系Ｓ⁰ の基本セル（ζ＝０）に配置されたダミー原子
（ｈ_dmy ；ζ＝０）は、リファレンス系Ｓ⁰ の周期性の
ためにイメージセル（ζ≠０）にも必然的に存在するこ
とになる。よって、ステップ［２−１３］で対応リスト
［Ｓ⁰ ←→Ｓ² ］に追加されるダミー原子として以下の
２種類のものが存在する。これを非周期系Ｓ² のダミー
原子（ｉ_dmy ）に関して説明すると、一つは、リファレ
ンス系Ｓ⁰ の基本セル（ζ＝０）に含まれているダミー
原子（ｈ_dmy ；ζ＝０）と対応リスト［Ｓ⁰ ←→Ｓ² ］
で対応付けられているダミー原子（ｉ^（ζ=0） _dmy ）で
あり、もう一つは、リファレンス系Ｓ⁰のイメージセル
（ζ≠０）に含まれているダミー原子（ｈ_{dmy ；}ζ≠
０）と対応リスト［Ｓ⁰ ←→Ｓ² ］で対応付けられてい
るダミー原子（ｉ^{（ζ≠０）} _dmy ）である。この２種類
の非周期系Ｓ² のダミー原子に関しては、それらのポテ
ンシャルパラメータを

【数２５】

【０１７０】と設定する。この２つの式が意味している
ところは、非周期系Ｓ² で本来存在していた原子（ｉ
_big ）の物理量をステップ［２−１３］以降ではダミー
原子（ｉ^{( ζ=0)} _dmy）の物理量として計算するために、
まず式（51）という設定を行う。そして、不必要である
にも拘らずリファレンス系Ｓ⁰ の周期性のために必然的
に設定されてしまうダミー原子（ｉ^{（ζ≠０）} _dmy ）に
関しては、式（52）によってそのポテンシャルパラメー
タを０とすることで物理的には存在していないことと同
義の状況を作り出している。次にこれらの原子の緩和量
に関しては、Δｒ_{（ ｉdmy ）α}の初期値を０と置く。Δ
ｒ（ｉ_big ）_αに関しては、ステップ［２−１３］以降
では不必要な量となるため値の設定変更は特に必要な
い。

【０１７１】以上の一連の処理は、テーラー展開の展開
点をずらすことで展開変数が小さくなるように設定し直
す処理に相当している。具体的には、原子（ｉ_big ）の
緩和量Δｒ（ｉ_big ）_αが大きくなったため、展開点を
ｒ⁰ （ｉ_big ）_αからｒ⁰ （ｉ_big ）_α＋Δｒ
（ｉ_big ）_αに取り直すために、リファレンス系Ｓ⁰ の
この位置（ｒ⁰ （ｈ_dmy ；ζ＝０）_α＝ｒ⁰ （ｉ_big ）
_α＋Δｒ（ｉ_big ）_α）にダミー原子（ｈ_dmy ；ζ＝
０）を仮想的に配置する。このダミー原子のポテンシャ
ルパラメータＡ⁰ （ｈ_dmy ；ζ＝０）は式（48）によっ
て０に設定されているため、ダミー原子（ｈ_dmy ；ζ＝
０）を加えたことによってリファレンス系Ｓ⁰ のポテン
シャルエネルギーや力が変化することは全くない。リフ
ァレンス系Ｓ⁰ に対して以上の処理を行った後、次いで
非周期系Ｓ² に対しては、ｒ² （ｉ_big ）_αの代わりに
ｒ² （ｉ^（ζ=0） _dmy ）_αで原子（ｉ_big ）の原子座標
を表すための処理を行う。即ち、原子（ｉ_big ）の座標
をｒ² （ｉ^（ζ=0） _dmy ）_αを用いて

【数２６】

【０１７２】と表し、ステップ［２−１２］以降のＭＭ
計算では原子（ｉ_big ）の緩和量としてΔｒ（ｉ_big ）
_αの代わりにΔｒ（ｉ^（ζ=0） _dmy ）_αを求める。

【０１７３】ダミー原子（ｉ_dmy ）を配置することによ
って非周期系Ｓ² に含まれる見かけ上の原子数は増えて
いるが、原子（ｉ_big ）の本来のポテンシャルパラメー
タであるＡ² （ｉ_big ）が式（50）によって０に変更さ
れ、かつダミー原子（ｉ^{（ζ ≠０）,} _dmy）のポテンシャ
ルパラメータも式（52）によって０に設定されているた
め、非周期系Ｓ² に含まれている物理的な原子数はステ
ップ［２−１３］前後で全く変化していない。

【０１７４】また、ポテンシャルエネルギーや力の計算
の際に考慮しなければならない原子の数（ペアの数）が
ダミー原子の設定によって急増してしまうように一見思
われるが、ダミー原子のほとんどはそのポテンシャルパ
ラメータＡが０になっているため、リファレンス系Ｓ⁰
計算部［１−４］や非周期系Ｓ² 計算部［１−５］で計
算される総和計算の際、ダミー原子に関する項は初めか
ら取り扱わなくてもよい場合が多い。従って、ダミー原
子が設定されたことによるポテンシャルエネルギーや力
の計算量の増加分はそれほど大きくはならない。

【０１７５】上述したステップ［２−１２］ではダミー
原子の設定によってリファレンス系Ｓ⁰ の結晶構造が見
かけ上変更されているため、ステップ［２−１３］終了
後はステップ［２−４］へ一旦戻ってｅ⁰(h0) とｆ⁰(h
0) _α の計算を再度実行し、その後通常のＭＭ計算に
よる繰り返し計算を引き続き実行する。

【０１７６】以上のステップ［２−１１］からステップ
［２−１３］までの処理を行うことによって、非常に大
きな格子緩和が生じるような非周期系Ｓ² に対しても十
分な計算精度でのＭＭ計算を実行することが可能とな
り、本実施形態による材料設計支援システムの汎用性が
格段にアップされるものである。

【０１７７】ここで、以下のことを注意しなければなら
ない。上述したステップ［２−１２］とステップ［２−
１３］での処理はこの２つのステップが初めて実行され
る場合の処理であり、ステップ［２−１２］とステップ
［２−１３］の２回目以降の実行ではその処理が以下の
ように変更される。ステップ［２−１１］において同じ
原子の緩和量Δｒ(i) _α が基準値Δｒ^max _whole より
も大きくなっていると２回以上判定された場合、この原
子に対応するダミー原子は既にリファレンス系Ｓ⁰ に設
定されている。よって、このような原子に対する２回目
以降のステップ［２−１２］及びステップ［２−１３］
ではダミー原子の設定及びポテンシャルパラメータの設
定・変更を行う必要はなく、式（49）のみを実行する。
即ち展開点の位置ｒ⁰ （ｈ_dmy ；ζ＝０）_αだけを再設
定し、Δｒ_{（ｉdmy ）α}の初期値を再度０に設定するだ
けでよい。

【０１７８】以上と同様な処理をＭＤ計算に対しても行
うことが可能であり、その場合はフローチャートが図７
から図１１に変更される。図１１のステップ［３−１
１］からステップ［３−１３］までの処理は、上述した
ステップ［２−１１］からステップ［２−１３］までの
処理と全く同じになっている。

【０１７９】この図１０と図１１に示したフローチャー
トは任意の大きさの格子緩和が生じる非周期系Ｓ² に対
する高精度なＭＭ／ＭＤ計算を実現するための手続きで
あるが、これと同様の機能を得るために初期設定段階に
おいて予めダミー原子を用意しておく方式を採用するこ
ともできる。即ち、非周期系Ｓ² の領域Ｒ_relax に含ま
れるＮ² 個の原子(i) 全てに対するダミー原子をもとも
とリファレンス系Ｓ⁰と非周期系Ｓ² に設定しておくこ
とで、図１０のステップ［２−１２］とステップ［２−
１３］（或いは図１１のステップ［３−１２］とステッ
プ［３−１３］）でダミー原子の新規設定という処理を
省略することも可能である。

【０１８０】以上ＭＭ計算とＭＤ計算に関する第１の実
施形態の変形例を、第２、第３の実施形態として説明し
てきた。ＭＭ計算に関する基本的なフローチャートは図
６に示されており、図８と図１０のフローチャートが変
形例となっている。説明を分かりやすいものとするため
にこの２つの変形例を別々に取り扱ったが、両方の機能
を備えた材料設計支援システムの構築も可能である。同
様にＭＤ計算に関しても、基本的なフローチャート図７
に対する変形例である図９と図１１の機能を兼ね備えた
材料設計支援システムの構築が可能である。

【０１８１】（第４の実施形態）以下、第４の実施形態
について説明する。

【０１８２】この実施形態においては、式（10）で与え
られる以外の関数形によって記述される原子間相互作用
の取り扱いについて説明する。以下説明を行う処理は、
上述してきた第１，第２，第３の実施形態全てに対して
適用が可能である。原子間の短距離反発力としては一般
にBorn-Mayerタイプの関数形

【数２７】

【０１８３】が用いられることが多い。ここで、Ｂ(i)
はポテンシャルパラメータであり、原子(i) の有効半径
に相当している。このようなポテンシャルパラメータＢ
(i) を含む関数形を、以下では

【数２８】

【０１８４】で表す。式（10）の関数ψはポテンシャル
パラメータを一切含んでいないのに対し、式（56）では
関数ψ＾がポテンシャルパラメータＢを含んでいること
を特徴としている。式（56）のような関数形を持つ原子
間ポテンシャルを本実施形態による材料設計支援システ
ムで取り扱う場合、以下のような処理を行う。

【０１８５】ポテンシャルパラメータＡに対する記述と
同様に、リファレンス系Ｓ⁰ に対するポテンシャルパラ
メータＢをＢ⁰(h0) と書き、非周期系Ｓ² に対するポテ
ンシャルパラメータＢをＢ²(i) と書く。そして、Ｂ
²(i) を

【数２９】

【０１８６】で表す。ここでＢ⁰(i)に対しては、対応リ
ストおよびリファレンス系Ｓ⁰ の周期性から、

【数３０】

【０１８７】が成り立っている。式（57）を用いてＢ
²(i) を表した場合、非周期系Ｓ² の不純物原子（σ）
のみが有限な値のΔＢ（σ）（≠０）を持ち、それ以外
の原子のΔＢ(i) は全て０になっている。

【０１８８】次に、非周期系Ｓ² 計算部［１−５］にお
いて、非周期系Ｓ² に対する物理量ΔＥ 及びΔＦ(i)
_α が以下のように計算される。まず、この２つの量は
それぞれ式（21）と式（22）の代わりに

【数３１】

【０１８９】という３つの項によって表される。それぞ
れの項は、先ずエネルギーに関して、

【数３２】

【０１９０】と表され、力に関しては

【数３３】

【０１９１】と表せる。ここで、式（61）（62）（65）
（66）は、それぞれ式（23）（24）（28）（29）と全く
同義であるが、式（61）と式（65）ではリファレンス系
Ｓ⁰ に対するポテンシャルパラメータＢ⁰ が、式（62）
と式（66）では非周期系Ｓ² に対するポテンシャルパラ
メータＢ² が計算に用いられている。

【０１９２】式（56）で与えられる関数形で表される原
子間ポテンシャルの取り扱いにおいて重要な点は、式
（59）と式（60）の第３項がΔＢに関するテーラー展開
の高次の項に相当していることである。この高次の項を
与えている式（63）及び式（67）では１次の項が書き下
されており、Ｏ［（ΔＢ）² ］は２次以上の任意の項ま
での高次項を表している。

【０１９３】以上のような方法によって非周期系Ｓ² に
対する物理量ΔＥ 及びΔＦ(i) _αを計算することで、
本実施形態による計算方法では式（56）で与えられる一
般的な関数形で表される原子間ポテンシャルを用いたＭ
Ｍ／ＭＤ計算を行うことが可能になり、Born-Mayerタイ
プやそれ以外の原子間相互作用を取り扱うことができ
る。

【０１９４】また、式（10）では記述されない関数形で
表される原子間相互作用に関しては、通常の方法によっ
てポテンシャルエネルギー及び力を求め、それらを本発
明によって計算したポテンシャルエネルギー或いは力と
足し合わせてＭＭ／ＭＤ計算を実行することも可能であ
る。

【０１９５】以上説明してきた第１〜第４の実施形態の
材料設計支援システムによると、従来法であるスーパー
セル法と比較して以下のような効果を得ることが可能で
ある。

【０１９６】（１）計算時間の短縮 （２）計算労力の削減 （３）完全な孤立形に対する計算を実現 （４）電気的中性が破れている系に対する補正なしの計
算を実現 以上４つの効果をＭＭ計算を例にとり、実際の計算結果
に基づいて具体的に説明していく。

【０１９７】まず、計算時間の短縮に関して説明する。
説明を分かりやすくするため、簡単な結晶系として代表
的なイオン結晶であるＮａＣｌを例にとって説明する。
ＮａＣｌはｆｃｃ構造を形成し、プラス電荷を持つＮａ
原子とマイナス電荷を持つＣｌ原子が交互に並んだ構造
になっている。いま、ＮａＣｌの完全結晶をリファレン
ス系Ｓ⁰ とし、リファレンス系Ｓ⁰ のＮａ原子が１個だ
けＣｌ原子と置き換わった非周期系Ｓ² を考える。この
ような非周期系Ｓ² に対し、本実施形態による計算方法
及び従来法であるスーパーセル法に基づいてＭＭ計算で
安定構造（緩和構造）を計算した結果が図１２である。
◇が本実施形態の計算方法による計算結果を、□がスー
パーセル法による計算結果を示す。

【０１９８】この計算例では、スーパーセル法でのユニ
ットセルを５１２原子に取っている（実線で表示）。計
算精度に関してこれとの比較を行うため、本実施形態に
よる計算方法では格子緩和を考慮する領域Ｒ_relax とし
てスーパーセル法のユニットセルの体積とほぼ等しい球
を設定している（点線で表示）。また、本実施形態によ
る計算方法でのリファレンス系Ｓ⁰ のユニットセルは、
立方晶としての最小単位胞である８原子からなるユニッ
トセルを考えている。図１２において中心に位置してい
る原子（黒く塗り潰した原子）がＮａからＣｌに置換さ
れた置換型不純物原子である。この図を見て分かる通り
置換型不純物原子を中心として格子緩和が発生してお
り、本実施形態による計算方法と従来法であるスーパー
セル法とで同じ計算結果が得られている。

【０１９９】ここで注意しなければならないのは、両者
の比較はスーパーセル法の基本セルに関してのみ行わな
ければならないと言うことである。何故なら、スーパー
セル法では周期性によってイメージセルにも不純物が含
まれているため、スーパーセル法による計算結果では図
１２に合計９つの不純物原子とその周囲での格子緩和が
表示されている。それに対し本実施形態による計算方法
では、完全に孤立した不純物系を扱うことが可能である
ため、結晶の全領域に渡って正しい計算結果を得ること
ができている。但し図１２に示した例では、格子緩和の
領域がそれほど大きくないためスーパーセル法で用いた
５１２原子のユニットセルは十分な大きさを持ってお
り、基本セルに関しては本発明の計算方法による計算結
果と同じ結果を得ることができている。

【０２００】基本セルに関して両者の結果が同じであっ
たと言うことは、本実施形態による計算方法ではテーラ
ー展開による近似計算を行っているにも拘らず、十分な
計算精度を達成し得ることを意味している。さらに重要
な点は、本実施形態による計算方法に従うと、従来法で
あるスーパーセル法と比較して約８０倍の計算時間の効
率化を図れる点である。このような大幅な計算時間の短
縮が可能である理由は以下の通りである。

【０２０１】ＭＭ／ＭＤ計算における計算時間は、主に
ポテンシャルエネルギーや力を計算する際に考慮しなけ
ればならない原子のペアの数で決まる。以下、力に関し
て説明を行うと、スーパーセル法では式（18）によって
各原子に働く力を計算する。式（18）では式（15）を用
いてｆ⁰(h0) _α が計算されるが、この計算ではユニッ
トセルに含まれる原子とそれらが相互作用を持つ結晶内
の全ての原子とのペアについて原子間ポテンシャルの一
次微分を計算しなければならない。そのペアの数はユニ
ットセルが大きくなった場合、飛躍的に増加する。その
ため、スーパーセル法では『孤立した不純物を含む結晶
系』を実現するためにユニットセルのサイズを大きくす
ると、計算時間が格段にアップしてしまう。

【０２０２】これを回避するため本実施形態では、長い
計算時間を要する式（15）の計算を小さいユニットセル
（例えば、結晶の最小単位胞）に対してのみ実行するこ
とで計算時間の短縮化を可能にしている。これを具体的
に実現しているのが式（22）によるΔＦ(i) _α の計算
である。この計算は、式（28）と式（29）で与えられる
ことを既に説明した。

【０２０３】まず、式（28）でｆ⁰(h0) _α を計算する
が、これはリファレンス系Ｓ⁰ に対する物理量であり、
リファレンス系Ｓ⁰ のユニットセルが小さいために計算
量はそれほど多くならない。

【０２０４】

【数３４】

【０２０５】以上のような理由から、従来法であるスー
パーセル法と比較して本実施形態の計算方法では考慮し
なければならないペアの数が格段に減っており、計算時
間の短縮が可能となっている。また、式（28）のｆ⁰(h
0) _α の計算では、クーロン相互作用に関してはエワ
ルド法（参考文献１）を用いることができ、従来法で用
いられている効率的な手法をそのまま応用することも可
能である。さらに、ｆ⁰(i)_α はリファレンス系Ｓ⁰ に
対する物理量であるので、ＭＭ計算においてリファレン
ス系Ｓ⁰ として平衡構造を用いる場合、初期設定段階に
おいて１度計算しておけばそれ以降のステップでは計算
が不必要となるため、さらに計算時間の短縮を図ること
ができる。

【０２０６】以上の力の計算に関する説明は、ポテンシ
ャルエネルギーΔＥ の計算に対しても同様に該当する
ものであり、以上のような理由によって本発明による計
算方法では従来法であるスーパーセル法と比較して格段
の計算時間短縮を図ることができる。

【０２０７】次に、本実施形態の２番目の効果である計
算労力の削減について説明する。これを説明するために
再度ＮａＣｌを例にとり、今回は２個のＮａ原子をＣｌ
原子と置き換えた系を想定する。この系の安定構造を本
実施形態による計算方法とスーパーセル法によってＭＭ
計算でシミュレ―トした結果が図１３である（黒く塗り
潰した２個の原子が不純物原子）。

【０２０８】この図１３では先に説明した図１２とは異
なり、スーパーセル法の基本セルに関しても特に基本セ
ルの端のところで２つの計算結果が一致しておらず、隣
接したイメージセルでは２つの計算結果のずれがさらに
目立っている。これは図１３の例では格子緩和の起きる
範囲が大きいため、スーパーセル法において仮定したユ
ニットセル（構成原子５１２個）が十分なサイズを持っ
ておらず、スーパーセル法で得られた計算結果は正しい
安定構造（緩和構造）を与えていないことを意味してい
る。スーパーセル法で得られた計算結果が正しいか否か
を確認するためには、ユニットセルのサイズを変えて計
算を何回か行う必要性がある。そのため、できるだけ妥
当なユニットセルのサイズを予測して設定し、さらに確
認のために計算を行うというように、スーパーセル法で
は計算労力だけでなく計算時間も余分に必要となる。

【０２０９】これに対し本実施形態による計算方法で
は、図８及び図９で与えられるフローチャートで説明し
たように、格子緩和を考慮する領域Ｒ_relax をシミュレ
―ション実行中に適宜変更することが可能である。その
ため、実際の格子緩和よりも小さいサイズの領域Ｒ
_relax を初めに設定してしまった場合でも、唯１回のシ
ミュレ―ションで正しい計算結果（ＭＭ計算の場合は緩
和構造）を得ることが可能であり、ユーザの労力及び計
算時間をも削減することができる。実際、図１３に示し
た本実施形態の計算方法による計算例では、格子緩和を
考慮する領域Ｒ_relaxが図１２に示した例よりも大きく
なっている。

【０２１０】次に、本実施形態による計算方法の３番目
の効果について説明する。従来法であるスーパーセル法
では、基本セルに配置された不純物が周期境界条件のた
めに周囲のイメージセルにも必然的に存在し、ユニット
セルを十分大きくとらない限り『孤立した不純物』とい
う状態を正しく実現することはできない。通常、これは
計算時間との兼ね合いから困難であり、近似的な系（不
純物同士の相互作用を完全には無視できない系）につい
て計算を行わざるを得ない場合が多い。この近似が悪か
った場合の例が、図１３に示されているスーパーセル法
による計算結果である。

【０２１１】これに対して本実施形態による計算方法で
は、計算対象である非周期系Ｓ² は周期性を持っておら
ず、スーパーセル法のようにイメージセルに含まれる
『不純物』との相互作用というものが非周期系Ｓ² では
もともと存在しないため、常に『完全に孤立した不純
物』に対するシミュレ―ションが可能になっている。

【０２１２】最後に、本実施形態の４番目の効果を説明
する。スーパーセル法では、不純物原子の設定によって
ユニットセル内の電荷の総和が０からずれた場合、０と
なるように不純物原子自体或いは周囲の原子に対して電
荷の補正を行わなければならないため、そのような場合
は正しい電荷分布を扱うことができないという問題点を
持っている。例えば、図１３に示したスーパーセル法に
よる計算では、２個のＮａ原子がＣｌ原子と置換されて
いるためにユニットセル当たりの電荷の総和が−４にな
っている。これを補正するために図１３の計算では、ユ
ニットセル内の全ての原子に対して電荷＋４／５１２が
一様に加算されている。このように電荷分布が補正され
た場合、系のポテンシャルエネルギーが変わる、平衡格
子定数や平衡原子座標が変わる等の不都合が生じる。

【０２１３】これに対し、本実施形態が計算対象として
いる非周期系Ｓ² では、この電気的中性を必ずしも満た
す必要があるわけではないため、図１２や図１３のよう
な電気的中性が破れている系に対しても電荷の補正を一
切行うことなく計算を実行することが可能である。

【０２１４】なお、本発明は上述した各実施形態に限定
されるものではない。実施形態に記載した各部の構成
は、ハードウェアで実現するのみならずソフトウェアで
実現することもできる。さらに、実施形態に記載した手
法は、コンピュータに実行させることのできるプログラ
ムとして、磁気ディスク（フロッピディスク，ハードデ
ィスクなど）、光ディスク（ＣＤ−ＲＯＭ，ＤＶＤな
ど）、半導体メモリなどの記録媒体に格納して頒布する
こともできる。その他、本発明の要旨を逸脱しない範囲
で、種々変形して実施することができる。

【０２１５】（第５の実施形態）周期系に対するＭＤ計
算 次に、本発明に係る第５の実施形態を図面を参照して説
明する。

【０２１６】この第５の実施形態は、請求項４に記載さ
れた発明に対応する。

【０２１７】まず、この実施形態を説明する前に、周期
的境界条件を課した無限系に対してＭＤ計算を行う場合
の力の計算について簡単に説明を行っておく。

【０２１８】原子（ｈζ）と原子（ｈ' ζ' ）に働く２
体の原子間相互作用をφ（ｈζ, ｈ' ζ' ）と書くと
（（ｈζ）はζ番目のユニットセルに含まれるｈ番目の
原子を表す）、着目しているあるユニットセル（ζ= ０
とする）に含まれる原子（ｈ０）に働く力のα成分は

【数３５】

【０２１９】

【数３６】

【０２２０】格子和の計算に対しては上述したように通
常エワルド法が適用され、その計算量を具体的な数値と
して表すことが難しいためここではＮ_all で表したが、
ユニットセル当たりの原子数Ｎと比較した場合、常にＮ
_all の方が大きな値（Ｎ_all＞Ｎ）＞となる。さらにＭ
Ｄ計算では系の時間変化を追うため通常数千から数万ス
テップ以上の繰り返し計算を行う必要があるので、力に
関するト−タルの計算量はＮ×Ｎ×Ｎ_all ×Ｎ_stepとな
る（係数３は省略）。

【０２２１】この計算量はユニットセル当たりの原子数
Ｎの二乗に比例するため、セルサイズを大きくすると計
算時問が飛躍的に長くなる。ＭＤ計算（スーパーセル
法）で通常用いられるセルサイズは原子数が数百から数
千個程度が主であるが、最近では超並列計算機による数
百万個以上の系に対する計算も行われている。このユニ
ットセル当たりの原子数Ｎは計算対象としている系や現
象に依存して決まるものであるが、周期系と非周期系に
対してそれぞれのＮの下限値に対する目安を次のように
示すことができる。

【０２２２】・周期系…統計力学的に許容される精度の
温度ゆらぎを与えるＮ_{Ｔ> ０} ・非周期系…周期的境界条件のために各ユニットセルに
人為的に配置された乱れ同士の相互作用を無視し得るセ
ルサイズを与えるＮ_SC まず周期系に関して説明すると、ＭＤ計算は原子・分子
の運動を有限温度（Ｔ＞０Ｋ）で解析するため、各原子
は温度効果による徴視的な振動を起こす。逆に言えばこ
の振動が系の温度を与えるため、各原子の振動（速度）
を用いてＭＤ計算の各時刻（繰り返し計算の１ステッ
プ）ごとに系の温度を計算することができる。この温度
の時間変化のゆらぎは、系の独立な原子の数（周期的境
界条件を課している場合はユニットセル当たりの原子数
Ｎ）に依存（１／ EMBED Equation.2 に比例）して決
まるものである。

【０２２３】Ｎが小さい場合、ある指定された温度Ｔを
実現するためには各原子が極端に激しく振動しなければ
ならず、その結果温度のゆらぎが１／ EMBED Equation.
2に比例して大きくなるため、このようなＭＤ計算は物
理的に正常でない系を取り扱っていることになる。

【０２２４】そのため周期系に対しては、統計力学的に
許容される精度での小さい温度ゆらぎを与え得る独立な
原子数Ｎ_{Ｔ> ０}以上のセルサイズを持つ系でのＭＤ計算
の実行が要求される。

【０２２５】次に非周期系に関しては、非周期系をス−
パ−セル法によって近似的に取り扱う場合、周期的境界
条件のために各ユニットセルに対して仮想的な乱れが配
置されることになる。そのため、現実の非周期系（乱れ
同士の相互作用のない系）を表すためにはユニットセル
のサイズを十分大きく設定しなければならない。この十
分なセルサイズを与えるＮ_SCは長距離相互作用を持つ系
においてより大きな値となる。

【０２２６】以上のＭＤ計算手法をふまえ、以下、第５
の実施形態を説明する。

【０２２７】図１４は周期系Ｓ¹ に対するＭＤ計算を実
行するための材料設計支援システムの構成を示したブロ
ック図であり、リファレンス系Ｓ⁰ 処理部［４−１］、
周期系Ｓ¹ 処理部［４−２］、原子間ポテンシャル処理
部［４−４］、リファレンス系Ｓ⁰ 計算部［４−５］、
周期系Ｓ¹ 計算部［４−６］、周期系ＭＤ計算部［４−
８］から構成される。第５の実施形態ではまず図１４の
各ブロックにおける処理動作を説明し、その後フロ−チ
ャ−トを用いて本発明による周期系Ｓ¹ に対するＭＤ計
算を実行するための処理の流れを説明する。

【０２２８】まず、図１４のリファレンス系Ｓ⁰ 処理部
［４−１］で行なう処理を説明する。リファレンス系Ｓ
⁰ 処理部［４−１］では、リファレンス系Ｓ⁰ の結晶構
造、ポテンシャルパラメ−タ等が保有・処理される。リ
ファレンス系Ｓ⁰ は周期系Ｓ¹ の絶対零度（Ｔ＝０Ｋ）
での構造に相当し、周期性を有する系である。またリフ
ァレンス系Ｓ⁰ の周期性の単位（ユニットセルＵ⁰ ）は
例えば結晶の最小単位胞で与えられ、後述するようにこ
のユニットセルＵ⁰ が周期系Ｓ¹ のユニットセルＵ¹ の
構成単位となる。

【０２２９】リファレンス系Ｓ⁰ 処理部［４−１］に
は、リファレンス系Ｓ⁰ の格子定数、ユニットセルＵ⁰
当たりの原子数Ｎ⁰ 、Ｎ⁰ 個の原子の原子座標ｒ⁰ （ｈ
０）_α、ポテンシャルパラメ−タＡ⁰ （ｈ０）、質量ｍ
⁰ （ｈ０）が入力され、それらの値が保有される。

【０２３０】ユニットセルＵ⁰ 当たりの原子数Ｎ⁰ は、
前述したＮ_{Ｔ= ０}に相当するものである。また添字（ｈ
ζ）はリファレンス系Ｓ⁰ のζ番目のユニットセルのｈ
番目の原子を表す。周期的境界条件の下でユニットセル
は全て等価であるため、ここでは０番目（ζ= ０）のユ
ニットセルに着目し、この基本セル（ζ= ０）に対して
原子座標やポテンシャルパラメ−タ等を指定している。

【０２３１】次に、ｒ⁰ （ｈ０）_αは基本セル（ζ=
０）に含まれるｈ番目の原子の座標のα成分を表す。任
意のユニットセルに含まれる原子の座標ｒ⁰ （ｈζ）_α
は、格子定数（基本並進ベクトルＡ_λα）と原子の相対
座標x _λ（ｈ）から

【数３７】

【０２３２】によって計算される。ここでλは３本の基
本並進ベクトルを区別するための添宇である。

【０２３３】リファレンス系Ｓ⁰ 処理部［４−１］で
は、原子座標ｒ⁰ （ｈ０）_αの代わりに基本並進ベクト
ルＡ_λαと相対座標x _λ（ｈ）を入力し、上記の式に従
ってｒ⁰ （ｈ０）_αを設定する方式を採用しても良い。
原子座標ｒ⁰ （ｈ０）_αポテンシャルパラメ−タＡ⁰
（ｈ０）、質量ｍ⁰ （ｈ０）に関しては、周期的境界条
件のため任意のζに対して

【数３８】

【０２３４】が成り立っている。ポテンシャルパラメ−
タＡ⁰ の取り扱いは図１４の原子間ポテンシャル処理部
［４−４］での処理動作とも関係するため、原子間ポテ
ンシャル処理部［４−４］の説明において詳述する。

【０２３５】次に、図１４の周期系Ｓ¹ 処理部［４−
２］について説明する。この周期系Ｓ¹ 処理部［４−
２］では、周期系Ｓ¹ の構造、ポテンシャルパラメ−タ
等が保有・処理される。周期系Ｓ¹ はリファレンス系Ｓ
⁰ の各原子に対して温度効果による徴視的な振動を考慮
した系であり、周期系Ｓ¹ のユニットセルＵ¹ はリファ
レンス系Ｓ⁰ のＭ個分のユニットセルＵ⁰ で構成され
る。

【０２３６】前述したように、このユニットセルＵ¹ に
含まれる原子数（Ｎ_{Ｔ> ０}）に相当）が少ないほどＭＤ
計算に要する計算時間は少なくて済むが、系の温度ゆら
ぎが大きくなるために統計力学的に要求される精度での
ＭＤ計算が実現できなくなる。そのためＭＤ計算では、
通常、ユニットセル当たりの原子数が数百個以上の系を
用いて計算が行われる。

【０２３７】周期系Ｓ¹ 処理部［４−２］には、周期系
Ｓ¹ の格子定数、ユニットセルＵ¹当たりの原子数Ｎ
¹ 、Ｎ¹ 個の原子の原子座標ｒ¹ （ｋ０）_α、速度ｖ¹
（ｋ０）_α、ポテンシャルパラメ−タＡ¹ （ｋ０）、質
量ｍ¹ （ｋ０）が入力され、それらの値が保有される。
ユニットセルＵ¹ 当たりの原子数Ｎ¹ は、前述したＮ_{Ｔ
> ０}に相当するものである。

【０２３８】このＮ¹ を直接指定する代わりに、ユ−ザ
−が要求する精度の温度ゆらぎΔＴを指定し、それに基
づいてＮ¹ を定める方式でもよい。

【０２３９】また添宇（ｋξ）は周期系Ｓ¹ のξ番目の
ユニットセルのｋ番目の原子を表す。周期的境界条件の
下でユニットセルは全て等価であるため、リファレンス
系Ｓ⁰ 処理部［４−１］における処理と同様に０番目
（ξ= ０）のユニットセルに着目し、この基本セル（ξ
= ０）に対して原子座標やポテンシャルパラメ−タ等を
指定している。

【０２４０】上述したように周期系Ｓ¹ のユニットセル
Ｕ¹ はリファレンス系Ｓ⁰ のＭ（=Ｍ_１×Ｍ_２×Ｍ_３）
個のユニットセルＵ⁰ で構成され、周期系Ｓ¹ の原子は
リファレンス系Ｓ⁰ の原子と１対１に対応付けられてい
る。この対応付けは、周期系Ｓ¹ の基本セル（ξ= ０）
の原子（ｋ０）とリファレンス系Ｓ⁰ の原子（ｈζ）に
関して例えばｋ= ｈ+ ζＮ⁰ （ｋ= １, …, Ｎ¹ ）とい
う関係式を用いて行うことができ、この情報が対応リス
ト［Ｓ⁰ ←→Ｓ¹ ］に保存される。

【０２４１】図１６にリファレンス系Ｓ⁰ （構成原子お
よびユニットセルを点線で表示）と周期系Ｓ¹ （構成原
子およびユニットセルを実線で表示）の対応関係の一例
を示した。図１６では二次元的な表記になっており、リ
ファレンス系Ｓ⁰ の９個（Ｍ= ９）のユニットセルＵ⁰
で周期系Ｓ¹ のユニットセルＵ¹ が構成されている場合
の例となっている。

【０２４２】また図１８に図１６に対する対応リスト
［Ｓ⁰ １←→Ｓ¹ ］を示す。図１６に示した例では、リ
ファレンス系Ｓ⁰ のユニットセルＵ⁰ 当たりの原子数Ｎ
⁰=４、周期系Ｓ¹ のユニットセルＵ¹ 当たりの原子数Ｎ
¹=３６となっている。周期系Ｓ¹ 処理部［４−２］で
は、周期系Ｓ¹ の格子定数、ユニットセルＵ¹ 当たりの
原子数Ｎ¹ 、Ｎ¹ 個の原子の原子座標ｒ¹ （ｋ０）_α、
ポテンシャルパラメ−タＡ¹ （ｋ０）、質量ｍ¹ （ｋ
０）の値をユ−ザ−が直接入力する代わりに、（Ｍ_１×
Ｍ_２×Ｍ_３）という情報を入力する方式でもよい。

【０２４３】そして、この情報に基づいてユニットセル
Ｕ¹ 当たりの原子数Ｎ¹ をＮ⁰ ×Ｍと計算し、対応リス
ト［Ｓ⁰ ←→Ｓ¹ ］を作成する。その際上述の関係式ｋ
= ｈ+ ζＮ⁰ において、例えばζ= （ζ_１, ζ_２,
ζ_３）, （ζ_１= １, …, Ｍ_１,ζ_２= １, …Ｍ_２, ζ
_３= １, …，Ｍ_３）とおく。さらに対応リスト［Ｓ⁰ ←
→Ｓ¹ ］とリファレンス系Ｓ⁰ の周期性から

【数３９】

【０２４４】のように、座標ｒ¹ （ｋ０）_αの初期値、
およびポテンシャルパラメ−タＡ¹ （ｋ０）と質量ｍ¹
（ｋ０）の値を設定する方式でもよい。

【０２４５】速度ｖ¹ （ｋ０）_αに関しては、例えば温
度ＴでのＭＤ計算を行う場合、この温度に対応した初速
度を設定する必要があるため、ユ−ザ−によって指定さ
れた温度Ｔを与える初速度を本発明による材料設計支援
システムが自動的に発生させる方式でもよいし、ユ−ザ
−が適宜設定する方式でもよい。

【０２４６】次に図１４の原子間ポテンシャル処理部
［４−４］では原子間相互作用を表す関数形が保持され
る。本実施例では原子(i) と原子（j ）に働く原子間相
互作用を

【数４０】

【０２４７】という関数形によって記述する。ここでＡ
(i) は、リファレンス系Ｓ⁰ 処理都［４−１］、周期系
Ｓ¹ 処理部［４−２］非周期系Ｓ² 処理部［４−３］で
述べたポテンシャルパラメ−タＡ⁰ （ｈζ）、Ａ¹ （ｋ
ξ）、Ａ²(i)に相当する。例えばφ（ｉ,j）がク−ロン
相互作用を表す場合、Ａ(i)=qi（qiは原子(i) の電
荷）、ψ（ｒ|i,j）= １／ｒ（i,j ） EMBED Equation.
2 （ｒ（i,j ）は原子(i)と原子（j ）の原子間距
離）となる。

【０２４８】式(79)のような関数形で記述される原子間
ポテンシャルとして、ク−ロン相互作用の他に、ファン
・デル・ワ−ルス力、多重極相互作用等がある。これら
の相互作用は通常のＭＤ計算において基本的に考慮され
るものであり、式(79)は原子間ポテンシャルの代表的な
関数形を表しているといえる。本発明の材料設計支援シ
ステムでは、例えばク−ロン相互作用＋ファン・デル・
ワ−ルス力というように、同時に複数の原子問相互作用
を扱うこともできる。その場合、原子(i) と原子（j ）
に働く原子間相互作用φ（i,j ）は

【数４１】

【０２４９】で記述される。原子間相互作用を１種類の
み考慮する場合も、複数種類考慮する場合も取り扱いは
同じであるため、上述したリファレンス系Ｓ⁰ 処理部
［４−１］周期系Ｓ¹ 処理部［４−２］および以下の説
明では原子間相互作用を識別するための添字μを省略す
る。また本実施例では原子間ポテンシャルに関する説明
を行なうが、本発明による材料設計支援システムでは分
子間ポテンシャルの取り扱いも同様に可能である。

【０２５０】本発明の材料設計支援システムでは周期系
Ｓ¹ に対するＭＤ計算を行うため、リファレンス系Ｓ⁰
と周期系Ｓ¹ のエネルギ−および各原子に働く力をぞれ
ぞれ計算する。以下この方法を説明する。

【０２５１】まず、図１４のリファレンス系Ｓ⁰ 計算部
［４−５］での処理動作を説明する。リファレンス系Ｓ
⁰ 計算部［４−５］では、リファレンス系Ｓ⁰ のポテン
シャルエネルギ−および各原子に働く力を計算する。そ
のために、まず２原子間の位置ベクトルのα成分ｒ⁰
（ｈζ, ｈ' ζ' ）_αを以下のように定義する。

【０２５２】

【数４２】

【０２５３】後述する周期系Ｓ¹ や非周期系Ｓ² に対す
る２原子間の位置べクトルも式(82)と同様に定義される
ものである。リファレンス系Ｓ⁰ 計算部［４−５］で
は、e⁰（ｈ０）とf⁰（ｈ０）_αで表される物理量を次式

【数４３】

【０２５４】原子間栢互作用が無視できるほど小さくな
る遠方の原子まで相互作用考慮することを意味してい
る。ク−ロン相互作用の場合、式(83)と式(84)の格子和
に対してエワルド法を用いて計算の効率化を図ることも
可能である。このエワルド法はク−ロン相互作用だけで
なく、関数形がｒ^−ｎで与えられる相互作用の格子和の
計算に対しても適用が可能である。またe⁰（ｈ０）とf⁰
（ｈ０）_αを用いれば、リファレンス系Ｓ⁰ のユニット
セルＵ⁰ 当たりのポテンシャルエネルギ−Ｅ⁰ _cellと原
子（ｈ０）に働く力のα成分F⁰（ｈ０）_αをそれぞれ

【数４４】

【０２５５】と計算することができる。

【０２５６】次に、図１４の周期系Ｓ¹ 計算部［４−
６］での動作処理を説明する。周期系Ｓ¹ 計算部［４−
６］では、周期系Ｓ¹ のユニットセルＵ¹ 当たりのポテ
ンシャルエネルギ− EMBED Equation.2 と原子（ｋ
０）に働く力のα成分F¹（ｋ０）_αが

【数４５】

【０２５７】のように計算される。ここでe¹（ｋ０）と
f¹（ｋ０）_αはそれぞれ

【数４６】

【０２５８】で与えられるものである。式(90)のe⁰（ｋ
０）と式(91)のf⁰（ｋ０）_αはリファレンス系Ｓ⁰ に対
する物理量であり、対応リスト［Ｓ⁰ ←→Ｓ¹ ］とリフ
ァレンス系Ｓ⁰ の周期性から

【数４７】

【０２５９】この領域は例えば半径Ｒ_cut の球として与
えられるものである。領域Ｒ_cut はユ−ザ−が指定する
方式でもよいし、本発明による材料設計支援システムが
適宜設定する方式でもよい。

【０２６０】なお、後で述べるようにユ−ザ−が必要と
する計算精度と計算時間の兼ね合いから領域Ｒ_cut を設
定することによって、周期系Ｓ¹ に対する効率的なＭＤ
計算を実現することができる。またＭＤ計算実行中は、
通常、系の全エネルギ−（運動エネルギ−＋ポテンシャ
ルエネルギ−）をモニタ−しエネルギ−の保存則を確認
するが、本発明の計算方法においても運動エネルギ−を

【数４８】

【０２６１】のように計算し、周期系Ｓ¹ のユニットセ
ルＵ¹ 当たりの全エネルギ−をＫ¹ _{cel l} ＋Ｅ¹ _cell EMBE
D Equation.2 と求めることができる。

【０２６２】次に、図１４の周期系ＭＤ計算部［４−
８］での処理動作を説明する。本実施例では、ＭＤ計算
のアルゴリズム（ニュ−トン方程式の差分化の方法）と
してベルレの方法を、またアンサンブルとしてはミクロ
カノニカル・アンサンブル（粒子数、体積、エネルギ−
一定）を例に取って説明を行うが、これら以外の方法の
適用も可能である。周期系ＭＤ計算部［４−８］では周
期系Ｓ¹ に対するＭＤ計算を行うため、基本セル（ξ=
０）に含まれるＮ¹ 個の原子（ｋ０）の時刻ｔ＋Δｔで
の位置と速度を

【数４９】

【０２６３】を用いて計算する。ここでF¹（ｋ０｜ｔ）
_αは、周期系Ｓ¹ 計算部［４−６］で式(89)に従って計
算される時刻ｔでの原子（ｋ０）に働く力のα成分であ
る。ΔｔはＭＤ計算における時間刻みであり、この値は
ユ−ザ−が任意に指定してもよいし、本発明の材料設計
支援システム側で設定する方式でもよい。

【０２６４】周期系ＭＤ計算部［４−８］では、式(96)
と式(97)を用いて周期系Ｓ¹ の構造の動的挙動を時刻ｔ
_end まで繰り返し計算する。この終了時刻ｔ_end はユ−
ザ−によって任意に指定される値である。

【０２６５】以上、本発明による周期系Ｓ¹ に対するＭ
Ｄ計算を実行する材料設計支援システムの構成を示した
ブロック図（図１４）における処理動作を各ブロック毎
に説明してきた。次に周期系Ｓ¹ に対するＭＤ計算を行
うための手順をフロ−チャ−トを用いて具体的に説明す
る。

【０２６６】図２０に周期系Ｓ¹ に対するＭＤ計算を行
うための手順を示した。

【０２６７】まず、ＭＤ計算のための諸情報（温度Ｔ、
時間刻みΔｔ、終了時刻ｔ_end 、Ｒ_cut 等）を入力する
（ステップ［５−１］）。これらの値の入力は、温度Ｔ
のようにこの後のステップでの処理に関係するものを除
けば、以下のステップ［５−２］やステップ［５−３］
の後に入力する方式でもよい。

【０２６８】次に、リファレンス系Ｓ⁰ 処理部［４−
１］にリファレンス系Ｓ⁰ に関する物理量を入力する
（ステップ［５−２］）。この物理量は、リファレンス
系Ｓ⁰ の格子定数、ユニットセルＵ⁰ 当たりの原子数Ｎ
⁰ 、Ｎ⁰ 個の原子の座標ｒ⁰ （ｈ０）_α、ポテンシャル
パラメ−タＡ⁰ （ｈ０）、質量ｍ⁰ （ｈ０）である。

【０２６９】次に、周期系Ｓ¹ 処理部［４−２］に周期
系Ｓ¹ に関する物理量を入力する（ステップ［５−
３］）。この物理量は、周期系Ｓ¹ の格子定数、ユニッ
トセルＵ¹ 当たりの原子数Ｎ¹ 、Ｎ¹ 個の原子の座標ｒ
¹ （ｋ０）_α、速度ｖ¹ （ｋ０）、ポテンシャルパラメ
−タＡ¹ （ｋ０）、質量ｍ¹ （ｋ０）である。

【０２７０】次いで、同じく周期系Ｓ¹ 処理部［４−
２］において対応リスト［Ｓ⁰ ←→Ｓ¹ ］を作成する
（ステップ［５−４］）。ブロック毎の説明で既に述べ
たように、周期系Ｓ¹ 処理部［４−２］に対する物理量
の入力は、リファレンス系Ｓ⁰ の物理量や温度Ｔなどに
基づいて本発明の材料設計支援システムが自動的に値
（座標と速度に関しては初期値）を設定する方式でもよ
い。

【０２７１】次いでリファレンス系Ｓ⁰ 計算部［４−
５］において、f⁰（ｈ０）_αを計算する（ステップ［５
−５］）。

【０２７２】以上までのステップがＭＤ計算の初期設定
に当たり、これ以降のステップがＭＤ計算に従う繰り返
し計算の処理動作となる。繰り返し計算の各ステップが
ＭＤ計算の時刻ｔに対応しており、例えば初期時刻をｔ
= ０とする。

【０２７３】まず周期系Ｓ¹ 計算部［４−６］において
F¹（ｋ０）_αを計算し（ステップ［５−６］）、次に周
期系ＭＤ計算部［４−８］でF¹（ｋ０）_αを用いてｒ¹
（ｋ０）_αとｖ¹ （ｋ０）_αを計算する（ステップ［５
−７］）。それを受けて、周期系Ｓ¹ 処理部［４−２］
においてｒ¹ （ｋ０）_αとｖ¹ （ｋ０）_αの値をそれぞ
れ更新する（ステップ［５−８］）。以上のステップ
［５−６］からステップ［５−８］までの処理を、ステ
ップ［５−９］で時刻ｔがｔ= ｔ_end になったと判定さ
れるまで繰り返し実行し、最後に時刻ｔ_end での周期系
Ｓ¹ の構造を出力あるいは表示する（ステップ［５−１
０］）。

【０２７４】本発明によるＭＤ計算では、繰り返し計算
実行中の任意の時刻において、ｒ¹（ｋ０）_α、ｖ¹
（ｋ０）_α、およびこれらの値を用いて計算した物理量
を出力あるいは表示することも可能である。またステッ
プ［５−６］では、ポテンシャルエネルギ− EMBED Equ
ation.2 と運動エネルギ− EMBED Equation.2 をユ
−ザ−の希望する任意の時刻において計算・出力するこ
とも可能である。ただしポテンシャルエネルギ−を計算
する場合は、ステップ［５−５］においてe⁰（ｈ０）の
計算を予め行っておく。

【０２７５】ここで周期系Ｓ¹ に対するＭＤ計算のフロ
−チャ−ト（図２０）に関する補足説明を行う。

【０２７６】図２０ではステップ［５−６］からステッ
プ［５−８］までの処理を周期系Ｓ¹ 計算部［４−
６］、周期系ＭＤ計算部［４−８］、周期系Ｓ¹ 処理部
［４−２］ごとの処理に分けた単純なステップ化で表し
たが、ベルレのアルゴリズムを用いる場合、式(97)によ
れば時刻ｔ+ Δｔの速度を求めるために時刻ｔと時刻ｔ
＋Δｔでの力が必要となる。

【０２７７】そのため実際の計算では、例えばベルレの
アルゴリズムを用いる場合、図２０のフロ−チャ−トを
図２１のように一部変更する。図Ｓではステップ［５−
１１］が新たに設けられ、このステップ［５−１１］で
周期系Ｓ¹ 計算部［４−６］において時刻ｔ= ０（初期
時刻を０とする）での力F¹（ｋ０｜ｔ（= ０））_αを求
める。そして、図２０のステップ［５−６］からステッ
プ［５−８］までを、図２１では次のような処理と置き
換える。

【０２７８】まず、周期系ＭＤ計算部［４−８］におい
て時刻ｔ十Δｔでの座標ｒ¹ （ｋ０｜ｔ＋Δｔ）_αを式
(96)を用いて計算し（ステップ［５−１２］）、これを
受けて周期系Ｓ¹ 処理部［４−２］において座標の値を
ｒ¹ （ｋ０｜ｔ）_αからｒ¹（ｋ０｜ｔ+ Δｔ）_αに更
新する（ステップ［５−１３］）。

【０２７９】次に、この更新された座標ｒ¹ （ｋ０｜ｔ
+ Δｔ）_αを用いて時刻ｔ+ Δｔでの力F¹（ｋ０｜ｔ+
Δｔ）_αを計算する（ステップ［５−１４］）。そし
て、このF¹（ｋ０｜ｔ+ Δｔ）_αと一つ前の繰り返し計
算のステップ［５−１４］（ｔ= ０の場合はステップ
［５−１１］）で既に計算済みのF¹（ｋ０｜ｔ）_αとを
用いて、周期系ＭＤ計算部［４−８］において時刻ｔ+
Δｔでの速度ｖ¹ （ｋ０｜ｔ+ Δｔ）_αを式(97)を用い
て計算し（ステップ［５−１５］）、それを受けて周期
系Ｓ¹ 処理部［４−２］において速度の値をｖ¹ （ｋ０
｜ｔ）_αからｖ¹ （ｋ０｜ｔ+ Δｔ）_αに更新する（ス
テップ［５−１６］）。以上のステップ［５−１１］と
ステップ［５−１２］〜ステップ［５−１６］は、用い
るアルゴリズム（例えぱベルレの方法）に応じて手順が
変化するものであり、本発明による材料設計支援システ
ムでは採用するアルゴリズムによって適宜フロ−チャ−
トを変更することが可能である。

【０２８０】（第６の実施形態）非周期系に対するＭＤ
計算 この第６の実施形態では、非周期系Ｓ² に対するＭＤ計
算を実行するための材料設計支援システムの一実施例を
説明する。

【０２８１】この実施形態は、請求項２に記載された発
明に対応するものである。

【０２８２】図１５は非周期系Ｓ² に対するＭＤ計算を
実行する材料設計支援システムの構成を示したブロック
図であり、リファレンス系Ｓ⁰ 処理部［４−１］、周期
系Ｓ¹ 処理部［４−２］、非周期系Ｓ² 処理部［４−
３］、原子間ポテンシャル処理部［４−４］、リファレ
ンス系Ｓ⁰ 計算部［４−５］、周期系Ｓ¹ 計算部［４−
６］、非周期系Ｓ² 計算部［４−７］、周期系ＭＤ計算
部［４−８］、非周期系ＭＤ計算部［４−９］から構成
される。

【０２８３】この図１５は図１４に対して非周期系Ｓ²
処理部［４−３］、非周期系Ｓ² 計算部［４−７］、非
周期系ＭＤ計算部［４−９］の３つを新たにつけ加えた
構成になっており、リファレンス系Ｓ⁰ 処理部［４−
１］、周期系Ｓ¹ 処理部［４−２］、原子間ポテンシャ
ル処理部［４−４］、リファレンス系Ｓ⁰ 計算部［４−
５］、周期系Ｓ¹ 計算部［４−６］、周期系ＭＤ計算部
［４−８］は図１４の各部と同しものであるため、ここ
では非周期系Ｓ² 処理部［４−３］、非周期系Ｓ² 計算
部［４−７］、非周期系ＭＤ計算部［４−９］について
説明を行った後、フロ−チャ−トを用いて非周期系Ｓ²
に対するＭＤ計算の手順を説明する。

【０２８４】まず、図１５の非周期系Ｓ² 処理部［４−
３］で行われる処理を説明する。非周期系Ｓ² 処理部
［４−３］では、非周期系Ｓ² の結晶構造、ポテンシャ
ルパラメ−タ等が保有・処理される。本実施形態では非
周期系Ｓ^２の構造を周期系Ｓ¹からのずれによって記述
する。

【０２８５】このずれは、局所的に乱れ（例えぱ不純
物）が存在することによって乱れを中心とするある領域
内で構造変形が起こる結果生じるものであり、以下、本
実施例では非周期系Ｓ² として局所的な不純物を含む系
を例にとって説明する。

【０２８６】このようなずれによる記述を行うため、対
応リスト［Ｓ¹ ←→Ｓ² ］を用いて周期系Ｓ¹ と非周期
系Ｓ² の原子を１対１に対応付ける。非周期系Ｓ² で
は、局所的な不純物の存在のために周期系Ｓ^１が有して
いた周期性が消滅する。そのため、周期性を表す（ｋ
ξ）という原子番号の表記方法は適切ではなく、非周期
系Ｓ² の構成原子に対しては通し番号(i) を設定する。
このような処理は例えばi=ｋ+ ξＮ¹ という関係式を用
いて行なうことが可能であり、この情報が対応リスト
［Ｓ¹ ←→Ｓ² ］に保存される。

【０２８７】図１７に、リファレンス系Ｓ⁰ （構成原子
とユニットセルを点線で表示）、周期系Ｓ¹ （構成原子
とユニットセルを実線で表示）、非周期系Ｓ² （構成原
子を灰色の丸で、後述する領域Ｒ_relax を一点鎖線で表
示）の対応関係の一例を示した。

【０２８８】図１７のリファレンス系Ｓ⁰ および周期系
Ｓ¹ は図１６で示した例と同じものであり、この図１６
の周期系Ｓ¹ に対して置換型不純物原子（黒丸で表示）
を１個導入した系が図１７の非周期系Ｓ² となってい
る。すなわち非周期系Ｓ² は、領域Ｒ_relax 内の原子は
灰色丸（不純物原子は黒丸）で、領域Ｒ_relax 外の原子
は白丸（実線）で表される原子配置を持つものである。

【０２８９】図１７に対する対応リスト［Ｓ¹ ←→Ｓ
² ］として、非周期系Ｓ² の２１個の構成原子の通し番
号(i) を登録した例を図１８に示した。この２１個の原
子は図１７の一点鎖線で囲まれた領域Ｒ_relax に合まれ
る原子に相当している。図１８に示した対応リスト［Ｓ
¹ ←→Ｓ² ］の例では、対応リスト［Ｓ⁰ ←→Ｓ¹ ］と
の対応関係を明らかにするために、非周期系Ｓ² の原子
と対応していない周期系Ｓ¹ の原子（例えば（ｋ= ５,
ξ= ０）や（ｋ= ６, ξ= ０））までリストアップされ
ているが、このような原子に関しては対応リスト［Ｓ¹
←→Ｓ² ］に登録しない方式でもよい。

【０２９０】また図１７に示した例では領域Ｒ_relax が
周期系Ｓ¹ のユニットセルＵ¹ よりも小さくなっている
が、領域Ｒ_relax がユニットセルＵ¹ よりも大きい場合
も同様な取り扱いが可能である。以下では説明を分かり
易いものとするため、非周期系Ｓ² の不純物原子を通し
番号(i) の代わりに不純物番号（σ）でも表す。

【０２９１】以下、非周期系Ｓ² 処理部［４−３］での
処理動作の流れをＭＤ計算の初期設定部分に関して説明
する。非周期系Ｓ² 処理部［４−３］には、非周期系Ｓ
² に含まれる不純物原子の個数Ｎ^imp 、不純物原子の番
号、不純物原子のポテンシャルパラメ−タＡ^imp （σ）
および質量ｍ^imp （σ）、格子緩和を考慮する領域Ｒ
_relax が入力される。

【０２９２】ここで不純物原子の番号とは、不純物原子
との置換えを行なう周期系Ｓ¹ の原子の番号である。す
なわち、不純物原子をＮ^imp 個発生させる場合は、{
（ｋ_１ξ_１）, （ｋ_２ξ_２）…，（ｋ_Ｎimp ξ
_Ｎimp ）} という組みを指定する。

【０２９３】図１７に示した例では不純物原子は１個
（Ｎ^imp=１）であり、{ （ｋ= ２; ξ= ０）} のみを指
定する。Ｎ^imp として任意の個数の不純物原子を取り扱
うことが可能であり、また、指定する不純物原子は周期
系Ｓ¹ の基本セル（ξ= ０）に含まれるものでもよい
し、周囲のユニットセル（ξ≠０）に含まれるものでも
よい。さらに、この不純物原子を指定する組み{ （ｋ_１
ξ_１）, （ｋ_２ξ_２）…，（ｋ_Ｎimp ξ_Ｎimp ）} は、
対応リスト［Ｓ¹ ←→Ｓ² ］によって不純物番号の組み
{ （σ_１）, （σ_２）…，（σ_Ｎimp ）} と対応付けら
れている。

【０２９４】図１８の対応リスト［Ｓ¹ ←→Ｓ² ］では
不純物原子は１個であり、非周期系Ｓ² の原子（i=２）
に不純物番号（σ= １）が割り振られている。次に、格
子緩和を考慮する領域Ｒ_relax について説明する。非周
期系Ｓ² では不純物を中心として格子緩和が生じ、各原
子の緩和量は不純物より遠方ではほとんど０になるた
め、不純物を中心とするある有限な領域内のみで格子緩
和を計算すればよい。この格子緩和を考慮する領域とし
て、例えば不純物原子の重心を中心とする半径Ｒ_relax
の球を仮定することができる。この領域は必ずしも球形
である必要はなく、木実施例では格子緩和を考慮する領
域を一般的にＲ_relax で表す。この領域Ｒ_relax は、ユ
−ザ−が任意に指定する方式でもよいし、本発明の材料
設計支援システムが適宜設定する方式でもよい。

【０２９５】また非周期系Ｓ⁰ 処理部［４−３］では、
領域Ｒ_relax に含まれる原子数Ｎ²（Ｎの二乗ではない
ことに注意。混乱を避けるため本明細書ではＮの二乗は
Ｎ×Ｎと表記する）を数え上げる。図１７に示した例で
はＮ²=２１となっており、常にＮ² ＞Ｎ^imp となる。補
足説明となるが、本発明が扱う非周期系Ｓ² と同し系を
ス−パ−セル法で取り扱う場合に必要となるユニットセ
ル当たりの原子数（前述したＮ_SCに相当）は、通常Ｎ_SC
> Ｎ² となる（ス−パ−セル法では乱れ同士の相互作用
を小さくするために、領域Ｒ_relax の周囲にさらにバッ
ファ−層のような領域を必要とするため）。

【０２９６】よって計算量を単純に比較しただけでも、
本発明による計算方法の方がス−パ−セル法よりも非周
期系に対する計算量が少なくなっていることが分かる。
次に非周期系Ｓ² 処理部［４−３］では、対応リスト
［Ｓ¹ ←→Ｓ² ］と周期系Ｓ¹の周期性から、次式

【数５０】

【０２９７】に基づいて座標ｒ²(i)_αと速度ｖ²(i)_αの
初期値を設定する。この初期値は必ずしも式(98)および
式(99)を用いて設定する必要はなく、ユ−ザ−が適宜設
定する方式でもよい。次に、非周期系Ｓ² に対するポテ
ンシャルパラメ−タと質量の取り扱いを説明する。上述
したように非周期系Ｓ² 処理部［４−３］には、Ｎ^imp
個の不純物原子（σ）のポテンシャルパラメ−タＡ^imp
（σ）と質量ｍ^imp （σ）が入力される。これを受けて
非周期系Ｓ² 処理部［４−３］では、領域Ｒ_{rela x} 、に
含まれるＮ² 個の原子(i) のポテンシャルパラメ−タＡ
²(i)および質量ｍ²(i)を

【数５１】

【０２９８】と設定する。

【０２９９】周期系Ｓ¹ 処理部［４−２］の説明で述べ
たように、周期系Ｓ¹ の原子は対応リスト［Ｓ⁰ ←→Ｓ
¹ ］で対応付けられたリファレンス系Ｓ⁰ の原子と同じ
ポテンシャルパラメ−タを持つ（式(77)参照）のに対
し、非周期系Ｓ² に関しては、不純物原子（σ）につい
ては対応リスト［Ｓ¹ ←→Ｓ² ］で対応付けられた周期
系Ｓ¹ の原子とは異なるポテンシャルパラメ−タを持つ
（式(101) 参照）。

【０３００】本発明による材料設計支援システムでは、
不純物原子として置換型、侵入型および欠陥の３種類を
取り扱うことが可能であり、それらは周期系Ｓ¹ と非周
期系Ｓ² のポテンシャルパラメ−タの関係によって図１
９のように指定される。図１９にはポテンシャルパラメ
−タＡと共に質量ｍの取り扱いに関する規則も示した
が、この質量ｍはポテンシャルパラメ−タＡと同様の取
り扱いがなされるものであり（上述のリファレンス系Ｓ
⁰ 処理都［４−１］、周期系Ｓ¹ 処理部［４−２］、お
よび非周期系Ｓ² 処理部［４−３］の説明参照）、質量
ｍをポテンシャルパラメ−タＡの一種とみなすこともで
きる。

【０３０１】図１９に示した３種類の不純物原子のう
ち、取扱いに特に注意を要するのが侵入型不純物原子で
ある。侵入型不純物原子はそれに対応する原子が周期系
Ｓ¹ に本来存在しないため、対応リスト［Ｓ¹ ←→Ｓ
² ］で１対１の対応付けを行うためには周期系Ｓ¹ （及
びリファレンス系Ｓ⁰ ）に予めダミ−原子を設定してお
く必要が生しる。このダミ−原子はポテンシャルパラメ
−タＡ¹ （及びＡ⁰ ）が０に設定されており、そのため
ダミ−原子の考慮によって周期系Ｓ¹ （及びリファレン
ス系Ｓ⁰ ）のエネルギ−や各原子に働く力が変化するよ
うな悪影響は一切生じない。

【０３０２】次に、図１５の非周期系Ｓ² 計算部［４−
７］での動作処理を説明する。非周期系Ｓ² 計算部［４
−７］では、非周期系Ｓ² の差分ポテンシャルエネルギ
−ΔE²および原子(i) に働く力のα成分F²(i) _αを、次
式

【数５２】

【０３０３】を用いて計算する。ここでF¹(i) _αは周期
系Ｓ¹ における力であり、対応リスト［Ｓ¹ ←→Ｓ² ］
と周期系Ｓ¹ の周期性から

【数５３】

【０３０４】が成り立ち、式(89)を用いて周期系Ｓ¹ 計
算部［４−６］で計算されるF¹（ｋ０）_αに相当する。
またΔE^{2（a ）}，ΔE^{2（b ）}，ΔF^{2（a ）}(i) _α，ΔF
^{2（b ）}(i) _αは、エネルギ−に関しては

【数５４】

【０３０５】と、力に関しては

【数５５】

【０３０６】と記述される量である。ここで式(105) の
e¹(i) と式(107) のf¹(i) _αは周期系Ｓ¹ に対する物理
量であり、対応リスト［Ｓ¹ ←→Ｓ² ］と周期系Ｓ¹ の
周期性から

【数５６】

【０３０７】ここで差分ポテンシャルエネルギ−ΔE²に
関する補足説明を行う。ΔE²は、次式

【数５７】

【０３０８】で定義される物理量であり、Ｅ¹ _WholeとＥ
² _Wholeはそれぞれ周期系Ｓ¹ と非周期系Ｓ² の系全体に
対するポテンシャルエネルギ−を表す。すなわちΔE
²は、周期系Ｓ¹ に不純物が局所的に導入されることに
よって生じるエネルギ−の変化分に相当する。非周期系
Ｓ² は周期性の単位（ユニットセル）を持たないため、
リファレンス系Ｓ⁰ や周期系Ｓ¹ に対するＥ⁰ _Cell やＥ
¹ _Cell のような量を定義することができない。そのた
め、本実施例では周期系Ｓ¹ からのポテンシャルエネル
ギ−変化分ΔＥ² によって非周期系Ｓ² のエネルギ−を
記述する。また運動エネルギ−の変化分ΔＫ² を

【数５８】

【０３０９】と計算することで、非周期系Ｓ² の全エネ
ルギ−を周期系Ｓ¹ からの変化分ΔＫ² ＋ΔE²として求
めることができる。次に式(106) と式(108) の計算方法
について補足説明を行う。

【０３１０】

【数５９】

【０３１１】次に、図１５の非周期系ＭＤ計算部［４−
９］での処理動作の説明を行う。第６の実施形態では、
第５の実施形態と同様にＭＤ計算のアルゴリズム（ニュ
−トン方程式の差分化の方法）としてベルレの方法を、
またアンサンブルとしてはミクロカノニカル・アンサン
ブル（粒子数、体積、エネルギ−一定）を例に取って説
明を行うが、これら以外の方法の適用も可能である。非
周期系ＭＤ計算部［４−９］では非周期系Ｓ² に対する
ＭＤ計算を行うため、領域Ｒ_relax に含まれるＮ² 個の
原子(i) の時刻ｔ+ Δｔでの位置と速度を

【数６０】

【０３１２】を用いて計算する。周期系ＭＤ計算部［４
−８］での扱いと同しく、ΔｔはＭＤ計算の時間刻みで
あり、F²（i ｜ｔ）_αは非周期系Ｓ² 計算部［４−７］
で式(103) を用いて計算される時刻ｔでの原子(i) に働
く力のα成分である。本発明による材料設計支援システ
ムで非周期系Ｓ² に対するＭＤ計算を行う場合、この非
周期系ＭＤ計算部［４−９］で式(114) と式(115) を用
いて非周期系Ｓ² の構造の時間変化を計算すると同時
に、上述した周期系ＭＤ計算部［４−８］においても周
期系Ｓ¹ に対するＭＤ計算を実行することになる。この
手順について、フロ−チャ−トを用いて以下説明を行
う。

【０３１３】図２２に非周期系Ｓ² に対するＭＤ計算の
手順を示した。ステップ［６−１］からステップ［６−
４］までは、先に説明した図２０のステップ［５−１］
からステップ［５−４］までと同じ処理を行う。次に、
非周期系Ｓ² 処理部［４−３］に非周期系Ｓ² に対する
物理量を入力する（ステップ［６−５］）。この物理量
は、非周期系Ｓ² に含まれる不純物原子の個数Ｎ^imp 、
Ｎ^imp 個の不純物原子（σ）の番号、不純物原子（σ）
のポテンシャルパラメ−タＡ^imp （σ）と質量ｍ^imp
（σ）、および格子綬和を考慮する領域Ｒ_relax であ
る。これらの物理量が入力されると、同じく非周期系Ｓ
² 処理部［４−３］において対応リスト［Ｓ¹ ←→Ｓ
² ］を作成し、非周期系Ｓ² の領域Ｒ_relax に含まれる
Ｎ² 個の原子(i) の座標ｒ²(i)_αおよび速度ｖ²(i)_αの
初期値、ポテンシャルパラメ−タＡ²(i)と質量ｍ²(i)の
値を設定する（ステップ［６−６］）。次いで、リファ
レンス系Ｓ⁰ 計算部［４−５］においてf⁰（ｈ０）_αを
計算する（ステップ［６−７］）。ここまでのステップ
がＭＤ計算の初期設定に当たり、これ以降のステップが
ＭＤ計算に従う繰り返し計算の処理動作となる。まず周
期系Ｓ¹ 計算部［４−６］においてF¹（ｋ０）_αを計算
し（ステップ［６−８］）、続いて周期系ＭＤ計算部
［４−８］でF¹（ｋ０）_αを用いてｒ¹ （ｋ０）_αとｖ
¹ （ｋ０）_αを計算する（ステップ［６−９］）。次
に、非周期系Ｓ² 計算部［４−７］においてF²(i) _αを
計算し（ステップ［６−１０］）、非周期系ＭＤ計算部
［４−９］でF²(i) _αを用いてｒ²(i)_αとｖ²(i)_αを計
算する（ステップ［６−１１］）。さらに、周期系Ｓ¹
処理部［４−２］においてｒ¹ （ｋ０）_αとｖ¹ （ｋ
０）_αの値を更新し（ステップ［６−１２］）、非周期
系Ｓ² 処理部［４−３］においてｒ²(i)_αとｖ²(i)_αの
値を更新する（ステップ［６−１３］）。以上のステッ
プ［６−８］からステップ［６−１３］までの処理をス
テップ［６−１４］で時刻ｔがｔ= ｔ_end になったと判
定されるまで繰り返し実行し、最後に時刻ｔ_end での非
周期系Ｓ² の構造を出力あるいは表示する（ステップ
［６−１５］）。本発明によるＭＤ計算では、繰り返し
計算実行中の任意の時刻において、ｒ²(i)_α、ｖ
²(i)_α、およびこれらの値を用いて計算した物理量を出
力あるいは表示することも可能である。また、ステップ
［６−１０］ではポテンシャルエネルギ−と運動エネル
ギ−に関してΔE²とΔＫ² をユ−ザ−の希望する任意の
時刻に計算・出力することも可能であり、ただしポテン
シャルエネルギ−を計算する場合は、ステップ［６−
７］でe⁰（ｈ０）を、ステップ［６−８］でe¹（ｋ０）
を予め計算しておく。上述したように本発明による材料
設計支援システムで非周期系Ｓ² に対するＭＤ計算を行
う場合、同時に周期系Ｓ¹ に対するＭＤ計算も実行する
ことになる。そのため周期系Ｓ¹ に対しても、任意の時
刻でのｒ¹ （ｋ０）_α、ｖ¹ （ｋ０）_α、およびこれら
の値を用いて計算した物理量の出力あるいは表示、同じ
く任意の時刻での EMBED Equation.2 および EMBED E
quation.2 の計算と出力（ステップ［６−８］）、時
刻ｔ_end での周期系Ｓ¹ の最終構造の出力あるいは表示
（ステップ［６−１５］）が可能である。この図２２の
フロ−チャ−トに関しては、用いるアルゴリズム（例え
ばベルレの方法）に応じて図２０に対する図２１のよう
な手順の変更が可能であり、本発明による材料設計支援
システムでは非周期系Ｓ² に対しても用いるアルゴリズ
ムに応じたＭＤ計算を適宣実行することが可能である。

【０３１４】（実施形態７の実施形態）以下、前記第
５、第６の実施形態の材料設計支援システムに関する変
形例を第７の実施形態として説明する。

【０３１５】まず、式(79)で与えられる関数形以外の原
子間相互作用の取り扱いについて説明する。

【０３１６】式(79)で記述されるｒ^−ｎ関数形（特にク
−ロン相互作用）は、長距離相互作用であるがために、
計算時間の増大、ス−パ−セル法を適用したことよる計
算誤差などの深刻な問題を抱えており、本発明による材
料設計支援システムではそれの解決を目的としている。

【０３１７】式(79)で記述されない原子問相互作用とし
ては、例えば短距離反発力を表すBorn-Mayerタイプ（指
数関数形）や共有結合を表すValence-Force Field モデ
ル（多体力）等があるが、これらは一般的に短距離相互
作用となっている。

【０３１８】この短距離相互作用に関しては上述した計
算時間の増大や計算誤差などの問題は発生しないため、
通常の方法でエネルギ−や力をリファレンス系Ｓ⁰ 、周
期系Ｓ¹ 、非周期系Ｓ² に対してそれぞれ計算すること
に問題はなく、本発明による材料設計支援システムでは
そのようにして計算した短距離相互作用によるエネルギ
−や力を本発明による計算方法で求めた長距離相互作用
によるエネルギ−や力と足し合わせることで、一般的な
原子間相互作用を用いたＭＤ計算を実行することが可能
である。

【０３１９】第６の実施形態の材料設計支援システムに
よって非周期系Ｓ² に対するＭＤ計算を行う場合、上述
したように同時に周期系Ｓ¹ に対してもＭＤ計算を行う
ことになる。材料設計等の目的で実際に非周期系Ｓ² に
対するＭＤ計算を行う場合、同じ母体（周期系Ｓ¹ ）に
対して不純物の個数や種類を変えて何度もＭＤ計算を実
行する必要が生じると予想される。そのような場合、周
期系Ｓ¹ に対しては全く同じ計算を幾度も繰り返すこと
になり経済的でない。これを解決するため、本発明によ
る材料設計支援システムでは周期系Ｓ¹ に対するＭＤ計
算を１度だけ行い、各時刻の周期系Ｓ¹ の物理量（原子
の位置、速度、力等）に関するデ−クベ−スを作成し、
そのデ−タベ−スを用いることで非周期系Ｓ² に対する
ＭＤ計算を効率的に実行することも可能である。このデ
−タベ−スは、周期系Ｓ¹ の種類および温度等のシミュ
レ−ション条件ごとに作成されるものである。

【０３２０】また、第１の実施形態によって非周期系Ｓ
² に対するＭＤ計算を行う場合には、データベースをリ
ファレンス系Ｓ⁰ の種類及び温度などのシュミレーショ
ン条件毎に作成することで効率化を図ることが可能にな
る。

【０３２１】ここで周期系Ｓ¹ と非周期系Ｓ² の座標と
速度に関して補足説明を行う。本発明の材料設計支援シ
ステムでは、周期系Ｓ¹ をリファレンス系Ｓ⁰ からのず
れで、非周期系Ｓ² を周期系Ｓ¹ からのずれで取り扱う
ことを基本概念としており、それを式で表すと、周期系
Ｓ¹ に関しては

【数６１】

【０３２２】と、非周期系Ｓ² に関しては

【数６２】

【０３２３】とそれぞれ記述することができる。ここ
で、Δｒ¹ （ｋ０）_αとΔｖ¹ （ｋ０）_αが温度効果に
よる微視的な振動（格子振動）によるリファレンス系Ｓ
⁰ からのずれに、Δｒ²(i)_αとΔｖ²(i)_αが局所的な不
純物の存在による周期系Ｓ¹ からのずれに相当してい
る。第５の実施形態と第６の実施形態ではｒ¹ （ｋ０）
_α，ｖ¹ （ｋ０）_α，ｒ²(i)_α，ｖ²(i)_αを直接扱う方
式を説明したが、まずΔｒ¹（ｋ０）_α，Δｖ¹ （ｋ
０）_α，Δｒ²(i)_α，Δｖ²(i)_αを計算した後、式(11
6) 〜式(119) を用いてｒ¹ （ｋ０）_α，ｖ¹ （ｋ０）
_α，ｒ²(i)_α，ｖ²(i)_αを求める方式でもよい。その場
合は、周期系Ｓ¹ と非周期系Ｓ² の力の計算として、F¹
（ｋ０）_α( 式(89)）およびF²(i) _α（式(103) ）の代
わりに

【数６３】

【０３２４】を求める。そして、周期系Ｓ¹ 計算部［４
−６］では式(96)と式(97)のｒ¹,ｖ¹,F¹をそれぞれΔｒ
¹,Δｖ¹,ΔF¹に、非周期系Ｓ² 計算部［４−７］では式
(114)と式(115) のｒ²,ｖ²,F²をそれぞれΔｒ²,Δｖ²,
ΔF²に置き換えて繰り返し計算を行う。式(117) と式(1
19) においてｖ⁰ はリファレンス系Ｓ⁰ の原子の速度を
表すが、リファレンス系Ｓ⁰ は絶対零度での系であるた
めｖ⁰ は一意的にすぺて０に設定される。

【０３２５】本発明による材料設計支援システムによれ
ば、非周期系Ｓ² に対して設定される格子緩和を考慮す
る領域Ｒ_relax をシミュレ−ションの途中で適宜変更す
ることが可能であり、それによって領域Ｒ_relax の最適
化を行い、非周期系Ｓ² の構造を効率的に計算すること
ができる。これを行うため、第２の実施形態で説明した
方法との組み合わせによるＭＤ計算も可能である。

【０３２６】本発明による材料設計支援システムでは、
周期系Ｓ¹ と比較して非周期系Ｓ²が大変形を起こすよ
うな場合でも高精度の計算が可能である。これは式(11
8) で示されるずれΔｒ²(i)_αが大きくなった場合に相
当し、このΔｒ²(i)_αがある程度以上に大きくなると、
式(106) および式(108) の計算精度が悪くなる可能性が
ある。この問題はずれΔｒ²(i)_αが小さくなるように周
期系Ｓ¹ 及びリファレンス系Ｓ⁰ にダミ−原子を設定す
ることで回避可能であり、本発明による材料設計支援シ
ステムでは非周期系Ｓ² が大変形を起こすような場合に
関しても精度の良い計算を行うことができる。

【０３２７】これを行うため、前述の実施形態３で説明
した方法との組み今わせによるＭＤ計算を行うことも可
能である。また、温度効果による格子振動が激しくなる
等の原因でリファレンス系Ｓ⁰ と周期系Ｓ¹ の構造差が
大きくなった場合（式(116)のずれΔｒ¹ （ｋ０）_αが
大きくなった場合に相当）、式(90)および式(91)の計算
精度が悪くなる可能性があるが、上記周期系Ｓ¹ に対す
る処理と同様にリファレンス系Ｓ⁰ にダミ−原子を設定
することでこの問題も回避可能であり、本発明による材
料設計支援システムは温度効果による格子振動が激しい
場合でも周期系Ｓ¹ に対するＭＤ計算を実行可能であ
る。

【０３２８】非周期系Ｓ² に対するＭＤ計算において、
リファレンス系Ｓ⁰ のユニットセルＵ⁰ 当たりの原子数
Ｎ⁰ と周期系Ｓ¹ のユニットセルＵ¹ 当たりの原子数Ｎ
¹ とが同程度である場合、周期系Ｓ¹ に対する計算を必
ずしも行う必要はなく、リファレンス系Ｓ⁰ と非周期系
Ｓ² のみを考慮してＭＤ計算を実行する方式でも良く、
これは第１の実施形態に相当する。

【０３２９】この実施形態ではＭＤ計算のみを説明した
が、本発明による材料設計支援システムではＭＭ（分子
力学）計算の実行も可能であり、その場合、原子の質量
をｍ= １、速度をｖ= ０とすることで、本実施例での説
明と同じ手続きによってＭＭ計算を行うことができる。
また、この処理によってリファレンス系Ｓ⁰ と周期系Ｓ
¹ は全く等価な系となるため、ＭＭ計算では周期系Ｓ¹
に関する計算を省略することも可能である。

【０３３０】以上説明してきた第５〜第７の実施形態に
かかる材料設計支援システムによると、従来法と比較し
て以下のような効果を得ることができる。

【０３３１】周期系Ｓ¹ に対して ・計算時間の短縮 非周期系Ｓ² に対して ・計算時間の短縮 ・完全に孤立した乱れを含む無限系に対する計算（ス−
パ−セル法でサイズ無限大のユニットセルを用いた場合
と等価）の実現 ・電気的中性が破れた系に対する正しい電荷分布での計
算の実現 以上の効果について以下順次説明する。

【０３３２】まず第５の実施形態による周期系Ｓ¹ に対
するＭＤ計算における計算時間の短縮について説明す
る。

【０３３３】

【数６４】

【０３３４】そして式(91)では、領域Ｒ_cut 外の原子
（ｋ' ξ' ）との相互作用としてリファレンス系Ｓ⁰ に
おける相互作用を近似的に用いる（これを与えるのが第
１項目のf⁰（ｋ０）_αである）。

【０３３５】これは、領域Ｒ_cut 内部の原子とは周期系
Ｓ¹ における正しい相互作用を考慮し、領域Ｒ_cut 外部
の原子とは温度効果による格子振動を無視した近似的な
相互作用を考慮することを意味する。

【０３３６】本発明による材料設計支援システムではこ
の領域Ｒ_cut を適切に設定することで、周期系Ｓ¹ に対
する効率的なＭＤ計算を次のように実現する。原子（ｋ
０）と相互作用している原子の数を領域Ｒ_cut 内部と外
部で比較した場合、領域Ｒ_{cu t} 内の原子数が少なくなる
ように領域Ｒ_cut を設定すれば、計算量の多い領域Ｒ
_cut 外の原子との相互作用にリファレンス系Ｓ⁰ の計算
結果を用いることができる。リファレンス系Ｓ⁰ に対す
る力の計算はユニットセルＵ⁰ が小さいためもともと計
算量が少ない上、ＭＤ計算の初期設定段階で一度計算を
行えば繰り返し計算中は再計算の必要がない（図２０参
照）。

【０３３７】そのためこのリファレンス系Ｓ⁰ に対する
計算結果を周期系Ｓ¹ の力の計算で利用することによっ
て、周期系Ｓ¹ に対するＭＤ計算の大幅な時間短縮を行
うことができる。

【０３３８】周期系Ｓ¹ に対するＭＤ計算に要する計算
量の見積は、従来法（スーパセル法）がＮ¹ ×Ｎ¹ ×Ｎ
_all ×Ｎ_stepであるのに対し、本発明による方法ではＮ
¹ ×Ｎ_cut ×Ｎ_stepであり（Ｎ_cut は領域Ｒ_cut に含ま
れる原子数）、Ｒ_cut の設定に応じて従来法よりも（Ｎ
¹ ×Ｎ_all ）／Ｎ_cut 倍の高速化が可能となる。

【０３３９】式(91)で行う近似が計算精度に与える影響
については、まず原子間の相互作用は近距離にある原子
からより大きな影響を受けるため、遠距離の原子との相
互作用のみに近似を導入する式(91)の取り扱いは妥当な
ものであると言える。

【０３４０】さらに式(91)はＲ_cut 外の原子との相互作
用を無視しているわけではなく、平衡位置に固定されて
いると仮定（温度効果による振動を無視）しているだけ
であるため、ＭＤ計算において通常計算される物理量
（温度、圧力、拡散係数等）は、長時問平均されたマク
ロな量であることから式(91)で扱う近似によって大きな
影響は一般的に受けないと考えられる。

【０３４１】以上、本発明による材料設計支援システム
では、ユ−ザ−が必要とする計算精度と計算時間の兼ね
合いからＲ_cut を適宜設定することによって、周期系Ｓ
¹ に対する効率的なＭＤ計算を実現することができる。
実際の材料設計では何種類もの系に対して条件を変えて
ＭＤ計算を何回も実行し、材料設計のための指針を得る
ことが重要となるため、そのような状況において本発明
が提供する材料設計支援システムは特に効果を発揮する
ものである。

【０３４２】次に第６の実施形態による非周期系Ｓ² に
対するＭＤ計算における効果を説明する。非周期系Ｓ²
に対するＭＤ計算を行う際、本発明ではリファレンス系
Ｓ⁰、周期系Ｓ¹ 、非周期系Ｓ² の３つの系を取り扱
う。そして、非周期系Ｓ² に対する力F²(i) _αの計算を
式(103) を用いて行う。

【０３４３】式(103) によれば、F²(i) _αは周期系Ｓ¹
の原子に働く力F¹(i) _αと非周期系Ｓ² に対する物理量
ΔF^2（ａ）(i) _αおよび△F^{2（b ）}(i) _αの３つの項で
与えられる。まずF¹(i) _αに関しては、上述した周期系
Ｓ¹ に対するＭＤ計算における取り扱いと同様に式(91)
において計算の効率化がなされている。

【０３４４】また式(107) で計算されるΔF^2（ａ）(i)
_αは、周期系Ｓ¹ に対する物理量f¹（ｋ０）_αとＮ^imp
（《Ｎ² ）個の不純物原子に関する総和で与えられるた
め、これも計算量が少ない項となっている。最後にΔF
^{2（b ）}(i) _αは、式(108) が差分量[ ψ'(ｒ² ｜i,j)
_α- ψ'(ｒ¹ ｜i,j)_α] に関する総和で与えられるもの
であり、この差分量は不純物から遠ざかるにつれて急速
に０になるため、式(108) に関する計算は収束が非常に
速いものになっている。

【０３４５】以上のような理由からF²(i) _αの計算の効
率化が行われており、非周期系Ｓ²に対するＭＤ計算の
時間短縮が可能となっている。

【０３４６】非周期系Ｓ² に対するＭＤ計算に要する計
算量の大まかな見積は、従来法（ス−パ−セル法）がＮ
_SC×Ｎ_SC×Ｎ_all ×Ｎ_stepであるのに対し、本発明によ
る方法は（Ｎ¹ ×Ｎ_cut ＋Ｎ² ×Ｎ'_cut）×Ｎ_stepとな
る（Ｎ'_cutは非周期系Ｓ² 計算部［４−７］に関する説
明で述べた領域Ｒ'_cutに含まれる原子数である）。

【０３４７】Ｎ_SC、Ｎ¹ 、Ｎ² はほぼ同じオ一ダ−であ
るため省略すると、本発明による材料設計支援システム
では領域Ｒ_cut およびＲ'_cutの設定に応して、ス−パ−
セル法よりも約（Ｎ_SC×Ｎ_all)／( Ｎ_cut ＋Ｎ'_cut）倍
の高速化が可能である。

【０３４８】そしてリファレンス系Ｓ⁰ と非周期系Ｓ²
の二つのみを取り扱う前記第１〜第４の実施形態と対比
して、第５〜第７の実施形態による材料設計支援システ
ムではリファレンス系Ｓ⁰ と非周期系Ｓ² に加えて周期
系Ｓ¹ を考慮することで、温度効果による格子振動によ
って周期系Ｓ¹ が大きなサイズのユニットセルＵ¹ （Ｎ
¹ 》Ｎ⁰ ）を有する場合でも非周期系Ｓ² に対する効率
的なＭＤ計算を実現するものである。

【０３４９】次に、上述した非周期系Ｓ² に対するＭＤ
計算での２番目、３番目の効果について説明する。この
「完全に孤立した乱れを含む無限系に対する計算」およ
び「電気的中性が破れた系に対する正しい電荷分布での
計算」は、いずれも非周期系Ｓ² の構造を周期系Ｓ¹ か
らのずれとして記述したことによって実現されたもので
ある。

【０３５０】本発明による材料設計支援システムでは非
周期系Ｓ² に対して周期的境界条件を課していないこと
から、ス−パ−セル法における異なるユニットセルに含
まれる乱れ同士の相互作用は一切存在せず、完全に孤立
した状態として乱れを取り扱うことができる。また、乱
れ（不純物）が電荷を有する場合、ス−パ−セル法では
ユニットセル当たりの電荷の総和を０にするため人為的
な電荷（通常一様電荷）を各構成原子に割り振る必要が
生じる。これは、ユニットセル当たりの総電荷が０から
ずれた系ではエネルギ−が発散してしまうため、それを
防ぐ目的でテクニカルに要求される制約であり、その結
果、ス−パ−セル法では正しい電荷分布での計算ができ
ないという問題を抱えている。それに対し、本発明では
非周期系Ｓ² に対して周期的境界条件を用いていないた
め、電気的中性が破れた非周期系Ｓ² に対しても正しい
電荷分布での計算が可能である。

【０３５１】以上本発明による材料設計支援システムに
よれば、周期系Ｓ¹ に対する高速なＭＤ計算および非周
期系Ｓ² に対する高精度かつ高速なＭＤ計算を実現する
ことができる。

【０３５２】

【発明の効果】以上説明したように本発明によれば、材
料の設計において、周期系や非周期系の構造安定性や動
的挙動などの解析を行うにあたり、高速でかつ精度の高
い計算（シュミレーション）を行うことができる。

【図面の簡単な説明】

【図１】第１の実施形態に係わる材料設計支援システム
の構成を示すブロック図。

【図２】リファレンス系Ｓ⁰ と非周期系Ｓ² との関係を
説明するための図。

【図３】図２に示したリファレンス系Ｓ⁰ と非周期系Ｓ
² の構成原子にそれぞれ原子番号（h ζ）及び通し番号
(i) を割り振って示す図。

【図４】図３のリファレンス系Ｓ⁰ と非周期系Ｓ² に対
応する対応リストの一例を示す図。

【図５】置換型・侵入型不純物原子及び欠陥のポテンシ
ャルパラメータＡをリファレンス系Ｓ⁰ と非周期系Ｓ²
に関して比較して示す図。

【図６】第１の実施形態で説明するＭＭ計算を実行する
ためのフローチャート。

【図７】第１の実施形態で説明するＭＤ計算を実行する
ためのフローチャート。

【図８】第２の実施形態で説明するＭＭ計算を実行する
ためのフローチャート。

【図９】第２の実施形態で説明するＭＤ計算を実行する
ためのフローチャート。

【図１０】第３の実施形態で説明するＭＭ計算を実行す
るためのフローチャート。

【図１１】第３の実施形態で説明するＭＤ計算を実行す
るためのフローチャート。

【図１２】置換型不純物を含むＮａＣｌのＭＭ計算によ
る安定構造（緩和構造）の計算例を示す図。

【図１３】置換型不純物を含むＮａＣｌのＭＭ計算によ
る安定構造（緩和構造）の計算例を示す図。

【図１４】第５の実施形態による周期系Ｓ¹ に対するＭ
Ｄ計算を実行する材料設計支援システムに係る一実施例
の構成を示すブロック図。

【図１５】第６の実施形態による非周期系Ｓ² に対する
ＭＤ計算を実行する材料設計支援システムに係る一実施
例の構成を示すブロック図。

【図１６】リファレンス系Ｓ⁰ と周期系Ｓ¹ の関係を説
明するための図。

【図１７】リファレンス系Ｓ⁰ 、周期系Ｓ¹ 、非周期系
Ｓ² の関係を説明するための図。

【図１８】対応リスト［Ｓ⁰ ←→Ｓ¹ ］および対応リス
ト［Ｓ¹ ←→Ｓ² ］の一例を示した図。

【図１９】第６の実施形態による非周期系Ｓ² に対する
ＭＤ計算で取り扱う置換型・侵入型不純物原子および欠
陥のポテンシャルパラメ−タＡおよび質量ｍを、周期系
Ｓ¹ と非周期系Ｓ² に関して比較した対応表。

【図２０】第５の実施形態による周期系Ｓ¹ に対するＭ
Ｄ計算を実行するためのフロ−チャ−ト。

【図２１】第５の実施形態による周期系Ｓ¹ に対するＭ
Ｄ計算のフロ−チャ−トの変形例。

【図２２】第６の実施形態による非周期系Ｓ² に対する
ＭＤ計算を実行するためのフロ−チャ−ト。

【符号の説明】

［１−１］…リファレンス系Ｓ⁰ 処理部 ［１−２］…非周期系Ｓ² 処理部 ［１−３］…原子間ポテンシャル処理部 ［１−４］…リファレンス系Ｓ⁰ 計算部 ［１−５］…非周期系Ｓ² 計算部 ［１−６］…ＭＭ計算部（第１の計算部） ［１−７］…ＭＤ計算部（第２の計算部）

───────────────────────────────────────────────────── フロントページの続き (54)【発明の名称】 材料設計支援システム、およびコンピュータシステムに系の安定構造や動的挙動を原子や分子レ ベルで解析させこれにより材料の設計支援を行うためのコンピュータプログラムを格納した記憶 媒体

Claims (10)

(57)【特許請求の範囲】
1. 【請求項１】 材料の設計支援に用いられるシステムで
   あって、局所的な乱れを有する系（非周期系Ｓ² ）の構
   造安定性や動的挙動を原子や分子レベルで解析する材料
   設計支援システムにおいて、 乱れを有さない系（リファレンス系Ｓ⁰ ）に対する物理
   量を処理するリファレンス系Ｓ⁰ 処理部と、 前記非周期系Ｓ² の構造をリファレンス系Ｓ⁰ からのず
   れで記述するために、前記非周期系Ｓ² に対する物理
   量、および前記リファレンス系Ｓ⁰ と前記非周期系Ｓ²
   の構成原子を１対１に対応付けるためのリストを処理す
   る非周期系Ｓ² 処理部と、 原子間または分子間相互作用を記述する関数形を処理す
   る原子間ポテンシャル処理部と、 前記リファレンス系Ｓ⁰ のエネルギ−および各原子また
   は各分子に働く力を計算するリファレンス系Ｓ⁰ 計算部
   と、 前記リファレンス系Ｓ⁰ 計算部で計算されたリファレン
   ス系Ｓ⁰ のエネルギ−および構成原子に働く力を用い
   て、非周期系Ｓ² のエネルギ−および構成原子に働く力
   を計算する非周期系Ｓ² 計算部と、 非周期系Ｓ² 計算部の計算結果を用いて、非周期系Ｓ²
   に対する安定構造あるいは動的挙動の計算（ＭＭ計算あ
   るいはＭＤ計算）を実行するＭＭ／ＭＤ計算部とを有す
   ることを特徴とする材料設計支援システム。
2. 【請求項２】 請求項１記載の材料設計支援システムに
   おいて、 前記非周期系Ｓ² の動的挙動の解析（ＭＤ計算）を有限
   温度（温度Ｔ＞０Ｋ）で行う場合に、 このシステムは、さらに、 前記非周期系Ｓ² の構造を、局所的な乱れを有さないが
   温度効果による原子の振動を考慮した系（周期系Ｓ¹ ）
   からのずれで記述する手段を有し、 この手段は、 周期系Ｓ¹ の構造をリファレンス系Ｓ⁰ （温度効果によ
   る原子の振動を考慮しない系）からのずれで記述するた
   めに、この周期系Ｓ¹ に対する物理量及びリファレンス
   系Ｓ⁰ と周期系Ｓ¹ の構成原子を１対１に対応付けるた
   めのリストを処理する周期系Ｓ¹ 処理部と、 前記リファレンス系Ｓ⁰ 計算部で計算されたリファレン
   ス系Ｓ⁰ のエネルギ及び構成原子に働く力を用いて、周
   期系Ｓ¹ のエネルギ及び構成原子に働く力を計算する周
   期系Ｓ¹ 計算部と、 非周期系Ｓ² の動的挙動を、周期系Ｓ¹ からのずれとし
   て解析するために、周期系ＭＤ計算を実行する周期系Ｍ
   Ｄ計算部とを有し、 前記非周期系ＭＭ／ＭＤ計算部は、前記周期系ＭＤ計算
   部及び非周期系Ｓ² 計算部の計算結果を用いて、非周期
   系Ｓ² に対する動的挙動の計算（ＭＤ計算）を実行する
   ことを特徴とする材料設計支援システム。
3. 【請求項３】 請求項１記載の材料設計支援システムに
   おいて、 非周期系Ｓ² 処理部は、 非周期系Ｓ² の格子緩和を考慮する領域をシュミレーシ
   ョン中に変更する手段を有することを特徴とする材料設
   計支援システム。
4. 【請求項４】 材料の設計支援に用いられるシステムで
   あって、局所的な乱れは有さないが、温度効果による原
   子の振動が存在する系（周期系Ｓ¹ ）の原子又は分子の
   動的挙動を解析する材料設計支援システムにおいて、 温度効果による原子の振動を考慮しない系（リファレン
   ス系Ｓ⁰ ）に対する物理量を処理するリファレンス系Ｓ
   ⁰ 処理部と、 前記周期系Ｓ¹ の構造をリファレンス系Ｓ^０からのずれ
   で記述するために、周期系Ｓ¹ に対する物理量およびリ
   ファレンス系Ｓ⁰ と周期系Ｓ¹ の構成原子を１対１に対
   応付けるためのリストを処理する周期系Ｓ¹ 処理部と、 リファレンス系Ｓ⁰ のエネルギ−および各原子または各
   分子に働く力を計算するリファレンス系Ｓ⁰ 計算部と、 前記リファレンス系Ｓ^０計算部で計算されたリファレン
   ス系Ｓ^０のエネルギ−および構成原子に働く力を用い
   て、周期系Ｓ¹ のエネルギ−および構成原子に働く力を
   計算する周期系Ｓ¹ 計算部と、 前記周期系Ｓ¹ 計算部の計算結果を用いて、周期系Ｓ¹
   に対する動的挙動の計算（ＭＤ計算）を実行する周期系
   ＭＤ計算部とを有することを特徴とする材料設計支援シ
   ステム。
5. 【請求項５】 請求項１あるいは４記載の材料設計支援
   システムにおいて、 リファレンス系Ｓ⁰ 処理部は、 リファレンス系Ｓ⁰ にダミー原子又はダミー分子を設定
   する手段を有することを特徴とする材料設計支援システ
   ム。
6. 【請求項６】 コンピュータシステムに、局所的な乱れ
   を有する結晶系（非周期系Ｓ² ）の安定構造性や動的挙
   動を原子や分子レベルで解析させこれにより材料の設計
   支援を行うためのコンピュータプログラムを格納した記
   憶媒体において、 乱れを有さない系（リファレンス系Ｓ⁰ ）に対する物理
   量を処理するリファレンス系Ｓ⁰ 処理手順と、 前記非周期系Ｓ² の構造をリファレンス系Ｓ⁰ からのず
   れで記述するために、前記非周期系Ｓ² に対する物理
   量、および前記リファレンス系Ｓ⁰ と前記非周期系Ｓ²
   の構成原子を１対１に対応付けるためのリストを処理す
   る非周期系Ｓ² 処理手順と、 原子間または分子間相互作用を記述する関数形を処理す
   る原子間ポテンシャル処理手順と、 前記リファレンス系Ｓ⁰ のエネルギ−および各原子また
   は各分子に働く力を計算するリファレンス系Ｓ⁰ 計算手
   順と、 前記周期系Ｓ¹ 計算部で計算された周期系Ｓ¹ のエネル
   ギ−および構成原子に働く力を用いて、非周期系Ｓ² の
   エネルギ−および構成原子に働く力を計算する非周期系
   Ｓ² 計算手順と、 非周期系Ｓ² 計算部の計算結果を用いて、非周期系Ｓ²
   に対する安定構造あるいは動的挙動の計算（ＭＭ計算あ
   るいはＭＤ計算）を実行するＭＭ／ＭＤ計算手順とを実
   行するための指令を格納することを特徴とする記憶媒
   体。
7. 【請求項７】 請求項６記載の記憶媒体において、 前記非周期系Ｓ² の動的挙動の解析（ＭＤ計算）を有限
   温度（温度Ｔ＞０Ｋ）で行う場合に、 この記憶媒体には、さらに、 前記非周期系Ｓ² の構造を、局所的な乱れを有さないが
   温度効果による原子の振動を考慮した系（周期系Ｓ¹ ）
   からのずれで記述する処理手順が格納され、 この処理手順は、 周期系Ｓ¹ の構造をリファレンス系Ｓ⁰ （温度効果によ
   る原子の振動を考慮しない系）からのずれで記述するた
   めに、この周期系Ｓ¹ に対する物理量及びリファレンス
   系Ｓ⁰ と周期系Ｓ¹ の構成原子を１対１に対応付けるた
   めのリストを処理する周期系Ｓ¹ 処理手順と、 前記リファレンス系Ｓ⁰ 計算部で計算されたリファレン
   ス系Ｓ⁰ のエネルギ及び構成原子に働く力を用いて、周
   期系Ｓ¹ のエネルギ及び構成原子に働く力を計算する周
   期系Ｓ¹ 計算手順と、 非周期系Ｓ² の動的挙動を、周期系Ｓ¹ からのずれとし
   て解析するために、周期系ＭＤ計算を実行する周期系Ｍ
   Ｄ計算手順とを有し、 前記非周期系ＭＭ／ＭＤ計算手順では、前記周期系ＭＤ
   計算手順及び非周期系Ｓ² 計算手順の計算結果を用い
   て、非周期系Ｓ² に対する動的挙動の計算（ＭＤ計算）
   を実行することを特徴とする記憶媒体。
8. 【請求項８】 請求項６記載の記憶媒体において、 非周期系Ｓ² 処理手順は、 非周期系Ｓ² の格子緩和を考慮する領域をシュミレーシ
   ョン中に変更する手順を含むことを特徴とする記憶媒
   体。
9. 【請求項９】 コンピュータシステムに、局所的な乱れ
   は有さないが、温度効果による原子の振動が存在する系
   （周期系Ｓ¹ ）での原子又は分子の動的挙動を解析させ
   これにより材料の設計支援を行うためのコンピュータプ
   ログラムを格納した記憶媒体において、 温度効果による原子の振動を考慮しない系（リファレン
   ス系Ｓ⁰ ）に対する物理量を処理するリファレンス系Ｓ
   ⁰ 処理手順と、 前記周期系Ｓ¹ の構造をリファレンス系Ｓ^０からのずれ
   で記述するために、周期系Ｓ¹ に対する物理量およびリ
   ファレンス系Ｓ⁰ と周期系Ｓ¹ の構成原子を１対１に対
   応付けるためのリストを処理する周期系Ｓ¹ 処理手順
   と、 リファレンス系Ｓ⁰ のエネルギ−および各原子または各
   分子に働く力を計算するリファレンス系Ｓ⁰ 計算手順
   と、 前記リファレンス系Ｓ^０計算手順で計算されたリファレ
   ンス系Ｓ^０のエネルギ−および構成原子に働く力を用い
   て、周期系Ｓ^１のエネルギ−および構成原子に働く力を
   計算する周期系ＭＤ計算手順と、 前記周期系Ｓ¹ 計算部の計算結果を用いて、周期系Ｓ¹
   に対する動的挙動の計算（ＭＤ計算）を実行する周期系
   ＭＤ計算手順とを有することを特徴とする記憶媒体。
10. 【請求項１０】 請求項６あるいは９記載の記憶媒体に
    おいて、 リファレンス系Ｓ^０処理手順は、 リファレンス系Ｓ^０にダミー原子又はダミー分子を設定
    する手順を含むことを特徴とする記憶媒体。