Lecture 7: Molecular Dynamics

9. Lecture 7: Molecular Dynamics#

9.1. Overview#

Questions

  • How can I use molecular dynamics (MD) to evolve a system over time?

  • How do I track system properties over a trajectory?

  • What information can MD give us about quantum emitters?

Objectives

  • Set up and run NVE and NVT molecular dynamics in ASE

  • Attach observers to log properties during a simulation

  • Analyse temperature and energy fluctuations

9.2. Lecture Slides#

The slides for this tutorial are embedded below. 📥 Download slides (.pptx)  |  Open in full screen

How to embed: Upload the .pptx to PowerPoint Online (SharePoint/OneDrive) → File → Share → Embed → copy the src="..." URL and replace SHAREPOINT_EMBED_URL_HERE above. The download link already points to slides/Tutorial9_MolecularDynamics.pptx on GitHub — just commit the file to that path.

9.3. What is Molecular Dynamics?#

Molecular dynamics (MD) simulates the time evolution of atoms by integrating Newton’s equations of motion:

\[m_i \ddot{\mathbf{r}}_i = \mathbf{F}_i = -\nabla_{\mathbf{r}_i} U(\{\mathbf{r}\})\]

where \(U\) is the potential energy surface provided by the calculator.

9.3.1. Why MD for quantum optics materials?#

  • Phonon-induced dephasing: thermal motion of atoms broadens optical lines and limits coherence times of quantum emitters. MD can sample these fluctuations.

  • Configurational disorder: in amorphous hosts, the emission wavelength of a defect depends on its local environment. MD explores this disorder.

  • Zero-point motion: even at 0 K, quantum zero-point fluctuations shift emission energies (though MD is classical — quantum corrections require additional treatment).

from ase.build import bulk
from ase.calculators.emt import EMT
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase.md.verlet import VelocityVerlet
from ase import units
import numpy as np
import matplotlib.pyplot as plt

# Set up FCC aluminium for a quick MD demo
al = bulk('Al', 'fcc', a=4.05) * (3, 3, 3)  # 3x3x3 supercell
al.calc = EMT()

# Initialise velocities at 300 K
MaxwellBoltzmannDistribution(al, temperature_K=300)

print(f"System: {len(al)} Al atoms")
print(f"Initial temperature: {al.get_temperature():.1f} K")
print(f"Initial kinetic energy: {al.get_kinetic_energy():.4f} eV")

9.4. NVE Ensemble (microcanonical)#

In the NVE ensemble, the number of atoms (N), volume (V), and total energy (E) are conserved. This uses the simple Velocity Verlet integrator.

from ase.io import Trajectory

# Set up Velocity Verlet integrator
timestep = 5 * units.fs   # 5 femtoseconds
dyn = VelocityVerlet(al, timestep=timestep)

# Storage
temperatures, times, pot_energies, kin_energies = [], [], [], []

def log_properties():
    T   = al.get_temperature()
    Epot = al.get_potential_energy()
    Ekin = al.get_kinetic_energy()
    t   = dyn.get_time() / units.fs
    temperatures.append(T)
    pot_energies.append(Epot)
    kin_energies.append(Ekin)
    times.append(t)

dyn.attach(log_properties, interval=1)

# Run 200 steps
dyn.run(200)

print(f"Simulation time: {times[-1]:.0f} fs")
print(f"Mean temperature: {np.mean(temperatures):.1f} ± {np.std(temperatures):.1f} K")
print(f"Energy conservation check:")
total_E = np.array(pot_energies) + np.array(kin_energies)
print(f"  ΔE_total / E_total = {(total_E.max()-total_E.min())/abs(total_E.mean()):.2e}")

# Plot the MD trajectory properties
fig, axes = plt.subplots(2, 2, figsize=(11, 7))

axes[0,0].plot(times, temperatures, 'steelblue', lw=0.8)
axes[0,0].axhline(np.mean(temperatures), color='red', lw=1.5, ls='--', label=f'Mean = {np.mean(temperatures):.0f} K')
axes[0,0].set_xlabel('Time (fs)'); axes[0,0].set_ylabel('Temperature (K)')
axes[0,0].set_title('Temperature'); axes[0,0].legend()

axes[0,1].plot(times, pot_energies, 'darkorange', lw=0.8, label='Potential')
axes[0,1].plot(times, kin_energies, 'seagreen', lw=0.8, label='Kinetic')
axes[0,1].set_xlabel('Time (fs)'); axes[0,1].set_ylabel('Energy (eV)')
axes[0,1].set_title('Potential and kinetic energy'); axes[0,1].legend()

total_E = np.array(pot_energies) + np.array(kin_energies)
axes[1,0].plot(times, total_E, 'purple', lw=0.8)
axes[1,0].set_xlabel('Time (fs)'); axes[1,0].set_ylabel('Total energy (eV)')
axes[1,0].set_title('Total energy (should be conserved)')

# Velocity autocorrelation (proxy for phonon spectrum)
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
# Reset and run fresh trajectory storing velocities
MaxwellBoltzmannDistribution(al, temperature_K=300)
dyn2 = VelocityVerlet(al, 5*units.fs)
velocities = []
def store_v():
    velocities.append(al.get_velocities().copy())
dyn2.attach(store_v, interval=1)
dyn2.run(200)

# Compute VACF for one atom
vacf = []
v0 = velocities[0][0]
for v in velocities:
    vacf.append(np.dot(v[0], v0))
vacf = np.array(vacf) / vacf[0]
axes[1,1].plot(np.array(range(len(vacf))) * 5, vacf, 'crimson', lw=0.8)
axes[1,1].axhline(0, color='gray', lw=0.8)
axes[1,1].set_xlabel('Time lag (fs)'); axes[1,1].set_ylabel('VACF (normalised)')
axes[1,1].set_title('Velocity autocorrelation (phonon sampling)')

plt.suptitle('NVE Molecular Dynamics — FCC Al (EMT)', fontsize=13)
plt.tight_layout()
plt.show()

9.5. NVT Ensemble (Langevin thermostat)#

In the NVT ensemble we control the temperature using a thermostat. The Langevin thermostat adds random forces and friction to maintain a target temperature — analogous to a molecule in a solvent bath.

This is useful for sampling thermally activated processes relevant to defect stability.

from ase.md.langevin import Langevin
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution

al2 = bulk('Al', 'fcc', a=4.05) * (3, 3, 3)
al2.calc = EMT()
MaxwellBoltzmannDistribution(al2, temperature_K=100)

# Langevin thermostat targeting 500 K
dyn_nvt = Langevin(al2,
                   timestep=5*units.fs,
                   temperature_K=500,
                   friction=0.01/units.fs)

temps_nvt, times_nvt = [], []
def log_nvt():
    temps_nvt.append(al2.get_temperature())
    times_nvt.append(dyn_nvt.get_time() / units.fs)

dyn_nvt.attach(log_nvt, interval=1)
dyn_nvt.run(300)

plt.figure(figsize=(8, 4))
plt.plot(times_nvt, temps_nvt, 'steelblue', lw=0.7, alpha=0.8, label='Instantaneous T')
plt.axhline(500, color='red', lw=2, ls='--', label='Target T = 500 K')
plt.axhline(np.mean(temps_nvt[-100:]), color='orange', lw=2, ls=':',
            label=f'Mean T (last 100 steps) = {np.mean(temps_nvt[-100:]):.0f} K')
plt.xlabel('Time (fs)'); plt.ylabel('Temperature (K)')
plt.title('NVT MD with Langevin thermostat')
plt.legend()
plt.tight_layout()
plt.show()

9.6. Key Points#

  • MD integrates Newton’s equations of motion using the calculator for forces

  • NVE (Velocity Verlet): conserves total energy — good for testing

  • NVT (Langevin/Nosé-Hoover): controls temperature — good for sampling

  • Observers (dyn.attach(func, interval=N)) log properties every N steps

  • For quantum optics: MD samples phonon-induced fluctuations that broaden emission lines

9.7. Exercise 7.1#

Run an NVT simulation of GaN using EMT at 300 K for 500 steps with a 2 fs timestep. Plot the temperature and compare the mean to 300 K. How many steps does it take to thermalise?

9.8. Exercise 7.2 (Research-level)#

The Huang-Rhys factor \(S\) quantifies electron-phonon coupling and determines the shape of a defect’s optical emission spectrum. It is defined as \(S = \sum_k \omega_k |\Delta Q_k|^2 / (2\hbar)\), where \(\Delta Q_k\) is the mass-weighted displacement between ground and excited state geometries. Look up how to compute \(S\) for the NV centre and discuss why a small \(S\) is desirable for quantum optics applications.

Contents