ASEのMDモジュールで遊ぶ

ASE (Atomic Simulation Environment) のase.mdモジュールを用いて，Ni結晶のMDシミュレーションをしてみます．

モデル

Ni (fcc) のバルクモデルを作成します．

import numpy as np
from ase.build import bulk
from ase.build.supercells import make_supercell
from ase.calculators.emt import EMT

# Niバルク
Ni_bulk= bulk('Ni', 'fcc', a=3.5, cubic=True) 

# 超格子
Ni_bulk = make_supercell(Ni_bulk, np.diag([3., 3., 3.]))

# EMTを用いる
Ni_bulk.set_calculator(EMT())

MDシミュレーションの設定

NVTシミュレーション (体積・温度一定) を行います．
はじめに系の温度を10Kに保ち，その後500Kに上げます．
また20ステップごとに温度の出力も行います．

from ase import units
from ase.md.nvtberendsen import NVTBerendsen
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution

dt = 2 * units.fs
temp0, nsteps0 = 10, 200
temp1, nsteps1 = 500, 400
taut = 20*units.fs
    
MaxwellBoltzmannDistribution(Ni_bulk, temp0*units.kB)
dyn = NVTBerendsen(Ni_bulk, dt, temp0, taut=taut, trajectory='md.traj')
def myprint():
    print(f'time={dyn.get_time() / units.fs: 5.0f} fs ' + \
          f'T={Ni_bulk.get_temperature(): 3.0f} K')
dyn.attach(myprint, interval=20)

MDシミュレーションを実行

dyn.run(nsteps0)

# 温度を上げる
dyn.set_temperature(temp1)
dyn.run(nsteps1)

結果を可視化する

温度・構造・動径分布関数 (RDF) の変化をmatplotlibを用いて可視化します．

%matplotlib widget 

import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from ase.visualize.plot import plot_atoms
from ase.io.trajectory import Trajectory
from ase.geometry.analysis import Analysis

traj =  Trajectory('md.traj')

fig, ax = plt.subplots(1, 3, figsize=(9,3), tight_layout=True)

t = np.arange(nsteps0+nsteps1+1) * dt
temp = [atoms.get_temperature() for atoms in traj]

nframes = 20

def update(iframe):
    idx = int((nsteps0+nsteps1)*iframe/nframes)
    
    ax[0].clear()
    ax[0].set_title('Temperature')
    ax[0].set_xlabel('time (fs)')
    ax[0].set_ylabel('T (K)')
    ax[0].plot(t, temp)
    ax[0].plot(t[idx], temp[idx], marker='X', markersize=10)
    
    ax[1].clear()
    ax[1].set_title('Structure')
    ax[1].axis('off')
    plot_atoms(traj[idx], ax=ax[1], rotation='45x,45y')
    
    distribution, distance = Analysis(traj[idx]).get_rdf(rmax=5., nbins=100, return_dists=True)[0]
    ax[2].clear()
    ax[2].set_title('RDF')
    ax[2].set_ylim((0,10))
    ax[2].set_xlabel('distance (A))')
    ax[2].set_ylabel('distribution')
    ax[2].plot(distance, distribution, color='darkblue')

ani = FuncAnimation(fig, update, np.arange(nframes), blit=True, interval=250.)
ani.save('ani.gif', writer="imagemagick")

参考

• ASE
  • Molecular Dynamicsのチュートリアル https://wiki.fysik.dtu.dk/ase/tutorials/md/md.html
  • mdモジュール https://wiki.fysik.dtu.dk/ase/ase/md.html#module-ase.md
• GitHub matplotlib/ipympl https://github.com/matplotlib/ipympl
• Qiita
  • とりあえずやってみる分子動力学シミュレーション
    https://qiita.com/tamaki_osamu/items/4bacfcc7dc78a5bfd94e
  • matplotlib でアニメーションを作る https://qiita.com/yubais/items/c95ba9ff1b23dd33fde2