Python
Julia
ASE
MateriAppsLive!

原子分子シミュレーションPythonパッケージAtomic Simulation Environment (ASE)を使ってみる その2:Quantum Espressoで第一原理計算

以前の記事

"原子分子シミュレーションPythonパッケージAtomic Simulation Environment (ASE)を使ってみる その1"

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)を使ってみる その1"

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.ba[:b]と書いていたものが、juliaでもa.bで良くなりました。

ますますPythonそのままに見えます。