とりあえずやってみる分子動力学シミュレーション

2020/03/29

この記事はQiita の記事を加筆修正したものです。


この記事では、分子動力学について簡単に説明して、簡単にシミュレーションを行ってみます。シミュレーションには、Atomic Simulation Environmentという Python モジュールを使用することにします。

分子動力学とは

分子動力学(Molecular Dynamics: MD)は、原子や分子間の相互作用(ポテンシャル)をもとに運動方程式を解いて行くことで、原子や分子の動的過程を得ます。

私たちの身の回りにあるものは、細かく見ていけば原子や分子の集合体です。分子動力学は、この原子や分子の動きを計算することで集合体としての振る舞いを理解しようとするものです。

とりあえずやってみる

難しい理論は置いておいて、とりあえずやってみましょう。ここでは、Atomic Simulation Environment(ASE)という便利なパッケージを使うことにします。

パッケージを入れる

Python が必要なので適宜環境を構築してください。Python 環境が整ったら、

pip install --upgrade --user ase

で入ります。もし、numpy, scipy, matplotlib がない場合はこの 3 つもインストールします。

pip install --upgrade --user numpy scipy matplotlib

これだけで準備は完了です。

銅(Cu)のシミュレーションをやってみる

銅原子(Cu)のシミュレーションをやってみましょう。 (ここで上げる内容はチュートリアルとほとんど同じです。)

初期配置を作る

シミュレーションを行うにはまずはじまり(初期配置)を設定してやる必要があります。ここでは、銅原子をFCC(面心立方格子)上に配置します。FCC とは結晶構造の一種で、銅の結晶はこの FCC の構造を持つことが知られています。

ASE を利用すれば簡単に FCC 上に銅原子を配置できます。

import matplotlib.pyplot as plt
from ase.visualize.plot import plot_atoms
from ase.lattice.cubic import FaceCenteredCubic

size = 3
atoms = FaceCenteredCubic(directions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]],
                          symbol="Cu",
                          size=(size, size, size),
                          pbc=True)
plot_atoms(atoms, rotation=('0x,0y,0z'))
plt.show()

これを実行すると、

が得られます。3×3×33\times3\times3の FCC が作れてますね!

初速度を決める

シミュレーションを行うには、初期配置だけでなく、初速度も決める必要があります。ここでは、温度300kB300k_Bのマクスウェル=ボルツマン分布に従い速度を決定します。

from ase import units
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
MaxwellBoltzmannDistribution(atoms, 300 * units.kB)

ASE はこれだけで初速度を設定してくれます。

ポテンシャルを決める

MD を行うには、銅原子間にどうような力が働くのか、すなわちポテンシャルが分かっている必要があります。このポテンシャルには、実験から求めた経験的なポテンシャルや量子力学を駆使し電子の振る舞いからポテンシャルを記述する非経験的なポテンシャルなど、研究用途に応じ実に様々なものが利用されています。

今回はEMT(effective medium theory)というポテンシャルを用いることにします。これも ASE で簡単に定義することができます。

from ase.calculators.emt import EMT
atoms.set_calculator(EMT())

はい、これでおわりです。

運動方程式を立てる、計算する

最後に、今まで決めてきた銅原子をどのように計算していくかを決めます、このとき、

  • アンサンブル(統計集団)
  • 微分方程式(運動方程式)の解法

を決める必要があります。アンサンブルとは、原子や分子の集団のことで集団に課せられる制約(例えば温度一定とか)により様々なアンサンブルに分けることができます。微分方程式の解法とは、コンピュータでいかに微分を解くかという数値計算の技法のことでこれにも様々な技法が存在します。

ここでは、最も簡単なミクロカノニカルアンサンブル(粒子数NN, 体積VV, エネルギーEEが一定)に従う運動方程式を、速度ベルレ法で解きます。

from ase import units
from ase.md.verlet import VelocityVerlet
dyn = VelocityVerlet(atoms, 5 * units.fs)

この時、時間刻みは55fs(フェムト秒)。ちなみに、11fs とは101510^{-15}s のことです。

MD を実行する。

これまでの準備をもとに 実施に 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

# 初期配置をつくる(FCC)
size = 3
atoms = FaceCenteredCubic(directions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]],
                          symbol="Cu",
                          size=(size, size, size),
                          pbc=True)
# ポテンシャルにEMT(effective medium theory)を使う
atoms.set_calculator(EMT())
# 300kbのマクスウェル=ボルツマン分布に従う運動量を設定する
MaxwellBoltzmannDistribution(atoms, 300 * units.kB)
# 速度ベルレ法でNVE一定のMD計算をする
dyn = VelocityVerlet(atoms, 5 * units.fs)


def printenergy(a=atoms):  # ポテンシャルエネルギー、運動エネルギーの出力
    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計算
dyn.attach(printenergy, interval=10)
dyn.run(1000)

上記を実行すると、次のような出力が得られます。

Energy per atom: Epot = -0.006eV  Ekin = 0.044eV (T=340K)  Etot = 0.038eV
Energy per atom: Epot = -0.006eV  Ekin = 0.044eV (T=340K)  Etot = 0.038eV
Energy per atom: Epot = 0.029eV  Ekin = 0.010eV (T= 76K)  Etot = 0.038eV
.....

Epotがポテンシャルエネルギー、Ekinが運動エネルギー、Etotが全エネルギーです。 グラフにしてみると、

が得られます。全エネルギーがよく保存していて、NVE 一定のシミュレーションが正しくできていることが分かります。

おわりに

今回のシミュレーションは、非常に小さく簡単な系に対して非常に短い時間行いました。

より実践的なシミュレーションを行うには様々な工夫が必要になります。例えば、今回は NVE 一定のシミュレーションを行いましたが、一般的に実験系に近い NVT 一定(カノニカルアンサンブル)のシミュレーションを行いたい場合の方が多いです。NVT 一定のシミュレーションを行うには、Nosé–Hoover thermostatLangevin dynamicsなどの方法を用いる必要があります。 あるいは、ポテンシャルに第一原理計算を利用する第一原理分子動力なども盛んに行われています。

実のところ ASE では大変ありがたいことにこれらも簡単に試すことができます。ぜひ皆さんも一度遊んでみてはどうでしょう。

logo

たまおさ

釣りとか登山とか好きです。(@tamaki_osamu)