So far, we have learned some fundamentals about how to write code and run an MD simulation to model the atomistic processes of materials. By running the simulation, we expect to generate a set of time-dependent atomic trajectories by solving Newton’s equations of motion for a system of particles. Next, it is important to understand these simulation results. Indeed, it is essential to extract meaningful physical properties from the simulation results. In this lecture, we will cover several fundamental post-analysis techniques:
- Validation of MD simulation by tracking the evolution of observables,
- Visualization of MD trajectories,
- Radial distribution function,
- Vibrational spectrum.
In the previous lecture and coding exercises, we have frequently mentioned the tracking of observables like total energy, volume, and temperature throughout an MD simulation to ensure the system reaches equilibrium and conserves energy (in appropriate ensembles).
- NVE: Total energy should remain constant.
- NVT: A fluctuation of temperature or kinetic energy around some values.
- NPT: Both temperature and volume fluctuate to maintain temperature and pressure.
In your MD simulation, it is advised to save these values at regular intervals. Many codes print out this information to some file. You must review these results to ensure your simulation results make sense!
import numpy as np
import matplotlib.pyplot as plt
# Load MD data (time, energy, volume, temperature)
time = np.loadtxt('time.dat') # Time data
energy = np.loadtxt('energy.dat') # Total energy
volume = np.loadtxt('volume.dat') # System volume
temperature = np.loadtxt('temperature.dat') # Temperature
# Plot evolution of observables
plt.figure(figsize=(12, 6))
plt.subplot(3, 1, 1)
plt.plot(time, energy, label='Total Energy', color='b')
plt.xlabel('Time (ps)')
plt.ylabel('Energy (eV)')
plt.title('Evolution of Total Energy')
plt.subplot(3, 1, 2)
plt.plot(time, volume, label='Volume', color='g')
plt.xlabel('Time (ps)')
plt.ylabel('Volume (ų)')
plt.title('Evolution of Volume')
plt.subplot(3, 1, 3)
plt.plot(time, temperature, label='Temperature', color='r')
plt.xlabel('Time (ps)')
plt.ylabel('Temperature (K)')
plt.title('Evolution of Temperature')Visualization allows researchers to visually inspect the dynamics of the system, spot abnormalities, and better understand atomic movements. Tools like OVITO are commonly used. Please refer to the OVITO paper to find more functions.
The Radial Distribution Function
Where:
-
$V$ is the volume of the system. -
$N$ is the number of particles. -
$r_{ij}$ is the distance between particles$i$ and$j$ . -
$\delta(r - r_{ij})$ is the Dirac delta function ensuring that only pairs with separation$r_{ij}$ equal to$r$ contribute. -
$4\pi r^2 \Delta r$ is the volume of a spherical shell at distance$r$ with thickness$\Delta r$ .
graph TD
A[Start] --> B[Compute the distance pairs]
B --> C[Group distances into bins]
C --> D[Update each bin's count]
D --> E[Normalization]
E --> F[Total number of pairs]
E --> G[Shell Volume 4πr²Δr]
E --> H[Particle density N/V]
G --> I[End]
def compute_rdf(positions, num_bins=100, r_max=10.0):
N = len(positions) # Number of atoms
rdf = np.zeros(num_bins)
dr = r_max / num_bins
for i in range(N):
for j in range(i+1, N):
r = np.linalg.norm(positions[i] - positions[j])
if r < r_max:
bin_index = int(r / dr)
rdf[bin_index] += 2 # Each pair counted twice
# Normalize RDF
r = np.linspace(0, r_max, num_bins)
rdf /= (4 * np.pi * r**2 * dr * N)
return r, rdf
# Example usage
positions = np.random.rand(100, 3) * 10 # Generate random positions
r, g_r = compute_rdf(positions)
# Plot RDF
plt.plot(r, g_r)
plt.xlabel('r (Angstrom)')
plt.ylabel('g(r)')
plt.title('Radial Distribution Function')
plt.show()RDF is commonly used to characterize the short-range order in liquids and gases. In a RDF, we are interested the location of Peaks and their spreads, as they indicate common interatomic distances (e.g., 1st and 2nd nearest-neighbor distance).
In addition to RDF, another commond characterization is to understand how particles in a system vibrate. In experiment, such information can be measured from Infrared (IR) or Raman Spectroscopy, and Inelastic Neutron/X-ray Scattering. An analogical measurement in MD is the Vibrational Density of States (VDOS).
In a MD simulation, the VDOS, called
Mathematically, the relationship is:
Where:
-
$C(\tau)$ is the velocity autocorrelation function. -
$\omega$ is the angular frequency. -
$\tau$ is the time lag.
The Fourier transform is performed over the time correlation
The VACF, called
-
$N$ is the total number of particles in the system. -
$\mathbf{v}_i(t)$ is the velocity of particle$i$ at time$t$ .
In practice, because simulations are finite and VACF data is computed over a limited time interval, we typically use the discrete Fourier transform (DFT) or fast Fourier transform (FFT) to compute the VDOS numerically:
from scipy.fft import fft
def compute_vacf(velocities):
vacf = np.correlate(velocities, velocities, mode='full')
return vacf[vacf.size // 2:] # Only take positive lag times
# Example: Simulate random velocities
velocities = np.random.randn(1000)
# Compute VACF
vacf = compute_vacf(velocities)
# Compute vibration spectrum via Fourier Transform
vibration_spectrum = np.abs(fft(vacf))
# Plot vibration spectrum
plt.plot(vibration_spectrum[:len(vibration_spectrum)//2])
plt.xlabel('Frequency (THz)')
plt.ylabel('Intensity')
plt.title('Vibration Spectrum')
plt.show()In a system of
In experiments, you can measure vibrations using IR, Raman, or Neutron scattering. How can we measure the vibrations from an MD trajectory? The key is to extract the vibration frequencies from each normal mode.
Let's imagine a simplest case of a single harmonic oscillator, the velocity of an atom is related to its position by:
For sinusoidal motion
This shows that the velocity oscillates at the same frequency
The power spectrum is given by the Fourier transform of the VACF:
Since these frequencies arise from the harmonic vibrations of atoms, the power spectrum
-
$D(\omega)$ counts the number of vibrational modes at frequency$\omega$ , -
$S(\omega)$ measures the velocity oscillations caused by these vibrational modes.
To express this mathematically, we use the relation:
This says that the VDOS is proportional to the power spectrum of the VACF. To get the VDOS from the power spectrum, we can normalize it so that:
where
Since
This form shows that the VDOS is directly obtained by taking the cosine transform of the VACF, which captures the oscillatory nature of the velocity correlations and connects them to the vibrational frequencies in the system.
If the system has only one normal mode, we can trivially find the frequency by computing the
To understand the concept of VACF, we can perform a simple numerical experiment to simulate two model systems based on an O2 diatomic molecule. In O2, the equilibrium distance is considered to be 1.16 Å, and the atoms are connected by a spring with
In the first model, we seek to simulate a simple harmonic vibration by creating the initial velocities along the spring. In the second model, we simulate the combination of vibration and rotation by creating the initial velocities not aligned with the spring.
Running the following MD code, you expect to simulate the motions and then compute the corresponding VACF and VDOS.
import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft
def harmonic_force(r1, r2, k, r_eq):
# Function to compute force due to harmonic potential
r12 = np.linalg.norm(r2 - r1)
force_mag = -k * (r12 - r_eq)
force = -force_mag * (r2 - r1) / r12
return force
def MD(r1, r2, v1, v2, N_steps):
velocities = np.zeros([N_steps, 6])
F12 = harmonic_force(r1, r2, k, r_eq)
for step in range(N_steps):
# Verlet update
r1 += v1 * dt + 0.5 * F12 / mass * dt ** 2
r2 += v2 * dt - 0.5 * F12 / mass * dt ** 2
F12_new = harmonic_force(r1, r2, k, r_eq)
v1 += 0.5 * (F12 + F12_new) * dt / mass
v2 -= 0.5 * (F12 + F12_new) * dt / mass
F12 = F12_new
velocities[step][:3] = v1
velocities[step][3:6] = v2
return velocities
if __name__ == "__main__":
# Parameters for the simulation
dt = 1e-15
data = [
('O$_2$ (vibration)', 1.16, 1180, 2.66e-26, 2000, False),
('O$_2$ (vibration + rotation)', 1.16, 1180, 2.66e-26, 5000, True),
#('H2', 0.74, 510, 1.67e-27),
#('N2', 1.10, 2294, 2.33e-26)
]
fig, axs = plt.subplots(2, len(data), figsize=(12, 6))
for i, (name, r, k, mass, N_steps, rotate) in enumerate(data):
print(name, r, k, mass)
r_eq = r * 1e-10 # equlibrium distance in m
# Initial positions
r1 = np.array([0.0, 0.0, 0.0])
r2 = np.array([0.0, 0.0, r_eq*1.2])
# Initialize velocities
v1 = np.zeros(3)
v2 = np.zeros(3)
if rotate:
v1[0] += 50
v2[0] -= 50
# MD simulation
velocities = MD(r1, r2, v1, v2, N_steps)
# Plot VACF
VACF = np.array([np.dot(velocities[0], velocities[t]) for t in range(N_steps)])
axs[0, i].plot(np.arange(N_steps)*dt*1e12, VACF)
axs[0, i].set_title(name)
axs[0, i].set_xlabel('Time (ps)')
axs[0, i].set_ylabel('VACF')
# Plot VDOS
VDOS = np.abs(fft(VACF))**2
freqs = np.fft.fftfreq(N_steps, dt) / 1e12
axs[1, i].plot(freqs[:N_steps//2], VDOS[:N_steps//2])
axs[1, i].set_xlabel('Frequency (THz)')
axs[1, i].set_ylabel('log-VDOS')
axs[1, i].set_xlim([0, 60])
axs[1, i].set_yscale('log')
plt.tight_layout()
plt.show()
For ideal harmonic oscillators, you should expect the VACF to maintain consistent behavior over time. Thanks to Fourier analysis, we can conveniently extract the vibrational frequencies for single and multiple harmonic oscillators. However, some numerical aspects must be noted. For instance, low-frequency modes cannot be sampled efficiently within a short VACF run time. The magnitude needs to be properly rescaled if needed.
Here is a report regarding the simulation of a more realistic system: VDOS study of LJ liquid and solids.
For realistic systems, LAMMPS allows the computation of both RDF and VACF. Therefore, one just needs to directly get them from the LAMMPS output and then plot them on your own.
- Interpretation of RDF: Peaks in RDF indicate atomic neighbor counts. Liquids and solids exhibit very different behaviors in their RDF.
- Interpretation of VDOS: Peaks in the VDOS correspond to characteristic vibrational modes of the system. Try to identify the low-frequency and high-frequency modes and link them to atomic motions from the MD trajectory.
- Practical issues in computing VDOS: Assuming there exists a very small frequency, how can you ensure that such a vibration mode can be reflected in your VACF?