Play with ASE MD module

Let's perform MD simulation of Ni crystal using ʻase.md` module of ASE (Atomic Simulation Environment).

model

Create a bulk model of 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 bulk
Ni_bulk= bulk('Ni', 'fcc', a=3.5, cubic=True) 

#Superlattice
Ni_bulk = make_supercell(Ni_bulk, np.diag([3., 3., 3.]))

#Use EMT
Ni_bulk.set_calculator(EMT())

MD simulation settings

Perform NVT simulation (constant volume / temperature). First keep the system temperature at 10K and then raise it to 500K. It also outputs the temperature every 20 steps.

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)

Run MD simulation

dyn.run(nsteps0)

#Raise the temperature
dyn.set_temperature(temp1)
dyn.run(nsteps1)

Visualize the results

Visualize changes in temperature, structure, and radial distribution function (RDF) using 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")

reference

• ASE --Molecular Dynamics Tutorial https://wiki.fysik.dtu.dk/ase/tutorials/md/md.html --md module https://wiki.fysik.dtu.dk/ase/ase/md.html#module-ase.md
• GitHub matplotlib/ipympl https://github.com/matplotlib/ipympl
• Qiita --Molecular dynamics simulation to try for the time being https://qiita.com/tamaki_osamu/items/4bacfcc7dc78a5bfd94e --Make animation with matplotlib https://qiita.com/yubais/items/c95ba9ff1b23dd33fde2