Jouez avec le module MD de ASE

Faisons une simulation MD du cristal Ni en utilisant le module ʻase.md` de ASE (Atomic Simulation Environment).

modèle

Créez un modèle en vrac de 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 en vrac
Ni_bulk= bulk('Ni', 'fcc', a=3.5, cubic=True) 

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

#Utiliser EMT
Ni_bulk.set_calculator(EMT())

Paramètres de simulation MD

Effectuer une simulation NVT (volume / température constants). Maintenez d'abord la température du système à 10K, puis augmentez-la à 500K. Il affiche également la température toutes les 20 étapes.

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)

Lancer la simulation MD

dyn.run(nsteps0)

#Élever la température
dyn.set_temperature(temp1)
dyn.run(nsteps1)

Visualisez les résultats

Visualisez les changements de température, de structure et de fonction de distribution radiale (RDF) à l'aide de 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")

référence

• ASE
• Didacticiel sur la dynamique moléculaire https://wiki.fysik.dtu.dk/ase/tutorials/md/md.html
• modulemd https://wiki.fysik.dtu.dk/ase/ase/md.html#module-ase.md
• GitHub matplotlib/ipympl https://github.com/matplotlib/ipympl
• Qiita
• Simulation de dynamique moléculaire à essayer pour le moment https://qiita.com/tamaki_osamu/items/4bacfcc7dc78a5bfd94e
• Créez une animation avec matplotlib https://qiita.com/yubais/items/c95ba9ff1b23dd33fde2