ボーズ分布
前回に続いてボーズ分布を可視化します。
※やっぱりPCサイトでしか表示されません。。
なんで急にボーズ分布なんだみたいな話は今度別で書く予定…書きたい…
関数以外は前回と同じなので関数のコードだけ載せます。
def plot_BE_distribution(x, beta=1, mu=-1):
f = 1.0 / (np.exp(beta * (x - mu)) - 1.0)
p = beta / (- np.log(1 - np.exp(beta * mu))) * f
source = plt.ColumnDataSource(data=dict(x=x, f=f, p=p))
y_range=(0, 1)
freq = plt.figure(
plot_height=200, plot_width=300, toolbar_location='below',
title='ボーズ分布関数', y_range=y_range)
loglog = plt.figure(
plot_height=200, plot_width=300, toolbar_location='below',
title='ボーズ分布関数（両対数）', x_axis_type='log', y_axis_type='log')
prob = plt.figure(
plot_height=200, plot_width=300, toolbar_location='below',
title='規格化されたボーズ分布関数', y_range=y_range)
freq.line('x', 'f', source=source)
loglog.line('x', 'f', source=source)
prob.line('x', 'p', source=source)
temp_slider = models.Slider(start=0.01, end=1.0, value=beta, step=0.01, title='inverse temperature')
pot_slider = models.Slider(start=-10, end=-0.01, value=mu, step=0.01, title='chemical potential')
callback = models.CustomJS(args=dict(source=source, temp=temp_slider, pot=pot_slider), code="""
const data = source.data;
const beta = temp.value;
const mu = pot.value;
const x = data['x'];
const f = data['f'];
const p = data['p'];
for (var i = 0; i < x.length; i++) {
f[i] = 1.0 / (Math.exp(beta * (x[i] - mu)) - 1.0);
p[i] = beta / (- Math.log(1.0 - Math.exp(beta * mu))) * f[i];
}
source.change.emit();
""")
temp_slider.js_on_change('value', callback)
pot_slider.js_on_change('value', callback)
layout = layouts.row(layouts.column(prob, freq, loglog), layouts.column(temp_slider, pot_slider))
return layout
ガンマ分布
いつの間にか社会人です。
最近業務でいくつかの確率分布を触っているのですが、パラメータが変わると分布がどう変わるのか、パッと知りたかったので
ガンマ分布を例に、PythonのBokehライブラリを使ってインタラクティブなグラフを作ってみました。
実際に動かせるので触ってみてください。
※PCサイトでないと表示されないみたいです。
ちなみにコードは以下のような感じです。JupyterNotebook上で実行してください。
最後のbokeh.plotting.show()の代わりにbokeh.embed.components()を使ってあげればscriptタグとdivタグを吐き出してくれます。
これに加えてブログのヘッダーにBokeh.jsを読み込むコードを追加すればOKです。詳しくは公式のドキュメントを見てください。
import numpy as np
import scipy.special as special
import bokeh
import bokeh.plotting as plt
import bokeh.models as models
import bokeh.layouts as layouts
import bokeh.embed as embed
print(bokeh.__version__)
def plot_gamma_distribution(x, k=2, beta=1):
gamma_y = beta**k / special.gamma(k) * x**(k-1) * np.exp(-beta * x)
power_y = x**(k-1)
exp_y = np.exp(-beta * x)
source = plt.ColumnDataSource(data=dict(x=x, gamma=gamma_y, power=power_y, exp=exp_y))
gamma = plt.figure(plot_height=200, plot_width=300, title='ガンマ分布', y_range=(0, 1))
power = plt.figure(plot_height=200, plot_width=300, title='べき函数部分')
exp = plt.figure(plot_height=200, plot_width=300, title='指数関数部分')
gamma.line('x', 'gamma', source=source)
power.line('x', 'power', source=source)
exp.line('x', 'exp', source=source)
shape_slider = models.Slider(start=0.1, end=10, value=k, step=0.01, title='shape parameter')
scale_slider = models.Slider(start=0.1, end=10, value=beta, step=0.01, title='scale parameter')
callback = models.CustomJS(args=dict(source=source, shape=shape_slider, scale=scale_slider), code="""
function Gamma(x) {
var p = [0.99999999999980993, 676.5203681218851, -1259.1392167224028,
771.32342877765313, -176.61502916214059, 12.507343278686905,
-0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7
];
var g = 7;
if (x < 0.5) {
return Math.PI / (Math.sin(Math.PI * x) * Gamma(1 - x));
}
x -= 1;
var a = p[0];
var t = x + g + 0.5;
for (var i = 1; i < p.length; i++) {
a += p[i] / (x + i);
}
return Math.sqrt(2 * Math.PI) * Math.pow(t, x + 0.5) * Math.exp(-t) * a;
}
const data = source.data;
const k = shape.value;
const beta = scale.value;
const x = data['x'];
const gamma = data['gamma'];
const power = data['power'];
const exp = data['exp'];
for (var i = 0; i < x.length; i++) {
gamma[i] = beta**k / Gamma(k) * x[i]**(k-1) * Math.exp(-beta * x[i]);
power[i] = x[i]**(k-1);
exp[i] = Math.exp(-beta * x[i]);
}
source.change.emit();
""")
shape_slider.js_on_change('value', callback)
scale_slider.js_on_change('value', callback)
layout = layouts.row(layouts.column(gamma, power, exp), layouts.column(shape_slider, scale_slider))
return layout
plt.output_notebook()
fig = plot_gamma_distribution(np.linspace(0, 10, 1001))
plt.show(fig)
JSでのガンマ関数の実装は
こちらのコードをそのままコピペしました。
pipenvで環境構築
気づけば卒業まであと半年もないです。
せっかく自分で開発したこのプログラム、来年以降も研究室で使ってもらえるように、ドキュメントを整備することになりました。
そこで手始めに、実行するための環境構築をしていきます。
自分で開発したPythonプログラムを他の環境でも使えるようにするためには、パッケージの依存関係を書き出す必要があります。
普通であればプロジェクトごとに仮想環境作ってpip freezeなりconda env exportなりで必要なパッケージを書き出すのですが、これまで環境を分けずに適当にいろんなパッケージを入れまくってたのでどれがこのプログラムに必要なのかわからなくなっていました…
調べてみると最近はpipenvなるものが流行ってるみたいなので、AnacondaとpipでごちゃごちゃしてしまったPython環境の再構築がてらこちらに移行してみます。
(Anaconda環境の削除方法はここでは書きません。)
(2018/11/14追記)
勘違いしていたのですが、pipenvではシステム上にないpythonバージョンのビルドはしてくれません。
そこはpyenvが便利なので、結論としては pyenv + pipenv でpythonの環境構築が完成するということになります
pyenvでビルドしたものは自動で探してくれるので、システム上にないバージョンが必要になった場合は、pyenvでビルドしてからpipenv installを実行すればOK。
参考：
pyenv と pipenv による python 環境 - Qiita
pipenv
pythonのバージョン管理とかパッケージ管理とか、色々あってどれ使えばいいかよくわからないですよね。
pyenv + virtualenv は導入までにやることが多くてめんどくさそうなので、
私はこれまでAnacondaを使っていました笑
Anacondaは導入は楽なんですけど、Anaconda Cloudを探しても見つからない場合はpipで入れるしかなく、Anacondaだけで完結しない上にはそれが危険だとか。
pipenvは導入も楽だし裏でpipを使うのでパッケージが衝突することもなさそうです。
公式ドキュメント
Pipenvと仮想環境 — pipenv 2018.10.14.dev0 ドキュメント
Pipenvの基本的な使い方 — pipenv 2018.10.14.dev0 ドキュメント
Pipenvの進んだ使い方 — pipenv 2018.10.14.dev0 ドキュメント
これらに従って使っていきましょう。
install
Mac
root権限を持つローカル環境を想定
$ brew install pyenv $ vim ~/.bashrc # 以下3行追記 export PYENV_ROOT="$HOME/.pyenv" export PATH="$PYENV_ROOT/bin:$PATH" eval "$(pyenv init -)" $ . ~/.bashrc $ brew install pipenv
linux
root権限を持たない、リモートマシンの一般ユーザーを想定
$ git clone git://github.com/yyuu/pyenv.git ~/.pyenv $ vim ~/.bashrc # 以下4行追記 export PYENV_ROOT="$HOME/.pyenv" export PATH="$PYENV_ROOT/bin:$PATH" eval "$(pyenv init -)" export PYTHONPATH= $ . ~/.bashrc $ pyenv install 3.6.7 $ pyenv global 3.6.7 $ pip install pipenv
エラーへの対処
私の場合ははじめにシステムデフォルトのpipを使って
$ pip install --user pipenv $ pipenv install
とした時に、
Traceback (most recent call last):
File "/usr/local/bin/pip", line 7, in <module>
from pip._internal import main
ImportError: No module named _internal
というようなエラーが発生しました。
(恐らくですが)
pip install pipenvで、pipenvに対する依存関係としてインストールされたpipコマンドと、参照されたpipモジュールのバージョンが食い違っていたことが原因です。
以下のように対処しました
- モジュール参照が混乱しないようにデフォルトのPYTHONPATHをクリア(
export PYTHONPATH=) - 適当なバージョンのpythonを自分所有のディレクトリにビルド(
pyenv install 3.6.7) - これをデフォルトに変更(
pyenv global 3.6.7)
create & remove environment
先に必要となるバージョンのpythonをビルドしておきましょう。
カレントディレクトリを対象としてPythonの仮想環境が構築されます。
恐らくディレクトリ単位で環境が管理されています。
デフォルトでは実体は~/.local/share/virtualenvs/に配置されます。
$ pyenv install <VER.> $ pipenv --python <VER.> $ pipenv --rm
activate & deactivate
$ pipenv shell (env_name)$ exit
install package
pipenv installはpip installと完全に互換性があるそうなので、使いやすそうです。
開発環境でのみ使うものは、--devをつけてインストールします。
pipenv install chainer pipenv install --dev jupyter
後処理として、Pipfile, Pipfile.lockを自動で更新してくれます。
Pipfileを使って環境構築
ローカルで作成したPipfileを使って、linuxのリモートサーバに必要な環境を構築してみます。
例：python2.7
on local
$ pyenv install 2.7 $ pipenv --python 2.7 $ pipenv install PKG1 PKG2 ... $ pipenv shell (env_name) $ ./test_script (env_name) $ exit $ git add . $ git commit $ git push
on remote
$ git clone https://github.com/hoge/hoge.git $ cd hoge/ $ pyenv install 2.7 $ pipenv install $ pipenv shell (env_name) $ ./test_script (env_name) $ exit
間違っている・勘違いしているなどあれば遠慮なく指摘してください
NEMDで熱伝導率計算 - その1
3ヶ月ほど放置してましたすみません
前回まではEMDで熱伝導率を計算することを目的としていくつか記事を書いていました
指導教官づてでLAMMPSの熱伝導率計算に詳しい方に聞いたのですが、
EMDはパラメータチューニングがかなりシビアらしいですね。
理論的背景も難解なのでEMDで熱伝導率を計算するのは諦めました。
はい、ということで今回からは
LAMMPSを使ってNEMD（非平衡分子動力学法）で熱伝導率計算する時の知見
をメモしていこうと思います
サンプル
何はともあれまずは公式ドキュメントとサンプルコードを確認しましょう。
以前にも書いた気がしますが、NEMDは適当に流しただけだったので再度見ていきます。
LAMMPSでの熱伝導率計算は大別すると4つの方法があり、
exapmle/KAPPA/ ではlangevinなどとタグ付けされた5つの例が確認できます。
まとめると以下のような感じ
- NEMD
- サーモスタットで温度を制御
langevin - 系に出入りするエネルギーを制御
heatehex - 原子組の運動エネルギーを交換してしまうことで温度勾配を作る
mp
- サーモスタットで温度を制御
- EMD
heatflux
NEMDは1〜3と番号を振ったので、それぞれの詳細は別記事に分けて書いていこうと思います
1〜3の共通点、相違点
実際に系に温度勾配がかかりますが、周期的境界条件を満たすために逆向きの温度勾配もかける必要があります。
熱流方向にn等分した領域を考えて、その中から高温領域hot・低温領域coldを用意します。
NEMDによる熱伝導率計算の流れは共通して以下のようになっています。
- まずは系全体が一定温度になるように熱平衡させる
hot-cold間に温度勾配を発生させて定常状態にさせる（各手法でここが異なる）hot、coldにサーモスタットを当てて一定温度になるように温度制御hot、coldに一定エネルギーを加える/取り除くhotの中で最低温度の原子とcoldの中で最高温度の原子のペアを作り、運動エネルギーを交換する
- 測定を開始して、累積熱流束や温度勾配の時間平均などを求める
- 最終的な熱伝導率を見積もる
ログファイルから必要な値をどうやって読み取るか、はそれぞれの記事で解説していきます。
なお、サンプルのログファイルとREADMEの更新時期がずれているため、
ログファイルから読み取れる値とREADMEに書いてある値が異なっていることに注意してください。
温度プロファイル
一様な物質なら、理想的な定常状態では温度プロファイルは直線的になっているはず。
サンプルコードでは、fix ave/chunkコマンドで熱流方向の温度プロファイルの時間平均を定期的にprofile.**というファイルに出力するようになっています。
- 測定開始直後
- 測定中
- 測定終了直前
の3点の温度プロファイルをそれぞれのサンプルで描画しました。
ちゃんと測定開始時点で定常状態になっていて、直線的な温度勾配になっていることが確認できますね。
注意点
サンプルは、
- 計算対象がレナードジョーンズ流体 (
pair_style lj/cut) - それに伴い、物理量を無次元で表す特殊な単位系を使用している (
units lj)
という特殊な例です。
私のように固体物性を調べたい場合にはそのまま流用できない部分もあると思います。
試行錯誤する中でつまづいたことなども次回以降の記事で書いていくので参考にしてください。
今回はここまで
LAMMPSで熱伝導率計算(実践) その3
タイトルが安定しないです
今日もEMDを使った熱伝導率計算、Green-Kubo公式とにらめっこしました
まずはGreen-Kubo公式の表式と、LAMMPSスクリプトとの対応から
: 熱伝導率
: ボルツマン定数
: 温度
: 計算モデルの体積
: 時間tにおける熱流束
: HCACF(Heat Current Auto-Correlation Function)
## 熱伝導率計算部分 # 熱流束を計算 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 # HCACFを計算 fix JJ all ave/correlate $s $p $d & c_flux[1] c_flux[2] c_flux[3] type auto & file profile.heatflux ave running # xyz各方向の熱伝導率を計算 # trap: HCACF時間積分に相当, scale: 係数部分に相当 variable scale equal ${convert}*$s*dt/$t/$t/vol/${kB} variable k11 equal trap(f_JJ[3])*${scale} variable k22 equal trap(f_JJ[4])*${scale} variable k33 equal trap(f_JJ[5])*${scale} variable kappa equal (v_k11+v_k22+v_k33)/3.0 thermo $d thermo_style custom step temp v_Jx v_Jy v_Jz v_k11 v_k22 v_k33 v_kappa run 1000000 # 1,000,000 steps x 0.8 fs = 800 ps
色々読んでみて、correlation lengthに相当する$pは短すぎるとよくないらしい
dt$s$pが数ps程度になるとよい？
profile.heatfluxにHCACFが出力されているという情報があったのでこう考えたが、間違っている可能性もあります
いまいちHCACFの正体とLAMMPSスクリプトとの対応がわからん。。
fix ave/correlate command — LAMMPS documentation
よく見てみると、trap()関数については以下の記述がありました
The trap function defined for equal-style variables can be used to perform a time integration of this vector of datums, using a trapezoidal rule. This is useful for calculating various quantities which can be derived from time correlation data. If a normalization factor is needed for the time integration, it can be included in the variable formula or via the prefactor keyword.
(自動翻訳)
トラップ用に定義された関数に等しいスタイル変数は、 台形法則を用いて、基準点のこのベクトルの時間積分を実行するために使用することができます。これは、時間相関データから導出できるさまざまな量を計算するのに便利です。時間積分に正規化係数が必要な場合は、変数式またはprefactorキーワードに含めることができます。
つまり、fix ave/correlateで計算したものがHCACFに相当し、時間積分がtrap()関数に相当するようですね
ところで、HCACFは時間経過とともに減衰して0に近づき、これが収束時間の目安になるらしいです
LAMMPSの出力からこの減衰を見るために以下の計算をしました
LAMMPSスクリプト
## 熱伝導率計算 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 fix 1 all nve 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 ${convert}*$s*dt/$t/$t/vol/${kB} variable k11 equal trap(f_JJ[3])*${scale} variable k22 equal trap(f_JJ[4])*${scale} variable k33 equal trap(f_JJ[5])*${scale} variable kappa equal (v_k11+v_k22+v_k33)/3.0 variable h11 equal sum(f_JJ[3]) variable h22 equal sum(f_JJ[4]) variable h33 equal sum(f_JJ[5]) variable HCACF equal (v_h11+v_h22+v_h33)/3.0 thermo $d thermo_style custom step temp v_h11 v_h22 v_h33 v_HCACF v_k11 v_k22 v_k33 v_kappa run 10000000 print "running average conductivity: ${kappa}"
これで標準出力にthermo_styleで指定した熱力学データが出力されます
v_HCACFは累積和になっているので、該当部分を切り出してthermo.datというファイルに保存し、gnuplotで以下のスクリプトを使ってグラフにしました
## normalizeするために最初のデータを参照 # 0 100 -0.014850832 -0.0038024632 -0.00058573814 29083.279 1906.6538 45.242732 10345.059 18.835168 1.2348039 0.029300495 6.6997575 # 10000 98.637033 -0.0018503786 0.017331996 -0.0073324707 12161.638 17045.085 9943.1863 13049.97 11.653242 17.012881 8.619362 12.428495 set logscale x set term png set output "hcacf_xyz.png" oldx = 0 oldy = 0 oldz = 0 plot "thermo.dat" using ($1*0.0008):(dx = ($6-oldx)/(12161.638-29083.279), oldx = $6, dx) w l title "HCACF x", \ "thermo.dat" using ($1*0.0008):(dy = ($7-oldy)/(17045.085-1906.6538), oldy= $7, dy) w l title "HCACF y", \ "thermo.dat" using ($1*0.0008):(dz = ($8-oldz)/(9943.1863-45.242732), oldz = $8, dz) w l title "HCACF z" set output "hcacf.png" old = 0 plot "thermo.dat" using ($1*0.0008):(dx = ($9-old)/(13049.97-10345.059), old = $9, dx) w l title "HCACF" set output "kappa.png" plot "thermo.dat" using ($1*0.0008):13 w l title "kappa"
HCACFは確かに減衰していますが、文献ではこんなに振動したり負の値になったりせずなめらかに減衰しているので何かがおかしいですね。。
確かめなければいけないパラメータがいくつもあるので骨が折れそうです
LAMMPSで熱伝導率計算(実践) - EMD
前回
GaNで熱伝導率計算をしてみよう
SWポテンシャルの下で構造最適化したw-GaNに対して、EMDで熱伝導率を解析してみます。
とりあえず実行
## 変数設定 variable p equal 200 # correlation length variable s equal 10 # sample interval variable d equal $p*$s # dump interval variable t equal 100 # temperature variable tstep equal 0.001 # timestep variable kB equal 1.3806504e-23 variable eV2J equal 1.6021766e-19 variable ps2s equal 1.0000000e-12 variable ang2m equal 1.0000000e-10 variable convert equal ${eV2J}*${eV2J}/${ps2s}/${ang2m} ## 基本設定 units metal dimension 3 boundary p p p atom_style atomic newton on timestep ${tstep} ## 構造データ読み込み、初期速度設定 read_data GaN.sw-optimized velocity all create $t 87287 mom yes rot yes dist gaussian ## ペアポテンシャル pair_style sw pair_coeff * * GaN.sw Ga N ## 熱平衡状態計算 fix 1 all nvt temp $t $t 0.01 drag 0.2 thermo 1000 run 10000 velocity all scale $t unfix 1 ## 熱伝導率計算：略 #=> kappa = -15.4177052995385 (W/mK)
🤔🤔🤔
(c.f.100KにおけるGaNの熱伝導率の実験値は100(W/mK)のオーダー)
とりあえずMDのtimestepを1fs~10fsで変えながら計算してみたところ、以下のようになりました(9,10fsの場合はtimestep大きすぎて系が爆発したので計算できず)
収束も何もあったもんじゃないですね笑
一応4fsではそれらしい値を示していますが偶然だと思います
もちろん他にもいじるべきパラメータがたくさんあるので、どこをどうしたらちゃんと計算できるのか、これから調べます
特に、Green-Kubo公式を解くために熱流束の自己相関関数のアンサンブル平均(もう既にナンノコッチャですが)を計算しますが、この計算におけるパラメータ(p,s,d)がどういう役割を果たすのかわかってないです
読み物
https://arxiv.org/pdf/1005.4478.pdf
熱伝導率を求めるGreen-Kubo公式の表式は以下の通り
係数はいいとして、重要なのは0~∞の時間積分と、被積分関数である熱流束の自己相関関数のアンサンブル平均です
シミュレーションでは有限時間しか計算できないので上の式を厳密には計算できず、なるべく長い時間(タイムステップ数)で計算することで真の値に収束させる方針になります
この論文では金属であるSiやGeを対象としているため、デバイ温度(この温度以上では電子による熱伝導への影響を無視できるという指標)以上の1000Kを設定温度にしていますが、GaNは半導体でありもともと格子振動だけが熱キャリアと考えられます
HCACFが素早く減衰するのは低波長フォノンの寄与が大きいことを示すらしいですが、存在するフォノンの波長は計算モデルの大きさに制限されます。論文中では計算モデルを大きくしても同じように減衰したことから、計算モデルの大きさとタイムステップ数はここで十分、としています
各アンサンブル平均におけるHCACFの標準偏差と平均からnumerical errorを見てみると、ある時点(tc)から急激に大きくなります。よってあまりに長時間計算してもこのnumerical errorが累積してしまいよくない。さらにこのtcは初期速度や計算モデルサイズにより変化し、決めうちでタイムステップ数を決めるのも難しい…
有限サイズ効果もあるため、計算モデルサイズを少しずつ大きくし、サチるところを探しています
読み取れた計算条件は
- タイムステップ: 0.8fs
- 平衡状態計算: 105 steps
- 熱伝導率計算: 2 x 106 steps
- サンプリング: every 1 step
- correlation length: 103
- 試行回数: 8
ひとまずここまで
進捗があり次第更新します
LAMMPSで構造最適化
原子間ポテンシャルに従ってMD計算をするだけでなく構造最適化もできるらしい
原子位置はMD計算中に緩和されるとしても、格子定数の最適化は別で行うしかなさそう？
キーとなるコマンド
minimize: 構造最適化を実行fix box/relax: 外圧を加え、計算モデルの大きさと形も最適化するように変更
Wurtzite-GaNを構造最適化(失敗)
から取得できるw-GaN構造を初期構造として使用
LAMMPSのpotentialsディレクトリにいくつかプリセットがあるのでその中でGaN.swを使用して最適化してみます
## 基本設定 units metal dimension 3 boundary p p p atom_style atomic newton on ## 計算モデル作成 variable a equal 3.216 variable c equal 5.240 variable b_x equal $a*(-0.5) variable b_y equal $a*sqrt(3.0)/2.0 variable 1_3 equal 1.0/3.0 variable 2_3 equal 2.0/3.0 lattice custom 1.0 & a1 $a 0.0 0.0 & a2 ${b_x} ${b_y} 0.0 & a3 0.0 0.0 $c & basis ${1_3} ${2_3} 0.0 & basis ${2_3} ${1_3} 0.5 & basis ${1_3} ${2_3} 0.375 & basis ${2_3} ${1_3} 0.875 region myreg block 0 5 0 5 0 5 create_box 2 myreg create_atoms 2 box & basis 1 1 & basis 2 1 & basis 3 2 & basis 4 2 ## 原子をGa,Nに設定 mass 1 69.72 mass 2 14.01 ## 原子間ポテンシャルをStillinger-Weberに設定 pair_style sw pair_coeff * * GaN.sw Ga N ## 構造最適化実行 timestep 0.0001 # 0.1fs thermo_style custom step etotal temp press lx ly lz vol thermo 1000 # aniso: xyzそれぞれに等しく外圧を加える&長さ比は保存しない # couple: 指定した軸の長さ比を保存 # 外圧0, 1回のiterationで変化する体積は最大1% fix 1 all box/relax aniso 0.0 couple xy vmax 0.01 # エネルギーによる収束条件なし, 原子間力他の収束条件を使用 # 最大で10,000step, 1,000,000eval minimize 0.0 1.0e-20 10000 1000000
実行ログより一部抜粋
LAMMPS (16 Mar 2018)
Lattice spacing in x,y,z = 4.824 2.78514 5.24
Created orthogonal box = (0 0 0) to (24.12 13.9257 26.2)
...
Step TotEng Temp Press Lx Ly Lz Volume
0 -2877.0956 0 73125.066 24.12 13.925688 26.2 8800.2553
1000 -3100.5642 0 212.02536 24.548929 14.173331 26.081804 9074.9054
2000 -3100.7013 0 0.40261942 24.562407 14.181112 26.076651 9083.0776
3000 -3100.7028 0 -0.00062158384 24.562064 14.180914 26.077135 9082.9924
4000 -3100.7028 0 -3.3787783e-06 24.562064 14.180914 26.077135 9082.9924
5000 -3100.7028 0 5.4301294e-10 24.562064 14.180914 26.077135 9082.9924
6000 -3100.7028 0 3.1547809e-09 24.562064 14.180914 26.077135 9082.9924
7000 -3100.7028 0 8.1352718e-10 24.562064 14.180914 26.077135 9082.9924
8000 -3100.7028 0 1.5203109e-09 24.562064 14.180914 26.077135 9082.9924
9000 -3100.7028 0 3.0510664e-09 24.562064 14.180914 26.077135 9082.9924
10000 -3100.7028 0 3.1167627e-09 24.562064 14.180914 26.077135 9082.9924
Loop time of 2.77113 on 144 procs for 10000 steps with 750 atoms
99.5% CPU use with 144 MPI tasks x no OpenMP threads
Minimization stats:
Stopping criterion = max iterations
Energy initial, next-to-last, final =
-2877.0955666 -3100.70275883 -3100.70275883
Force two-norm initial, final = 727.066 4.09501e-11
Force max component initial, final = 450.072 3.9289e-11
Final line search alpha, max atom move = 1 3.9289e-11
Iterations, force evaluations = 10000 19987
なお、LAMMPSでは実行時に計算モデルが直方体となるように自動で変換されます
今回は初期構造として底面の辺=a, 高さ=cの菱形柱の単位格子を設定したが、
x=3a/2, y=√3a/2, z=cの直方体に変換されてます
最適化前後の構造を可視化したところ、以下のようになりました(OVITOで可視化)
不自然な曲がり方をしている。。
上の曲がった構造を周期的境界条件(PBC)に従って拡張すると原子同士が重なったりぽっかり穴が開いたり……もしかしてPBCが動作してない？
うまくいかない原因は…？
サンプルコードがどこかにないかよく探しました
なんだかんだでこれが一番大事で近道だったりします！！
今回はLAMMPSのexamplesディレクトリの中から以下の条件に当てはまるサンプルを探しました。
minimizeコマンドを使用pair_style swコマンドを使用
最終的な目的は違いますが、neb/in.neb.sivacが当てはまったのでこれを参考にしました
minimizeと合わせるとPBCが機能しない？
=>サンプルコードを動かしたところ、物質はGaNではなくSiだが構造が崩れることなくちゃんと動作しているので違う初期状態で速度を持ってしまっている？
=>velocityコマンドで明示的に速度0にしたが変わらずfix box/relaxコマンドで与えた外圧が悪さしている？
=>minimizeコマンドだけでも構造が崩れるので違うStillinger-WeberポテンシャルではWurtziteを再現できない？
=>むしろWurtziteのような4配位構造を再現するためのポテンシャルtimestepやmaxstep、
minimizeコマンドの収束条件などが不適？ =>いくつか試したが変わらず初期構造がオカシイ
=>自分で直方体の単位格子を取り直し、初期構造を設定したところちゃんと動作した！！
Wurtzite-GaNを構造最適化(今度こそ成功)
LAMMPSスクリプト
# 基本設定 units metal atom_style atomic boundary p p p # 単位格子 variable a equal 3.216 variable b equal $a*sqrt(3.0) variable c equal 5.240 variable 1_2 equal 1.0/2.0 variable 2_3 equal 2.0/3.0 variable 1_6 equal 1.0/6.0 variable 3_8 equal 3.0/8.0 variable 7_8 equal 7.0/8.0 lattice custom 1.0 & a1 $a 0.0 0.0 & a2 0.0 $b 0.0 & a3 0.0 0.0 $c & basis 0.0 0.0 0.0 & basis ${1_2} ${1_2} 0.0 & basis 0.0 ${2_3} ${1_2} & basis ${1_2} ${1_6} ${1_2} & basis 0.0 0.0 ${3_8} & basis ${1_2} ${1_2} ${3_8} & basis 0.0 ${2_3} ${7_8} & basis ${1_2} ${1_6} ${7_8} region myreg block 0 5 0 5 0 5 create_box 2 myreg create_atoms 2 box & basis 1 1 & basis 2 1 & basis 3 1 & basis 4 1 & basis 5 2 & basis 6 2 & basis 7 2 & basis 8 2 mass 1 69.72 mass 2 14.01 # ペアポテンシャルを設定 pair_style sw pair_coeff * * GaN.sw Ga N # 構造最適化 thermo_style custom step etotal temp press lx ly lz vol thermo 10 fix 1 all box/relax aniso 0.0 vmax 0.01 minimize 0.0 1.0e-20 1000 10000
最適化前後の構造は以下のようになりました(格子定数が変わっただけなのでぱっと見は変化なし
SWポテンシャルの下で構造最適化されたw-GaNの格子定数はそれぞれ以下のようになりました
a: 3.216 -> 3.190597
c: 5.240 -> 5.2102232
以上です！！
