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())

Ni_bulk.png

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")

anim_Ni.gif

reference

Recommended Posts

Play with ASE MD module
Play with PIR sensor module [DSUN-PIR]
Play with Prophet
Play around with the pythonista3 ui module
Play with PyTorch
Play with 2016-Python
Play with CentOS 8
Play with Pyramid
Play with Fathom
Play with Othello (Reversi)
Let's play with 4D 4th
Play with reinforcement learning with MuZero
Play with push notifications with imap4lib
Play around with Linux partitions
Let's play with Amedas data-Part 4
Play with Jupyter Notebook (IPython Notebook)
[Python] Play with Discord's Webhook.
Try MD simulation with ANN potential using aenet and ASE
Play RocketChat with API / Python
Let's play with Amedas data-Part 3
Let's play with Amedas data-Part 2
Play with A3RT (Text Suggest)
Play with numerical calculation of magnetohydrodynamics
Play with a turtle with turtle graphics (Part 1)
Play with Poincare series and SymPy
Let's play with Excel with Python [Beginner]
Sort names with Python's random module
Play with Pythonista UI implementation [Action implementation]
Play around with Linux partitions ~ Continued ~
Detects people with motion sensor module
Spark play with WSL anaconda jupyter (2)
Play with Turtle on Google Colab
Play with demons because it's setsubun