2020-08-28

EXAFS解析の"Thorough Search"法について。

www.jstage.jst.go.jp

Curve fitting (CF)だと最小のR因子を与えるパラメータ（の組）しか最終結果が得られないが、Thorough search (TS)だと（大雑把に言って）他の極小値もわかる、というもの。
具体的には、R因子に上限（例えばR<0.05等）を設定して、それを満たすパラメータ群の中で、一つのパラメータに注目したヒストグラム（ある幅を持たせた値の中に、条件を満たしたパラメータが何個存在するかをプロットしたもの。R因子の値そのものは関係がない点に注意。）のピークが、CFで言うところの極小値に対応する（つまり、ピークを与えるパラメータが解析結果の候補となる）というもの。
複数の候補が得られるメリットとしては、他の実験事実との整合が付く「より正しい構造パラメータ」を抽出出来るという点である。スペクトルと幾何構造が一対一対応している保証がないので、一体何対応なのかがわかっておくと、そこから条件をかけて絞ることができる。

しかし、「R因子」を見る代わりに「ヒストグラム」に置き換えることは正当化されるのだろうか？少し考察する。

ヒストグラムのピークがR因子の極小値に対応して良いのか？

R因子がパラメータに対して連続関数であれば（そしてこれは通常満たされると期待して良い）、R因子の極小値においてパラメータによる微分はゼロになる。すなわち、極小値近傍でパラメータを微小変化させてもR因子は変化しないから、極小値を与えるパラメータ近傍もR因子に対する条件を満たしており、ヒストグラムを取った時に極小値近傍のパラメータは全部カウントされることになる。なので、一見良さそうに見える。
しかし、重要なのがパラメータの数である。パラメータの数が少ないとこの方法は使えない。例えばパラメータ１個では、条件を満たすか満たさないかだけでカウントするから、ゼロか定数かになってピーク構造が現れない。パラメータの数が多くなって条件を満たすことが難しくなった時に、極小値に対応したピークがヒストグラムに現れるようになる。
実際、論文のFig.2で示されている、パラメータが２個の場合の模式図では、TSは単に中心を与えていて、最小値はおろか極小値からもズレている。一方、Fig.4辺りで示されているPt L3-edge EXAFSの解析では、パラメータは4個であり、結果もCFとよく対応している。

ヒストグラムのピークの幅は誤差として評価して良いのか？

前述したことを踏まえると、「ピークが狭ければ狭いほど、『ヒストグラムのピーク』と『R因子の極小値』の対応が良くなる」ということは言えると思う。論文のFig.2のヒストグラムのピークは明らかにlocal minimumを与えるピークより広く、対応が悪いことがわかる。
ただし、R因子が元々パラメータにどれくらい強く依存しているかを見ていないので、そもそも依存性がそんなに強くないパラメータでは元々幅が広いだろうから、「幅が広いからダメ」と直ぐに決め付けることは出来ない。
また、「スペクトルと構造モデルの誤差」そのものはR因子に含まれているので、ヒストグラムのピーク幅が統計的な意味での誤差には対応しないように思う。

ヒストグラムにすることのメリットについて。

極小値を与えるパラメータ（の組み）の候補が知りたいのであれば、R因子を小さい順に並べて、下からいくつか選べば済むように直観的には感じられる。しかし、前述した通り、極小値近傍のパラメータも小さいR因子を与えるため、並び替えるだけでは極小値が単一なのか複数あるのかがわからないのである。重要なのは、小さいR因子を与えるパラメータ群が近いのか遠いのかを判断することである。その意味において、ヒストグラムは「ピーク構造」という形で各パラメータ群を類別することができる。

以上のように、パラメータの多い複雑な系に対して上手い手法であることがわかる。
一つ注意としては、結局はパラメータフィッティングであるため、TSであれCFと同じようにフィッティングパラメータの数はデータ点の数によって制限を受けているという点である。そのため、複雑な系に対して有効であるが、複雑な系であればパラメータの数は増えるので、やはりパラメータの数をいかに減らすかというところに職人芸が必要であることには変わらない。

2020-06-10

運動量演算子のエルミート性。

量子力学数学

いつも忘れる。

参考サイト。
位置演算子と運動量演算子はエルミート演算子であることの証明 - 三浦と窮理とブログ

結局、運動量演算子が微分演算を含み、期待値（つまり積分）の中で部分積分を使うと符号と作用する関数を反転できる、ということによる。
ただし、無限遠で部分積分がゼロになることが暗に仮定されている。これは、直観的に言えば「平面波に対しては運動量演算子のエルミート性はどうなの？」となりそうだが、平面波はめちゃくちゃ特別でそもそも部分積分が要らないというチートな性質があるため、ここの一般的な導出には乗らないのである。

2020-05-28

単一三角格子の固有値問題。

数学固体物理

単一三角格子において、各サイトへのホッピングを $t$ で表すと、ハミルトニアンは
$\displaystyle H = \begin{pmatrix} 0&t&t\\t&0&t\\t&t&0 \end{pmatrix}$
と書ける。
このハミルトニアンの固有値と固有ベクトルを求める。

$\displaystyle \begin{vmatrix} -\varepsilon&t&t\\t&-\varepsilon&t\\t&t&-\varepsilon \end{vmatrix} = 0 \\ \displaystyle \therefore ( - \varepsilon^3 + 2 t^3 ) - ( - 3 t^2 \varepsilon) = 0 \\ \displaystyle \Rightarrow \varepsilon^3 - 3 t^2 \varepsilon - 2 t^3 = 0$

三次方程式の解き方を忘れたので、以下のサイトを参照。
三次方程式の解き方と例題３問 | 高校数学の美しい物語
 方程式の有理数解 | 高校数学の美しい物語
要するに、 $\varepsilon^3$ の係数は $1$ なので、解が有理数であるならば、定数次項である $-2t^3$ の約数が一つの解になっているはずである。
実際、 $\varepsilon = - t$ で確かに方程式がゼロになる。
よって、 $(\varepsilon + t)$ で約分して行けば、
$\displaystyle (\varepsilon + t )^2 (\varepsilon - 2t ) = 0 \\ \displaystyle \therefore \varepsilon = -t, 2t$

$\varepsilon = 2t$

$\displaystyle \begin{pmatrix} -2t&t&t\\t&-2t&t\\t&t&-2t \end{pmatrix} \begin{pmatrix} c_1 \\ c_2 \\ c_3 \end{pmatrix} = 0 \\ \displaystyle \therefore c_1 = c_2 = c_3 \\ \displaystyle \begin{pmatrix} c_1 \\ c_2 \\ c_3 \end{pmatrix} = c_1 \begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix}$

$\varepsilon = -t$

$\displaystyle \begin{pmatrix} t&t&t\\t&t&t\\t&t&t \end{pmatrix} \begin{pmatrix} c_1 \\ c_2 \\ c_3 \end{pmatrix} = 0 \\ \displaystyle \therefore c_1 = - c_2 - c_3 \\ \displaystyle \begin{pmatrix} c_1 \\ c_2 \\ c_3 \end{pmatrix} = c_2 \begin{pmatrix} -1 \\ 1 \\ 0 \end{pmatrix} + c_3 \begin{pmatrix} -1 \\ 0 \\ 1 \end{pmatrix}$

2020-03-26

等方的な関数のフーリエ変換。

数学

1次元

$r = |x|, q = |k_x|$

$\displaystyle f(k_x) = \int^\infty_{-\infty} dx\, f( r ) e^{ - i k_x x } = \int^\infty_0 dx\, f( r ) e^{ - i k_x x } + \int^0_{-\infty} dx\, f( r ) e^{ - i k_x x } \\ \displaystyle \qquad = \int^\infty_0 dr\, f( r ) e^{ - i k_x r } + \int^\infty_0 dr\, f( r ) e^{ i k_x r } = 2 \int^\infty_0 dr\, f( r ) \cos{( k_x r )} \\ \displaystyle \qquad = 2 \int^\infty_0 dr\, f( r ) \cos{( q r )} = f( q )$

2次元

$r = \sqrt{ x^2 + y^2 }, q = \sqrt{ k_x^2 + k_y^2 }$

$\displaystyle f( \vec k ) = \int^\infty_{-\infty} dx\, \int^\infty_{-\infty} dy\, f( r ) e^{ - i \vec k \cdot \vec r } = \int^\infty_0 r dr\, f( r ) \int^{2\pi}_0 d \theta \, e^{ - i q r \cos( \theta - \theta_k ) } \\ \displaystyle \int^{2\pi}_0 d \theta \, e^{ - i q r \cos( \theta - \theta_k ) } = \int^{2\pi}_0 d \theta \, e^{ - i q r \cos \theta } = 2\pi J_0( q r ) \\ \displaystyle \therefore f( \vec k ) = \int^\infty_{-\infty} dx\, \int^\infty_{-\infty} dy\, f( r ) e^{ - i \vec k \cdot \vec r } = 2 \pi \int^\infty_0 dr\, f( r ) \left( r J_0( q r ) \right) = f( q )$

3次元

$r = \sqrt{ x^2 + y^2 + z^2 }, q = \sqrt{ k_x^2 + k_y^2 + k_z^2 }$

$\displaystyle f( \vec k ) = \int^\infty_{-\infty} dx\, \int^\infty_{-\infty} dy\, \int^\infty_{-\infty} dz\, f( r ) e^{ - i \vec k \cdot \vec r } \\ \displaystyle \qquad = \int^\infty_0 r^2 dr\, f( r ) \int^{\pi}_0 \sin\theta d \theta \, e^{ - i q r \cos\theta } \int^{2\pi}_0 d\varphi \\ \displaystyle \qquad = 2\pi \int^\infty_0 r^2 dr\, f( r ) \int^1_{-1} du \, e^{ - i q r u } \\ \displaystyle \qquad = 4\pi \int^\infty_0 dr\, r^2 f( r ) \frac{ \sin( qr ) }{ qr } \\ \displaystyle \qquad = 4\pi \int^\infty_0 dr\, f( r ) \left( r^2 \frac{ \sin( qr ) }{ qr } \right) \\ \displaystyle \qquad = 4\pi \int^\infty_0 dr\, f( r ) \left( r^2 j_0( qr ) \right) \\ \displaystyle \qquad = f( q )$

2020-03-25

Gaussianの半値幅。

数学

規格化されたGaussianは次の形で与えられる。
$\displaystyle g( x ) = \frac{ 1 }{ \sigma \sqrt{ 2 \pi } } e^{ - \frac{ ( x - \mu )^2 }{ 2 \sigma^2 } }$

規格化されている（積分値が1になっている）ことは、次のようにして確かめられる。
$\displaystyle \int^\infty_{-\infty} dx \, e^{ - a x^2 } = \sqrt{ \frac{ \pi }{ a } } \\ \displaystyle \therefore \int^\infty_{-\infty} dx \, g( x ) = \frac{ 1 }{ \sigma \sqrt{ 2 \pi } } \sqrt{ \frac{ \pi }{ \frac{ 1 }{ 2 \sigma^2 } } } = 1$

$\sigma$ を用いた表記は、統計学的な処理を行う際に便利である（と想像される）が、直感的にどれくらいの幅を持つのかが分かりにくい。
そこで、半値半幅 $\Gamma$ を用いた表記を考えることにする。
$\Gamma$ の定義は、次で与えられる。
$\displaystyle \frac{ g( \Gamma + \mu ) }{ \max{ g } } = \frac{ g( \Gamma + \mu ) }{ g( \mu ) } = e^{ - \frac{ \Gamma^2 }{ 2 \sigma^2 } } = \frac{1}{2}$

したがって、 $\sigma \rightarrow \Gamma$ の変換は、
$\displaystyle - \frac{ \Gamma^2 }{ 2 \sigma^2 } = \ln{ \frac{1}{2} } = - \ln{2} \\ \displaystyle \Gamma^2 = 2 \sigma^2 \ln{2} \\ \displaystyle \therefore \Gamma = \sigma \sqrt{ 2 \ln{2} } \Leftrightarrow \sigma = \frac{ \Gamma }{ \sqrt{ 2 \ln{2} } }$
数値的に見積もると、 $\Gamma \approx 1.2 \sigma$ となる。半値半幅 $\Gamma$ は $\sigma$ よりも気持ちちょっと大きい程度である。

元のGaussianの式を $\Gamma$ で表すと、
$\displaystyle g( x ) = \frac{1}{ \Gamma } \sqrt{ \frac{ \ln{2} }{ \pi } } e^{ - \frac{ \ln{2} }{ \Gamma^2 } ( x - \mu )^2 }$

（おそらくLorentz関数でもそうだと思うが、）Gaussian型のピーク形状の場合、半値全幅 $2\Gamma$ （文献によっては、 $\Gamma$ を半値全幅の表記に用いる場合もあるため、注意が必要）を分解能として扱うため、Gaussian broadening等で理論値に実験誤差を入れ込む場合には、半値幅の用いた表記の方が便利である。

2020-02-28

NiOの反強磁性結晶構造における対称性をspglibを使って求める。

通常の結晶構造であれば、データベースか何かから引っ張って来れるが、それをベースにした超周期構造を作ろうと思うと、対称性（空間群）が変わるため、手でやるには歯が立たない。
そこで、入力した原子座標に対して対称性を見つけてくれる spglib を使って、反強磁性でのNiOの対称性を決定することを試みる。

spglibはC言語、FORTRAN, pythonに対応しているが、いずれにおいても「実行ファイル（今風に言えばアプリケーション）」ではなく、（~libとあるように）ライブラリーとして提供されている。そのためGUIのようなものに結晶構造を入れるようなものではないため、多少プログラミングに慣れている必要がある。
Spglib — spglib 1.12.0 documentation

個人的には、pythonが簡便でオススメである。なので、ここではpythonでの例を紹介する。
Spglib for Python — spglib 1.12.0 documentation
pythonバージョンのインストールは、pipですぐ入る。

pip install spglib

これを使って、まずは通常のNiOの結晶構造で試してみる。

import numpy as np
import spglib

# lattice constant
a = 4.19

# structural information ( [lattice constants], [coordination], [atomic numbers] )
# Caution: it should be a tuple
nio = \
    ( [ [ a, 0, 0 ], [ 0, a, 0 ], [ 0, 0, a ] ],
      [ [ 0,   0,   0   ],  # Ni
        [ 0,   0.5, 0.5 ],  # Ni
        [ 0.5, 0,   0.5 ],  # Ni
        [ 0.5, 0.5, 0   ],  # Ni
        [ 0.5, 0,   0   ],  # O
        [ 0,   0.5, 0   ],  # O
        [ 0,   0,   0.5 ],  # O
        [ 0.5, 0.5, 0.5 ]   # O
      ],
      [28,] * 4 + [8,] * 4
    )

# just space group 
print( spglib.get_spacegroup( nio ) )
#Fm-3m (225)

# primitive cell
print( spglib.find_primitive( nio ) )
#(array([[0.   , 2.095, 2.095],
#       [2.095, 0.   , 2.095],
#       [2.095, 2.095, 0.   ]]), array([[0. , 0. , 0. ],
#       [0.5, 0.5, 0.5]]), array([28,  8], dtype=int32))

$2 \times 2 \times 2$ のsuper cell でもちゃんと元の構造の対称性を返してくれる。

# some coordinates are just added 1 to discribe them in the other unit cell
# therefore, lattice constans and coordinates should be multipled and divided by 2, respectively
nio_sc \
    = [ [ ( a, 0, 0 ), ( 0, a, 0 ), ( 0, 0, a ) ],
        [ [ 0,   0,   0   ],  # Ni 0, 0, 0
          [ 0,   0.5, 0.5 ],
          [ 0.5, 0,   0.5 ],
          [ 0.5, 0.5, 0   ],
          [ 1,   0,   0   ],  # Ni 1, 0, 0
          [ 1,   0.5, 0.5 ],
          [ 1.5, 0,   0.5 ],
          [ 1.5, 0.5, 0   ],
          [ 0,   1,   0   ],  # Ni 0, 1, 0
          [ 0,   1.5, 0.5 ],
          [ 0.5, 1,   0.5 ],
          [ 0.5, 1.5, 0   ],
          [ 0,   0,   1   ],  # Ni 0, 0, 1
          [ 0,   0.5, 1.5 ],
          [ 0.5, 0,   1.5 ],
          [ 0.5, 0.5, 1   ],
          [ 1,   1,   0   ],  # Ni 1, 1, 0
          [ 1,   1.5, 0.5 ],
          [ 1.5, 1,   0.5 ],
          [ 1.5, 1.5, 0   ],
          [ 1,   0,   1   ],  # Ni 1, 0, 1
          [ 1,   0.5, 1.5 ],
          [ 1.5, 0,   1.5 ],
          [ 1.5, 0.5, 1   ],
          [ 0,   1,   1   ],  # Ni 0, 1, 1
          [ 0,   1.5, 1.5 ],
          [ 0.5, 1,   1.5 ],
          [ 0.5, 1.5, 1   ],
          [ 1,   1,   1   ],  # Ni 1, 1, 1
          [ 1,   1.5, 1.5 ],
          [ 1.5, 1,   1.5 ],
          [ 1.5, 1.5, 1   ],
          [ 0.5, 0,   0   ],  # O  0, 0, 0
          [ 0,   0.5, 0   ],
          [ 0,   0,   0.5 ],
          [ 0.5, 0.5, 0.5 ],
          [ 1.5, 0,   0   ],  # O  1, 0, 0
          [ 1,   0.5, 0   ],
          [ 1,   0,   0.5 ],
          [ 1.5, 0.5, 0.5 ],
          [ 0.5, 1,   0   ],  # O  0, 1, 0
          [ 0,   1.5, 0   ],
          [ 0,   1,   0.5 ],
          [ 0.5, 1.5, 0.5 ],
          [ 0.5, 0,   1   ],  # O  0, 0, 1
          [ 0,   0.5, 1   ],
          [ 0,   0,   1.5 ],
          [ 0.5, 0.5, 1.5 ],
          [ 1.5, 1,   0   ],  # O  1, 1, 0
          [ 1,   1.5, 0   ],
          [ 1,   1,   0.5 ],
          [ 1.5, 1.5, 0.5 ],
          [ 1.5, 0,   1   ],  # O  1, 0, 1
          [ 1,   0.5, 1   ],
          [ 1,   0,   1.5 ],
          [ 1.5, 0.5, 1.5 ],
          [ 0.5, 1,   1   ],  # O  0, 1, 1
          [ 0,   1.5, 1   ],
          [ 0,   1,   1.5 ],
          [ 0.5, 1.5, 1.5 ],
          [ 1.5, 1,   1   ],  # O  1, 1, 1
          [ 1,   1.5, 1   ],
          [ 1,   1,   1.5 ],
          [ 1.5, 1.5, 1.5 ] ],
        [28,] * 32 + [8,] * 32
      ]

# check whether there is dupulication
for i, p1 in enumerate( nio_sc[1] ):
    for p2 in nio_sc[1][ i+1 : ]:
        if p1 == p2:
            print( p1 )

# lattice constants should be twice in 2*2*2 super cell
nio_sc[0] = [ list( a_ ) for a_ in np.array( nio_sc[0] ) * 2 ]

# coordinate should be divided by 2 
nio_sc[1] = [ list( p_ ) for p_ in np.array( nio_sc[1] ) / 2 ]

nio_sc = tuple( nio_sc )

print( spglib.get_spacegroup( nio_sc ) )
print( spglib.find_primitive( nio_sc ) )

座標の重複があると、'None'が返って来るため、チェックする必要がある。目で見て一生懸命探すのには限界がある。

最後に反強磁性だが、格子定数と座標はsuper cellのままでよく、原子番号で反強磁性秩序の情報を入れる。仮に磁化軸方向と平行な原子をNi(28)として、反平行なものを仮にPd(46)とする。酸素はそのままで良い。
NiOの反強磁性磁化軸は[111]方向なので、それを考慮して原子番号を対応させていく。

nio_af \
    = [ [ ( a, 0, 0 ), ( 0, a, 0 ), ( 0, 0, a ) ],
        [ [ 0,   0,   0   ],  # Ni 0, 0, 0
          [ 0,   0.5, 0.5 ],
          [ 0.5, 0,   0.5 ],
          [ 0.5, 0.5, 0   ],
          [ 1,   0,   0   ],  # Ni 1, 0, 0
          [ 1,   0.5, 0.5 ],
          [ 1.5, 0,   0.5 ],
          [ 1.5, 0.5, 0   ],
          [ 0,   1,   0   ],  # Ni 0, 1, 0
          [ 0,   1.5, 0.5 ],
          [ 0.5, 1,   0.5 ],
          [ 0.5, 1.5, 0   ],
          [ 0,   0,   1   ],  # Ni 0, 0, 1
          [ 0,   0.5, 1.5 ],
          [ 0.5, 0,   1.5 ],
          [ 0.5, 0.5, 1   ],
          [ 1,   1,   0   ],  # Ni 1, 1, 0
          [ 1,   1.5, 0.5 ],
          [ 1.5, 1,   0.5 ],
          [ 1.5, 1.5, 0   ],
          [ 1,   0,   1   ],  # Ni 1, 0, 1
          [ 1,   0.5, 1.5 ],
          [ 1.5, 0,   1.5 ],
          [ 1.5, 0.5, 1   ],
          [ 0,   1,   1   ],  # Ni 0, 1, 1
          [ 0,   1.5, 1.5 ],
          [ 0.5, 1,   1.5 ],
          [ 0.5, 1.5, 1   ],
          [ 1,   1,   1   ],  # Ni 1, 1, 1
          [ 1,   1.5, 1.5 ],
          [ 1.5, 1,   1.5 ],
          [ 1.5, 1.5, 1   ],
          [ 0.5, 0,   0   ],  # O  0, 0, 0
          [ 0,   0.5, 0   ],
          [ 0,   0,   0.5 ],
          [ 0.5, 0.5, 0.5 ],
          [ 1.5, 0,   0   ],  # O  1, 0, 0
          [ 1,   0.5, 0   ],
          [ 1,   0,   0.5 ],
          [ 1.5, 0.5, 0.5 ],
          [ 0.5, 1,   0   ],  # O  0, 1, 0
          [ 0,   1.5, 0   ],
          [ 0,   1,   0.5 ],
          [ 0.5, 1.5, 0.5 ],
          [ 0.5, 0,   1   ],  # O  0, 0, 1
          [ 0,   0.5, 1   ],
          [ 0,   0,   1.5 ],
          [ 0.5, 0.5, 1.5 ],
          [ 1.5, 1,   0   ],  # O  1, 1, 0
          [ 1,   1.5, 0   ],
          [ 1,   1,   0.5 ],
          [ 1.5, 1.5, 0.5 ],
          [ 1.5, 0,   1   ],  # O  1, 0, 1
          [ 1,   0.5, 1   ],
          [ 1,   0,   1.5 ],
          [ 1.5, 0.5, 1.5 ],
          [ 0.5, 1,   1   ],  # O  0, 1, 1
          [ 0,   1.5, 1   ],
          [ 0,   1,   1.5 ],
          [ 0.5, 1.5, 1.5 ],
          [ 1.5, 1,   1   ],  # O  1, 1, 1
          [ 1,   1.5, 1   ],
          [ 1,   1,   1.5 ],
          [ 1.5, 1.5, 1.5 ] ],
        [ 28, 46, 46, 46 ]
      + [ 46, 28, 28, 28 ] * 3
      + [ 28, 46, 46, 46 ] * 3
      + [ 46, 28, 28, 28 ]
      + [ 8, ] * 32
      ]

# check whether there is dupulication
for i, p1 in enumerate( nio_af[1] ):
    for p2 in nio_af[1][ i+1 : ]:
        if p1 == p2:
            print( p1 )

# lattice constants should be twice in 2*2*2 super cell
nio_af[0] = [ list( a_ ) for a_ in np.array( nio_af[0] ) * 2 ]

# coordinate should be divided by 2 
nio_af[1] = [ list( p_ ) for p_ in np.array( nio_af[1] ) / 2 ]

nio_af = tuple( nio_af )

print( spglib.get_spacegroup( nio_af ) )
#R-3m (166)

print( spglib.find_primitive( nio_af ) )
#(array([[ 1.48138871,  0.85528017,  4.83819526],
#       [-1.48138871,  0.85528017,  4.83819526],
#       [ 0.        , -1.71056034,  4.83819526]]), array([[-3.70074342e-17,  0.00000000e+00,  0.00000000e+00],
#       [ 5.00000000e-01,  5.00000000e-01,  5.00000000e-01],
#       [ 7.50000000e-01,  7.50000000e-01,  7.50000000e-01],
#       [ 2.50000000e-01,  2.50000000e-01,  2.50000000e-01]]), array([28, 46,  8,  8], dtype=int32))

対称性がFm-3mからR-3mに落ちているのが分かる。

ここで、結晶軸がベクトル表記のままだと結晶構造作成プログラムでは不便な場合があるので、格子定数の形式に変換すると、

prim_nio_af = spglib.find_primitive( nio_af )
axes = prim_nio_af[0]

abc = [ np.linalg.norm( axis ) for axis in axes ]
print( abc )
#[5.131681011130759, 5.131681011130759, 5.131681011130759]

# alpha beta gamma
abg =[]
abg.append( np.degrees( np.arccos( np.dot( axes[1], axes[2] ) / np.linalg.norm(axes[1]) / np.linalg.norm(axes[2]) ) ) )
abg.append( np.degrees( np.arccos( np.dot( axes[2], axes[0] ) / np.linalg.norm(axes[2]) / np.linalg.norm(axes[0]) ) ) )
abg.append( np.degrees( np.arccos( np.dot( axes[0], axes[1] ) / np.linalg.norm(axes[0]) / np.linalg.norm(axes[1]) ) ) )
print( abg )
#[33.55730976192072, 33.55730976192072, 33.55730976192072]

これで、第一原理計算プログラムで反強磁性が計算できるようになった。

2020-02-09

gfortran での変数の次元の制約。

プログラミング

何も考えずにプログラムを作っていたら、配列の次元が7個を超えるとコンパイルできないらしい。

program test_array7

real( kind( 0d0 ) ) :: r( 1, 2, 3, 4, 5, 6, 7, 8 )

end program test_array7

これを gfortran でコンパイルすると

test_array7.f90:3:46:

 real( kind( 0d0 ) ) :: r( 1, 2, 3, 4, 5, 6, 7, 8 )
                                              1
Error: Array specification at (1) has more than 7 dimensions

添字が11個ある変数の計算しようとしていたので、これはちょっと困った。。。
構造体とかを使って、何とか逃げるしかないか。。。
ちなみに構造体なら7個以上の要素があっても大丈夫であった。

program test_type7                                                              
 
implicit none
 
type :: test
    real( kind( 0d0 ) ) :: r1
    real( kind( 0d0 ) ) :: r2
    real( kind( 0d0 ) ) :: r3
    real( kind( 0d0 ) ) :: r4
    real( kind( 0d0 ) ) :: r5
    real( kind( 0d0 ) ) :: r6
    real( kind( 0d0 ) ) :: r7
    real( kind( 0d0 ) ) :: r8
end type
 
type( test ) :: t
 
end program test_type7

nano_exit

基礎的なことこそ、簡単な例が必要だと思うのです。