特許 7381055 結晶構造算出方法、結晶構造算出プログラムおよび結晶構造算出装置

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2023-11-07

(45)【発行日】2023-11-15

(54)【発明の名称】結晶構造算出方法、結晶構造算出プログラムおよび結晶構造算出装置

(51)【国際特許分類】

   G01N 23/20 20180101AFI20231108BHJP

   G16C 60/00 20190101ALI20231108BHJP

【ＦＩ】

G01N23/20

G16C60/00

【請求項の数】 7

(21)【出願番号】P 2019150697

(22)【出願日】2019-08-20

(65)【公開番号】P2021032617

(43)【公開日】2021-03-01

【審査請求日】2022-07-22

【国等の委託研究の成果に係る記載事項】（出願人による申告）平成３０年度よりの、国立研究開発法人科学技術振興機構、戦略的創造研究推進事業「機械学習と計算科学による塗布型電子材料の構造・機能予測」委託研究、産業技術力強化法第１７条の適用を受ける特許出願

(73)【特許権者】

【識別番号】304036754

【氏名又は名称】国立大学法人山形大学

(74)【代理人】

【識別番号】110000383

【氏名又は名称】弁理士法人エビス国際特許事務所

(72)【発明者】

【氏名】松井 弘之

【審査官】嶋田 行志

(56)【参考文献】

【文献】特開２０００－０１１０２５（ＪＰ，Ａ）

【文献】特表２０１３－５０９１２６（ＪＰ，Ａ）

【文献】RandSpg:An open-source program for generatinf atomistic crystal structures with specific spacegroups，Computer Physics Communications，Elsevier B.V.，2016年12月21日，doi: 10.17632/v2nzgmpd37.1

【文献】Virtual hexagonal and multi-scale operator for fuzzy rank order testure classification using one-dimensional generalised Fourier analysis，IEEE Conference Proceedings，Vol. 2016, No. ICPR，2016年，pp. 2018-20240，doi: 10.1109/ICPR.2016.7899933

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｇ０１Ｎ ２３／００－Ｇ０１Ｎ ２３／２２７６

Ｇ１６Ｃ １０／００－Ｇ１６Ｃ ９９／００

ＪＳＴＰｌｕｓ／ＪＳＴ７５８０／ＪＳＴＣｈｉｎａ（ＪＤｒｅａｍＩＩＩ）

ＡＣＳ ＰＵＢＬＩＣＡＴＩＯＮＳ

Ｓｃｉｅｎｃｅ Ｄｉｒｅｃｔ

Ｓｃｉｔａｔｉｏｎ

ＡＰＳ Ｊｏｕｒｎａｌｓ

(57)【特許請求の範囲】

【請求項1】

単位格子の基本格子ベクトルを示す格子ベクトル変数と前記単位格子内の各原子の位置ベクトルを示す原子位置変数とを含む目的物質の初期結晶構造を、前記目的物質の前記格子ベクトル変数と前記原子位置変数との直積集合の範囲で順次変更し、
変更された前記格子ベクトル変数と前記原子位置変数の組の自由エネルギーを算出し、前記直積集合のうち下記式を満たす連続的な部分集合を定義域とする自由エネルギーの連続関数において、算出された前記自由エネルギーが極小を示すか否かに基づいて前記目的物質の前記格子ベクトル変数と前記原子位置変数の組を算出する結晶構造算出方法。

【数1】

ここで、Ｕ_ｓｔは結晶の格子エネルギー、ベクトルａ、ｂおよびｃは格子ベクトル変数、ベクトルｒ_ｉ（ｉ＝１、２、・・・、Ｎ）は原子位置変数、Ｎは単位格子内の原子の数を示す。

【請求項2】

前記直積集合の範囲で前記格子ベクトル変数を独立変数として順次変更し、
変更された前記格子ベクトル変数に対して前記格子エネルギーが極小を示す前記原子位置変数を従属変数として算出することにより、前記部分集合に含まれる前記格子ベクトル変数と前記原子位置変数の組を探索し、
前記部分集合に含まれる前記格子ベクトル変数と前記原子位置変数の組の前記自由エネルギーを順次算出して、前記自由エネルギーが極小を示す前記格子ベクトル変数と前記原子位置変数の組を算出する請求項１に記載の結晶構造算出方法。

【請求項3】

前記直積集合の範囲で前記格子ベクトル変数および前記原子位置変数を独立変数として順次変更し、
変更された前記格子ベクトル変数および前記原子位置変数の組が、ラグランジュの未定乗数法に基づいて、前記部分集合に含まれ且つ前記部分集合の中で極小の前記自由エネルギーを有するか否かを判定する請求項１に記載の結晶構造算出方法。

【請求項4】

前記自由エネルギーは、量子化学計算または古典力場計算から擬調和振動子近似によって算出されるフォノンの固有振動数に基づいて算出される請求項１～３のいずれか一項に記載の結晶構造算出方法。

【請求項5】

過去に得られた結晶構造情報に含まれる前記格子ベクトル変数および前記原子位置変数の値に基づいて前記自由エネルギーを算出するための経験式を予め算出し、前記経験式に基づいて前記自由エネルギーを算出される請求項１～３のいずれか一項に記載の結晶構造算出方法。

【請求項6】

請求項１～５のいずれか一項に記載の結晶構造算出方法の各ステップをコンピュータに実行させるための結晶構造算出プログラム。

【請求項7】

単位格子の基本格子ベクトルを示す格子ベクトル変数と前記単位格子内の各原子の位置ベクトルを示す原子位置変数とを含む目的物質の初期結晶構造を、前記格子ベクトル変数と前記原子位置変数との直積集合の範囲で順次変更する変数変更部と、
変更された前記格子ベクトル変数と前記原子位置変数の組の自由エネルギーを算出し、前記直積集合のうち下記式を満たす連続的な部分集合を定義域とする自由エネルギーの連続関数において、算出された前記自由エネルギーが極小を示すか否かに基づいて前記目的物質の前記格子ベクトル変数と前記原子位置変数の組を算出する極小条件判定部とを備える結晶構造算出装置。

【数1】

ここで、Ｕ_ｓｔは結晶の格子エネルギー、ベクトルａ、ｂおよびｃは格子ベクトル変数、ベクトルｒ_ｉ（ｉ＝１、２、・・・、Ｎ）は原子位置変数、Ｎは単位格子内の原子の数を示す。

【発明の詳細な説明】

【技術分野】

【0001】

この発明は、結晶構造算出方法、結晶構造算出プログラムおよび結晶構造算出装置に係り、特に、格子エネルギーおよび自由エネルギーに基づいて目的物質の結晶構造を算出する結晶構造算出方法、結晶構造算出プログラムおよび結晶構造算出装置に関する。

【背景技術】

【0002】

物質の結晶構造は、その物質の性質を決定づける主要因の一つである。例えば、有機半導体の移動度や医薬品の溶解性は物質の結晶構造に大きく依存する。この物質の結晶構造は、一般的に、単結晶Ｘ線回折測定によって決定されており、新規物質を合成する前にその結晶構造を予測することは困難であった。新規物質の合成前に計算機上で結晶構造を予測することが可能となれば、新規物質開発を大幅に効率化できる。

【0003】

そこで、結晶構造を計算機上で予測する技術として、例えば、非特許文献１には、格子エネルギーを最小とする結晶構造を探索する格子エネルギー最小化（ＬＥＭ；Ｌａｔｔｉｃｅ Ｅｎｅｒｇｙ Ｍｉｎｉｍｉｚａｔｉｏｎ）法が提案されている。このＬＥＭ法は、原子が静止しているものとみなして原子の運動エネルギーを無視した上で、原子の位置エネルギーの総和（格子エネルギーＵ_ｓｔ）が最小または極小となる結晶構造を探索するものである。

【先行技術文献】

【非特許文献】

【0004】

【文献】Ｒｅｐｏｒｔ ｏｎ ｔｈｅ ｓｉｘｔｈ ｂｌｉｎｄ ｔｅｓｔ ｏｆ ｏｒｇａｎｉｃ ｃｒｙｓｔａｌ ｓｔｒｕｃｔｕｒｅ ｐｒｅｄｉｃｔｉｏｎ ｍｅｔｈｏｄｓ、Ａｎｔｈｏｎｙ Ｍ． Ｒｅｉｌｌｙ ｅｔ ａｌ．、Ａｃｔａ Ｃｒｙｓｔ．（２０１６）．Ｂ７２、Ｐ．４３９－Ｐ．４５９

【発明の概要】

【発明が解決しようとする課題】

【0005】

しかしながら、ＬＥＭ法は、原子の熱振動や零点振動を無視しているため、熱膨張や温度による構造転移などの熱的現象を再現することができず、現実の結晶構造との間に差異が生じるおそれがある。また、現実の実験では数個の結晶構造（多形）しか得られない物質をＬＥＭ法でシミュレーションすると、数百個以上の結晶構造（多形）が得られる場合があり、現実の実験結果からかけ離れた結晶構造が算出されるといった問題もあった。

【0006】

この発明は、このような従来の問題点を解消するためになされたもので、目的物質の結晶構造を高精度に算出する結晶構造算出方法、結晶構造算出プログラムおよび結晶構造算出装置を提供することを目的とする。

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

【0007】

この発明に係る結晶構造算出方法は、単位格子の基本格子ベクトルを示す格子ベクトル変数と単位格子内の各原子の位置ベクトルを示す原子位置変数とを含む目的物質の初期結晶構造を、目的物質の格子ベクトル変数と原子位置変数との直積集合の範囲で順次変更し、変更された格子ベクトル変数と原子位置変数の組の自由エネルギーを算出し、直積集合のうち下記式（１）を満たす連続的な部分集合を定義域とする自由エネルギーの連続関数において、算出された自由エネルギーが極小を示すか否かに基づいて目的物質の格子ベクトル変数と原子位置変数の組を算出するものである。

【0008】

ここで、直積集合の範囲で格子ベクトル変数を独立変数として順次変更し、変更された格子ベクトル変数に対して格子エネルギーが極小を示す原子位置変数を従属変数として算出することにより、部分集合に含まれる格子ベクトル変数と原子位置変数の組を探索し、部分集合に含まれる格子ベクトル変数と原子位置変数の組の自由エネルギーを順次算出して、自由エネルギーが極小を示す格子ベクトル変数と原子位置変数の組を算出することが好ましい。

【0009】

また、直積集合の範囲で格子ベクトル変数および原子位置変数を独立変数として順次変更し、変更された格子ベクトル変数および原子位置変数の組が、ラグランジュの未定乗数法に基づいて、部分集合に含まれ且つ部分集合の中で極小の自由エネルギーを有するか否かを判定することができる。

【0010】

また、自由エネルギーは、量子化学計算または古典力場計算から擬調和振動子近似によって算出されるフォノンの固有振動数に基づいて算出することができる。

【0011】

また、過去に得られた結晶構造情報に含まれる格子ベクトル変数および原子位置変数の値に基づいて自由エネルギーを算出するための経験式を予め算出し、経験式に基づいて自由エネルギーを算出することもできる。

【0012】

この発明に係る結晶構造算出プログラムは、上記のいずれかに記載の結晶構造算出方法の各ステップをコンピュータに実行させるものである。

【0013】

この発明に係る結晶構造算出装置は、単位格子の基本格子ベクトルを示す格子ベクトル変数と単位格子内の各原子の位置ベクトルを示す原子位置変数とを含む目的物質の初期結晶構造を、格子ベクトル変数と原子位置変数との直積集合の範囲で順次変更する変数変更部と、変更された格子ベクトル変数と原子位置変数の組の自由エネルギーを算出し、直積集合のうち下記式（１）を満たす連続的な部分集合を定義域とする自由エネルギーの連続関数において、算出された自由エネルギーが極小を示すか否かに基づいて目的物質の格子ベクトル変数と原子位置変数の組を算出する極小条件判定部とを備えるものである。

【発明の効果】

【0014】

この発明によれば、格子ベクトル変数と原子位置変数の直積集合のうち下記式（１）を満たす連続的な部分集合を定義域とする自由エネルギーの連続関数において、格子ベクトル変数と原子位置変数を順次変更して算出される自由エネルギーが極小を示すか否かに基づいて目的物質の格子ベクトル変数と原子位置変数の組を算出するので、目的物質の結晶構造を高精度に算出する結晶構造算出方法、結晶構造算出プログラムおよび結晶構造算出装置を提供することが可能となる。

【図面の簡単な説明】

【0015】

【図1】この発明の実施の形態１に係る結晶構造算出装置の構成を示すブロック図である。

【図2】実施の形態１に係る結晶構造算出方法を示すフローチャートである。

【図3】実施の形態１において目的物質の格子ベクトル変数と原子位置変数の組を探索する様子を示す図である。

【図4】格子エネルギーの極小点を探索する様子を示す図である。

【図5】従来のＬＥＭ法で目的物質の格子ベクトル変数と原子位置変数の組を探索する様子を示す図である。

【図6】結晶構造候補を算出する様子を示す図である。

【図7】自由エネルギーの連続関数において極小点を示す図である。

【図8】実施の形態２に係る結晶構造算出装置の構成を示すブロック図である。

【図9】実施の形態２に係る結晶構造算出方法を示すフローチャートである。

【図10】実施の形態２において目的物質の格子ベクトル変数と原子位置変数の組を探索する様子を示す図である。

【図11】実施の形態３に係る結晶構造算出装置の構成の要部を示すブロック図である。

【図12】実施の形態４に係る結晶構造算出装置の構成の要部を示すブロック図である。

【図13】実施の形態５に係る結晶構造算出装置の構成の要部を示すブロック図である。

【図14】実施の形態１の結晶構造算出方法を用いて温度に対する格子定数の変化を算出した結果を示すグラフである。

【図15】実施の形態１の結晶構造算出方法を用いて格子定数に対する自由エネルギーの変化を算出した結果を示すグラフである。

【発明を実施するための形態】

【0016】

以下、この発明の実施の形態を添付図面に基づいて説明する。
実施の形態１
図１に、この発明の実施の形態１に係る結晶構造算出装置１の構成を示す。結晶構造算出装置１は、初期構造生成部２を有し、この初期構造生成部２に格子エネルギー算出部３、極小条件判定部４、自由エネルギー算出部５、極小条件判定部６および出力部７が順次接続されている。また、極小条件判定部４が原子位置変数変更部８に接続されると共に極小条件判定部６が格子ベクトル変数変更部９に接続され、原子位置変数変更部８および格子ベクトル変数変更部９がそれぞれ格子エネルギー算出部３に接続されている。そして、原子位置変数変更部８および格子ベクトル変数変更部９により、変数変更部１０が形成されている。
また、初期構造生成部２、格子エネルギー算出部３、極小条件判定部４、自由エネルギー算出部５、極小条件判定部６、原子位置変数変更部８および格子ベクトル変数変更部９に制御部１１が接続され、この制御部１１に操作部１２と格納部１３が接続されている。また、結晶構造算出装置１の出力部７に表示部１４が接続されている。

【0017】

初期構造生成部２は、使用者から入力される目的物質の分子構造情報に基づいて、単位格子の基本格子ベクトルを示す格子ベクトル変数と単位格子内の各原子の位置ベクトルを示す原子位置変数とを含む目的物質の初期結晶構造を生成する。
ここで、分子構造情報は、目的物質の分子構造に関する情報で、例えば、目的物質の分子に含まれる原子の種類、原子の数および原子間の化学結合の種類などが挙げられる。

【0018】

原子位置変数変更部８は、初期構造生成部２で生成された目的物質の初期結晶構造において、目的物質の格子ベクトル変数と原子位置変数との直積集合の範囲で格子ベクトル変数を固定しつつ原子位置変数を順次変更する。ここで、直積集合は、目的物質の格子ベクトル変数と原子位置変数の全ての組み合わせを示し、分子構造情報に基づいて定義することができる。

【0019】

格子エネルギー算出部３は、使用者から入力された目的物質の格子エネルギー算出情報に基づいて、初期構造生成部２、原子位置変数変更部８または格子ベクトル変数変更部９で設定される原子位置変数と格子ベクトル変数の組を用いて格子エネルギーを算出する。
ここで、格子エネルギー算出情報は、格子エネルギーを算出するための情報で、例えば、分子構造情報および力場パラメータなどが挙げられ、さらに環境圧力（結晶が生成される環境の圧力）を含むことが好ましい。なお、力場パラメータとしては、例えば、格子エネルギーにおける共有結合エネルギーをｆ（ｒ）＝Ａ（ｒ－ｒ_０）^２で表したとき（Ａはばね定数、ｒは原子間の距離、ｒ_０は平衡状態の原子間の距離）のＡおよびｒ_０などが挙げられる。

【0020】

極小条件判定部４は、格子ベクトル変数を固定しつつ原子位置変数を順次変更して算出される格子エネルギーの連続関数において、格子エネルギー算出部３で算出された格子エネルギーが極小を示すか否かを判定する。これにより、目的物質の格子ベクトル変数と原子位置変数の直積集合のうち、下記式（１）を満たす格子ベクトル変数と原子位置変数の組が連続的に続く部分集合に、原子位置変数変更部８で設定された格子ベクトル変数と原子位置変数の組が含まれるか否かを判定することになる。

【0021】

【数1】

【0022】

ここで、Ｕ_ｓｔは結晶の格子エネルギー、ベクトルａ、ｂおよびｃは格子ベクトル変数、ベクトルｒ_ｉ（ｉ＝１、２、・・・、Ｎ）は原子位置変数、Ｎは単位格子内の原子の数を示す。

【0023】

格子ベクトル変数変更部９は、極小条件判定部４で格子エネルギーが極小を示すと判定された格子ベクトル変数と原子位置変数の組について、格子ベクトル変数の値を直積集合の範囲で順次変更する。
このようにして、格子ベクトル変数変更部９が、目的物質の格子ベクトル変数と原子位置変数の直積集合の範囲で格子ベクトル変数を独立変数として順次変更し、極小条件判定部４が、変更された格子ベクトル変数に対して格子エネルギーが極小を示す原子位置変数を従属変数として算出することにより、部分集合に含まれる格子ベクトル変数と原子位置変数の組が探索されることになる。

【0024】

自由エネルギー算出部５は、使用者から入力された目的物質の自由エネルギー算出情報に基づいて、極小条件判定部４で格子エネルギーが極小を示すと判定された、すなわち部分集合に含まれると判定された格子ベクトル変数と原子位置変数の組で自由エネルギーを算出する。
ここで、自由エネルギー算出情報は、自由エネルギーを算出するための情報で、例えば、分子構造情報、力場パラメータおよび環境温度（結晶が生成される環境の温度）などが挙げられ、さらに環境圧力（結晶が生成される環境の圧力）を含むことが好ましい。

【0025】

極小条件判定部６は、部分集合に含まれる格子ベクトル変数と原子位置変数の組から算出される自由エネルギーの連続関数において、自由エネルギー算出部５で順次算出される自由エネルギーが極小を示すか否かに基づいて目的物質の結晶構造を示す格子ベクトル変数と原子位置変数の組を算出する。

【0026】

出力部７は、極小条件判定部６で判定された目的物質の結晶構造を表示部１４に出力する。
制御部１１は、使用者により操作部１２から入力される指示に基づいて結晶構造算出装置１の各部の制御を行う。

【0027】

操作部１２は、使用者からの指令を入力するためのもので、ボタン、タッチパネル、キーボード、マウス、トラックボール等から形成することができる。
格納部１３は、動作プログラム等を格納するもので、ハードディスク、フレキシブルディスク、ＭＯ、ＭＴ、ＲＡＭ、ＣＤ－ＲＯＭ、ＤＶＤ－ＲＯＭ、ＳＤカード、ＣＦカード、ＵＳＢメモリ等の記録メディア、またはサーバ等を用いることができる。

【0028】

なお、初期構造生成部２、格子エネルギー算出部３、極小条件判定部４、自由エネルギー算出部５、極小条件判定部６、原子位置変数変更部８、格子ベクトル変数変更部９および変数変更部１０は、ＣＰＵと、ＣＰＵに各種の処理を行わせるための動作プログラムから構成されるが、それらをデジタル回路で構成してもよい。
表示部１４は、ディスプレイ装置を含み、出力部７から出力される目的物質の結晶構造を表示するものである。

【0029】

次に、図２に示すフローチャートを参照して、実施の形態１の結晶構造算出方法について説明する。なお、本発明の結晶構造算出プログラムは、結晶構造算出装置１が有するコンピュータに以下の結晶構造算出方法を実行させるプログラムである。

【0030】

まず、初期構造生成部２に使用者から目的物質の分子構造情報が入力されると、初期構造生成部２は、ステップＳ１で、分子構造情報に基づいて、格子ベクトル変数と原子位置変数の組を選択して目的物質の初期結晶構造を生成する。一般的に、結晶構造は、単位格子の３つの基本格子ベクトルのｘ、ｙ、ｚ成分を示す格子ベクトル変数（自由度が６）と、単位格子内の各原子の位置ベクトルのｘ、ｙ、ｚ成分を示す原子位置変数（自由度が３Ｎ）とで表すことができ、初期構造生成部２は、格子ベクトル変数と原子位置変数の組を選択することで目的物質の初期結晶構造を生成することができる。
なお、格子ベクトル変数は、回転の自由度を除外しているが、回転の自由度を含めてもよい。また、原子位置変数においてＮは単位格子内の原子の個数を示す。また、分子構造情報は予め保存してもよい。

【0031】

ここで、初期構造生成部２は、分子構造情報に基づいて算出される目的物質の格子ベクトル変数と原子ベクトル変数の直積集合の範囲で、格子ベクトル変数と原子位置変数の組を選択する。
例えば、初期構造生成部２は、図３に示すように、目的物質の格子ベクトル変数と原子ベクトル変数の直積集合Ｔ１の範囲で格子ベクトル変数と原子位置変数の組Ｖ１をランダムに選択することができる。ここで、等高線Ｌは、直積集合Ｔ１における格子エネルギーの分布を示す。初期構造生成部２は、選択した格子ベクトル変数と原子位置変数の組Ｖ１を格子エネルギー算出部３に出力する。
なお、初期構造生成部２は、分子構造情報に基づいて選択範囲などを限定して格子ベクトル変数と原子位置変数の組Ｖ１を選択することが好ましい。

【0032】

初期構造生成部２で選択された格子ベクトル変数と原子位置変数の組Ｖ１が格子エネルギー算出部３に入力されると、格子エネルギー算出部３は、ステップＳ２で、使用者から入力された目的物質の格子エネルギー算出情報に基づいて、格子ベクトル変数と原子位置変数の組Ｖ１を用いて格子エネルギーを算出する。算出された格子エネルギーは、格子エネルギー算出部３から極小条件判定部４に出力される。格子エネルギーとしては、例えば、化学結合した２原子間の結合距離エネルギー、化学結合した３原子間の結合角エネルギー、化学結合した４原子間の二面角エネルギーおよび化学結合していない２原子の分子間力エネルギーなどが挙げられる。
なお、格子エネルギー算出情報は、使用者が分子構造情報と共に入力してもよく、予め保存してもよい。

【0033】

格子エネルギー算出部３で算出された格子エネルギーが極小条件判定部４に入力されると、極小条件判定部４は、ステップＳ３で、格子ベクトル変数を固定しつつ原子位置変数を順次変更して算出される格子エネルギーの連続関数において、格子エネルギー算出部３で算出された格子エネルギーが極小を示すか否かを判定する。
ここで、格子ベクトル変数と原子位置変数の組Ｖ１から算出された格子エネルギーは最初に入力されたものであるため極小を示すか否かを判定できず、極小条件判定部４は、極小を示さないと判定して、その判定結果を原子位置変数変更部８に出力する。

【0034】

極小条件判定部４の判定結果が原子位置変数変更部８に入力されると、ステップＳ４に進んで、原子位置変数変更部８が、格子ベクトル変数と原子位置変数の組Ｖ１から原子位置変数のみを直積集合Ｔ１の範囲で変更して格子ベクトル変数と原子位置変数の組Ｖ２とし、この格子ベクトル変数と原子位置変数の組Ｖ２を格子エネルギー算出部３に出力する。

【0035】

原子位置変数変更部８で生成された格子ベクトル変数と原子位置変数の組Ｖ２が格子エネルギー算出部３に入力されると、格子エネルギー算出部３は、ステップＳ２に戻り、格子ベクトル変数と原子位置変数の組Ｖ２を用いて格子エネルギーを算出する。
このようにして、原子位置変数変更部８が、格子ベクトル変数を固定しつつ原子位置変数の値を順次変更し、その格子ベクトル変数と原子位置変数の組Ｖ２，Ｖ３およびＶ４の格子エネルギーが格子エネルギー算出部３から極小条件判定部４に入力される。

【0036】

極小条件判定部４は、格子エネルギー算出部３から順次入力された、これまでの格子エネルギーに基づいて、ステップＳ３で、格子エネルギー算出部３で算出された格子エネルギーが格子エネルギーの連続関数において極小を示すか否かを判定する。

【0037】

例えば、図４に示すように、極小条件判定部４は、これまでに格子エネルギー算出部３から格子エネルギーＵ１～Ｕ４が順次入力された場合には、その格子エネルギーＵ１～Ｕ４において最も小さな値を示す格子エネルギーＵ３を極小点として判定する。ここで、格子エネルギーＵ１～Ｕ４は、原子位置変数を一定の間隔で変えて算出する必要はなく、原子位置変数を異なる間隔で変えて算出することができる。例えば、格子エネルギーＵ２は格子エネルギーＵ１に対して原子位置変数の値を大きく空けて算出し、格子エネルギーＵ３およびＵ４は原子位置変数の値を小さく変えて算出することができる。
なお、極小条件判定部４は、格子エネルギー算出部３から順次入力される格子エネルギーの微分、すなわち上記式（１）を算出し、その微分がゼロを示す格子エネルギーＵ３を極小点として判定することもできる。

【0038】

このようにして、極小条件判定部４は、図３に示すように、直積集合Ｔ１のうち、上記式（１）を満たす格子ベクトル変数と原子位置変数の組が連続的に続く部分集合Ｔ２に、格子ベクトル変数と原子位置変数の組Ｖ１～Ｖ４が含まれるか否かを判定する。なお、極小条件判定部４は、従来の格子エネルギー最小化（ＬＥＭ；Ｌａｔｔｉｃｅ Ｅｎｅｒｇｙ Ｍｉｎｉｍｉｚａｔｉｏｎ）法と同様にして、原子の熱振動および零点振動を無視して格子エネルギーの極小点を探索する。
これにより、極小条件判定部４は、格子ベクトル変数を固定したときに、原子振動の中心を表す原子位置変数を探索、すなわち自由エネルギーが算出可能な部分集合Ｔ２に含まれる原子位置変数を探索することができる。

【0039】

そして、極小条件判定部４において格子エネルギーが極小を示すと判定されるまで、原子位置変数変更部８において原子位置変数が順次変更される。極小条件判定部４は、例えば格子ベクトル変数と原子位置変数の組Ｖ３の格子エネルギーが極小を示すと判定すると、格子ベクトル変数と原子位置変数の組Ｖ３を自由エネルギー算出部５に出力する。

【0040】

これにより、自由エネルギー算出部５が、ステップＳ５で、使用者から入力された目的物質の自由エネルギー算出情報に基づいて、格子ベクトル変数と原子位置変数の組Ｖ３を用いて自由エネルギーを算出する。なお、自由エネルギー算出情報は、使用者が分子構造情報と共に入力してもよく、予め保存してもよい。
ここで、自由エネルギー算出部５は、ヘルムホルツの自由エネルギーＦ＝Ｕ_ｓｔ＋Ｕ_ｐｈ－ＴＳ_ｐｈまたはギブズの自由エネルギーＧ＝Ｆ＋ＰＶに基づいて算出することができる。ここで、Ｕ_ｐｈはフォノン（原子振動）の内部エネルギー、Ｔは絶対温度、Ｓ_ｐｈはフォノンのエントロピー、Ｐは圧力、Ｖは体積である。擬調和振動子近似に基づけば、Ｕ_ｐｈは、フォノンの固有振動の組｛ν_ｉ｝を用いて、Ｕ_ｐｈ＝Σ_ｉ（ｈν_ｉ／２）ｃｏｔｈ（ｈν_ｉ／２ｋ_ＢＴ）と表され、Ｓ_ｐｈはＳ_ｐｈ＝Σ_ｉ｛－ｋ・ｌｎ［２ｓｉｎｈ（ｈν_ｉ／２ｋ_ＢＴ）］＋（ｈν_ｉ／２Ｔ）ｃｏｔｈ（ｈν_ｉ／２ｋ_ＢＴ）｝と表される。

【0041】

自由エネルギー算出部５は、算出された自由エネルギーを極小条件判定部６に出力する。自由エネルギー算出部５で算出された自由エネルギーが極小条件判定部６に入力されると、極小条件判定部６は、ステップＳ６で、部分集合Ｔ２に含まれる格子ベクトル変数と原子位置変数の組で算出される自由エネルギーの連続関数、すなわち部分集合Ｔ２を定義域とする自由エネルギーの連続関数において、自由エネルギー算出部５で算出された自由エネルギーが極小を示すか否かを判定する。
極小条件判定部６は、自由エネルギー算出部５で算出された自由エネルギーが極小を示さないと判定した場合には、その判定結果を格子ベクトル変数変更部９に出力する。

【0042】

極小条件判定部６の判定結果が格子ベクトル変数変更部９に入力されると、ステップＳ７に進んで、格子ベクトル変数変更部９が、格子ベクトル変数と原子位置変数の組Ｖ３から格子ベクトル変数のみを直積集合Ｔ１の範囲で変更して格子ベクトル変数と原子位置変数の組Ｖ５とし、この格子ベクトル変数と原子位置変数の組Ｖ５を格子エネルギー算出部３に出力する。

【0043】

このようにして、格子ベクトル変数変更部９で生成された格子ベクトル変数と原子位置変数の組Ｖ５が格子エネルギー算出部３に入力されると、格子エネルギー算出部３は、ステップＳ２に戻り、格子ベクトル変数と原子位置変数の組Ｖ５を用いて格子エネルギーを算出する。続いて、極小条件判定部４において、上記と同様に、格子エネルギー算出部３で算出された格子エネルギーが極小を示すと判定されるまで、格子ベクトル変数を固定しつつ原子位置変数を順次変更する。
すなわち、直積集合Ｔ１の範囲で格子ベクトル変数を独立変数として順次変更し、変更された格子ベクトル変数に対して格子エネルギーが極小を示す原子位置変数を従属変数として算出する。これにより、部分集合Ｔ２に含まれる格子ベクトル変数と原子位置変数の組を確実に算出することができる。

【0044】

例えば、極小条件判定部４において格子ベクトル変数と原子位置変数の組Ｖ６の格子エネルギーが極小を示すと判定されると、自由エネルギー算出部５が、格子ベクトル変数と原子位置変数の組Ｖ６の自由エネルギーを算出する。
極小条件判定部６は、自由エネルギー算出部５から順次入力される格子ベクトル変数と原子位置変数の組Ｖ３およびＶ６の自由エネルギーに基づいて、ステップＳ６で、自由エネルギー算出部５で算出された格子ベクトル変数と原子位置変数の組Ｖ６の自由エネルギーが自由エネルギーの連続関数において極小を示すか否かを判定する。

【0045】

極小条件判定部６は、例えば、格子ベクトル変数と原子位置変数の組Ｖ６で自由エネルギーが極小を示すと判定した場合には、その格子ベクトル変数と原子位置変数の組Ｖ６を最適化結晶構造とする。
なお、極小条件判定部６は、極小条件判定部４と同様に、部分集合Ｔ２に含まれる格子ベクトル変数と原子位置変数の組を一定の間隔で変える必要はなく、部分集合Ｔ２に含まれる格子ベクトル変数と原子位置変数の組を異なる間隔で変えて自由エネルギーの極小を判定することができる。また、極小条件判定部６は、自由エネルギー算出部５から順次入力される自由エネルギーの微分がゼロを示す自由エネルギーを極小点として判定することもできる。

【0046】

従来、ＬＥＭ法では、図５に示すように、格子エネルギーの等高線Ｌにおいて値が極小を示す離散的な格子ベクトル変数と原子位置変数の組Ｖ１’、Ｖ２’およびＶ３’を目的物質の結晶構造として算出する。ここで、ＬＥＭ法は、原子が静止していると見なしているため、温度の概念が存在せず、概ね０Ｋにおける結晶構造しか予測できない（零点振動が無視されているため０Ｋでもない）といった問題がある。このとき、格子ベクトル変数と原子位置変数の組Ｖ１’、Ｖ２’およびＶ３’に対して自由エネルギーを算出して目的物質の結晶構造を決定した場合でも、格子ベクトル変数と原子位置変数の組Ｖ１’、Ｖ２’およびＶ３’が離散的であるため、例えば熱膨張のような連続的な構造変化を再現することはできず、過剰数の多形の問題も解決することができなかった。

【0047】

本発明は、自由エネルギーの連続関数において極小を示す格子ベクトル変数と原子位置変数の組Ｖ６を探索する自由エネルギー最小化（ＦＥＭ；Ｆｒｅｅ Ｅｎｅｒｇｙ Ｍｉｎｉｍｉｚａｔｉｏｎ）法である。ここで、熱力学第二法則によれば、結晶が等温等圧条件下で熱平衡状態であることと、結晶のギブズの自由エネルギーが極小であることは等価である。従って、ギブズの自由エネルギーが極小となる結晶構造を探索すれば、熱振動および零点振動を考慮したより現実に近い結晶構造を算出することができる。
具体的には、連続的な温度変化におけるフォノン（原子振動）の大きさを考慮して目的物質の結晶構造を高精度に算出することができる。また、結晶構造の予測では、室温での結晶構造が重要であるため、０Ｋではなく、有限温度における結晶構造の算出は極めて重要となる。本発明では、熱膨張の再現および温度による構造転移などを高精度に再現することができる。さらに、本発明では、有限温度で安定して存在可能な結晶構造のみが得られるため、過剰数の結晶構造（多形）が算出されることを抑制することができる。

【0048】

このような結晶構造の算出方法は並列的にまたは連続的に行われる。すなわち、初期構造生成部２が、図３に示すように、直積集合Ｔ１の範囲で新たな格子ベクトル変数と原子位置変数の組Ｖ１ａをランダムに選択し、上記のステップＳ２～Ｓ７を繰り返すことにより、極小条件判定部６が、部分集合Ｔ２において自由エネルギーが極小を示す格子ベクトル変数と原子位置変数の組Ｖ２ａを最適化結晶構造として順次算出することができる。
例えば、図６に示すように、初期構造生成部２において１０００個の初期結晶構造がランダムに生成されると、極小条件判定部６において自由エネルギーが極小を示す１０００個の最適化結晶構造が算出されることになる。

【0049】

そして、極小条件判定部６は、１０００個の最適化結晶構造のうち同じ結晶構造を有するものを１つにまとめて、例えば９つの目的物質の結晶構造候補Ｒ１，Ｒ２・・・Ｒ９を算出する。これにより、図７に示すように、部分集合Ｔ２を定義域とする自由エネルギーの連続関数ｆにおいて極小点を示す結晶構造候補Ｒ１～Ｒ９をそれぞれ算出することができる。このとき、結晶構造候補Ｒ１～Ｒ９の数が多い場合には、極小条件判定部６は、算出された結晶構造候補Ｒ１～Ｒ９のうち自由エネルギーが所定の値以下のもの、例えば結晶構造候補Ｒ４，Ｒ５およびＲ６を目的物質の結晶構造として算出する。
このように、極小条件判定部６は、結晶構造候補Ｒ１～Ｒ９から自由エネルギーの値が大きいものを除くため、目的物質の結晶構造Ｒ４，Ｒ５およびＲ６をより高精度に算出することができる。

【0050】

極小条件判定部６は、算出した目的物質の結晶構造Ｒ４，Ｒ５およびＲ６を出力部７を介して表示部１４に出力し、目的物質の結晶構造Ｒ４，Ｒ５およびＲ６が表示部１４に表示される。これにより、例えば、表示部１４に表示された目的物質の結晶構造Ｒ４，Ｒ５およびＲ６に基づいて、実験で得られたＸ線解析データなどを高精度に解析することができる。

【0051】

本実施の形態によれば、格子ベクトル変数と原子位置変数の組の自由エネルギーを算出し、直積集合Ｔ１のうち上記式（１）を満たす連続的な部分集合Ｔ２を定義域とする自由エネルギーの連続関数ｆにおいて、算出された自由エネルギーが極小を示すか否かを判定するため、目的物質の結晶構造を高精度に算出することができる。
具体的には、直積集合Ｔ１の範囲で格子ベクトル変数を独立変数として順次変更し、変更された格子ベクトル変数に対して格子エネルギーが極小を示す原子位置変数を従属変数として算出するため、部分集合Ｔ２に含まれる格子ベクトル変数と原子位置変数の組Ｖ３およびＶ６を確実に探索することができる。さらに、部分集合Ｔ２に含まれる格子ベクトル変数と原子位置変数の組の自由エネルギーを順次算出して、その自由エネルギーが極小を示す格子ベクトル変数と原子位置変数の組Ｖ６を算出するため、目的物質の結晶構造をより高精度に算出することができる。

【0052】

実施の形態２
上記の実施の形態１では、極小条件判定部６が、部分集合Ｔ２に含まれる格子ベクトル変数と原子位置変数の組Ｖ３およびＶ６の自由エネルギーを順次算出してその自由エネルギーが極小を示す格子ベクトル変数と原子位置変数の組Ｖ６を算出したが、部分集合Ｔ２を定義域とする自由エネルギーの連続関数ｆにおいて格子ベクトル変数と原子位置変数の組の自由エネルギーが極小を示すか否かに基づいて目的物質の結晶構造を算出すればよく、これに限られるものではない。

【0053】

例えば、図８に示すように、実施の形態１の変数変更部１０に換えて変数変更部２２を配置すると共に極小条件判定部４および６に換えて極小条件判定部２３を配置し、さらに格子エネルギー算出部３および自由エネルギー算出部５を除いてずれ方向算出部２１を新たに配置することができる。ここで、ずれ方向算出部２１、変数変更部２２および極小条件判定部２３は、初期構造生成部２と出力部７の間に順次接続して配置されている。

【0054】

ずれ方向算出部２１は、格子エネルギーおよび自由エネルギーに基づいて格子ベクトル変数と原子位置変数の組を変更するずれ方向を算出する。具体的には、格子ベクトル変数と原子位置変数の組、格子エネルギー算出情報および自由エネルギー算出情報を用いて下記式（２）で表されるラグランジュの未定乗数法に基づいてラグランジュ関数Ｌを算出し、格子ベクトル変数と原子位置変数の値を複数方向に少しだけ変更したときに、ラグランジュ関数Ｌが最も小さくなるずれ方向を算出する。

【0055】

【数2】

【0056】

ここで、Ｌはラグランジュ関数、Ｇは目的物質の自由エネルギー、ベクトルａ、ｂおよびｃは格子ベクトル変数、ベクトルｒ_ｉ（ｉ＝１、２、・・・、Ｎ）は原子位置変数、ｘ_ｉ、ｙ_ｉ、ｚ_ｉはベクトルｒ_ｉのｘ、ｙ、ｚ成分、Ｎは目的物質の単位格子内の原子の数、λはラグランジュの未定乗数、Ｕ_ｓｔは目的物質の格子エネルギーを示す。

【0057】

変数変更部２２は、ずれ方向算出部２１で算出されたずれ方向に基づいて格子ベクトル変数と原子位置変数の値を、格子ベクトル変数と原子位置変数の直積集合Ｔ１の範囲で変更する。

【0058】

極小条件判定部２３は、変数変更部２２で変更された格子ベクトル変数と原子位置変数の組が、上記式（２）に基づいて、部分集合Ｔ２に含まれ且つ部分集合Ｔ２の中で極小の自由エネルギーを有するか否かを判定する。すなわち、極小条件判定部２３は、上記式（２）のラグランジュ関数Ｌが停留点（極小点、極大点および変曲点など）を示す格子ベクトル変数と原子位置変数の組を探索する。そして、極小条件判定部２３は、部分集合Ｔ２に含まれ且つ部分集合Ｔ２の中で極小の自由エネルギーを有する格子ベクトル変数と原子位置変数の組を目的物質の結晶構造として算出する。

【0059】

次に、図９に示すフローチャートを参照して、実施の形態２の結晶構造算出方法について説明する。
まず、実施の形態１と同様に、初期構造生成部２に使用者から目的物質の分子構造情報が入力されると、初期構造生成部２は、ステップＳ１で、分子構造情報に基づいて、格子ベクトル変数と原子位置変数の組を選択して目的物質の初期結晶構造を生成する。例えば、初期構造生成部２は、図１０に示すように、目的物質の格子ベクトル変数と原子ベクトル変数の直積集合Ｔ１の範囲で格子ベクトル変数と原子位置変数の組Ｖ１をランダムに選択する。初期構造生成部２は、選択した格子ベクトル変数と原子位置変数の組Ｖ１をずれ方向算出部２１に出力する。

【0060】

初期構造生成部２で選択された格子ベクトル変数と原子位置変数の組Ｖ１がずれ方向算出部２１に入力されると、ステップＳ２１で、ずれ方向算出部２１が、上記式（２）に基づいて、格子ベクトル変数と原子位置変数の組Ｖ１を変更するずれ方向を算出する。具体的には、ずれ方向算出部２１は、格子ベクトル変数と原子位置変数の組Ｖ１に対して複数方向に少しだけ値を変更してラグランジュ関数Ｌを算出し、そのラグランジュ関数Ｌが最も小さくなる方向をずれ方向Ｎ１に決定する。ずれ方向算出部２１は、算出したずれ方向Ｎ１を変数変更部２２に出力する。

【0061】

ずれ方向算出部２１で算出されたずれ方向Ｎ１が変数変更部２２に入力されると、変数変更部２２は、ステップＳ２２で、格子ベクトル変数と原子位置変数の組Ｖ１を直積集合Ｔ１の範囲でずれ方向Ｎ１に向かって変更する。例えば、変数変更部２２は、格子ベクトル変数と原子位置変数の組Ｖ１を格子ベクトル変数と原子位置変数の組Ｖ２に大きく変更する。変数変更部２２は、変更した格子ベクトル変数と原子位置変数の組Ｖ２を極小条件判定部２３に出力する。
このように、変数変更部２２は、ラグランジュ関数Ｌを用いて算出されたずれ方向Ｎ１に基づいて格子ベクトル変数と原子位置変数の組Ｖ１を変更するため、格子エネルギーおよび自由エネルギーがより低い格子ベクトル変数と原子位置変数の組Ｖ２を効率よく探索することができる。

【0062】

変数変更部２２で変更された格子ベクトル変数と原子位置変数の組Ｖ２が極小条件判定部２３に入力されると、極小条件判定部２３は、ステップＳ２３で、格子ベクトル変数と原子位置変数の組Ｖ２を用いて上記式（２）に基づいて、ラグランジュ関数Ｌを算出し、そのラグランジュ関数Ｌが停留点を示すか否かを判定する。これにより、格子ベクトル変数と原子位置変数の組Ｖ２が、部分集合Ｔ２に含まれ且つ部分集合Ｔ２の中で極小の自由エネルギーを有するか否かを判定することになる。
極小条件判定部２３は、格子ベクトル変数と原子位置変数の組Ｖ２がラグランジュ関数Ｌの停留点を示さないと判定した場合には、その判定結果をずれ方向算出部２１に出力する。

【0063】

極小条件判定部２３の判定結果がずれ方向算出部２１に入力されると、ステップＳ２１に戻り、ずれ方向算出部２１が、上記式（２）に基づいて、格子ベクトル変数と原子位置変数の組Ｖ２を変更するずれ方向Ｎ２を算出する。
このようにして、極小条件判定部２３で格子ベクトル変数と原子位置変数の組を用いて算出されるラグランジュ関数Ｌの値がこれまでに算出されたラグランジュ関数Ｌにおいて停留点を示すと判定されるまで、変数変更部２２が、ずれ方向Ｎ２およびＮ３に向かって格子ベクトル変数と原子位置変数の組Ｖ２およびＶ３を順次変更する。なお、極小条件判定部２３は、格子ベクトル変数と原子位置変数の組を一定の間隔で変える必要はなく、格子ベクトル変数と原子位置変数の組を異なる間隔で変えてラグランジュ関数Ｌの停留点を判定することができる。

【0064】

そして、格子ベクトル変数と原子位置変数の組Ｖ４のラグランジュ関数Ｌの値が、これまでに算出されたラグランジュ関数Ｌにおいて停留点を示すと判定された場合には、極小条件判定部２３は、その格子ベクトル変数と原子位置変数の組Ｖ４を最適化結晶構造とする。
このように、極小条件判定部２３が、ラグランジュの未定乗数法に基づいて判定するため、部分集合Ｔ２を定義域とする自由エネルギーの連続関数ｆにおいて極小を示す格子ベクトル変数と原子位置変数の組Ｖ４を効率よく探索することができる。

【0065】

このような結晶構造の算出方法は、実施の形態１と同様に、並列的にまたは連続的に行われる。そして、極小条件判定部２３が、最適化結晶構造を順次判定し、その複数の最適化結晶構造の中から目的物質の結晶構造が抽出される。

【0066】

本実施の形態によれば、極小条件判定部２３が、ラグランジュの未定乗数法に基づいて、格子ベクトル変数および原子位置変数の組Ｖ１～Ｖ４が部分集合Ｔ２に含まれ且つ部分集合Ｔ２の中で極小の自由エネルギーを有するか否かを判定するため、その条件を満たす格子ベクトル変数と原子位置変数の組Ｖ４を効率よく探索することができる。

【0067】

実施の形態３
上記の実施の形態１および２において、自由エネルギーは、量子化学計算または古典力場計算に基づいて算出することが好ましい。
例えば、図１１に示すように、実施の形態１にパラメータ保存部３１を新たに配置することができる。

【0068】

パラメータ保存部３１は、自由エネルギー算出部５に接続され、量子化学計算または古典力場計算に基づいて算出された力場パラメータを保存する。パラメータ保存部３１は、例えば、目的物質の力場パラメータおよび目的物質に類似する物質の力場パラメータなどを保存することができる。ここで、力場パラメータは、分子構造情報に基づいて量子化学計算により算出することができる。また、力場パラメータは、過去に得られた結晶構造情報に基づいて古典力場計算により算出することもできる。

【0069】

これにより、自由エネルギー算出部５は、極小条件判定部４から部分集合Ｔ２に含まれる格子ベクトル変数と原子位置変数の組が入力されると、パラメータ保存部３１に保存された目的物質の力場パラメータまたは目的物質に類似する物質の力場パラメータを探索する。続いて、自由エネルギー算出部５は、パラメータ保存部３１に保存された力場パラメータおよび極小条件判定部４から入力された格子ベクトル変数と原子位置変数の組から、擬調和振動子近似によってフォノンの固有振動数を算出し、そのフォノンの固有振動数に基づいて格子ベクトル変数と原子位置変数の組の自由エネルギーを算出する。
そして、算出された自由エネルギーは自由エネルギー算出部５から極小条件判定部６に出力され、その自由エネルギーが部分集合Ｔ２を定義域とする自由エネルギーの連続関数ｆにおいて極小を示すか否かが判定される。

【0070】

このように、パラメータ保存部３１が量子化学計算または古典力場計算に基づいて算出された力場パラメータを保存するため、自由エネルギー算出部５が、パラメータ保存部３１に保存された力場パラメータに基づいて自由エネルギーを正確に算出することができる。

【0071】

本実施の形態によれば、自由エネルギー算出部５が、量子化学計算または古典力場計算から擬調和振動子近似に基づいてフォノンの固有振動数を算出するため、そのフォノンの固有振動数に基づいて自由エネルギーを正確に算出することができる。

【0072】

実施の形態４
上記の実施の形態３では、自由エネルギー算出部５は、量子化学計算または古典力場計算から算出されるフォノンの固有振動数に基づいて自由エネルギーを算出したが、自由エネルギーを算出することができればよく、これに限られるものではない。
例えば、図１２に示すように、実施の形態３のパラメータ保存部３１に換えて経験式保存部４１を配置することができる。

【0073】

経験式保存部４１は、過去に得られた結晶構造に含まれる格子ベクトル変数および原子位置変数の値に基づいて自由エネルギーを算出する経験式を保存するものである。経験式としては、例えば、目的物質に類似する結晶構造について格子ベクトル変数および原子位置変数の値から自由エネルギーの値を過去に算出したときの式などが挙げられる。

【0074】

これにより、自由エネルギー算出部５は、極小条件判定部４から部分集合Ｔ２に含まれる格子ベクトル変数と原子位置変数の組が入力されると、経験式保存部４１に保存された経験式に基づいて、格子ベクトル変数と原子位置変数の組の自由エネルギーを算出する。
そして、算出された自由エネルギーは自由エネルギー算出部５から極小条件判定部６に出力され、その自由エネルギーが部分集合Ｔ２を定義域とする自由エネルギーの連続関数ｆにおいて極小を示すか否かが判定される。

【0075】

このように、経験式保存部４１が自由エネルギーを算出するための経験式を保存するため、自由エネルギー算出部５が、経験式保存部４１に保存された経験式に基づいて自由エネルギーを速やかに算出することができる。

【0076】

本実施の形態によれば、自由エネルギー算出部５が、過去に得られた結晶構造情報に含まれる格子ベクトル変数および原子位置変数の値に基づいて経験式を算出するため、その経験式に基づいて自由エネルギーを速やかに算出することができる。

【0077】

実施の形態５
上記の実施の形態１～４において、極小条件判定部は、目的物質のＸ線解析データが既に得られている場合には、そのＸ線解析データを参照して目的物質の格子ベクトル変数と原子位置変数の組を算出することが好ましい。
例えば、図１３に示すように、実施の形態１にＸ線データ保存部５１を新たに配置することができる。

【0078】

Ｘ線データ保存部５１は、極小条件判定部６に接続され、実験により得られた目的物質のＸ線解析データを保存する。Ｘ線解析データとしては、例えば、Ｘ線解析で得られた目的物質の格子ベクトル変数および原子位置変数の値などが挙げられる。

【0079】

これにより、極小条件判定部６は、自由エネルギー算出部５から自由エネルギーが入力されると、Ｘ線データ保存部５１から目的物質のＸ線解析データを取得し、Ｘ線解析データを参照して目的物質の格子ベクトル変数と原子位置変数の組を算出する。具体的には、極小条件判定部６は、自由エネルギー算出部５から入力された自由エネルギーの格子ベクトル変数および原子位置変数の値と、Ｘ線データ保存部５１に保存されたＸ線解析データの格子ベクトル変数および原子位置変数の値との一致度を算出する。続いて、極小条件判定部６は、自由エネルギー算出部５から入力された自由エネルギーに一致度を加えた値がその連続関数ｆ、すなわち一致度を加えた自由エネルギーの連続関数ｆにおいて極小を示すか否かを判定する。このとき、極小条件判定部６は、自由エネルギーと一致度に所定の重み付けをして加えることが好ましい。

【0080】

そして、極小条件判定部６は、極小を示すと判定された自由エネルギーの格子ベクトル変数と原子位置変数の組を目的物質の結晶構造として算出する。これにより、目的物質の格子ベクトル変数と原子位置変数の組をより高精度に算出することができる。

【0081】

本実施の形態によれば、極小条件判定部６が、自由エネルギーだけでなく、Ｘ線解析データを参照して格子ベクトル変数と原子位置変数の組を算出するため、目的物質の結晶構造を示す格子ベクトル変数と原子位置変数の組をより高精度に算出することができる。また、極小条件判定部６は、直積集合Ｔ２の範囲で変更された格子ベクトル変数と原子位置変数の組とＸ線解析データとの一致度を算出するため、一致度と自由エネルギーに基づいて目的物質の結晶構造を示す格子ベクトル変数と原子位置変数の組をさらに高精度に算出することができる。

【0082】

なお、上記の実施の形態１～５において、極小条件判定部は、部分集合Ｔ２に含まれる格子ベクトル変数と原子位置変数の組を順番に変える必要はなく、格子ベクトル変数と原子位置変数の組を様々変更して連続関数ｆの極小点を探索することができる。

【実施例】

【0083】

（実施例１）
実施の形態１の結晶構造算出方法を用いてアルゴンの結晶構造を算出した。結晶構造としては、０Ｋから９０Ｋにおける格子定数の変化を算出した。その結果を図１４に示す。
また、０Ｋから９０Ｋのそれぞれの温度について格子定数に対する自由エネルギーの変化を算出した。その結果を図１５に示す。
（比較例１）
ＬＥＭ法を用いてアルゴンの結晶構造を算出した。結晶構造としては格子定数を算出した。その結果を図１４に示す。
また、０Ｋから９０Ｋのそれぞれの温度について格子定数に対する格子エネルギーの変化を算出した。その結果を図１５に示す。

【0084】

図１４に示す結果から、実施例１は、アルゴンの熱膨張の実験値（Ｄｏｂｂｓ＆Ｊｏｎｅｓ，１９５７）に対する格子定数の差が０．７％以下であるのに対して、比較例１は１．６％であり、実験値に近い値を算出できることがわかった。また、比較例１は温度の軸が存在しないのに対して、実施例１は昇温に伴った熱膨張の割合が実験値と良く一致することがわかった。さらに、実施例１は、８７Ｋ以上で自由エネルギーの極小点が消失して無限に膨張している。この温度は、アルゴンの融点（８３．８Ｋ）と沸点（８７．３Ｋ）の中間値であり、沸点以上では結晶状態が安定でないことを再現できることがわかった。このように、ＦＥＭ法を用いることによって、熱的現象である熱膨張を再現できると共に、有限温度で安定でない結晶構造を適切に除外することができた。

【0085】

図１５に示す結果から、実施例１は、温度が上昇するほど自由エネルギーの曲率が小さくなり、原子振動の大きさが大きくなることがわかった。このとき、温度が０Ｋから８０Ｋでは自由エネルギーの極小点を算出できるが、温度が９０Ｋでは極小点を算出できなかった。これにより、有限温度で不安定な構造を除外することができ、過剰数の結晶構造（多形）を抑制することができた。

【符号の説明】

【0086】

１ 結晶構造算出装置、２ 初期構造生成部、３ 格子エネルギー算出部、４，６，２３ 極小条件判定部、７ 出力部、８ 原子位置変数変更部、９ 格子ベクトル変数変更部、１０，２２ 変数変更部、１１ 制御部、１２ 操作部、１３ 格納部、２１ ずれ方向算出部、３１ パラメータ保存部、４１ 経験式保存部、５１ Ｘ線データ保存部、Ｔ１ 直積集合、Ｔ２ 部分集合、Ｖ１～Ｖ６，Ｖ１ａ，Ｖ２ａ，Ｖ１’～Ｖ３’ 格子ベクトル変数と原子位置変数の組、Ｌ 等高線、Ｕ１～Ｕ９ 格子エネルギー、Ｒ１～Ｒ９ 結晶構造候補、ｆ 自由エネルギーの連続関数、Ｎ１～Ｎ３ ずれ方向。

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】

【図14】

【図15】