Let's perform MD simulation of Ni crystal using ʻase.md` module of ASE (Atomic Simulation Environment).
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())
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)
dyn.run(nsteps0)
#Raise the temperature
dyn.set_temperature(temp1)
dyn.run(nsteps1)
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")
md
module https://wiki.fysik.dtu.dk/ase/ase/md.html#module-ase.mdRecommended Posts