(57)【要約】【課題】これまでの方法では、与えたモデルのうちどれがより良いかの評価材料は得られるものの、真の解（結晶構造と磁気構造）である保証が出来ないという問題があった。【解決手段】そこで、本発明では、構造解析において、大域的最適化手法である内点法により半正定値計画法を解く方法と、多項式緩和と呼ばれる方法を組み合わせることで、大域的収束が保証されている状況で最適化問題を解き、大域的最適解が実際に得られた計算結果から、実験値とよく合う構造パラメータを計算結果から判断することが可能な構造解析装置を提供する。【選択図】図３
2020071133.tif 000002
【特許請求の範囲】
【請求項１】
物質に放射線を照射して回折・散乱スペクトルを得る構造解析装置であって、
前記回折・散乱スペクトルのピーク位置・ピーク強度を取得する反射強度取得部と、
前記ピーク位置・ピーク強度の解析により得られる２次計画問題の目的関数・制約条件を取得して、緩和により半正定値計画問題を作成する変数制約条件取得・半正定値計画緩和実行部と、
前記半正定値計画問題を大域的最適化手法により解くことで、大域的最適解を得る大域的最適化処理実行部と、
前記大域的最適解から複数の最適な解を算定する半正定値計画法実行結果精査・構造パラメータ取得部と、
計算結果から測定により得られた回折・散乱スペクトルのピーク位置・ピーク強度を精密化する構造パラメータ変数精密化実行部と
を含む構造解析装置。
【請求項２】
前記大域的最適化処理実行部は、前記大域的最適化手法として主双対内点法を用いて、前記２次計画問題から緩和によって作成された前記半正定値計画問題を解くために、前記主双対内点法により大域的最適化を実行する、
請求項１に記載の構造解析装置。
【請求項３】
前記半正定値計画法実行結果精査・構造パラメータ取得部は、半正定値計画問題（ＳＤＰ）ソルバーから出力される双対ギャップの値（Δｆ）が、計算機イプシロンの値から得られる値と比較して小さいことを判断基準とする、
請求項１又は２に記載の構造解析装置。
【請求項４】
前記半正定値計画法実行結果精査・構造パラメータ取得部は、前記半正定値計画問題の実施結果から得られた構造モデルの強度と実測の強度との差から求まる値を目的関数（ｆ(gl)）としたときに、前記目的関数を最小とする複数の構造モデルの内、最小のものを選択することを判断基準とする、
請求項１又は２に記載の構造解析装置。
【請求項５】
前記半正定値計画法実行結果精査・構造パラメータ取得部は、最適解行列のランク（階数）を判断基準とする、請求項１又は２に記載の構造解析装置。
【請求項６】
前記反射強度取得部は、前記物質ｓに放射線を放射して得られる前記回折・散乱スペクトルの前記ピーク位置・ピーク強度のデータを回折装置から取り込む機能を有する、
請求項１から５のいずれか１項に記載の構造解析装置。
【請求項７】
前記変数制約条件取得・半正定値計画緩和実行部は、前記回折・散乱スペクトルの前記ピーク位置・ピーク強度の集合、および、求める物質構造に関わる制約条件より、求めたい構造パラメータを変数とする前記２次計画問題の前記目的関数・制約条件を取得する機能を有する、請求項１から６のいずれか１項に記載の構造解析装置。
【請求項８】
前記構造パラメータ変数精密化実行部は、前記半正定値計画問題の実施結果と実測値との差から定まる複数のＲ因子・カイ二乗値の内、最小のものを選択することで、前記回折・散乱スペクトルのピーク位置・ピーク強度に基づき局所最適化を用いて精密化する、請求項１から７のいずれか１項に記載の構造解析装置。
【請求項９】
コンピュータを、請求項１から８のいずれか１項に記載の構造解析装置として動作させるためのプログラムであって、
前記コンピュータを前記構造解析装置が備える各部として機能させるためのプログラム。
【請求項１０】
情報処理装置によって実行される、物質に放射線を照射して回折・散乱スペクトルを得るための構造解析方法であって、
前記回折・散乱スペクトルのピーク位置・ピーク強度を取得する反射強度取得ステップと
前記ピーク位置・ピーク強度の解析により得られる２次計画問題の目的関数・制約条件を取得して、緩和により半正定値計画問題を作成する変数制約条件取得・半正定値計画緩和実行ステップと、
前記半正定値計画問題を大域的最適化手法により解くことで、大域的最適解を得る大域的最適化処理実行ステップと、
前記大域的最適解から複数の最適な解を算定する半正定値計画法実行結果精査・構造パラメータ取得ステップと、
計算結果から測定により得られた回折・散乱スペクトルのピーク位置・ピーク強度を精密化する構造パラメータ変数精密化実行ステップと
を含む構造解析方法。
【請求項１１】
前記大域的最適化処理実行ステップは、前記大域的最適化手法として主双対内点法を用いて、前記２次計画問題から緩和によって作成された前記半正定値計画問題を解くために、前記主双対内点法により大域的最適化を実行する、
請求項１０に記載の構造解析方法。
【請求項１２】
前記半正定値計画法実行結果精査・構造パラメータ取得ステップは、半正定値計画問題（ＳＤＰ）ソルバーから出力される双対ギャップの値（Δｆ）が、計算機イプシロンの値から得られる値と比較して小さいことを判断基準とする、
請求項１０又は１１に記載の構造解析方法。
【請求項１３】
前記半正定値計画法実行結果精査・構造パラメータ取得ステップは、前記半正定値計画問題の実施結果から得られた構造モデルの強度と実測の強度との差から求まる値を目的関数（ｆ(gl)）としたときに、前記目的関数を最小とする複数の構造モデルの内、最小のものを選択することを判断基準とする、
請求項１０又は１１に記載の構造解析方法。
【請求項１４】
前記半正定値計画法実行結果精査・構造パラメータ取得ステップは、最適解行列のランク（階数）が１であることを判断基準とする、請求項１０又は１１に記載の構造解析方法。
【請求項１５】
前記反射強度取得ステップは、前記物質ｓに放射線を放射して得られる前記回折・散乱スペクトルの前記ピーク位置・ピーク強度のデータを回折装置から取り込む、
請求項１０から１４のいずれか１項に記載の構造解析方法。
【請求項１６】
前記変数制約条件取得・半正定値計画緩和実行ステップは、前記回折・散乱スペクトルの前記ピーク位置・ピーク強度の集合、および、求める物質構造に関わる制約条件より、求めたい構造パラメータを変数とする前記２次計画問題の前記目的関数・制約条件を取得する機能を有する、請求項１０から１５のいずれか１項に記載の構造解析方法。
【請求項１７】
前記構造パラメータ変数精密化実行ステップは、前記半正定値計画問題の実施結果と実測値との差から定まる複数のＲ因子・カイ二乗値の内、最小のものを選択することで、前記回折・散乱スペクトルのピーク位置・ピーク強度に基づき局所最適化を用いて精密化する、請求項１０から１６のいずれか１項に記載の構造解析方法。
【請求項１８】
コンピュータに、請求項１０から１７のいずれか１項に記載の構造解析方法を実行させ、前記コンピュータを前記情報処理装置として動作させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、物質に放射線を照射したときに得られる回折スペクトルから該物質の結晶構造又は磁気構造を解析する技術に関する。具体的には、放射線を照射した物質から得られた値に基づいて２次計画問題の目的関数・制約条件を取得して、当該制約条件を緩和することで得られた半正定値計画問題を解くことで、大域的最適解を得ることができる構造解析装置、方法及びプログラムに関する。
【背景技術】
【０００２】
構造解析装置に関する従来技術として、パソコン（ＰＣ）、ワークステーション（ＷＳ）等を用いた、結晶の構造解析を自動的に行う、放射線（Ｘ線）解析装置がある。この種の解析装置は、結晶の試料をＸ線回折装置にセットするだけで、自動的に格子定数、空間群、回折強度等を測定し、これらの測定データに基づいて構造解析が行なわれ、結合距離、結合角等をグラフィック画面上に出力することができる。
【０００３】
一般に、物質中の原子や電子の配置（結晶構造）はＸ線回折などの実験データを解析することにより決定される。例えば、従来の微構造解析方法（特許文献１）では、磁場発生部に対して、微結晶を懸濁させた試料を回転させることで三次元配向を実現し、試料（微結晶）を擬単結晶化することができ、試料に対して、回転を継続しながらＸ線を照射し、その試料を透過して回折したＸ線をＸ線検出部が検出することで、擬単結晶のＸ線回折像を得るものである。
【０００４】
その際、Ｘ線遮蔽装置により、試料の回転方向の一部である特定部位が所望の方向を向いていないときはＸ線の照射を遮蔽し、特定部位が所望の方向を向いたときにのみＸ線の照射を許容している。これにより、前記特定部位が所望の方向を向いた状態の試料を透過して回折したＸ線を、Ｘ線検出部で断続的に検出することができるため、擬単結晶化した試料を回転させながら良好なＸ線回折像を得ることができる。
【０００５】
また、物質中の原子や電子の配置（結晶構造）や磁気モーメントの配列（磁気構造）は中性子回折などの実験データを解析することにより決定される（特許文献２）。これらの解析は、非線形パラメータによって表現された関数を、実験データにフィッティングすることで行われている。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】国際公開第２０１３／１１８７６１号パンフレット
【特許文献２】特開平１０−２４７５９９号公報
【非特許文献】
【０００７】
【非特許文献１】L. Vandenberghe ら、SIAM Review 38, 49 (1996).
【非特許文献２】M. X. Goemansら, Journal of the ACM, 42 (1995), pp. 1115−1145.
【非特許文献３】C. Bachocら, Journal of the American Mathematical Society 21, (2008), pp. 909−924.
【非特許文献４】A. Esmaeiliら、 J. Appl. Cryst. 50, 651 (2017).
【非特許文献５】K. Tomiyasuら、 J. Phys. Soc. Jpn. 81, 034709 (2012).
【非特許文献６】Makoto Yamashitaら, "A high-performance software package for semidefinite programs: SDPA 7," Research Report B-460 Dept. of Mathematical and Computing Science, Tokyo Institute of Technology, Tokyo, Japan, September, (2010).
【非特許文献７】R.H Tutuncu, K.C. Tohら, "Solving semidefinite-quadratic-linear programs using SDPT3," Mathematical Programming Ser. B, 95 (2003), pp. 189--217.
【非特許文献８】A. S. Willsら、J. Phys. Condens. Matter Lett. 18, L37 (2006).
【非特許文献９】K. Matsuhiraら、J. Phys. Soc. Jpn. 80, 094701 (2011).
【非特許文献１０】J. Rodr´iguez-Carvajal、Phys. B (Amsterdam, Neth.) 192, 55 (1993).
【発明の概要】
【発明が解決しようとする課題】
【０００８】
しかしながら、これまでの方法では、与えたモデルのうちどれがより良いかの評価材料は得られるものの、真の解（結晶構造と磁気構造）である保証が出来ないという問題があった。そのため、暫定的に得た解が間違いであること、もっともらしい解であるという傍証を揃えるために大変な時間と労力と費用コストを要すること、解の候補が多すぎて暫定解すら得られないことなどが生じて、機能性材料や製薬に関する広範な産業や理工学分野の発展の足枷となってきた。
【０００９】
このような解析は、高次元のパラメータ空間で行われる非線形最適化問題の一種とみなすことができ、上述した解析上の問題は、(i) 非線形最適化で大域的な最適解を得ることが困難であること、(ii) ローカルミニマムと呼ばれる局所的な意味での最適解に陥る可能性があること、(iii) 回折データから得られる構造因子強度は、知りたい構造パラメータの２次式となるため、結晶構造が持つ対称性やその他の状況次第で複数の解が存在し得ることなどに起因する。
【００１０】
このため、シミュレーティドアニーリング法、モンテカルロ法、遺伝的アルゴリズムなどの大域的最適化法が使用されることがある。しかし、いずれも大域的最適解に行き着くことが保証されている手法ではなく、より慎重な解析を実施したとしても、解析結果を真の解だと保証するものではなかった。
【００１１】
一般に、非線形計画問題の大域的最適化を短時間で実施するための方法として、半正定値計画緩和法を挙げることができる。半正定値計画緩和法は、双対定理や相補性定理や内点法の方法論によって、解析結果を大域的最適解だと保証することができる（非特許文献１）。ただし、全ての非線形計画問題について大域的最適解が得られたことの保証が実施できるわけではなく、主に２次計画問題と呼ばれるものが対象となる。
【００１２】
また、数学分野において、大域的最適値の下界や上界を求めることや、近似解を短時間で求めるために緩和が使用されている(非特許文献２、非特許文献３)。しかしながら、一般には、２次計画問題に緩和を適用しても、実際には大域的最適解の保証までは行えなくなる状況も発生する。同手法において、大域的最適解の保証が実施可能となる問題の種類を明確化することは容易な問題ではなく、実際に計算してみるまでは分からないことも多いといった課題が生じ得る。
【００１３】
図１は、物質中の原子や電子の配置（結晶構造）や磁気モーメントの配列（磁気構造）は中性子回折などの実験データから、局所的範囲の計算結果により、非線形パラメータによって表現された関数ｆのＳＤＲ（半定値緩和、semidefinite relaxation）を実行する前の状態の様子を示す。図１に示す例では、真の解とは異なる解に収束してしまい実測値と計算結果が異なり、計算結果の価値が全くないものとなる。この様な解とは異なる解に収束してしまうのを防止するために大域的範囲での解を求める必要がある。
【００１４】
また、図２はＳＤＲを実行した結果、大域的収束が保証されている状況で最適化問題を解くことで、関数ｆが大域的に最小化された様子を示す。しかしながら、大域的最適解から対応する構造パラメータを求めることと、真の解である保証が必要となる。
【００１５】
このような課題の解決のために、本発明者が誠意検討した結果、
(因子１)  大域的収束が保証されている状況で最適化問題を解くこと、
(因子２)  大域的最適解であることが計算結果から判断できること、
(因子３)  実験値とよく合う構造が複数ある場合に、そのことを計算結果から判断でき、計算結果から複数の解を取得できること、
の３つの因子を解決できれば構造解析で得られた解の保証を可能にできることを見出し、従前知られていなかった、半正定値計画緩和法が物質構造解析の大域的最適解の保証に有効に働くという事実を見出した。
【課題を解決するための手段】
【００１６】
そこで、本発明では、構造解析において、大域的最適化手法である内点法により半正定値計画法を解く方法と、多項式緩和と呼ばれる方法を組み合わせることで、大域的収束が保証されている状況で最適化問題を解き（因子１の解決）、大域的最適解が実際に得られた計算結果（因子２の解決）から、実験値とよく合う構造パラメータを計算結果から判断（因子３の解決）することが可能な構造解析装置を提供する。
【００１７】
本発明に係る構造解析装置の１つの実施形態として、
物質に放射線を照射して回折・散乱スペクトルを得る構造解析装置は、
前記回折・散乱スペクトルのピーク位置・ピーク強度を取得する反射強度取得部と、
前記ピーク位置・ピーク強度の解析により得られる２次計画問題の目的関数・制約条件を取得して、緩和により半正定値計画問題を作成する変数制約条件取得・半正定値計画緩和実行部と、
前記半正定値計画問題を大域的最適化手法により解くことで、大域的最適解を得る大域的最適化処理実行部と、
前記大域的最適解から複数の最適な解を算定する半正定値計画法実行結果精査・構造パラメータ取得部と、
計算結果から測定により得られた・散乱スペクトルのピーク位置・ピーク強度を精密化する構造パラメータ変数精密化実行部と
を含むことを特徴とする。
【００１８】
本発明に係る構造解析装置の好ましい実施形態として、
前記大域的最適化処理実行部は、前記大域的最適化手法として主双対内点法を用いて、前記２次計画問題から緩和によって作成された前記半正定値計画問題を解くために、前記主双対内点法により大域的最適化を実行することを特徴とする。
【００１９】
本発明に係る構造解析装置の好ましい実施形態として、
前記半正定値計画法実行結果精査・構造パラメータ取得部は、半正定値計画問題（ＳＤＰ）ソルバーから出力される双対ギャップの値（Δｆ）が、計算機イプシロンの値から得られる値と比較して小さいことを判断基準とすることを特徴とする。
【００２０】
本発明に係る構造解析装置の好ましい実施形態として、前記半正定値計画法実行結果精査・構造パラメータ取得部は、前記半正定値計画問題の実施結果から得られた構造モデルの強度と実測の強度との差から求まる値を目的関数（ｆ(gl)）としたときに、前記目的関数を最小とする複数の構造モデルの内、最小のものを選択することを判断基準とすることを特徴とする。
【００２１】
本発明に係る構造解析装置の好ましい実施形態として、前記半正定値計画法実行結果精査・構造パラメータ取得部は、最適解行列のランク（階数）を判断基準とすることを特徴とする。
【００２２】
本発明に係る構造解析装置の好ましい実施形態として、前記反射強度取得部は、前記物質ｓに放射線を放射して得られる前記回折・散乱スペクトルの前記ピーク位置・ピーク強度のデータを回折装置から取り込む機能を有することを特徴とする。
【００２３】
本発明に係る構造解析装置の好ましい実施形態として、前記変数制約条件取得・半正定値計画緩和実行部は、前記回折・散乱スペクトルの前記ピーク位置・ピーク強度の集合、および、求める物質構造に関わる制約条件より、求めたい構造パラメータを変数とする前記２次計画問題の前記目的関数・制約条件を取得する機能を有することを特徴とする。
【００２４】
本発明に係る構造解析装置の好ましい実施形態として、前記構造パラメータ変数精密化実行部は、前記半正定値計画問題の実施結果と実測値との差から定まる複数のＲ因子・カイ二乗値の内、最小のものを選択することで、前記回折・散乱スペクトルのピーク位置・ピーク強度に基づき局所最適化を用いて精密化することを特徴とする。
【００２５】
本発明に係る構造解析のためのプログラムの１つの実施形態として、前記プログラムは、コンピュータを前記構造解析装置として動作させ、前記構造解析装置が備える各部として機能させることを特徴とする。
【００２６】
本発明に係る構造解析方法の１つの実施形態として、前記方法は、前記構造解析装置の各部に相当するステップを含み、情報処理装置によって実行されることを特徴とする。
【００２７】
本発明に係る構造解析のためのプログラムの別の実施形態として、コンピュータに、前記構造解析方法を実行させ、前記コンピュータを前記情報処理装置として動作させることを特徴とする。
【発明の効果】
【００２８】
本発明に係る構造解析装置等は、
(因子１)  大域的収束が保証されている状況で最適化問題を解くこと、
(因子２)  大域的最適解であることが計算結果から判断できること、
(因子３)  実験値とよく合う構造が複数ある場合に、そのことを計算結果から判断でき、計算結果から複数の解を取得できることといった３つの因子を解決することにより、結晶構造又は磁気構造といった構造解析において、保証された真の解を提供することができる。
【図面の簡単な説明】
【００２９】
【図１】局所的範囲における最適値に収束してしまう非線形最適化プロセスの模式図である。
【図２】大域的収束が保証されている状況で最適化問題を解くことで大域的範囲における最適値に収束する非線形最適化プロセスの模式図である。
【図３】本発明の１つの実施形態に係る構造解析装置の構成を示す概念図である。
【図４】半正定値計画緩和（ＳＤＲ）と組み合わせた双対定理に基づいた非線形最適化プロセスの模式図である。
【図５】実測値と計算結果とが合致しているかどうか比較するための相関図である。
【発明を実施するための形態】
【００３０】
以下に図面を参照して、本発明の一実施形態について説明する。なお、実施の形態を説明するための全ての図において、同じものには原則として同一の符号を付し、その繰り返しの説明は省略する。本発明の個々の実施形態は、独立したものではなく、それぞれ組み合わせて適宜実施することができる。
【００３１】
本発明は、各原子座標または各サイトに割り当てられたパラメータを求める構造解析において、大域的最適化手法である内点法により半正定値計画法を解く方法と、多項式緩和と呼ばれる方法を組み合わせることで、大域的収束が保証されている状況で最適化問題を解くため（因子１の解決）の大域的最適化処理実行部（３）を、大域的最適解が実際に得られた計算結果（因子２の解決）から、実験値とよく合う構造パラメータを計算結果から判断し取得（因子３の解決）する、半正定値計画法実行結果精査・構造パラメータ取得部（４）を有する構造解析情報処理装置である。
【００３２】
情報処理装置（例えば、コンピュータ）上で実現した、本発明の１つの実施形態に係る構造解析装置を図３に示す。本発明の１つの実施形態に係る構造解析装置は、反射強度取得部（１）、変数制約条件取得部・半正定値計画緩和処理実行部（２）、大域的最適化処理実行部（３）、半正定値計画法実行結果精査・構造パラメータ取得部（４）、構造パラメータ精密化実行部（５）から構成される。
【００３３】
情報処理装置によって実行される構造解析方法は、まず、反射強度取得部によって、前記回折・散乱スペクトルのピーク位置・ピーク強度を取得し、次に、変数制約条件取得・半正定値計画緩和実行部によって、ピーク位置・ピーク強度の解析により得られる２次計画問題の目的関数・制約条件を取得して、緩和により半正定値計画問題を作成する。構造解析方法は、さらに、大域的最適化処理実行部によって、半正定値計画問題を大域的最適化手法により解くことで、大域的最適解を得て、半正定値計画法実行結果精査・構造パラメータ取得部によって、大域的最適解から（複数解があれば複数の）最適な解を算定することができる。最終的に、構造パラメータ変数精密化実行部によって、計算結果から測定により得られた回折・散乱スペクトルのピーク位置・ピーク強度を精密化することができる。
【００３４】
本発明の１つの実施形態に係る構造解析装置は、コンピュータを構造解析装置として動作させるためのプログラムによって実現することができ、プログラムは、コンピュータを構造解析装置が備える各部として機能させることができる。例えば、構造解析に係る各種プログラムは、記憶装置に記憶され、これらのプログラムは、制御部（ＣＰＵ、メモリ等）により必要に応じて読み出されてＲＡＭに移され、ＣＰＵに読み出されて実行される。このように、コンピュータに構造解析に係る各種プログラムを実行させることで、コンピュータを構造解析装置として機能させることができる。また、構造解析装置の少なくとも一部の機能を、特定用途向け集積回路（ＡＳＩＣ）等でハードウェア的に実装することも可能である。
【００３５】
本発明の１つの実施形態に係る構造解析装置において実現される大域的収束が保証されている状況について、関数ｆを用いた一例を図４に示す。関数ｆが大域的に最小化された ｆ(gl) が、大域的収束が保証されている状況で得られたもので、真の最小値であることを保証するため、半正定値計画問題におけるＡk、Ｃ、ｂkが（式１）と（式２）の間で共通であるときの、（式１）（主問題：primal problem）と（式２）（双対問題：dual problem）のその双対ギャップΔｆが実質的にゼロ（計算機イプシロンの平方根と同じくらい小さい）になり、全体での最適化が達成されることが保証される様子が図４に示される。
【００３６】
【数１】
2020071133.tif 000003
【数２】
2020071133.tif 000004
ただし、（式１）の最適化変数 Ｘ＝（Ｘij）は、ｎ次対称行列で、制約条件の
【数３】
2020071133.tif 000005
は、半正定値行列であることを意味する。また、（式２）の最適化変数 Ｙ＝（Ｙij） は、ｎ次対称行列で、制約条件の
【数４】
2020071133.tif 000006
は、半正定値行列であるという意味である。さらに、実験データとよく合う構造パラメータの一意性は、得られた Ｘ＝（Ｘij）の行列としてのランクが１であることから保証される。
【００３７】
以下に各部の詳細について順に説明する。
【００３８】
＜反射強度取得部＞
物質に放射線を放射して得られる回折スペクトルのピーク位置・ピーク強度の集合を反射強度と呼ぶ。反射強度は、回折装置により測定され、また、多くの回折装置には、実測の反射強度、および、計算から導き出される反射強度の推定誤差を出力するシステムが付属する。更に、粉末結晶の回折パターンにおけるピークサーチなど、装置が出力した実験データを処理して強度抽出を実施するソフトウェアも種々存在する（非特許文献４）。
【００３９】
本発明における反射強度取得部は、回折装置から出力された反射強度を入力として読み込む機能、さらに、実測と計算の反射強度の推定誤差も読み込み、解析に用いる機能、さらにまた、実験データから強度抽出を実施するソフトウェアの全て、または、一部を含む。反射強度取得部で得られた情報は、次の方法により、２次計画問題として取り扱われる。
【００４０】
一般に、Ｎ個の観測強度 I1obs,...,INobs が得られたとき、これらの反射強度と解析で求めたいパラメータ x = (x1, …, xn) の間に以下の関係（式３）が成立している。
【数５】
2020071133.tif 000007
【００４１】
【数６】
2020071133.tif 000008
は、観測強度は誤差を含むものの、計算強度からあまり遠くない値であることを指す。例えば、磁気構造解析の反射強度と磁気モーメントの関係式、構造解析の反射強度と各座標点における散乱長密度の関係式、構造解析の反射強度と構造パラメータである、各原子に固有の数値としての原子散乱長、または各原子座標の占有率、または原子散乱長と各原子座標の占有率との積の関係式を（式３）によって表現することができる。
【００４２】
つまり、強度の測定値 Ikobs は、解析で求めたい変数 x の各成分の２次以下の多項式で近似できる状況になる。
【００４３】
＜変数制約条件取得部＞
本発明における変数制約条件取得部では、大域的最適解の取得を実施する半正定値計画問題において適用する制約条件について、強度情報以外の、ソフトウェアユーザが指定する制約条件（文献値、実験条件、理論限界などから推測可能な占有率、散乱長、磁気モーメントなどを表す変数を x として、その上限などの取り得る値の範囲）を考慮することで、例えば、 x の多項式で書くことが可能な制約条件を、ここで付け加えることができる。
【００４４】
＜半正定値計画緩和処理実行部＞
半正定値計画緩和処理実行部では、大域的最適解の取得を困難にする制約条件を取り除き（最適化分野で「緩和」と呼ばれる）、次の処理で大域的最適化を実施する半正定値計画問題の準備を行う。半正定値計画問題とは、上述した（式１）もしくは（式２）の形に記述することが可能な最適化問題を指す(非特許文献５)。
【００４５】
（式２）では yk も最適化変数である。対称行列の内積は、以下（式４）で定義される:
2020071133.tif 000009
【００４６】
具体的なｎ次対称行列Ａk、ベクトルＣの成分とｂkの値を指定すれば、半正定値計画法ソルバー（以下、ＳＤＰソルバー）により、（式１）、（式２）の形の最適化問題を解くことができる。
【００４７】
＜大域的最適化処理実行部＞
本発明における大域的最適化処理実行部は、２次計画問題の制約条件の緩和によって得られた半正定置計画法の大域的最適解を求める。その目的のために、ＳＤＰソルバーを使用する。ＳＤＰソルバーは、半正定値計画問題（ＳＤＰ：semidefinite programming）を解くために開発されているソフトの総称で、数多くソフトが存在する。内点法により効率的に大域的最適解を見つけることができる。
【００４８】
ＳＤＰソルバーとして、例えば、ＳＤＰＡ７（非特許文献６）、ＭａｔｌａｂのソルバーＳＤＰＴ（非特許文献７）などを挙げることができる。
【００４９】
＜半正定値計画法実行結果精査・構造パラメータ取得部＞
本発明における半正定値計画法実行結果精査・構造パラメータ取得部は、大域的最適化処理実行部の実施結果をソフトウェアユーザが理解可能な形に解釈するため、３つの方法として、双対ギャップ、実測と計算結果の差、最適解行列のランクにより、大域的収束および解の一意性を確認、精査する手段を提供する。
【００５０】
（１）双対ギャップ
半正定値計画法の実施結果として大域的最適解が得られていることを双対定理より、双対ギャップがほぼ０に等しいことでチェックする。大域的最適化処理実行部では、大域的最適化手法により半正定値計画問題を解くことで、大域的最適解を得る。大域的最適解への収束が行われたかどうかは、パラメータＡk、Ｃ成分とｂkを共通とする（式１）、（式２）の最適値が一致することで確認できる。
【００５１】
一般に、（式１）、（式２）の最適値の差は双対ギャップと呼ばれ、大域的収束が得られたかどうかの保証は、この双対ギャップの計算値がほぼ０かどうかを見て行うことができる。ほぼ０かどうかの判断基準として、計算機イプシロンから求まる値が目安となる。計算機イプシロンの値は計算機により異なるが、１０-5〜１０-7程度の範囲に収束した場合には「ゼロ」と判断される。この判断は双対定理と呼ばれる、数学分野で知られた定理に基づく理論的に保証されたものである。
【００５２】
（２）ある構造モデルが実験データから棄却されることの証明
（式３）に基づいた構造モデルが実験データから棄却されることを証明するため、実測値と計算値の差を表現する式として、Ｎ個の強度の差を足し合わせた式（式５）、またはＮ個の強度の差の最大値を示す式（式６）を最小化する目的関数として選択する。目的関数は緩和の前後で共通であるが、緩和により制約条件を満たす変数の範囲は広がるので、緩和後に得られた半正定値計画法の大域的最適値は、緩和前の最適化問題の最適値より真に小さいことは常に成立する。
【００５３】
このことから、ＳＤＰソルバーで最適値があまり小さくならなかった構造モデルは、実験データに合わないモデルであることが計算結果より証明できたことになる。実験データに反しない構造モデルでは、半正定値計画問題の結果が、大域的に最小化された大域的最適値では０（ゼロ）に近い正の値になる。
【数７】
2020071133.tif 000010
【数８】
2020071133.tif 000011
【００５４】
（３）最適解行列のランクによる解の一意性の確認
（式５）、（式６）を目的関数としたとき、実験データに反しない構造モデルでは、半正定値計画問題（式１）の大域的最適値は、０に近い正の値になる。さらに、半正定値計画法の大域的最適解 X(opt) は、（式７）のように様々なランクの大域的最適解の重ね合わせとして得られることが成立している。
【００５５】
すなわち、緩和の状況から、半正定値計画法のランク１の大域的最適解 Xi がそれぞれ構造パラメータの大域的最適解 xi に対応しており、 Xi から xi の算出は容易にできる状況である。
【数９】
2020071133.tif 000012
ただし、
Xi : ランク１の半正定値対称行列により与えられる半正定値計画問題の大域的最適解、
Zi : ランク＞１の半正定値対称行列でさらに小さいランクの大域的最適解である半正定値対称行列の和として書けないもの
【００５６】
さらに、（式７）の半正定値計画法の大域的最適解 X(opt) のランクが１である場合、数値的には、１つの固有値を除いて固有値が非常に０に近い値になる状況を指すので、該当の構造モデルを採用した場合の構造パラメータの大域的最適解 x(opt) は唯一つであることが計算結果から示される。 X(opt) のランクが＞１であれば、それぞれの Xi に対応して複数の大域的最適解 xi が求まることを示す。
【００５７】
なお、 X のランク計算は特異値分解などを用いて行うことができる。 X から各 xi を求める計算は、ニュートン法やグレブナ―基底などを用いた代数的方程式の求解を実施する方法により実施することもできる。
【００５８】
大域的最適解として得られた構造パラメータ（複数候補が挙がっている場合もある）に、マルカート法や共役勾配法など局所最適化の手法を適用して、構造パラメータを精密化する。本発明における構造パラメータ精密化実行部は、前記大域的最適化と代数計算により得られた複数解を、構造解析分野で使用される目的関数を有する最適化問題に局所最適化を用いて、強度情報に合うよう構造パラメータを精密化する。
【００５９】
目的関数が（式５）、（式６）の最小化問題は、計算効率の面で他の目的関数より優る。そのため、半正定値計画問題に帰着させる緩和処理では、（式５）、（式６）の様な目的関数が採用されている。この最適値に基づき構造モデルの棄却を行った結果は、半正定値計画法の大域的収束性により理論的に保証されている。
【００６０】
他方で、例えば以下に定義されるＲ因子（式８）やカイ二乗値（式９）は、構造解析において広く使用されている。
【数１０】
2020071133.tif 000013
【数１１】
2020071133.tif 000014
【００６１】
本発明によって、ソフトウェアユーザに理解容易な、構造モデルの適合度を表す指標数値を提供することができる。
【実施例１】
【００６２】
ここでは、フラストレーションと５ｄ電子の特殊性を持ち、スピントロニクス素子の構成要素となりうる３次元Ｗｅｙｌ半金属ないしはその近似材料と期待されるパイロクロアＮｄ2Ｉｒ2Ｏ7を例にとる。パイロクロアＮｄ2Ｉｒ2Ｏ7は、磁気伝播ベクトル
【数１２】
2020071133.tif 000015
によって記述されるオールインオールアウトタイプの磁気構造（正四面体上のスピンが全て内向きまたは外向きに整列している）を示すことが理論的、実験的に示唆されている（非特許文献５、非特許文献８）。しかし、中性子回折実験において Kmag = (0, 0, 0) の磁気構造の磁気反射は核反射と重なり、非偏極回折ではそれらの分離にあいまいさが避けられない。さらに、一般的に、許容可能な他の磁気構造解が存在しないかどうかの不確実性も避けられない。
【００６３】
本発明では、この材料の偏極中性子回折結果と磁気構造の数学的証明の併用により、本発明の効果を示すが、これに限定されるものではない。以下、各部における実施内容について記述する。
【００６４】
＜反射強度取得部＞
偏極中性子散乱実験は、米国オークリッジ国立研究所（ORNL）のHFIR（High Flux Isotope Reactor）に設置された熱中性子３軸分光器ＨＢ１で実施した。 表１に、測定した積分磁気散乱強度を示した。表中 Imag,obs と Imag,cal の Imag は積分磁気散乱強度を、添え字 obs は測定値であることを、添え字 cal は計算値であることを表す。
表１   測定された積分磁気散乱強度と計算された積分磁気散乱強度 
【表１】
2020071133.tif 000016
【００６５】
＜変数制約条件取得・半正定値計画緩和処理実行部＞
空間群
【数１３】
2020071133.tif 000017
における座標 Kmag = (0, 0, 0) の任意の磁気構造は、ユニットセル内の８個の磁性原子と３個のデカルト方向（Ｘ、Ｙ、又はＺ）によって生成され、
【数１４】
2020071133.tif 000018
による narb = 24 個の変数で記述される。ここでは、磁気対称性を既約表現（ＩＲ）解析で分類する（添え字arb は、arbitrary で任意を意味する；ＮｄとＩｒはＮｄ2Ｉｒ2Ｏ7 の磁性原子を表す）。可能な磁気構造は、 ΓNd＝Γ3＋２Γ5＋３Γ7＋６Γ9とΓIr＝Γ3＋２Γ5＋３Γ7＋６Γ9 に分類される。
【００６６】
Γの係数は基底ベクトルの数を表し、各Γに対する変数の総数は
【数１５】
2020071133.tif 000019
に減少することが判る。さらに、Γ9は、６つの強磁性基底ベクトルを含む。 しかし、Ｎｄ2Ｉｒ2Ｏ7は、中性子回折法の検出限界を大きく下回る自発磁化（〜 １０-4 μＢ /式単位）のみを示す（非特許文献９）。したがって、中性子回折の解析では、強磁性基底ベクトルは除外され、ｎΓ9 が１２からさらに６へ減少する。
【００６７】
したがって、Ｎｄ2Ｉｒ2Ｏ7における制約条件の一部削除により本磁気構造解析のために使用される変数は、基底ベクトルの係数として、
【数１６】
2020071133.tif 000020
と定義される。
【００６８】
表２に、London大学の Dr. A. S. Wills が作成した群論的考察による磁気構造モデル構築ソフトＳＡＲＡｈ（非特許文献１０）により得られる磁気構造を構成する基底ベクトル
【数１７】
2020071133.tif 000021
の既約表現（ＩＲ）を示した。
表２  既約表現（ＩＲ）
【表２】
2020071133.tif 000022
Variablesの欄中の記号「―」は、強磁性基底ベクトルが除外されていることを示す。
【００６９】
＜大域的最適化処理実行部＞
測定データを元にＳＤＲ法を適用し、半正定値計画問題を作成する。作成した半正定値計画問題をＳＤＰソルバーで解いて出力される値を示す。本発明では、ソルバーＳＤＰＡ（非特許文献６）を使用したがこれに限定されるものではない。
【００７０】
ＳＤＰソルバーは、半正定置計画法の大域的最適解 X を以下の様に与えた。
【数１８】
2020071133.tif 000023
に大域的に収束している（optは、最適化されたことを意味する）。ここで、 xΓs （xΓ3, xΓ5, xΓ7, xΓ9）はボーア磁子 μB 単位であり、負号は強磁性の最近接 Nd-Ir 相関を示す。
【００７１】
＜半正定値計画法実行結果精査・構造パラメータ取得部＞
（１）ＳＤＰソルバーから出力される双対ギャップ Δf は、 Γ3, Γ5, Γ7, Γ9 の順に、２．１×１０-7、１．４×１０-6，３．１×１０-6、および８．０×１０-7であった。すべてのギャップ値は、１．×１０-5の値以下であり、実質的にゼロであり、大域最適解への収束を保証した。
【００７２】
（２）目的関数（式５）の最小値
【数１９】
2020071133.tif 000024
は、Γ3, Γ5, Γ7, Γ9 の順に、４．５、１９、 ３２、２７であった。この結果は、Γ3において f(gl) が最小値をとることを示す。より明確にするため、Imag,cal の結果を表１に与える。
【００７３】
また、各Γの Imag,cal と Imag,obs の合致度合いを示す相関図を図５に示す。
【００７４】
対角線は Imag,obs と Imag,cal が等価である位置を示す。Γ5, Γ7, Γ9 における計算値と実測値の相関は見られないが、 Γ3 でのみ相関が見られる。 また、弱い強磁性の存在の報告（非特許文献８）にもかかわらず、 Γ9 が解として否定されるのは、弱い強磁性がバルク磁気構造と対称性的に分離されていることを示す。例えば、表面または界面に保護された特性として分離されている場合である。本発明は、この予想を微視的かつ数学的に実証し、将来的には磁壁スピン／エレクトロニクスの応用可能性を指摘するものである。
【００７５】
（３）上記の大域的最適解
【数２０】
2020071133.tif 000025
において、 Γ3 ではランク１の最適解が得られ、（式３）の大域的最適解が ±xΓ3(opt1) のみであることが直ちに保証される。
【００７６】
他方、 Γ5 ではランク２、 Γ9 ではランク３の解が得られている。このとき双対問題の最適解 Y(opt) ((式２)の半正定値対称行列)を用いると、線形方程式
【数２１】
2020071133.tif 000026
と強度情報以外に構造パラメータ x に与えた制約条件が、大域的最適解全体を定める代数方程式を与える。そこでこの代数方程式を解くと、同じ目的関数値を持つ、
【数２２】
2020071133.tif 000027
と
【数２３】
2020071133.tif 000028
の全てが磁気構造パラメータの大域最適解となり、無数の解の存在を示している。さらに、ＳＤＲ法は、Γ7 でゼロ xΓ7(opt1) の結果を与えている。これは、Γ7 では、xΓ7 の値が非ゼロであるとき、
【数２４】
2020071133.tif 000029
となる強度 I において
【数２５】
2020071133.tif 000030
が大きく増加するためである。
【００７７】
一方、ＳＤＲ法を用いなければ、Γ5、Γ7、Γ9でより良い解が存在しないことの確認に多大な時間を要する。また、ＳＤＲ法を用いたΓ5、Γ7、Γ9の事例のように、代数方程式も得られないので、複数の最適解候補を列挙することは本発明のように短い時間では難しい。このように、ＳＤＲ法により困難を容易に克服することができる。
【００７８】
＜構造パラメータ精密化実行部＞
以上の構造パラメータ取得部において、Γ3の一意性が数学的に証明された上で、比較のため、
【数２６】
2020071133.tif 000031
(XΓ3(opt) の固有ベクトル)として得られた XΓ3(opt) を初期値に用いた従来の最小二乗法を行うと、|mNd| = 1.22(5)μBと|mIr| = 0.14(5)μB を得た。これらの値は、報告されている値、|mNd| = 1.27(1)μB 、 |mIr| = 0.34(1)μB に近い。 
【００７９】
また、構造解析によく使われるＲ因子
【数２７】
2020071133.tif 000032
（式７）は、Γ3, Γ5, Γ7, Γ9 に対し、それぞれ、 ０．１３、０．６３、 １．０と ０．８１であった。
さらに、カイ二乗値
【数２８】
2020071133.tif 000033
（式９）は、Γ3, Γ5, Γ7, Γ9 に対し、それぞれ、２．６、１４、 ２５、２０となった。
【００８０】
本発明において、ミラー指数 (hmag, kmag, lmag) と磁気散乱強度 Imag,obs を測定することにより、ＳＤＲ法は、大域最適解として数学的に証明された磁気構造をミラー指数 (hmag, kmag, lmag) と磁気散乱強度 Imag,obs に最も合う様に決めることができる。
【００８１】
その性質としては、第一に、kmag =（0,0,0）型を扱え、多くの応用磁石材料を含む研究開発において非常に高い需要が期待されることを示している。第二に、この方法は偏極した実験に限定されず、kmag≠（0,0,0）または磁気モーメントが検出可能な kmag =（0,0,0）型の磁気構造に対する非偏極実験データの解析にも適用できる。得られたＳＤＲ解を磁気構造パラメータの初期値として用いる通常の最小二乗法を組み合わせることにより、磁気構造と結晶構造の両方を同時に精密化することもできる。
【００８２】
第三に、ＳＤＲ法は何千もの自由パラメータを扱うことができる。したがって、ターゲット物質が複雑な組成で記述されていても、使用可能な空間群が不正確である場合でも、ＩＲ分析が使用できない場合でも、原理的には、多くの反射点を提供する高度な回折計を組み合わせることによって、任意のタイプの磁気構造を解決することができる。また、原子核・電子・磁気モーメント等の散乱長一般の密度分布を決定することにも適用できる。
【実施例２】
【００８３】
実施例１は磁気構造解析であるが、結晶構造解析（式１０）も数学的に同じ構造を有している。特に、磁気反射強度
【数２９】
2020071133.tif 000034
は磁気モーメント
【数３０】
2020071133.tif 000035
の２次多項式であるのと同様に、結晶構造解析で用いられる
【数３１】
2020071133.tif 000036
は原子種によって決まる散乱長fj, または各サイトごとの占有率oj（０以上１以下の値を取る）の2次多項式になる。これらの結晶構造解析の変数についても大域的最適化の需要は存在し、実施例１と同様に大域的最適化・解の保証を実施することができる。以下に具体例を示す。
【００８４】
（結晶構造解析）
【数３２】
2020071133.tif 000037
s: 反射強度のスケールを表すスケールパラメータ、
L: 強度補正のためのローレンツ因子（Qnucや回折角θに依存）
【数３３】
2020071133.tif 000038
: ミラー指数hnuc knuc lnucに該当する逆格子ベクトル、
rj: 原子の座標、
fj: rjに位置する原子の散乱長、
oj: rjに位置する原子の占有率。
【００８５】
ローレンツ因子 L は既定の式による計算値として算出可能である。そこで、原子座標 rj は既知、もしくは原子座標rj の候補がある有限個の範囲に絞られている状況を考える。ここでは以下のいずれかの方法により、温度因子 Ti に由来する不定性を排除する:
ローレンツ因子などの補正を実施した強度情報をフーリエ変換した結果として得られるパターソン関数では、原子座標の各差ベクトルri−rjにピークが存在する。パターソン関数には異なる差ベクトルri−rj, rk−rlがほぼ一致するという状況によるピークオーバーラップが発生するが、該当する座標
【数３４】
2020071133.tif 000039
におけるピーク強度はオーバーラップピークの積分強度の和なので、積分強度Sxの値は以下のように原子散乱長に由来する定数fjと占有率ojの積fjojの２次多項式になる。
（X線散乱因子）
【数３５】
2020071133.tif 000040
ただしajk,cjは原子種から定まる既知の定数。
（中性子散乱長 fj の場合）
【数３６】
2020071133.tif 000041
特に、観測値の回折像からのSxの値の導出は、理論的に厳密な形で実施可能である。
【００８６】
さらに、化学式から得られる以下の値
【数３７】
2020071133.tif 000042
の値）の情報を用いるとスケール因子sの値も推定できる。
【００８７】
結果として得られている強度情報は変数Xj=fjojの２次多項式である。また化学式は変数Xjの１次式となる。そこで実施例１と全く同様に半正定値計画緩和を行い、大域的最適化と解の保証を行うことができる。
【数３８】
2020071133.tif 000043
の取り得る値は原子種によって既知の値に限定されるので、各サイトの原子種と占有率の可能な組合せを全て導出できる。
【００８８】
上述のSxの値を用いる方法を、無機結晶構造データベースの複数の結晶構造およびそのシミュレーション回折像に適用したところ、3572中3019個の結晶構造においては半正定値計画法の結果得られた行列のランクが1となった。この場合は直ちに、導出された各変数Xjの値が唯一の大域的最適解であることが結論された。残りの810個のうち、ランク２，３，４の結晶構造はそれぞれ217,115,45個で、5個以上は176個であった。
【００８９】
例えばNaClのように単位胞内に原子座標の数が二つしかない場合、原子座標r1,r2に載る原子を交換しても得られる結晶構造は同一となるため、問題の性質上必ず、半正定値計画法の結果得られる行列のランクは２となる。しかし、この状況では、S0,Sr1-r2の値から変数 X1, X2の大域的最適値を導出するには単に２次方程式を解けばよいので、ここまで述べた方法も特に必要ない。以下では同様にランクは２となる (Tl2O3).89の事例を示す。
【００９０】
((Tl2O3).89の場合）この場合、最適化変数の数=原子座標の数は4つである。散乱長・占有率の割り当てを行った原子座標は以下の4つである。
(0     0     0),  (1/2   1/2   1/2),  (3/4  3/4  3/4),   (1/4  1/4  1/4)
【００９１】
以下はX線の散乱因子を用いた場合の数値を示す。簡単のため異常分散の寄与は考慮していない。中性子線の場合も同様になる。強度データおよび化学式から得られる制約条件を元にＳＤＲ法を適用し、半正定値計画問題を作成した。作成した半正定値計画問題をＳＤＰソルバーで解いたところ、以下の固有値を持つ大域的最適解X (opt)が得られた。ほぼランク2の行列とみなせる。
[Singular values    ] 1.53109 0.219057 1.18718e-09 2.70513e-10
【００９２】
このことから、上記の原子座標に該当するfjojの値は、以下の式によって与えられる。
【数３９】
2020071133.tif 000044
ただし、右辺の二つのベクトルは、最初の二つに固有値に対応する固有ベクトルから得られるもので、λはこの時点では未定の実数である。
【数４０】
2020071133.tif 000045
であることから以下の２次方程式が得られる。
【数４１】
2020071133.tif 000046
【００９３】
結果として、λ=1.08,0.904が得られるので、「ランク1」の制約条件を除く前の二次計画問題には、以下の2つの大域的最適解が存在したことが分かる:
【数４２】
2020071133.tif 000047
【産業上の利用可能性】
【００９４】
本発明に係る構造解析装置等は、物質中の原子や電子の配置（結晶構造）をＸ線回折などによって解析することや、磁気モーメントの配列（磁気構造）を中性子回折などによって解析することに利用することができる。
【符号の説明】
【００９５】
１  反射強度取得部
２  変数制約条件取得・半正定値計画緩和実行部
３  大域的最適化処理実行部
４  半正定値計画法実行結果精査・構造パラメータ取得部
５  構造パラメータ変数精密化実行部
【図１】
2020071133.tif 000048
【図２】
2020071133.tif 000049
【図３】
2020071133.tif 000050
【図４】
2020071133.tif 000051
【図５】
2020071133.tif 000052