VASPなど第一原理計算プログラムで計算した全エネルギー及び原子間力を元に、それを再現するような原子間ポテンシャルを作成します
作成したポテンシャル(正確に言えば最適化したパラメータ)を用いて、フォノン分散を計算したり、後述のMD計算プログラムLAMMPSで使用してMD計算することも可能です
また、このプログラムを転用してボルン有効電荷(と分極)の予測を行うことも可能です

学習データ

単一の.xyzファイルかつ以下の条件を満たすように用意してください

forceを必ず含めること
全エネルギーはcohesive_energy属性として追加すること
同一の構造としてトレーニングデータをまとめるためにconfig_type属性を追加すること

使い方

主に使用するのは実行ファイルhdnnpyで、argparserを使用しているので./hdnnpy -hでオプションなどを確認できます

usage: hdnnpy [-h] {training,param_search,sym_func,test,phonon,optimize} ...

High Dimensional Neural Network Potential

positional arguments:
  {training,param_search,sym_func,test,phonon,optimize}
    training            see `training -h`
    param_search        see `training -h`
    sym_func            see `sf -h`
    test                see `test -h`
    phonon              see `phonon -h`
    optimize            see `optimize -h`

optional arguments:
  -h, --help            show this help message and exit

settings.py.sample、phonopy_settings.py.sampleはそれぞれリネームしてpython スクリプトとして使用してください
それぞれ全体的な設定、フォノン分散計算時の設定を記述します。

$ mv settings.py{.sample,}
$ mv phonopy_settings.py{.sample,}

settings.py

## import 略
# matplotlib全体の設定
mpl.use('Agg')
mpl.rc('font', size=20)

# コマンドライン引数を取得(変更不要)
args = get_parser()

# ファイル・ディレクトリ設定
file = DictAsAttributes(
    xyz_file='GaN/debug_test/GaN.xyz',
    config=['CrystalGa2N2', 'CrystalGa16N16'],
    out_dir='output',
    test_dir='test',
)

# MPI設定(今の所GPUは対応してないので変更不要)
mpi = DictAsAttributes(
    comm=MPI.COMM_WORLD,
    rank=MPI.COMM_WORLD.Get_rank(),
    size=MPI.COMM_WORLD.Get_size(),
    gpu=-1,
    chainer_comm=chainermn.create_communicator('naive', MPI.COMM_WORLD),
)

# 対称関数(Symmetry Functions)のハイパーパラメータ設定
sym_func = DictAsAttributes(
    Rc=[5.0],
    eta=[0.01, 0.1, 1.0],
    Rs=[2.0, 3.2, 3.8],
    lambda_=[-1, 1],
    zeta=[1, 2, 4],
)

# 機械学習のハイパーパラメータ設定
model = DictAsAttributes(
    epoch=1000,
    batch_size=100,
    preproc=None,
    input_size=20,
    init_lr=5.0e-3,
    final_lr=1.0e-6,
    lr_decay=1.0e-6,
    mixing_beta=1.0,
    l1_norm=0.0e-4,
    l2_norm=0.0e-4,
    layer=[
        {'node': 30, 'activation': 'tanh'},
        {'node': 30, 'activation': 'tanh'},
        {'node': 1, 'activation': 'identity'},
    ],
    metrics='validation/main/tot_RMSE',
)

# ハイパーパラメータ探索(ベイズ最適化)に使用するパラメータ設定
skopt = DictAsAttributes(
    kfold=2,
    init_num=5,
    max_num=10,
    space=[
        Real(1.0e-4, 1.0e-2, prior='log-uniform', name='init_lr'),
        Real(1.0e-6, 1.0e-4, prior='log-uniform', name='final_lr'),
        Real(1.0e-6, 1.0e-3, prior='log-uniform', name='lr_decay'),
    ],
    acq_func='LCB',
    callback=[
        SamePointStopper(),
        DeltaYStopper(1.0e-3, n_best=3),
    ],
)

training mode

$ mpirun -np 8 ./hdnnpy training

chainerを使用してトレーニングを実行終了すると、output/<datetime>/ディレクトリに計算結果result.yamlなどが出力されます
lammps.nnpは後述のLAMMPSで使用するパラメータセットです(.txtフォーマット)

param_search mode

$ mpirun -np 8 ./hdnnpy param_search

scikit-optimizeを使用してハイパーパラメータ探索を実行
さらに、実行後最良のパラメータを使用してトレーニングを実行
パラメータ探索の過程は出力されず、skopt_result.yamlのみ出力される

sym_func mode

$ mpirun -np 8 ./hdnnpy sym_func

対称関数計算のみを実行

phonon mode

$ ./hdnnpy phonon

作成したHDNNPを用いて、phonopy_settings.pyの設定のもとでtest_dir(in settings.py)にあるPOSCAR構造のフォノン分散図を計算する

f:id:justice_vsbr:20180802005440p:plain

optimize, test mode

おそらく不要

HDNNP-LAMMPS

有名なMD計算プログラムLAMMPSに、HDNNPを組み込んでみた
ソースコードは同じようにGithubで公開してます

github.com

こちらはC++で書かれたLAMMPSの拡張コードとなっており、LAMMPSのsrc/ディレクトリにコピーするかリンクを作成し、コンパイルすることで使えるようになります

仕様

MKL及びC++ヘッダライブラリEigenを使用するので、よしなに設定しておいてください(投げやり

使い方

$ cp *.h path/to/lammps/src/
$ cp *.cpp path/to/lammps/src/
$ cd path/to/lammps/src/
$ make mpi  # どのconfigでコンパイルするかは各自環境に合わせて
$ lmp_mpi  # 生成される実行ファイル名もconfigにより異なる

LAMMPSスクリプト中では

pair_style nnp
pair_coeff * * lammps.nnp A B

pair_coeffコマンドの第3引数にはHDNNPプログラムで出力されるような形式に従って記述したパラメータファイルを、第4引数以降には各atom typeの原子記号を指定してください

~~………これをREADMEにすればよかったのでは？~~

2018-08-01

LAMMPSで熱伝導率計算

計算物理 LAMMPS MD

LAMMPS

LAMMPS Molecular Dynamics Simulator

C++で書かれた高速で動作する分子動力学法(molecular dynamics: MD)計算パッケージ
詳しい説明は省きます(自分も理解してない部分が多い)

熱伝導率解析

MDで熱伝導率を解析する手法は主に２つに分類されます

非平衡MD(non-equilibrium MD: NEMD):
系に高温領域と低温領域を設け、その間に発生する熱流束と温度勾配から、フーリエ則を用いて熱伝導率を解析
平衡MD(equilibrium MD: EMD):
グリーン–久保公式 - Wikipedia
Green-Kubo公式は平衡状態における物理量の"揺らぎ"から輸送係数(今回は熱伝導率)を求めるもの

正直に言うとこの"平衡"と"非平衡"の違い、完全には理解してません。。
ただ、熱力学的平衡 - Wikipediaによると

熱力学的に非平衡 (non-equilibrium) であるとは、上記の熱的、力学的、化学的平衡のいずれかが満たされていない状態であり、系に物質またはエネルギーの正味の流れ、あるいは相転移などが生じる。またこのような非平衡状態は不安定であるため別の状態へ転移するが、転移速度が極めて遅いために不安定な状態が維持される場合、この状態を準安定状態という。

とある。
NEMDはシミュレーション中高温領域で加熱&低温領域で放熱してエネルギーの流れがあるという点で非平衡、EMDは断熱状態であるという点で平衡、という理解でいいんでしょうかね？

サンプルコードとその解説

LAMMPSのexample/KAPPA/ディレクトリには熱伝導率解析を行うためのサンプルコードとその実行ログが入ってます

6. How-to discussions — LAMMPS documentation

READMEをざっと訳すと、

このディレクトリにはレナードジョーンズ流体に対して5つの方法(NEMD4つ、EMD1つ)で熱伝導率を解析するためのスクリプトが入っている
ただし、あくまでサンプルであって、十分な平衡状態・精度が保証されているものではない
固体に対しても同様に適用可能
in.langevin: ランジュバン動力学を用いて高温/低温領域の温度を制御(NEMD)
in.heat: HEXアルゴリズムを用いて高温/低温領域の熱を制御(NEMD)
in.ehex: eHEXアルゴリズムを用いて(ry
in.mp: Muller-Plathe法を用いて運動エネルギーが最大/最小の粒子同士のエネルギーを交換(NEMD)
in.heatflux: Green-Kubo公式を用いる(EMD)
熱伝導率を計算するためには、NEMDの場合計算後にログファイルから熱流束と温度勾配を計算する必要があり、EMDの場合はスクリプトの中でそのまま計算、表示することができる

LAMMPSはコマンドが大量にあり、スクリプトをパッと見てもさっぱりわからないので、上にある5つの手法それぞれのサンプルスクリプトを可能な範囲で解説していきます。間違っていたら指摘してください。。

in.langevin

## 変数設定
# 単位格子を各軸方向に何倍するか
variable        x equal 10
variable        y equal 10
variable        z equal 20

# レナードジョーンズ流体のパラメータ(単位系はlj)
variable        rho equal 0.6      # 密度
variable        t equal 1.35       # 初期温度
variable        rc equal 2.5       # カットオフ半径
variable        tlo equal 1.0      # 低温領域温度
variable        thi equal 1.70     # 高温領域温度


## 初期設定
units           lj                 # 単位系
atom_style      atomic             # 原子が持つ属性情報を指定(atomicはデフォルト)

lattice         fcc ${rho}         # 単位格子を指定
region          box block 0 $x 0 $y 0 $z   # 単位格子をx,y,z倍した領域を作成 & "box"とタグ付け
create_box      1 box              # "box"を計算モデルに設定
create_atoms    1 box              # 計算モデルにlatticeで指定した構造に従い原子を配置
mass            1 1.0              # type1の原子の質量を指定

velocity        all create $t 87287  # 温度t&シード値87287の下でランダムに速度を設定

pair_style      lj/cut ${rc}       # 原子間ポテンシャルとしてレナードジョーンズポテンシャルを指定
pair_coeff      1 1 1.0 1.0        # ポテンシャルに用いるパラメータ

neighbor        0.3 bin            # 近接原子リストを作成するアルゴリズムを指定
neigh_modify    delay 0 every 1    # 近接原子リストをいつ更新するかを指定


## 高温/低温領域作成(この段階ではまだ温度設定はしない)
# z方向に熱流束を発生させるため、xy平面方向に広がるような領域"hot", "cold"を作成
# それぞれの領域の温度を計算して"Thot", "Tcold"とタグ付け
region          hot block INF INF INF INF 0 1
region          cold block  INF INF INF INF 10 11
compute         Thot all temp/region hot
compute         Tcold all temp/region cold


## 第一段階：全体の温度を$tで均一にして平衡状態になるまでMD
# NVTアンサンブル, 能勢-フーバーサーモスタット, 温度はtで一定
fix             1 all nvt temp $t $t 0.5
thermo          100                # 熱動力データ(温度や圧力など)を100ステップごとにロギング
run             1000               # 1000ステップMD計算を実施

velocity        all scale $t       # 現在の温度を計算し、温度tになるように速度をスケーリング

unfix           1                  # 上のfixコマンドの設定を破棄


## 第二段階：用意した高温/低温領域に温度を設定して再度平衡状態になるまでMD
fix             1 all nve          # NVEアンサンブルはミクロカノニカルアンサンブルであり、全体が断熱状態となる
# 高温/低温領域にランジュバンサーモスタットを設定, tally yesで熱浴と授受したエネルギーを監視
fix             hot all langevin ${thi} ${thi} 1.0 59804 tally yes
fix             cold all langevin ${tlo} ${tlo} 1.0 287859 tally yes
# fixとcomputeをタグを使って紐付け(hot-Thot, cold-Tcold)
fix_modify      hot temp Thot
fix_modify      cold temp Tcold

variable        tdiff equal c_Thot-c_Tcold  # c_XXは、"compute XX"で計算した値に置換される
# ロギングする数値を設定 f_XX, v_XXは同様にfixコマンド、variableコマンドで計算した値に置換される
thermo_style    custom step temp c_Thot c_Tcold f_hot f_cold v_tdiff
thermo          1000
run             10000


## 第三段階：熱伝導率計算
compute         ke all ke/atom     # 運動エネルギー(kinetic energy)を計算
variable        temp atom c_ke/1.5 # atom styleを指定することで、ログファイルに出力する度に原子ごとの平均を計算してくれる

# 高温/低温領域をリセット
fix             hot all langevin ${thi} ${thi} 1.0 59804 tally yes
fix             cold all langevin ${tlo} ${tlo} 1.0 287859 tally yes
fix_modify      hot temp Thot
fix_modify      cold temp Tcold

# 高温/低温領域の実際の温度差(v_tdiff)の累積時間平均を計算
# この設定では、10, 20, ..., 1000ステップの平均値を1000ステップ目における値とする
fix             ave all ave/time 10 100 1000 v_tdiff ave running
thermo_style    custom step temp c_Thot c_Tcold f_hot f_cold v_tdiff f_ave

# z軸に沿って原子を20個(=1/0.05)のグループに分割
compute         layers all chunk/atom bin/1d z lower 0.05 units reduced
# グループごとにv_tempを計算、平均をとってprofile.langevinに出力
fix             2 all ave/chunk 10 100 1000 layers v_temp file profile.langevin

run             20000

in.heat

NEMDなので流れはlangevinと同じ
主な違いはfixコマンドにある
tdfifの時間平均の仕方がlangevinと違う理由は不明

## 第二段階
# 毎ステップ, それぞれのregionに対して100の熱量(単位はunitsコマンドに従う)を加算/減算
fix             hot all heat 1 100.0 region hot
fix             cold all heat 1 -100.0 region cold

## 第三段階
fix             ave all ave/time 1 1 1000 v_tdiff ave running start 13000

in.ehex

in.heatとほぼ同じで、fix heatコマンドがfix ehexコマンドになっているだけ

in.mp

第二段階からの流れが少し違う少し複雑なアルゴリズムになっているので、詳細は以下のリンク参照

fix thermal/conductivity command — LAMMPS documentation

## 第二段階
compute         ke all ke/atom
variable        temp atom c_ke/1.5

fix             1 all nve

compute         layers all chunk/atom bin/1d z lower 0.05 units reduced
fix             2 all ave/chunk 10 100 1000 layers v_temp file profile.mp

# モデルをz方向に20分割し、Muller-Platheアルゴリズムに従って10ステップごとに1つの原子ペアのz方向の運動エネルギーを交換
fix             3 all thermal/conductivity 10 z 20

variable        tdiff equal f_2[11][3]-f_2[1][3]
thermo_style    custom step temp epair etotal f_3 v_tdiff

thermo          1000
run             20000              # 他のNEMDより多めのステップを取る


## 第三段階
# リセット
fix             3 all thermal/conductivity 10 z 20

fix             ave all ave/time 1 1 1000 v_tdiff ave running
thermo_style    custom step temp epair etotal f_3 v_tdiff f_ave

run             20000

in.heatflux

以下のリンクも参照

compute heat/flux command — LAMMPS documentation

## 変数設定
# x,y,z,rho,t,rc,tlo,thiはNEMDと同じなので略
# Green-Kubo公式では熱流束
variable        p equal 200        # 熱流束の自己相関を計算するときの長さ
variable        s equal 10         # 熱流束をサンプリングする間隔
variable        d equal $p*$s      # 出力間隔


## 初期設定
# NEMDと同じなので略


## 第一段階：熱平衡状態になるまで緩和
# NEMDと同じなので略


## 第二段階：熱伝導率計算

reset_timestep  0                  # タイムステップカウントをリセット

# KE:運動エネルギー、PE:ポテンシャルエネルギー、Stress:圧力を計算
# 上記3つを使って熱流束fluxを計算
compute         myKE all ke/atom
compute         myPE all pe/atom
compute         myStress all stress/atom NULL virial
compute         flux all heat/flux myKE myPE myStress
variable        Jx equal c_flux[1]/vol
variable        Jy equal c_flux[2]/vol
variable        Jz equal c_flux[3]/vol

fix             1 all nve
# 計算した熱流束を元にその時間相関を計算し、その平均を取る
# type autoを指定することで自己相関を計算
fix             JJ all ave/correlate $s $p $d &
                c_flux[1] c_flux[2] c_flux[3] type auto &
                file profile.heatflux ave running

# 熱伝導率を計算するための変数設定
variable        scale equal $s*dt/$t/$t/vol
variable        k11 equal trap(f_JJ[3])*${scale}
variable        k22 equal trap(f_JJ[4])*${scale}
variable        k33 equal trap(f_JJ[5])*${scale}

thermo          $d
thermo_style    custom step temp v_Jx v_Jy v_Jz v_k11 v_k22 v_k33

run             100000             # 精度よく求めるためには長時間のシミュレーションをする必要がある

variable        kappa equal (v_k11+v_k22+v_k33)/3.0
print           "running average conductivity: ${kappa}"

注意事項

格子定数がそのポテンシャルに対して最適かどうか

上記のNEMD/EMD計算中は計算モデルの大きさは固定されていますが、
原子間ポテンシャルに対して大きすぎたり小さすぎたりすると、圧力値がおかしなことになり負の熱伝導率が出力されたりします(EMDで確認)
LAMMPSで構造最適化も可能なので、事前に行っておく必要があります

justice-vsbr.hatenablog.com

熱伝導率計算に入る前に平衡状態になっているかどうか確認すること

サンプルでは流体を取り扱っている(流体と固体どう区別するのかは謎)ので少ないタイムステップで平衡状態になるが、そのまま流用すると平衡状態にならないまま次に進んでしまうことがあります
サンプルスクリプトに準備段階におけるロギングも仕込まれているので、これを元に温度やエネルギーの変化をよく見ておく必要があります

単位系を確認