原子分子シミュレーションPythonパッケージAtomic Simulation Environment (ASE)を使ってみる その３：第一原理分子動力学計算

この記事は、
"原子分子シミュレーションPythonパッケージAtomic Simulation Environment (ASE)を使ってみる　その１"
https://qiita.com/cometscome_phys/items/cd1f4d5f025872dfaae5

"原子分子シミュレーションPythonパッケージAtomic Simulation Environment (ASE)を使ってみる　その２：Quantum Espressoで第一原理計算"
https://qiita.com/cometscome_phys/items/4af134de6d959a7718b9
の続編の記事となります。
その２で導入したMateriApps LIVE!は導入ずみであるとします。

前回はASEを使って、Quantum Espressoで第一原理計算をしました。
今回はASEを使って、Quantum Espressoを用いて第一原理分子動力学計算(第一原理MD)をしようと思います。

分子動力学法

分子動力学法というのは、原子同士に働く力に従って原子を動かしていく手法のことです。
それぞれの原子間にどのような力が働くかは、これまで様々な研究がなされてきました。
近年では、第一原理分子動力学法というものが登場しています。これは、力というものはエネルギーの微分であるために、第一原理計算である原子配置におけるエネルギーを計算できれば、力も計算できる、ということを利用したものです。
ASEでは、ある計算手法がある原子配置を考えた時にエネルギーと力を出すことができるのであれば、それを利用して分子動力学計算ができます。
つまり、第一原理計算ソフトウエアでエネルギーと力を計算できれば、第一原理計算によって得られる原子間の力による分子動力学計算ができます。
これを実際にやってみましょう。

EMTを使った分子動力学計算

ASEのページにもありますが、Pythonで実装された簡単なソルバーEMTを用いて分子動力学(MD)シミュレーションしてみましょう。
https://wiki.fysik.dtu.dk/ase/tutorials/md/md.html
このあとQuantum Espressoを用いて第一原理MDをする関係上、ここで紹介されているチュートリアルのコードよりも計算時間が短くなるようなコードにしておきます。

コードは以下の通りです。チュートリアルのコードとほぼ同じですが、コメントを日本語にしておきました。

md.py

"""一定エネルギーでのMDシミュレーション"""

from ase.lattice.cubic import FaceCenteredCubic
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase.md.verlet import VelocityVerlet
from ase import units


from ase.calculators.emt import EMT
size = 2

# 結晶のセットアップ
atoms = FaceCenteredCubic(directions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]],
                          symbol='Cu',
                          size=(size, size, size),
                          pbc=True)

# Effective Medium Theoryを用いて、原子間の力を評価することとする。
atoms.set_calculator(EMT())

# 温度300Kの運動量分布をセット
MaxwellBoltzmannDistribution(atoms, 300 * units.kB)

# VelocityVerlet アルゴリズムを使って、一定エネルギーでのMD計算をする
dyn = VelocityVerlet(atoms, 5 * units.fs)  # 5 fs time step.


def printenergy(a):
    """ポテンシャル、運動エネルギー、全エネルギーをプリントするFunction"""
    epot = a.get_potential_energy() / len(a)
    ekin = a.get_kinetic_energy() / len(a)
    print('Energy per atom: Epot = %.3feV  Ekin = %.3feV (T=%3.0fK)  '
          'Etot = %.3feV' % (epot, ekin, ekin / (1.5 * units.kB), epot + ekin))

# ここでMD計算をする。
printenergy(atoms)
for i in range(20):
    dyn.run(1)
    printenergy(atoms)

このコードでは、銅原子のMDシミュレーションをしています。
最初、温度が300系の時に原子がもつ速度の分布をマクスウェルボルツマン分布に従ってMaxwellBoltzmannDistributionで生成し、その後はエネルギーが一定になるような時間発展を使ってMDをしています。チュートリアルでは3x3x3のスーパーセルの計算でしたが、2x2x2としています。出てくるエネルギーは、1ステップごとに表示しています。なお、チュートリアルでは10ステップごとに表示していました。

MD_EMTというディレクトリを作って、その中に上のコードをmd.pyとして保存し、それをpython md.pyで実行します。

得られる結果は、

user@malive:~/MD_EMT$ python md.py 
Energy per atom: Epot = -0.006eV  Ekin = 0.044eV (T=340K)  Etot = 0.038eV
Energy per atom: Epot = -0.004eV  Ekin = 0.043eV (T=330K)  Etot = 0.038eV
Energy per atom: Epot = -0.001eV  Ekin = 0.039eV (T=303K)  Etot = 0.038eV
Energy per atom: Epot = 0.004eV  Ekin = 0.034eV (T=262K)  Etot = 0.038eV
Energy per atom: Epot = 0.011eV  Ekin = 0.027eV (T=213K)  Etot = 0.038eV
Energy per atom: Epot = 0.018eV  Ekin = 0.021eV (T=162K)  Etot = 0.038eV
Energy per atom: Epot = 0.023eV  Ekin = 0.015eV (T=116K)  Etot = 0.038eV
Energy per atom: Epot = 0.028eV  Ekin = 0.011eV (T= 82K)  Etot = 0.039eV
Energy per atom: Epot = 0.030eV  Ekin = 0.008eV (T= 64K)  Etot = 0.039eV
Energy per atom: Epot = 0.031eV  Ekin = 0.008eV (T= 61K)  Etot = 0.039eV
Energy per atom: Epot = 0.029eV  Ekin = 0.009eV (T= 73K)  Etot = 0.038eV
Energy per atom: Epot = 0.026eV  Ekin = 0.012eV (T= 93K)  Etot = 0.038eV
Energy per atom: Epot = 0.023eV  Ekin = 0.015eV (T=119K)  Etot = 0.038eV
Energy per atom: Epot = 0.020eV  Ekin = 0.019eV (T=144K)  Etot = 0.038eV
Energy per atom: Epot = 0.017eV  Ekin = 0.021eV (T=165K)  Etot = 0.038eV
Energy per atom: Epot = 0.015eV  Ekin = 0.023eV (T=180K)  Etot = 0.038eV
Energy per atom: Epot = 0.014eV  Ekin = 0.024eV (T=188K)  Etot = 0.038eV
Energy per atom: Epot = 0.014eV  Ekin = 0.025eV (T=191K)  Etot = 0.038eV
Energy per atom: Epot = 0.014eV  Ekin = 0.024eV (T=189K)  Etot = 0.038eV
Energy per atom: Epot = 0.014eV  Ekin = 0.024eV (T=185K)  Etot = 0.038eV
Energy per atom: Epot = 0.015eV  Ekin = 0.023eV (T=180K)  Etot = 0.038eV

このようになります。Epotがポテンシャルエネルギー、Ekinが運動エネルギーを表しています。今回は全エネルギーが保存する時間発展をしていますので、二つの和であるEtotは常に一定になっています。

Quantum Espressoを使った分子動力学計算

次に、Quantum Espressoでエネルギーと力を評価してMDをやってみます。
上のコードでは、計算をEMTでやる、という指示は

# Effective Medium Theoryを用いて、原子間の力を評価することとする。
atoms.set_calculator(EMT())

の部分です。ですので、ここをQuantum Espressoに変えれば、第一原理MDになります。
前回の記事を参考にすると、

# Quantum Espressoを用いて、原子間の力を評価することとする。
pseudopotentials = {'Cu': 'Cu.pbe-kjpaw.UPF '}
calc = Espresso(pseudopotentials=pseudopotentials,
                tprnfor=True,ecutwfc = 20.0,occupations = 'smearing',smearing = 'gauss',degauss = 0.01)
atoms.set_calculator(calc)

となります。なお、ecutwfcはQuantum Espressoのパラメータで、平面波のカットオフエネルギーです。occupationsなどもパラメータで、金属の場合は設定した方が良い、とありました。

残りのコードは同じです。
ですので、全体としては、

md.py

"""一定エネルギーでのMDシミュレーション"""

from ase.lattice.cubic import FaceCenteredCubic
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase.md.verlet import VelocityVerlet
from ase import units
from ase.calculators.espresso import Espresso
from ase.calculators.emt import EMT

size = 2

# 結晶のセットアップ
atoms = FaceCenteredCubic(directions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]],
                          symbol='Cu',
                          size=(size, size, size),
                          pbc=True)


# Quantum Espressoを用いて、原子間の力を評価することとする。
pseudopotentials = {'Cu': 'Cu.pbe-kjpaw.UPF '}
calc = Espresso(pseudopotentials=pseudopotentials,
                tprnfor=True,ecutwfc = 20.0,occupations = 'smearing',smearing = 'gauss',degauss = 0.01)
atoms.set_calculator(calc)

# 温度300Kの運動量分布をセット
MaxwellBoltzmannDistribution(atoms, 300 * units.kB)

# VelocityVerlet アルゴリズムを使って、一定エネルギーでのMD計算をする
dyn = VelocityVerlet(atoms, 5 * units.fs)  # 5 fs time step.


def printenergy(a):
    """ポテンシャル、運動エネルギー、全エネルギーをプリントするFunction"""
    epot = a.get_potential_energy() / len(a)
    ekin = a.get_kinetic_energy() / len(a)
    print('Energy per atom: Epot = %.3feV  Ekin = %.3feV (T=%3.0fK)  '
          'Etot = %.3feV' % (epot, ekin, ekin / (1.5 * units.kB), epot + ekin))

# ここでMD計算をする。
printenergy(atoms)
for i in range(20):
    dyn.run(1)
    printenergy(atoms)

となります。エネルギーと力を出すのに異なるソフトウェアを使っていても同じように書けるのは、ASEの利点です。
これを実行すると、

user@malive:~/MD_QE$ python md.py 
Energy per atom: Epot = -2899.786eV  Ekin = 0.034eV (T=266K)  Etot = -2899.751eV
Energy per atom: Epot = -2899.785eV  Ekin = 0.034eV (T=261K)  Etot = -2899.751eV
Energy per atom: Epot = -2899.783eV  Ekin = 0.032eV (T=245K)  Etot = -2899.751eV
Energy per atom: Epot = -2899.780eV  Ekin = 0.029eV (T=221K)  Etot = -2899.751eV
Energy per atom: Epot = -2899.776eV  Ekin = 0.024eV (T=189K)  Etot = -2899.751eV
Energy per atom: Epot = -2899.771eV  Ekin = 0.020eV (T=154K)  Etot = -2899.751eV
Energy per atom: Epot = -2899.767eV  Ekin = 0.016eV (T=120K)  Etot = -2899.751eV
Energy per atom: Epot = -2899.763eV  Ekin = 0.012eV (T= 92K)  Etot = -2899.751eV
Energy per atom: Epot = -2899.760eV  Ekin = 0.009eV (T= 73K)  Etot = -2899.751eV
Energy per atom: Epot = -2899.759eV  Ekin = 0.008eV (T= 64K)  Etot = -2899.751eV
Energy per atom: Epot = -2899.759eV  Ekin = 0.008eV (T= 65K)  Etot = -2899.751eV
Energy per atom: Epot = -2899.761eV  Ekin = 0.010eV (T= 75K)  Etot = -2899.751eV
Energy per atom: Epot = -2899.763eV  Ekin = 0.012eV (T= 90K)  Etot = -2899.751eV
Energy per atom: Epot = -2899.765eV  Ekin = 0.014eV (T=106K)  Etot = -2899.751eV
Energy per atom: Epot = -2899.767eV  Ekin = 0.016eV (T=121K)  Etot = -2899.751eV
Energy per atom: Epot = -2899.768eV  Ekin = 0.017eV (T=132K)  Etot = -2899.751eV
Energy per atom: Epot = -2899.769eV  Ekin = 0.018eV (T=138K)  Etot = -2899.751eV
Energy per atom: Epot = -2899.769eV  Ekin = 0.018eV (T=139K)  Etot = -2899.751eV
Energy per atom: Epot = -2899.769eV  Ekin = 0.017eV (T=135K)  Etot = -2899.751eV
Energy per atom: Epot = -2899.768eV  Ekin = 0.017eV (T=129K)  Etot = -2899.751eV
Energy per atom: Epot = -2899.767eV  Ekin = 0.016eV (T=121K)  Etot = -2899.751eV

となります。
EMTと違ってEpotの大きさが大きく異なっています。しかし、エネルギーの大きさには意味はありません。また、使う擬ポテンシャルの種類が異なると異なる値が出ることも知られています。重要なのは、エネルギーの差分です。この計算でも、Etotがちゃんと保存していることがわかります。
これで、第一原理MDを実行することができました。

なお、第一原理MDは毎回毎回第一原理計算をするため、非常に重たいです。

Julia言語での実装

もはやPythonとほぼそのままに見えるのでJuliaでやる意義がわからなくなってきますが、PyCallの使い方の例としては役に立つと思います。

Quantum Espressoで第一原理MDをするコードは、

md.jl

#一定エネルギーでのMDシミュレーション
using PyCall
using Printf
const FaceCenteredCubic =  pyimport("ase.lattice.cubic").FaceCenteredCubic
const MaxwellBoltzmannDistribution = pyimport("ase.md.velocitydistribution").MaxwellBoltzmannDistribution
const VelocityVerlet = pyimport("ase.md.verlet").VelocityVerlet
const units = pyimport("ase").units
const Espresso = pyimport("ase.calculators.espresso").Espresso
const EMT = pyimport("ase.calculators.emt").EMT

size = 2

# 結晶のセットアップ
atoms = FaceCenteredCubic(directions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]],
                          symbol="Cu",
                          size=(size, size, size),
                          pbc=true)


# Quantum Espressoを用いて、原子間の力を評価することとする。
pseudopotentials = Dict("Cu" =>"Cu.pbe-kjpaw.UPF")
calc = Espresso(pseudopotentials=pseudopotentials,
                tprnfor=true,ecutwfc = 20.0,occupations = "smearing",smearing = "gauss",degauss = 0.01)
atoms.set_calculator(calc)

# 温度300Kの運動量分布をセット
MaxwellBoltzmannDistribution(atoms, 300 * units.kB)

# VelocityVerlet アルゴリズムを使って、一定エネルギーでのMD計算をする
dyn = VelocityVerlet(atoms, 5 * units.fs)  # 5 fs time step.


function printenergy(a)
    #ポテンシャル、運動エネルギー、全エネルギーをプリントするFunction
    epot = a.get_potential_energy() / length(a)
    ekin = a.get_kinetic_energy() / length(a)
    @printf("Energy per atom: Epot = %.3feV  Ekin = %.3feV (T=%3.0fK)   Etot = %.3feV \n", epot, ekin, ekin / (1.5 * units.kB), epot + ekin)
end

# ここでMD計算をする。
printenergy(atoms)
for i=0:20-1
    dyn.run(1)
    printenergy(atoms)
end

となります。