Simulating the movement of a group of atoms

This code is simulating the movement of a group of atoms in a material over time. The simulation has several parameters that control how the atoms behave, such as the number of atoms in the material, the temperature of the material, and the time step for the simulation.

The code first generates a lattice, or a regular arrangement, of atoms using the simulation parameters. It then calculates the forces between each pair of atoms using the Lennard-Jones potential, which is a mathematical equation that describes the interaction between atoms based on their distance from each other. The code then uses these forces to update the positions of the atoms over time, running the simulation for a specified number of time steps. Finally, the code generates a 3D scatter plot of the positions of the atoms at the final time step, which is displayed using the matplotlib library.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Dec 22 00:08:33 2022

@author: ramnot
"""

import numpy as np
import matplotlib.pyplot as plt

# Set up the simulation parameters
num_atoms = 100  # Number of atoms in the material
lattice_constant = 3.0  # Lattice constant of the material
temperature = 300  # Temperature of the material

# Generate a lattice of atoms using the simulation parameters
lattice = np.random.uniform(-lattice_constant/2, lattice_constant/2, size=(num_atoms, 3))

# Calculate the forces between the atoms using the Lennard-Jones potential
forces = np.zeros((num_atoms, 3))
for i in range(num_atoms):
    for j in range(i+1, num_atoms):
        r = lattice[i] - lattice[j]
        r_mag = np.linalg.norm(r)
        f = 24 * r * (2 * (r_mag**(-12)) - (r_mag**(-6))) / r_mag
        forces[i] += f
        forces[j] -= f

# Use the forces to update the positions of the atoms over time
time_step = 0.01  # Time step for the simulation
num_steps = 10000  # Number of time steps to run the simulation
positions = np.zeros((num_steps+1, num_atoms, 3))
positions[0] = lattice
for t in range(1, num_steps+1):
    lattice += forces * time_step / temperature  # Update the positions of the atoms
    positions[t] = lattice  # Save the positions for later

# Generate an image of the atomic structure at the final time step
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(positions[-1,:,0], positions[-1,:,1], positions[-1,:,2])
plt.show()