結晶の次元数

pymatgenのサブモジュールを眺めていたらpymatgen.analysis.dimensionalityという結晶の次元数を調べるモジュールを見つけた。 pymatgen.org

関数の使い方よりもそもそも結晶の次元数とはどのように定義されるか、そしてどのようにそれを計算するかに興味があったので参照されている論文を読んだ。

実装のもとになっている論文

読んだ時のメモを残しておく。

[1] Computational identification of promising thermoelectric materials among known quasi-2D binary compounds https://pubs.rsc.org/en/content/articlelanding/2016/ta/c6ta04121c#!divAbstract

  • quasi-2D material の中からZTが高い物質をhigh-throughput計算で探す
  • quasi-2D materialの見つけるアルゴリズム: 様々な(hkl)面に安定なスラブを作ってみて表面を作ることができなければバルク、表面を1つ作れればquasi-2D、2つ以上作れればquasi-1D
  • 原子間にボンドが存在するかは原子半径から判定する

[2] Data Mining for New Two- and One-Dimensional Weakly Bonded Solids and Lattice-Commensurate Heterostructures. https://www.ncbi.nlm.nih.gov/pubmed/28191965

  • 原子間にボンドが存在するかはJMolが使うパラメターで判定する
  • 結晶中の原子数最大のクラスターの原子数をmax1とする
  • 2x2x2のsupercell中の原子数最大のクラスターの原子数をmax2とする
  • log_2(max2/max1)が結晶の実質の次元

[3] Definition of a scoring parameter to identify low-dimensional materials components. https://arxiv.org/pdf/1808.02114.pdf

  • 結晶の実質次元を求める既存手法ではunit cellのとり方がまずいときや、隣接しないセル間にボンドがあるときにうまく判定ができなかった。
  • crystal netが与えられたとき、あるサイトと同クラスターに属する他のセルの同サイトの集合を考える。ユークリッド空間上でこの集合の張る部分空間のrankをサイトの次元と定義する。
  • 上の意味の次元は幅優先探索で求めることができる。
  • ボンドを結ぶ基準を連続的に変化させて、ロバストに次元を判定する。

結晶の次元数

[3]による結晶の次元数の定義が一番妥当だと思うし、pymatgen.analysis.dimensionalityの実装は主に[3]がもとになっている。

次元数を求めるアルゴリズムは[3]に擬似コードが載っているのでそれを参照。 擬似コードよりもpymatgen.analysis.dimensionality.get_dimensionality_larsenの実装を読むほうが分かりやすいかもしれない。

ところで

dimensionとdimensionalityは意味が違うらしい。

wikidiff.com

周期表を描く(pymatgen+matplotlib)

pymatgen.util.plotting.periodic_table_heatmapを使うと周期表の形で元素ごとにある実数値をヒートマップで可視化することができる。

一方、例えば元素ごとにイオン半径を可視化したいときなどは対応する関数がpymatgenにないので自分でいろいろと書かないといけない。 それでもこの場合もperiodic_table_heatmapの実装を真似れば同じように可視化することができる。

本当はカラーバーをつけたいのだがmatplotlibに力尽きたのでその部分はまたいつか考える。

gist.github.com

f:id:lan496:20190202122108p:plain

参照リンク

qiita.com

pymatgen.org

整数行列のスミス標準形・エルミート標準形を計算する

整数行列のスミス標準形とエルミート標準形とそれぞれの変換行列が欲しい状況に遭遇したのだが、標準的なpythonのライブラリには実装されていなかったので自分で書いた。

github.com

スミス標準形の実装には以下のブログを大いに参考にした(というかほぼそのまま):
www.dlfer.xyz


スミス標準形は大まかには、行列の左上から右下にかけて値を確定させていき、対角線上にまだ値が確定していない範囲の整数の最大公約数が来るように行と列を入れ替えていくことで得られる。

ユニモジュラ行列による基本変形によって最大公約数が求められるのはユークリッドの互除法が動くのと同じ理屈になっている。


また、エルミート標準形を求めるときはスミス標準形を求める手続きのなかで行(列)を動かす操作を取り除いて同様にユークリッドの互除法をすればよい。


ところでwikipedia を読んで初めて知ったのだが、整数行列AをH=UAと分解するのをrow-style Hermite normal form、H=ARと分解するのをcolumn-style Hermite normal form と呼ぶらしい。

en.wikipedia.org

pymatgenで状態図を描く

pymatgenで状態図を描く方法の備忘録。

ついでにMaterials Project Rest を使う方法も調べた。

Materials Project API key の登録

まず、Materials Project Rest を使うにはユーザー登録してAPIキーを作成しなければならない。

登録は以下のリンクから行える: https://www.materialsproject.org/dashboard

コマンドラインで以下を実行すると.pmgrc.yamlに自分のAPIキーを設定できて、後でいちいちAPIキーをベタ書きしなくてよくなる

pmg config -a PMG_MAPI_KEY your_api_key

状態図を描く

pymatgen.analysis.phase_diagramを使うと、基準となる物質(単体でも化合物でもよい)の組成に関してギブス自由エネルギー密度のconvex hullを描くことができる(すごい)。 凸包の計算自体はscipy.spatial.ConvexHullでやっているみたい。

gist.github.com

参照リンク

graphillion.graphset.graphsの引数

graphillionでグラフセットを作成する方法を公式ドキュメントで調べたときのメモ。

graphillionのバージョンはv1.01。

複雑な条件でグラフセットを作成するにはgraphillion.graphset.graphsを使えば良さそう。 graphillion.graphset.graphsで指定できる引数は以下の通り: graphillion/graphset.py at master · takemaru/graphillion · GitHub

graphs(vertex_groups=None, degree_constraints=None, num_edges=None, no_loop=False, graphset=None, linear_constraints=None)
    Returns a GraphSet with graphs under given constraints.
    
    This is the base method for specific graph classes, e.g.,
    paths and trees.
    
    This method can be parallelized with OpenMP by specifying the
    environmental variable `OMP_NUM_THREADS`:
    
      `$ OMP_NUM_THREADS=4 python your_graphillion_script.py`
    
    Args:
      vertex_groups: Optional.  A nested list.  Vertices in an
        inner list are connected while those in different inner
        lists are disconnected.  For `[[1, 5], [3]]`, 1 and 5 are
        connected, while they are not connected with 3.
    
      degree_constraints: Optional.  A dict with a vertex and a
        range or int.  The degree of a vertex is restricted by the
        range.  For `{1: 2, 6: range(2)}`, the degree of vertex 1
        is 2 and that of 6 is less than 2, while others are not
        cared.
    
      num_edges: Optional.  A range or int.  This argument
        specifies the number of edges used in graphs to be stored.
        For `range(5)`, less than 5 edges can be used.
    
      no_loop: Optional.  True or False.  This argument specifies
        if loop is not allowed.
    
      graphset: Optional.  A GraphSet object.  Graphs to be stored
        are selected from this object.
    
      linear_constraints: Optional.  A list of linear constraints.
        A linear constraint consists of weighted edges and
        lower/upper bounds.  An edge weight is a positive or
        negative number, which defaults to 1.  Weights of the edges
        that are not included in the constraint are zeros.  For
        instance, `linear_constraints=[([(1, 2, 0.6), (2, 5),
        (3, 6, 1.2)], (1.5, 2.0))]`, feasible graph weights are
        between 1.5 and 2.0, e.g., `[(1, 2), (2, 3), (3, 6)]` or
        `[(1, 2), (2, 5), (5, 6)]`.
        See graphillion/test/graphset.py in detail.
    
    Returns:
      A new GraphSet object.

vertex_groups, degree_constraints, no_loopはdocstringに書いてある通りの意味。

num_edgesはgraphillion.graphset.cliquesの実装を読めば分かるが、列挙されるグラフの辺の数(または範囲)。

graphsetにGraphSet オブジェクトを指定するとユニバースではなくこのgraphsetからグラフが列挙される(?)

linear_constraintsは $$ l \leq \sum_{e \in E} w_e \leq u $$ の形の線形制約を複数与えることができる。

pymatgenで結晶の体積を予測する

pymatgenにはpymatgen.analysis.structure_prediction.volume_predictor.DLSVolumePredictorというクラスがあって、ボンド長のデータベースに基づいて与えられた結晶の体積を予測することができる(すごい)。

DLSVolumePredictor を使えるようにする

pymatgen.analysis.structure_prediction.volume_predictor.DLSVolumePredictorをimportしようとするとDLS_bond_params.yamlがないと怒られた。 多分私のインストールの仕方が悪いのだと思うが、とりあえず https://github.com/materialsproject/pymatgen/blob/master/pymatgen/analysis/structure_prediction/DLS_bond_params.yaml からyamlファイルをダウンロードしてpymatgen/analysis/structure_prediction以下に配置すると動くようになる。

アルゴリズム

具体的なアルゴリズム実装にしか説明がないみたい。 概要を自分用にまとめておく。

まず、結晶中の原子種Aのボンド長bp(M)を $$ bp(M) = r(M) + k(M) * (standard\,deviation\,of\,electronnegativity) $$ で予測する。 ここで、r(M)とk(M)はフィッティングパラメター。

各サイトの近傍原子までの距離と予測したボンド長bp(M)から期待される距離の比を取り、この比の最小値が1になるように結晶を拡大する。 icsd_vol=Trueにすると、ICSD用に更に体積を1.05倍する。

cubic CaTiO3 でDLSVolumePredictorを試してみた。 material projectのデータによると体積は58.440Å^{3}なのに対して、DLSVolumePredictorは49.506Å^{3}と予測している。 お手軽なわりにかなりいい結果だと思う。

gist6dc418f28adafd0bf7e8deb9c40f4e60