Edited at

原子分子の動きをパソコンでみてみよう:Atomic Simulation Environmentで遊んでみた


ASE (Atomic Simulation Environment)とは

量子化学計算や、古典分子動力学法をpython上で実行できるライブラリです。以下のサイトを参照していただければ、どのようなコンセプトで作られたか理解できると思います。

https://wiki.fysik.dtu.dk/ase/about.html


インストールのしかた

以下のURLを参考にしました:

https://wiki.fysik.dtu.dk/ase/install.html

私は、git clone -b 3.18.1 https://gitlab.com/ase/ase.gitで問題なくできました。


古典分子動力学をやってみる


まず初めに物理系を作成する

テスト系として窒素分子(N2)を作成してみます。ASEは基本的にase.Atomsに物理系の座標情報を入れていきます。

from ase import Atoms

d = 1.1
molecule = Atoms('2N', positions=[(0., 0., 0.), (0., 0., d)])

moleculeにちゃんと座標データが入っているか確認してみましょう。そのためには、ase.ioを使います:

from ase import io

write('N2.pdb',molecule)

.pdbと出力ファイルを指定するとPDBファイル形式を取り扱える。これは生物系の人には便利)

すると以下のような構造が得られます(PyMOLか何かで確認)。確かに窒素の2原子分子が作成されているようです。

aaa.png


分子動力学計算の設定

以下の設定で計算を走らせます

- 積分器:Velocity Verlet

- ポテンシャル:Effective medium theory (EMT)、(https://wiki.fysik.dtu.dk/ase/ase/calculators/emt.html#module-ase.calculators.emt)

- 時間刻み:0.01 fs

- 20 stepsごとに構造の保存

from ase.md.verlet import VelocityVerlet

from ase import units
from ase.calculators.emt import EMT
molecule.set_calculator(EMT())
dyn = VelocityVerlet(molecule, dt=0.01 * units.fs)
for i in range(100):
pot = molecule.get_potential_energy()
kin = molecule.get_kinetic_energy()
print('%2d: %.5f eV, %.5f eV, %.5f eV' % (i, pot + kin, pot, kin))
H.append(pot + kin)
K.append(kin)
U.append(pot)
dyn.run(steps=20)
write(f'traj.pdb',molecule,append=True)

# エネルギーの図作成
import matplotlib.pyplot as plt
%matplotlib inline
plt.plot(H, c='black')
plt.plot(K, c='red')
plt.plot(U, c='blue')
plt.savefig('energy.png')


得られたトラジェクトリのエネルギーの確認

下図をみていただくと、たしかに、全エネルギーは保存しており、NVEアンサンブルになっているようです(黒線が全エネルギー、赤線が運動エネルギー、青線がポテンシャルエネルギー)。

energy.png


感想

かなり簡単に計算結果を得ることができました。オブジェクト指向で書かれているのでソフトウェアの設計の勉強にもなりました。このような形でとある過去の遺産をラップできたらいいのになあ。とりあえず基本的な使い方はりかいしましたが、グラフェンの計算など気になることがあるので、もう少し遊んでみようと思います。