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}と予測している。 お手軽なわりにかなりいい結果だと思う。
pymatgenでVASPの入力ファイルを作る
pymatgenを使うとVASPの入力ファイルを簡単に作ることができるが、公式ドキュメントにはまとまった記述がなかったので自分用に備忘録を残しておく。 正直、pymatgenのAPIドキュメントを読めばこの記事を読む必要はないと思う。
POTCAR
http://pymatgen.org/installation.html#potcar-setup に従ってpymatgenが読むPOTCARの配置場所を設定する。
まず、例えばVASPに付属しているPOTCARのディレクトリ(<EXTRACTED_VASP_POTCAR>)を<MY_PSP>以下にpymatgenが読み込める形式でコピーする:
pmg config -p <EXTRACTED_VASP_POTCAR> <MY_PSP>
次に<MY_PSP>をpymatgenが読めるようにする:
pmg config --add PMG_VASP_PSP_DIR <MY_PSP>
実は上のコマンドを実行することによって、ホームディレクトリ直下に以下の内容の.pmgrc.yaml
を作成されている。
以上の設定でpymatgenからPOTCARを読み込めるようになる。
POTCARはpymatgen.io.vasp.inputs.Potcar
で作ることができる。
POSCAR
pymatgen.core.Structure
からpymatgen.io.vasp.inputs.Poscar
で作ることができる。
既存のPOSCARからpymatgenのPoscarオブジェクトを作りたい場合はpymatgen.io.vasp.inputs.Poscar.from_file
を使えば良い。
一般にPOTCARにおける原子種の並びとPOSCARにおける原子種の並びは同じになっていなければならない(http://cms.mpi.univie.ac.at/wiki/index.php/POTCAR)。
そのため、pymatgen.io.vasp.inputs.Potcar
に渡す原子種(symbols)とpymatgen.core.Structure
の原子種(species)は同じになっていなければならない。
具体的には
potcar = pymatgen.io.vasp.inputs.Potcar(symbols) structure = pymatgen.core.Structure(lattice, species, coords) poscar = pymatgen.io.vasp.inputs.Poscar(structure)
みたいなときに
symbols == [a[0] for a in itertools.groupby(species)]
が成り立っていなければならない。
INCAR
INCARタグをキーとするdictからpymatgen.io.vasp.inputs.Incar
で作ることができる。
既存のINCARからpymatgenのIncarオブジェクトを作りたい場合はpymatgen.io.vasp.inputs.Incar.from_file
を使えば良い。
KPOINTS
APIドキュメントを読めば分かるので省略。
例
VASPの公式チュートリアルのfcc Si を計算する場合の例を示す。 交換相関エネルギー汎関数にはPBEを用いた。
import numpy as np from pymatgen.core import Structure, Lattice from pymatgen.io.vasp.inputs import Incar, Kpoints, Poscar, Potcar, VaspInput def fccSi(a): cell = a * np.array([[0., 0.5, 0.5], [0.5, 0., 0.5], [0.5, 0.5, 0.]]) coords = np.array([[0., 0., 0.]]) struct = Structure(Lattice(cell), ['Si'], coords) return struct def write_input(a, output_dir='.'): # poscar structure = fccSi(a) comment = None poscar = Poscar(structure, comment) # incar params = { 'System': 'fcc Si', # not read WAVECAR, and start from superposition of atomic charge density 'ISTART': 0, 'ICHARG': 2, 'ENCUT': 240, 'ISMEAR': 0, 'SIGMA': 0.1 } incar = Incar(params) # kpoints kpoints = Kpoints(num_kpts=0, style="Monkhorst", kpts=((11, 11, 11), ), kpts_shift=(0, 0, 0)) # potcar potcar = Potcar(['Si'], 'PBE') # write input files vaspinput = VaspInput(incar, kpoints, poscar, potcar) vaspinput.write_input(output_dir) if __name__ == '__main__': for a in np.linspace(3.5, 4.3, num=9): output_dir = "fccSi_{}".format(int(10 * a)) write_input(a, output_dir)
周期カーネルに対するFeature Map
ガウス過程回帰の課題のひとつとしてスケーラビリティがあります。
データ数をNとしてfitting の時間計算量、predict の時間計算量がであることから大量のデータを扱いたいときは何らかの近似が必要です。
そのような近似手法のひとつとして、カーネルを有限次元の特徴ベクトルの内積で近似してベイズ線形回帰に帰着するものがあります。
特徴量ベクトルの次元をDとして、ベイズ線形回帰のfitting の時間計算量は 、predict の時間計算量はなので特徴量ベクトルの次元がデータ数よりも小さければより高速に計算ができることになります。
与えられたカーネルに対応する特徴量ベクトルを作る方法としてRandom Fourier Features
があります。*1
Gaussian カーネルやLaplacian カーネルについてはRandom Feature Map で特徴量ベクトルを得られますが、周期カーネルにこの方法を適用するのは難しそうでした。
そこでいろいろ調べたところ、周期カーネルを近似する論文*2を見つけたので概要と実装を備忘録としてまとめておきます。
(追記: 一部の数式が上手く表示されないようなので画像で置き換えました。何もわからない。)
Random Fourier Features
そもそもRandom Fourier Features がどのようなアルゴリズムなのか簡単に述べておきます。
上のstationaryなカーネルがpositive definite であることと、ある有限非負測度のフーリエ変換で書けることが同値である(Bochner の定理)ことを使います。
Bocher の定理より、適当に定数倍したpositive definite なカーネルはある確率密度関数のフーリエ変換で書けます。
\begin{align}
k(\mathbf{x} - \mathbf{y})
&= \int_{\mathbb{R}^{D}} p(\mathbf{\omega}) e^{i (\mathbf{x} - \mathbf{y})^{T} \mathbf{\omega}} d\mathbf{\omega} \\
&= \mathbb{E}_{\mathbf{\omega}, b} \left[ z_{\mathbf{\omega}}(\mathbf{x}) z_{\mathbf{\omega}}(\mathbf{y})\right]
\quad (\mbox{where} \quad z_{\mathbf{\omega}}(\mathbf{x}) := \sqrt{2} \cos(\mathbf{x}^T \mathbf{\omega} + b))
\end{align}
最後の式において、はに従い、bは]の一様分布に従います。
とbをそれぞれD個ずつサンプリングして
\begin{align}
\mathbf{z}(\mathbf{x}) := \sqrt{\frac{2}{D}} \left( \cos(\mathbf{x}^T \mathbf{\omega}_{1} + b_{1}), \cdots, \cos(\mathbf{x}^T \mathbf{\omega}_{D} + b_{D})\right)^{T}
\end{align}
とすると、カーネルをモンテカルロ積分の要領で近似できます。
\begin{align}
k(\mathbf{x} - \mathbf{y})
&\simeq \frac{1}{D} \sum_{j = 1}^{D} z_{\mathbf{\omega}_{j}}(\mathbf{x}) z_{\mathbf{\omega}_{j}}(\mathbf{y}) \\
&= \mathbf{z}(\mathbf{x})^{T} \mathbf{z}(\mathbf{y})
\end{align}
例えばGaussian カーネルは正規分布をフーリエ変換で書けるのでこの方法で対応する特徴量ベクトルを得ることができます。
これは既にsklearn に実装されています。*3
Fourier Features for Periodic Kernels
次に周期カーネルを考えます。
具体的には次のようなカーネルです。*4
\begin{align}
k_{perSE}(t, t') = k_{perSE}(t, t')
&= \exp \left( -2 \left( \frac{\sin (\omega_{0} (t - t'))}{l} \right)^{2} \right) \\
&= \exp \left( \frac{\cos \omega_{0} (t - t') - 1}{l^{2}} \right)
\end{align}
ここでで、このカーネルは周期Tの関数を説明したい時などに使われます。
Random Fourier Features で特徴量ベクトルを作るなら、周期カーネルを逆フーリエ変換して得られる確率密度関数に従って周波数をサンプリングすることになりますが、どうみても簡単にサンプリングができる確率分布を得られそうにありません。
そこで、周期カーネルをフーリエ級数展開して近似することを考えます。
は偶関数であることと、第一種変形ベッセル関数の積分表示
*5を使うと、は次のように計算できます。
フーリエ級数の第K項以降を切り捨てて近似します。
よって、周期カーネルを2K次元の特徴量ベクトルで近似できることがわかりました。
実装とテスト
論文の公式実装が見当たらなかったので自分で書きました。*6
論文ではランダムに2000個の入力に対して、周期T=2の周期カーネルでグラム行列を計算したものと、第C項までを使って近似で得られる行列のフロベニウスノルムでの誤差を測っていたので、自分の実装で検証してみました。
特徴量ベクトルの次元とハイパーパラメターlength_scale(l)を変えた時の誤差は下図のとおりです。
図の左上部分の誤差が元論文と異なりますが、傾向としては同じ結果が得られました。
下図はlength_scale(l)=0.2の場合にカーネルの値をプロットしたものです。
フーリエ級数展開の16次まで取り込めば、元の周期カーネルをほぼ再現できていることがわかります。
*1:https://people.eecs.berkeley.edu/~brecht/papers/07.rah.rec.nips.pdf
*2:http://web.it.usyd.edu.au/~framos/Publications_files/Tompkins-Ramos.pdf
*3:sklearn.kernel_approximation.RBFSampler — scikit-learn 0.19.1 documentation
*4:sklearn.gaussian_process.kernels.ExpSineSquared — scikit-learn 0.19.1 documentation
*5:Modified Bessel Function of the First Kind -- from Wolfram MathWorld
化合物の構造を表現する記述子
無機化合物の化学式とユニットセルの情報及び原子の座標から物性値(結合エネルギー・バンドギャップ)を予測するKaggleコンペに参加していた。
Nomad2018 Predicting Transparent Conductors | Kaggle
feature engineering の過程で分子および無機結晶の構造から特徴量抽出する手法に関する論文をいろいろ読んだので、以下にまとめておく。 基本的に記述子の作り方の部分しか読んでいないので、詳しくは論文を参照。
分子に対する手法
Coulomb Matrix (CM)
- [1109.2618] Fast and Accurate Modeling of Molecular Atomization Energies with Machine Learning
次のような"Coulomb" Matrix を考える。
ここでは原子Iの形式電荷、は原子Iと原子Jの距離。
分子を構成する原子Iと原子Jの間の反発力を表している。 対角成分については、原子の自己エネルギーを多項式でフィッティングしたもの。
Coulomb Matrix の固有値を絶対値の大きい順で並べ、適当に0でpadding して次元を揃えたものを記述子として使う(分子ごとに原子数が異なるので記述子の次元を揃える必要がある)。
分子の回転・並進移動の元で記述子は不変。
Bag of Bonds (BoB)
- http://pubs.acs.org/doi/abs/10.1021/acs.jpclett.5b00831
- ボンドを単語だと思ってbag of words する。
- 分子を構成するボンドを種類ごとにカウントして、それぞれの個数を成分とするベクトルを記述子とする。
Bonds, angles, machine learning (BAML)
Extended Connectivity Fingerprints (ECFP)
Molecular atomic radial angular distribution (MARAD)
Histogram of distances, angles and dihedral angles (HDAD)
- [1702.05532] Machine learning prediction errors better than DFT accuracy
- 原子間距離のヒストグラムを記述子にする(Histogram of disrances, HD)。
- HDに加えて、結合角のヒストグラムも記述子にする(Histogram of distances and angles, HDA)。
- HDAに加えて、二面角*1のヒストグラムも記述子にする(Histogram of distances and angles, HDAD)。
- ヒストグラムのbinは一つの区間に極値が一つだけ含まれるように選ぶ。
Fourier series of radial distance distributions
- [1307.2918] Fourier series of atomic radial distribution functions: A molecular fingerprint for machine learning models of quantum chemical properties
- 原子の動径分布関数をフーリエ変換する。
- 動径分布関数を考えるので回転のもとで不変。
- フーリエ変換するので並進移動のもとで不変。
Molecular Graph (MG)
- [1704.01212] Neural Message Passing for Quantum Chemistry
- [1603.00856] Molecular Graph Convolutions: Moving Beyond Fingerprints
無機結晶に対する手法
分子に対する手法ほど進展していない印象。 分子の場合は並進移動と回転に対する対称性を考えればよかったが、結晶の場合は更にユニットセルの選び方に関して不変であるべきなのでより難しそう。 crystal net *2を一意に表す方法があればいいのだが、そもそもグラフ同型判定問題はNPに属しているという厳しさがある。
Partial radial distribution function (PRDF)
- [1307.1266] How to represent crystal structures for machine learning: towards fast prediction of electronic properties
- 元素Aを中心とした元素Bの動径分布関数。配位環境が異なるとPRDFも異なるので、等価でないサイトを占める元素A全体で平均を取る。
- (cutoff 半径とビン幅の決め方がはっきりしない)
Atom centered symmetry functions
- https://journals.aps.org/prl/pdf/10.1103/PhysRevLett.98.146401
- http://aip.scitation.org/doi/full/10.1063/1.4966192
- まず次のようなカットオフ関数を定義する。
- 各原子に対して次のような動径分布関数を考える。 をいろいろ変えて記述子とする。
- 各原子に対して次のような角度分布関数を考える。 をいろいろ変えて記述子とする。
Bond orientational order parameter (BOP, BOOP)
- Phys. Rev. B 28, 784 (1983) - Bond-orientational order in liquids and glasses
- 原子の分布関数を動径方向に周辺化して球面調和関数で展開したときのスペクトル。
- pymatgen.analysis.local_env で簡単に計算できる。
Bispectrum
- [0910.1019] Gaussian Approximation Potentials: the accuracy of quantum mechanics, without the electrons
- [1209.3140] On representing chemical environments
- 原子の分布関数を三次元球面に埋め込んで、4次元の球面調和関数で展開する。
F-Fingerprint
- How to quantify energy landscapes of solids. - PubMed - NCBI
- 化合物に含まれる原子A, Bに対して次の関数を考える。
- ここで、は原子と原子の距離、はユニットセル内の原子A,Bの個数、はユニットセルの体積、はこの関数をヒストグラムにするビン幅、はgaussian-smeared delta function
Ewald Sum Matrix
- [1503.07406] Crystal Structure Representations for Machine Learning Models of Formation Energies
- Coulomb Matrix における をサイトIを占めるすべての原子とサイトJを占めるすべての原子との間の静電ポテンシャルの和に修正。
- Ewald法を用いて静電ポテンシャルの総和を計算する。
- pymatgen にEwald法を用いて結晶の静電ポテンシャルの総和を計算する関数があるが、計算する式が少し違う。ただEwald Sum Matrix を実際に実装するときに大いに参考にはなる。
Extended Coulomb-like Matrix
Sine Matrix
Crystal Graph (CG)
- [1710.10324] Crystal Graph Convolutional Neural Networks for Accurate and Interpretable Prediction of Material Properties
- primitive cell 内の原子とボンドを特徴量ベクトルで表現して、Graph Convolution
- 畳み込み層のあとでプーリング層に入れて、primitive cell内の原子数によらない特徴量ベクトルを得る。
マルチカノニカル法で個数推定
はじめに
去年の11,12月にマルコフ連鎖モンテカルロ法(MCMC)の集中講義を受けてきました。 講義の中で紹介されたマルチカノニカル法を応用してレアイベントをサンプリングする手法が特に面白かったので、復習も兼ねて実装してみました。
マルチカノニカル法
まずレアイベントサンプリング用のマルチカノニカル法を説明します。
状態の確率分布がが与えられているときに、統計量の分布が知りたいとします。 単純にをサンプリングしてもいいのですが、統計量の周辺分布に極めて確率の小さい区間があるとその部分をうまくサンプリングすることができません。
この問題を回避するために適当な重みを用意して、の代わりに重み に従ってをサンプリングすることを考えます。 これは普通のメトロポリス・ヘイスティングス法で実現できます。 すなわち、状態から状態 に遷移する確率を
とすればいいわけです。
重みとして、求めたい周辺分布に近い分布の逆数を用いるのがマルチカノニカル法です。
重みに従ってをサンプリングしたのときの統計量 の周辺分布は
となって、この分布はほぼ平坦であることがわかります。 平坦な周辺分布からサンプリングすることで、普通にサンプリングすると確率0と推定されてしまうようなレアイベントにもMCMCで訪れられるのがポイントです。
Wang-Landau 法
残る問題は、求めたい周辺分布に近い分布をどのように得るかです。
周辺分布を大まかに求める方法としてWang-Landau 法というものがあります。 詳細はwikiの擬似コードを読むのが一番分かりやすいと思います(https://en.wikipedia.org/wiki/Wang_and_Landau_algorithm)。
大まかに言うと、のヒストグラムが平坦になるように出現回数の低いところの重みを逐次大きくしていき、ある程度平坦になったらより細かく重みを調整していく感じです。
演習問題
具体例として、次のような数え上げ問題を考えます。
L×Lの正方格子において、角と向かい側の角とをつなぐ自己交差のない経路の個数を求めよ。
文章だと分かりづらいですが、要は↓の動画で扱われている問題です。 www.youtube.com
一辺の長さがのとき、辺は全部で本あるので、各辺について使うか否かを考えると辺の配置は全部で通りです。探索する配置の数が指数関数で増大するので、愚直に全通りを試していてはすぐに計算が終わらなくなってしまいます。
ただ、この問題には効率的なアルゴリズムがあり、L=26までの厳密解が求まっているようです(http://oeis.org/A007764?language=japanese)。 なので個数推定をする意義はあまりないのですが、演習問題を解いているということでそこら辺は気にせずいきます。
この問題をマルチカノニカル法で扱うために、正方格子の各辺について使うか否か決めたもの (以下これを「状態 」と呼ぶ)に対してエネルギーを次式で定めます。
ここで、であることと状態が左上と右下の頂点をつなぐ自己交差のない経路であることが同値になっています。 先のマルチカノニカル法を用いてエネルギーの分布を求めて、に状態の総数をかければ今求めたい経路の数を推定することができます。
具体的な実装はこちら。 github.com
結果
以上の方法で自己交差のない経路の数を100回推定し、その平均と不偏分散の平方根を真値と比較したものが下図です。 有効数字2桁くらいの精度が得ることができました。
ちなみにL=6が手元のMacBook Air で4時間くらいかかっています。 L=6の状態の数がで、その中で自己交差のない経路になっているものはおよそ分の1しかないことを考えると、これほどレアな事象をそれなりの精度でサンプリングできるのはすごいと思えます。
L | 推定値 | 誤差 | 真値 |
---|---|---|---|
3 | 185.072 | 2.6228 | 184 |
4 | 8241.03 | 275.432 | 8512 |
5 | 1.2242e+06 | 75237.9 | 1262816 |
6 | 5.58901e+08 | 6.03206e+07 | 575780564 |
参考文献
http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0125062 https://sites.google.com/site/iwanamidatascience/vol2/number-place http://www.springer.com/us/book/9783642031625
AOJ 1320 City Merger
問題文(http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=1320)
AOJ2022(http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2022)とほとんど同じ問題。他の都市名の部分文字列になっている都市名を予め取り除いてからbitDPをする。
#include <bits/stdc++.h> using namespace std; #define REP(i,n) for(int i=0;i<(int)(n);i++) #define ALL(x) (x).begin(), (x).end() typedef long long ll; typedef long double ld; const int INF = 1e9; const ld EPS = 1e-8; bool is_contained(string s, string t) { for(int i = 0; i + s.length() - 1 < t.length(); ++i) { if(s == t.substr(i, s.length())) return true; } return false; } int d[16][16]; int memo[16][1<<16]; int dp(int n, int k, int S) { //cout << k << " " << S << endl; if(memo[k][S] != INF) return memo[k][S]; if(S == (1<<k)) return d[k][k]; int res = INF; REP(i, n){ if(i == k) continue; if(S & (1<<i)){ res = min(res, d[k][i] + dp(n, i, S - (1<<k))); } } //cout << k << " " << S << ": " << res << endl; return memo[k][S] = res; } int main(){ int n; while(cin >> n && n) { vector<string> cities(n); REP(i,n) cin >> cities[i]; vector<string> tmp; REP(i,n){ bool flag = true; REP(j,n){ if(i == j) continue; if(is_contained(cities[i], cities[j])){ flag = false; break; } } if(flag) tmp.push_back(cities[i]); } cities = tmp; n = cities.size(); REP(i,16)REP(j,16) d[i][j] = INF; REP(i,n)REP(j,n){ string s = cities[i]; string t = cities[j]; bool flag = true; for(int k = max(0, (int)(s.length() - t.length())); k < s.length(); ++k){ if(s.substr(k) == t.substr(0, s.length() - k)){ d[i][j] = k; flag = false; break; } } if(flag) d[i][j] = s.length(); } REP(i,n) d[i][i] = cities[i].length(); int res = INF; REP(i,16)REP(j,1<<16) memo[i][j] = INF; REP(i,n){ int tmp = dp(n, i, (1<<n) - 1); res = min(res, tmp); } cout << res << endl; } return 0; }
AOJ 2534 Dictionary
問題文(http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2534)
まず、与えられた文字列の組(string_i, string_j)(i < j)が辞書式順序を満たすようにアルファベットの間に順序を入れていく。アルファベットx, yに順序が入っていることは有向グラフに(x, y)の辺を張ることで表現する。文字列の長さが違うときは'@'(='a' - 1)で埋めておいた。
このグラフの構築段階で(i,j)に辺があるのに(j,i)にも辺が張られているば矛盾しているのでfalseを返す。また、'@'は最小であるはずなのでそれが満たされていない時もfalseを返す。
'@'を始点として辺の張られているアルファベットに深さ優先探索で移動していくとき、もしも矛盾がなければ高々深さは26にしかならないから、それ以上の深さを探索しようとしていたらfalseを返す。
#include <bits/stdc++.h> using namespace std; #define REP(i,n) for(int i=0;i<(int)(n);i++) bool dfs(int d[][27], int k, int cnt) { if(cnt > 30) return false; bool flag = true; REP(j,27) if(d[k][j] == 1) flag = false; if(flag) return true; bool res = true; REP(j,27){ if(d[k][j] == 1){ res = res && dfs(d, j, cnt + 1); } } return res; } bool solve(vector<string> &s) { int n = s.size(); int d[27][27] = {0}; REP(i,n)for(int j = i + 1; j < n; ++j) { int pos = 0; while(pos < s[i].length() && s[i][pos] == s[j][pos]) pos++; if(pos == s[i].length()) continue; if(d[s[j][pos]-'a'+1][s[i][pos]-'a'+1] == 1) return false; d[s[i][pos]-'a'+1][s[j][pos]-'a'+1] = 1; } REP(i,27) if(d[i][0] == 1) return false; return dfs(d, 0, 0); } int main(){ int n; while(cin >> n && n) { vector<string> s(n); REP(i,n) cin >> s[i]; int lmx = 0; REP(i,n) lmx = max(lmx, (int)s[i].length()); REP(i,n)REP(j,lmx-s[i].length()) s[i] += 'a' - 1; if(solve(s)) cout << "yes" << endl; else cout << "no" << endl; } return 0; }