以前の記事
"原子分子シミュレーションPythonパッケージAtomic Simulation Environment (ASE)を使ってみる その１"
https://qiita.com/cometscome_phys/items/cd1f4d5f025872dfaae5
の続きとなります。
今回は、ASE上で第一原理計算をしてみましょう。
第一原理計算ソフトウエアとしては、フリーソフトウエアであるQuantumEspressoを使うことにします。
第一原理計算をして、NaClの格子定数を求めてみます。
Quantum Espressoのインストール
Quantum Espressoは
https://www.quantum-espresso.org
からダウンロードできます。そして、インストールも指示に従えばできると思います。
しかし、この記事では、インストールは行いません。以下で紹介するMateriApps LIVE!を使うことで、インストール作業を省略してしまいましょう。
MateriApps LIVE!
何らかのシミュレーションをするとき、最近は様々なフリーソフトウエアがあるおかげで、誰でも気軽に試すことができるようになってきました。
しかし、様々なフリーソフトウエアを試してみる時にまず最初の障害となるのが、そのソフトウエアのインストール作業です。インストールをするのに失敗して、結局投げ出してしまう、ということはよくあることだと思います。
例えば、第一原理計算ソフトウエアにはQuantum EspressoやABINIT、Open MXなど様々なフリーに使えるソフトウェアがあります。しかし、これらがどんな感じなのかを試そうとした時にそれぞれのインストールに四苦八苦して結局ソフトを試すことができない、ということはあるあるだと思います。
そんな時に便利なものがあります。それは、MateriApps LIVE!です。
https://cmsi.github.io/MateriAppsLive/
このサイトの説明には、
MateriApps LIVE!は、手持ちのノートPCなどを用いて、気軽に計算物質科学アプリケーションを試せる環境を提供します。MateriAppsアプリケーション、OS (Debian GNU/Linux)、エディタ、可視化ツールなど、チュートリアルを始めるのに必要な環境は全て1本のUSBメモリに収められています。
とあります。
これはVirtualBoxという仮想PCソフトウエアのディスクイメージで、MacでもWindowsでもLinuxでも、VirtualBoxをインストールすれば、このディスクイメージを読むことができます。
VirtualBoxでは、Linuxが立ち上がります。
そのLinuxの中に、
ABINIT, AkaiKKR, ALPS, CP2K, DCore, DDMRG, DSQSS, ERmod, feram, Gromacs, HPhi, LAMMPS, mVMC, OpenMX, Quantum ESPRESSO, RESPACK, SALMON, SMASH, TRIQS/CTHYB, TRIQS/DFT tools, xTAPP
が入っています。つまり、すでに有名なソフトウエアがインストール済みです。
ソフトウエアの詳細については
https://github.com/cmsi/MateriAppsLive/wiki/GettingStarted
こちらをみるとわかりやすいと思います。
最初にすること
VirtualBoxのインストールと、ディスクイメージの起動まではできたとします。ディスクイメージはここでは64ビット版のamd64の方を使ったとします。
ユーザー名はuser、パスワードはlive
でログインすることができます。
デフォルトではキーボードが英語配列になっていますので、日本のキーボードを使っている方は、
左下のツバメみたいなアイコンをクリックして、System ToolsからLXTerminalを選んでターミナルを出し、
setxkbmap -layout jp
と入れれば日本語キーボード配列になります。
このままだと起動するたびに英語キーボード配列になるので、立ち上げたターミナルで
echo setxkbmap -layout jp >> .bashrc
と打って.bashrcに設定を書き込んでおきましょう。
あとは、Wikiを見るなり何なりで、色々試すことができます。
#ASEのインストール
以前の記事
"原子分子シミュレーションPythonパッケージAtomic Simulation Environment (ASE)を使ってみる その１"
https://qiita.com/cometscome_phys/items/cd1f4d5f025872dfaae5
で取り上げたAtomic Simulation Environment (ASE)をインストールしてみましょう。
このディスクイメージにはPythonが入っていますので、
pip install ase
とすればASEを簡単に入れることができます。このPythonのバージョンは2.7です。
ここで注意ですが、以下のサンプルコードはPython3系では動きませんでした。Python2系では動きました。
追記:Python3で動かすためには、site-package/are/espresso.pyの292行目のmapをlist()で囲めば良いようです。
サンプルコードの実行
NaCl結晶の格子定数を求めてみましょう。
実験的な格子定数は
http://crystalbase.co.jp/index/b/nacl.html
で、関連する日本語論文は
https://www.jstage.jst.go.jp/article/jsssj/28/3/28_3_144/_pdf
にあります。
https://wiki.fysik.dtu.dk/ase/ase/calculators/espresso.html
にあるサンプルコードはそのままでは動きません。
rocksalt.set_calculator(calc)
が足りません。また、擬ポテンシャルをちゃんと持っているものを指定しなければなりません。
ですので、
from ase.build import bulk
from ase.calculators.espresso import Espresso
from ase.constraints import UnitCellFilter
from ase.optimize import LBFGS
pseudopotentials = {'Na': 'Na.pbe-spn-kjpaw_psl.1.0.0.UPF',
'Cl': 'Cl.pbe-n-rrkjus_psl.1.0.0.UPF'}
rocksalt = bulk('NaCl', crystalstructure='rocksalt', a=6.0)
calc = Espresso(pseudopotentials=pseudopotentials,
tstress=True, tprnfor=True, kpts=(3, 3, 3))
rocksalt.set_calculator(calc)
ucf = UnitCellFilter(rocksalt)
opt = LBFGS(ucf)
opt.run(fmax=0.005)
# cubic lattic constant
print((8*rocksalt.get_volume()/len(rocksalt))**(1.0/3.0))
とします。
ここで、MateriApps LIVE!でのQuantum Espressoの擬ポテンシャルは
/usr/share/espresso/pseudo
にありますので、
cd /usr/share/espresso/pseudo
sudo wget https://www.quantum-espresso.org/upf_files/Na.pbe-spn-kjpaw_psl.1.0.0.UPF
sudo wget https://www.quantum-espresso.org/upf_files/Cl.pbe-n-rrkjus_psl.1.0.0.UPF
として、NaとClの擬ポテンシャルをダウンロードします。なお、他の擬ポテンシャルは
https://www.quantum-espresso.org/pseudopotentials
にあります。
コードでは、格子定数を最初は6として、設定したあと、構造最適化をしています。最後に、構造最適化して出てきた格子定数をプリントしています。
サンプルコードをファイルに保存して、実行すると、
user@malive:~$ python test.py
Step Time Energy fmax
LBFGS: 0 06:20:31 -1960.772207 0.5120
LBFGS: 1 06:20:36 -1960.780670 0.4990
LBFGS: 2 06:20:41 -1960.819410 0.3478
LBFGS: 3 06:20:46 -1960.836529 0.1518
LBFGS: 4 06:20:50 -1960.831494 0.0331
LBFGS: 5 06:20:55 -1960.834038 0.0063
LBFGS: 6 06:21:00 -1960.834046 0.0009
5.659067832564933
こんな感じで出力されて、第一原理計算を用いてNaClの結晶構造の最適化を行うことができました。
Julia言語での実装
最後に、Julia言語での実装も書いておきます。
MateriApps LIVE!のディスクイメージへのインストールは、
wget https://julialang-s3.julialang.org/bin/linux/x64/1.1/julia-1.1.0-linux-x86_64.tar.gz
tar xzvf julia-1.1.0-linux-x86_64.tar.gz
となります。
~/julia-1.1.0/bin/julia
で立ち上げて、さっきまで使っていたPythonを使うように
ENV["PYTHON"]="/usr/bin/python"
とし、
]を押してパッケージモードにして、
add PyCall
としてPyCallをインストールします。
そして、コードは、
using PyCall
const bulk = pyimport("ase.build").bulk
const Espresso = pyimport("ase.calculators.espresso").Espresso
const UnitCellFilter = pyimport("ase.constraints").UnitCellFilter
const LBFGS = pyimport("ase.optimize").LBFGS
pseudopotentials =Dict("Na"=> "Na.pbe-spn-kjpaw_psl.1.0.0.UPF","Cl" => "Cl.pbe-n-rrkjus_psl.1.0.0.UPF")
rocksalt = bulk("NaCl", crystalstructure="rocksalt", a=6.0)
calc = Espresso(pseudopotentials=pseudopotentials,
tstress=true, tprnfor=true, kpts=(3, 3, 3))
rocksalt.set_calculator(calc)
ucf = UnitCellFilter(rocksalt)
opt = LBFGS(ucf)
opt.run(fmax=0.005)
# cubic lattic constant
println((8*rocksalt.get_volume()/length(rocksalt))^(1.0/3.0))
となります。
なお、PyCallのバージョンが上がったためか、これまではPythonのa.bはa[:b]と書いていたものが、juliaでもa.bで良くなりました。
ますますPythonそのままに見えます。