Help us understand the problem. What is going on with this article?

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

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

分子動力学とは

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

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

とりあえずやってみる

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

パッケージを入れる

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

pip install --upgrade --user ase

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

pip install --upgrade --user numpy scipy matplotlib

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

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

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

初期配置を作る

シミュレーションを行うには、まず銅原子の初期配置が必要です。ここでは、銅原子を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()

これを実行すると、
Figure_1.png
が得られます。$3\times3\times3$のFCCが作れてますね!

ポテンシャルを決める

MDを行うには、銅原子間のポテンシャルが分かっている必要があります。今回はEMT(effective medium theory)ポテンシャルを使います。

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

初速度を決める

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

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

運動方程式を解く

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

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

この時、時間刻みは$5$fs(フェムト秒)。ちなみに、$1$fsとは$10^{-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が全エネルギーです。
グラフにしてみると、

Figure_2.png

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

おわりに

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

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

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

(気づいたらASEの宣伝になってた...)

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした