Spielen Sie mit dem MD-Modul von ASE

Lassen Sie uns eine MD-Simulation von Ni-Kristallen mit dem Modul ase.md von ASE (Atomic Simulation Environment) durchführen.

Modell-

Erstellen Sie ein Massenmodell von 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 Masse
Ni_bulk= bulk('Ni', 'fcc', a=3.5, cubic=True) 

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

#Verwenden Sie EMT
Ni_bulk.set_calculator(EMT())

MD-Simulationseinstellungen

Führen Sie eine NVT-Simulation durch (konstantes Volumen / konstante Temperatur). Halten Sie zuerst die Temperatur des Systems bei 10 K und erhöhen Sie sie dann auf 500 K. Außerdem wird die Temperatur alle 20 Schritte ausgegeben.

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)

Führen Sie die MD-Simulation aus

dyn.run(nsteps0)

#Erhöhen Sie die Temperatur
dyn.set_temperature(temp1)
dyn.run(nsteps1)

Visualisieren Sie die Ergebnisse

Visualisieren Sie Änderungen in Temperatur, Struktur und radialer Verteilungsfunktion (RDF) mithilfe von 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")

Referenz

• ASE
• Tutorial zur Molekulardynamik https://wiki.fysik.dtu.dk/ase/tutorials/md/md.html --md Modul https://wiki.fysik.dtu.dk/ase/ase/md.html#module-ase.md
• GitHub matplotlib/ipympl https://github.com/matplotlib/ipympl
• Qiita
• Molekulardynamiksimulation vorerst versuchen https://qiita.com/tamaki_osamu/items/4bacfcc7dc78a5bfd94e
• Erstellen Sie eine Animation mit matplotlib https://qiita.com/yubais/items/c95ba9ff1b23dd33fde2