Statistical mechanics
Develop a framework to simulate and visualize molecular collisions, and extract thermodynamic insights utilizing Python
Think about somebody always throwing a ball at your head. Relying upon the dimensions and mass of the ball, chances are you’ll really feel something between delicate annoyance and insufferable agony. And but, there are numerous such “balls” colliding with, not solely your head, however your entire physique each immediate and you don’t really feel a factor. Nitrogen, oxygen, and water molecules within the air are in fixed random movement round you, even when the air seems to be stationary. The results of the collisions of those molecules with you is that air exerts a sure stress on you (termed atmospheric stress). Since we’re “used to” the atmospheric stress, it doesn’t really feel like something out of the peculiar, however in case you are ever in an airplane or a hyperbaric chamber, the place the stress deviates from the atmospheric stress, chances are you’ll discover totally different sensations, like ear-popping, which can be your physique’s responses to adjustments in stress.
A macroscopic understanding of stress, and the way it varies with the amount {that a} fuel occupies was given by Robert Boyle within the seventeenth century, that’s generally termed as Boyle’s regulation. Nevertheless, it was solely a century later {that a} qualitative molecular image was offered by Daniel Bernoulli to clarify stress, and in addition relate it to the temperature and kinetic vitality of molecules. This “kinetic principle,” nonetheless, didn’t discover a lot acceptance till, within the nineteenth century, James Clerk Maxwell, constructing on the work of Rudolf Clausius, formalized it by way of a statistical regulation and Ludwig Boltzmann illuminated its relationship with entropy, resulting in the formulation of the Maxwell-Boltzmann distribution of velocities of molecules in a fuel. This ultimately led to the event of an in depth theoretical remedy of gases from the molecular standpoint, which was in keeping with the well-understood macroscopic view.
On this article, we are going to attempt to hyperlink the molecular image of gases with its macroscopic properties — that’s, stress, quantity, and temperature — by performing a sequence of numerical simulations utilizing Python.
Kinetic principle of gases
The kinetic principle of gases is a mannequin that describes a fuel by way of molecules which can be in fixed, random movement.
- Molecules journey in straight traces till they’re interrupted by elastic collisions with different molecules or the partitions of the container. Elasticity implies that no kinetic vitality is misplaced in collisions.
- Molecules don’t “work together” with one another, that’s, they don’t exert any engaging or repulsive intermolecular forces.
- The quantity occupied by the molecules themselves is taken into account to be negligible in comparison with the amount occupied by the fuel.
If molecules following these guidelines are sampled and their velocities measured, they are going to conform to the Maxwell-Boltzmann distribution. The temperature of the fuel is definitely a property that we outline based mostly on the form of the distribution. The upper the typical kinetic vitality is (or the flatter the distribution is), the upper the temperature of the fuel is, and vice versa. It doesn’t make sense to outline a temperature of a person molecule; reasonably, it’s an averaged property of a set of molecules.
Observe that the kinetic principle of gases is a mannequin, it’s not essentially a real reflection of actuality (as is the case with all scientific fashions). Molecular collisions aren’t strictly elastic, and molecules do have short-range interactions. Nevertheless, a mannequin with these comparatively easy assumptions is ready to replicate the properties of a great fuel very nicely, which is a really helpful macroscopic mannequin to foretell the habits of gases underneath varied circumstances.
Excellent fuel regulation
An excellent fuel is a hypothetical fuel whose properties are associated to one another via a easy equation-of-state, known as the best fuel regulation.
Right here P refers back to the stress of the fuel, V refers back to the quantity of the container, T is the temperature, and n is the variety of moles of the fuel. If any of those three portions are held fixed, the fourth doesn’t change as nicely, and that is mirrored via R, which is the common fuel fixed. The perfect fuel regulation is an correct mannequin for a fuel that’s at a excessive temperature or low stress, since a lot of the assumptions of the kinetic principle of gases are true in these regimes.
We’ll take a look at two relationships based mostly on the best fuel regulation with our mannequin.
- At a continuing temperature and variety of moles, the stress is inversely proportional to the amount.
- At a continuing quantity and variety of moles, the stress is immediately proportional to the temperature.
The previous is known as Boyle’s regulation and the latter is known as Homosexual-Lussac’s regulation, each named after researchers that first found the relationships between these variables via experiments.
To simulate the movement of molecules in accordance with the kinetic principle of gases, we have to arrange a dynamic n-body simulation with collisions. Observe that the target of this text is to supply a simulation for pedagogical functions. Due to this fact, the code is ready as much as maximize understanding, not execution velocity.
Defining molecular properties
We begin by importing all of the required modules.
import os
import sys
import time
import numpy as np
import scipy as sci
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
First, we are going to outline a category known as Molecule. Objects of this class will retailer properties such because the plenty, positions and velocities of molecules within the simulation. We additionally outline an attribute known as coloration, that could be used to differentiate between totally different gases and set the colour of the factors within the animation.
class Molecule(object):
#Constructor
def __init__(self,molecule,place,velocity,mass,coloration="black"):
self.molecule=molecule
self.place=np.array([x_i for x_i in position])
self.velocity=np.array([v_i for v_i in velocity])
self.mass=mass
self.coloration=coloration#Setters for place, velocity, mass and coloration
def set_position(self,place):
self.place=np.array([x_i for x_i in position])
def set_velocity(self,velocity):
self.velocity=np.array([v_i for v_i in velocity])
def set_color(self,coloration):
self.coloration=coloration
def set_mass(self,mass):
self.mass=mass
#Getters for place, velocity, mass and coloration
def get_position(self):
return self.place
def get_velocity(self):
return self.velocity
def get_color(self):
return self.coloration
def get_mass(self):
return self.mass
Including molecules to a field
Subsequent, we create a Simulation class and outline key enter parameters, comparable to the scale of the simulation field, that in the end management the amount of the fuel. We additionally initialize variables for bookkeeping, comparable to for the wall momentum, that may enable us to calculate stress.
class Simulation(object):
#Constructor
def __init__(self,title,box_dim,t_step,particle_radius):
#Set simulation inputs
self.title=title #Identify of the simulation
self.box_dim=[x for x in box_dim] #Dimensions of the field
self.t_step=t_step #Timestep
self.particle_radius=particle_radius #Radius of the particles#Calculate quantity and variety of dimensions
self.V=np.prod(self.box_dim) #Space/Quantity of the field
self.dim=int(len(box_dim)) #Variety of dimensions (2D or 3D)
#Initialize paramters
self.molecules=[] #Create empty listing to retailer objects of sophistication Molecule
self.n_molecules=0 #Create variable to retailer variety of molecules
self.wall_collisions=0 #Create variable to retailer variety of wall collisions
self.wall_momentum=0 #Create variable to retailer internet momentum exchanged with wall
So as to add molecules, we first must initialize their positions and velocities. For each, we are able to outline capabilities that create arrays of values based mostly on a distribution enter by the person. Observe that it doesn’t matter what distributions the initialized positions and velocities observe (so long as there may be nothing unphysical, like a molecule outdoors the field). We will see that the velocities ultimately observe the identical distribution.
#Nonetheless in Simulation class
def _generate_initial_positions(self,n,dist="uniform"):
#Uniform distribution
if dist=="uniform":
_pos=np.random.uniform(low=[0]*self.dim,excessive=self.box_dim,measurement=(n,self.dim))#Retailer positions in non permanent variable
self._positions=_pos
def _generate_initial_velocities(self,n,v_mean,v_std,dist="regular"):
#Regular distribution with imply v_mean and std v_std
if dist=="regular":
self.v_mean=v_mean
self.v_std=v_std
_vel=np.random.regular(loc=v_mean,scale=v_std,measurement=(n,self.dim))
#Uniform distribution with decrease sure v_mean and better sure v_std
if dist=="uniform":
self.v_mean=v_mean
self.v_std=v_std
_vel=np.random.uniform(low=v_mean,excessive=v_std,measurement=(n,self.dim))
#All velocities equal to v_mean
if dist=="equal":
self.v_mean=v_mean
self.v_std=v_std
_vel=v_mean*np.ones((n,self.dim))
#Randomly change velocities to unfavorable with chance 0.5
for i in vary(_vel.form[0]):
for j in vary(_vel.form[1]):
if np.random.uniform() > 0.5:
_vel[i,j]=-_vel[i,j]
#Retailer velocities in non permanent variable
self._velocities=_vel
It’s time to add the molecules to the field. We name the 2 beforehand outlined capabilities to generate arrays of preliminary positions and velocities, create objects of the Molecule class, and assign them to the objects.
#Nonetheless in Simulation class
def add_molecules(self,molecule,n,v_mean,v_std,pos_dist="uniform",v_dist="regular",coloration="black"):
#Generate preliminary positions and velocities
self._generate_initial_positions(n,dist=pos_dist)
self._generate_initial_velocities(n,v_mean,v_std,dist=v_dist)#Initialize objects of sophistication Molecule in an inventory (set mass to 1 as default)
_add_list=[Molecule(molecule,position=self._positions[i,:],velocity=self._velocities[i,:],coloration=coloration,mass=1) for i in vary(n)]
self.molecules.lengthen(_add_list)
self.n_molecules+=len(_add_list)
Lastly, we write a operate for normal bookkeeping, that’s, to make matrices to retailer the place and velocity vectors in addition to their magnitudes. We additionally make a distance matrix, that shops the gap between each two molecules. It will come in useful to detect collisions.
#Nonetheless in Simulation class
def make_matrices(self):
#Make empty matrices to retailer positions, velocities, colours, and much
self.positions=np.zeros((self.n_molecules,self.dim))
self.velocities=np.zeros((self.n_molecules,self.dim))
self.colours=np.zeros(self.n_molecules,dtype="object")
self.plenty=np.zeros(self.n_molecules)#Iterate over molecules, get their properties and assign to matrices
for i,m in enumerate(self.molecules):
self.positions[i,:]=m.get_position()
self.velocities[i,:]=m.get_velocity()
self.colours[i]=m.get_color()
self.plenty[i]=m.get_mass()
#Make vectors with magnitudes of positions and velocities
self.positions_norm=np.linalg.norm(self.positions,axis=1)
self.velocities_norm=np.linalg.norm(self.velocities,axis=1)
#Make distance matrix
self.distance_matrix=np.zeros((self.n_molecules,self.n_molecules))
for i in vary(self.distance_matrix.form[0]):
for j in vary(self.distance_matrix.form[1]):
self.distance_matrix[i,j]=np.linalg.norm(self.positions[i,:]-self.positions[j,:])
#Set diagonal entries (distance with itself) to a excessive worth
#to stop them for showing within the subsequent distance filter
np.fill_diagonal(self.distance_matrix,1e5)
Detecting and dealing with collisions
Right here, we come to the meat of the simulation, that’s, modeling the dynamics of the molecules. Collision physics are outlined for a 2D airplane, however extension to a 3D field is feasible with minor adjustments. We outline a operate that goes via a three-item guidelines:
- Examine if there are molecule-molecule collisions and if sure, replace velocities accordingly
- Replace positions of the molecules based mostly on their velocities
- Examine if there are molecule-wall collisions and if sure, replace velocities and wall momentum accordingly
To filter molecules which can be shut sufficient to collide, we scour the gap matrix and return the pairs of indices of molecules which have distances lower than the sum of their radii. It’s potential for molecules to be inside this cutoff, and but be departing away from one another. So, we apply one other criterion (see Determine 1) to analyze whether or not molecules are approaching one another. If they’re, we replace their velocities based mostly on the equations given in Determine 1. We arrive at these equations by conserving the kinetic energies and the linear momenta (alongside the collision axis) of the molecules.
Updating the positions of the molecules is sort of easy. We add the product of the speed vector and the timestep to the earlier place to get the brand new place. To establish collisions with the partitions, we merely test if the brand new place of a molecule exceeds both the decrease or higher bounds of the scale of the field. If it does, we modify the signal of the speed regular to the wall, and set the brand new place based mostly on this velocity. Additional, we add twice the magnitude of this velocity to the variable monitoring the momentum exchanged with the wall.
#Nonetheless in Simulation class
def update_positions(self):
#1: Examine molecule collisions#Discover molecule pairs that may collide
collision_pairs=np.argwhere(self.distance_matrix < 2*self.particle_radius)
#If collision pairs exist
if len(collision_pairs):
#Undergo pairs and take away repeats of indices
#(for eg., solely take into account (1,2), take away (2,1))
pair_list=[]
for pair in collision_pairs:
add_pair=True
for p in pair_list:
if set(p)==set(pair):
add_pair=False
break
if add_pair:
pair_list.append(pair)
#For each remaining pair, get the molecules, positions, and velocities
for pair in pair_list:
m_1=self.molecules[pair[0]]
m_2=self.molecules[pair[1]]
x_1=m_1.get_position()
x_2=m_2.get_position()
u_1=m_1.get_velocity()
u_2=m_2.get_velocity()
#Examine if molecules are approaching or departing
approach_sign=np.signal(np.dot(u_1-u_2,x_2-x_1))
#If molecules are approaching
if approach_sign == 1:
#Get plenty
ms_1=m_1.get_mass()
ms_2=m_2.get_mass()
#Calculate closing velocities
v_1=u_1 - 2*ms_2/(ms_1 + ms_2) * (np.dot(u_1-u_2,x_1-x_2)/np.linalg.norm(x_1-x_2)**2) * (x_1 - x_2)
v_2=u_2 - 2*ms_1/(ms_1 + ms_2) * (np.dot(u_2-u_1,x_2-x_1)/np.linalg.norm(x_2-x_1)**2) * (x_2 - x_1)
#Replace velocities of the molecule objects
m_1.set_velocity(v_1)
m_2.set_velocity(v_2)
#2: Replace positions
#Iterate over all of the molecule objects
for i,m in enumerate(self.molecules):
#Get the place, velocity, and mass
_x=m.get_position()
_v=m.get_velocity()
_m=m.get_mass()
#Calculate new place
_x_new=_x + _v * self.t_step
#3: Examine wall collisions
#Examine collisions with the highest and proper partitions
_wall_diff=_x_new - np.array(self.box_dim)
#If wall collisions current
if _wall_diff[_wall_diff>=0].form[0] > 0:
#Increment collision counter
self.wall_collisions+=1
#Examine whether or not collision in x or y route
_coll_ind=np.argwhere(_wall_diff>0)
#For part(s) to be mirrored
for c in _coll_ind:
#Mirror velocity
_v[c]=-_v[c]
#Increment wall momentum
self.wall_momentum+=2*_m*np.abs(_v[c])
#Replace velocity
m.set_velocity(_v)
#Replace place based mostly on new velocity
_x_new=_x + _v * self.t_step
#Examine collisions with the underside and left partitions
if _x_new[_x_new<=0].form[0] > 0:
#Increment collision counter
self.wall_collisions+=1
#Examine whether or not collision in x or y route
_coll_ind=np.argwhere(_x_new<0)
#For part(s) to be mirrored
for c in _coll_ind:
#Mirror velocity
_v[c]=-_v[c]
#Increment wall momentum
self.wall_momentum+=2*_m*np.abs(_v[c])
#Replace velocity
m.set_velocity(_v)
#Replace place based mostly on new velocity
_x_new=_x + _v * self.t_step
#Replace place of the molecule object
m.set_position(_x_new)
#Assemble matrices with up to date positions and velocities
self.make_matrices()
That’s it, the arduous half is completed! Lastly, we have to write a operate to run the simulation. This entails calling the capabilities that we have now outlined beforehand in a loop that runs for a specified variety of iterations, calculated based mostly on the required simulation time and time step.
#Nonetheless in Simulation class
def safe_division(self,n,d):
if d==0:
return 0
else:
return n/ddef run_simulation(self,max_time):
#Print "Beginning simulation"
print("Beginning simulation...")
#Make matrices
self.make_matrices()
#Calculate variety of iterations
self.max_time=max_time
self.n_iters=int(np.flooring(self.max_time/self.t_step))
#Make tensors to retailer positions and velocities of all molecules at every timestep f
self.x_dynamics=np.zeros(((self.n_molecules,self.dim,self.n_iters)))
self.v_dynamics=np.zeros((self.n_molecules,self.n_iters))
#In every iteration
for i in vary(self.n_iters):
#Save positions and velocities to the outlined tensors
self.x_dynamics[:,:,i]=self.positions
self.v_dynamics[:,i]=self.velocities_norm
#Calculate rms velocity
self.v_rms=np.sum(np.sqrt(self.velocities_norm**2))/self.velocities_norm.form[0]
#Print present iteration data
_P=self.safe_division(self.wall_momentum,i*self.t_step*np.sum(self.box_dim))
print("Iteration:{0:d}tTime:{1:.2E}tV_RMS:{2:.2E}tWall Strain:{3}".format(i,i*self.t_step,self.v_rms,_P))
#Name the update_positions operate to deal with collisions and replace positions
self.update_positions()
#Caclulate closing stress
self.P=self.wall_momentum/(self.n_iters*self.t_step*np.sum(self.box_dim))
print("Common stress on wall: {0}".format(self.P))
return self.P
That concludes the code involving the physics of the simulation. Working a simulation isn’t any enjoyable, nonetheless, when you can not visualize it. Let’s make the most of matplotlib to create an animation of the simulation.
Animating the simulation field
We begin by writing a operate to create a determine with a side ratio that’s in keeping with the offered dimensions of the field.
#Nonetheless in Simulation class
def create_2D_box(self):
fig=plt.determine(figsize=(10,10*self.box_dim[1]/self.box_dim[0]),dpi=300)
return fig
We’ll use the Animation module in matplotlib to make the animation. To make the most of that, we have to outline a operate that takes the iteration quantity (body) as enter and creates a plot. This operate is then offered as an argument to the FuncAnimation operate within the Animation module.
#Nonetheless in Simulation class
def show_molecules(self,i):
#Clear axes
plt.cla()#Plot a line displaying the trajectory of a single molecule
plt.plot(self.x_dynamics[0,0,:i+1],self.x_dynamics[0,1,:i+1],coloration="crimson",linewidth=1.,linestyle="-")
#Plot a single molecule in crimson that's being tracked
plt.scatter(self.x_dynamics[0,0,i],self.x_dynamics[0,1,i],coloration="crimson",s=20)
#Plot the remainder of the molecules
plt.scatter(self.x_dynamics[1:,0,i],self.x_dynamics[1:,1,i],coloration=self.colours[1:],s=20)
#Take away ticks on the plot
plt.xticks([])
plt.yticks([])
#Set margins to 0
plt.margins(0)
#Set the boundaries of the field in accordance with the field dimensions
plt.xlim([0,self.box_dim[0]])
plt.ylim([0,self.box_dim[1]])
def make_animation(self,filename="KTG_animation.mp4"):
#Name the operate to create the determine
fig=self.create_2D_box()
#Create the animation
anim=FuncAnimation(fig,self.show_molecules,frames=self.n_iters,interval=50,blit=False)
#Save animation as a file
anim.save(filename,author="ffmpeg")
There may be one other animation that we are able to make, displaying histograms of velocities each iteration. It will enable us to look at the convergence of the speed distribution to a Maxwell-Boltzmann distribution.
#Nonetheless in Simulation class
def plot_hist(self,i):
#Clear axes
plt.cla()#Make histogram
plt.hist(self.v_dynamics[:,i],density=True,coloration="plum",edgecolor="black")
#Outline axis limits
plt.xlim([0,3])
plt.ylim([0,3])
def make_velocity_histogram_animation(self,filename="KTG_histogram.mp4"):
#Create empty determine
fig=plt.determine(figsize=(5,5),dpi=500)
#Create animation
anim_hist=FuncAnimation(fig,self.plot_hist,frames=self.n_iters,interval=50,blit=False)
#Save animation
anim_hist.save(filename,author="ffmpeg")
It’s time to reap the rewards of our arduous work. The trouble spent in writing the code in an object-oriented vogue will repay now, since we now have a generalized solver and may run simulations with totally different enter parameters with only a few traces of code.
if __name__=="__main__":
#Create simulation object and outline enter parameters
sim=Simulation(title="kinetic_theory_simulation",box_dim=[1.0,1.0],t_step=1e-2,particle_radius=1e-2)#Add N2 molecules to the field
sim.add_molecules("N2",n=100,v_mean=1.0,v_std=0.2,v_dist="regular")
#Run the simulation and retailer the stress output in P
P=sim.run_simulation(15)
#Make the field animation
sim.make_animation()
#Make the histogram animation
sim.make_velocity_histogram_animation()
Animations of the simulation field and velocity histogram similar to the above simulation are proven beneath. Within the first graphic, the actions of the molecules within the field, together with collisions, are proven, and the trajectory of 1 chosen molecule is highlighted in crimson, for illustrative functions. Within the second graphic, the histogram of velocities of all of the molecules within the field is proven at every iteration, and it’s clear from the animation that the preliminary distribution, which is Gaussian (as specified), adjustments to a distribution that has a narrower left tail and broader proper tail, mimicking the traits of a Maxwell-Boltzmann distribution. Extra rigorous statistical checks can be utilized to assist this quantitatively.
Extracting thermodynamic insights
We return to the best fuel regulation that relates varied thermodynamic variables to one another. As talked about beforehand, we take a look at two relationships — stress towards quantity and stress towards temperature. We preserve the variety of molecules within the field fixed for all the next simulations. The three variables — stress, quantity, and temperature — are calculated as follows: the stress is the online momentum exchanged with the partitions throughout your complete simulation divided by the product of the overall simulation time and perimeter of the field. The quantity is outlined because the product of the size and breadth of the field (technically it’s the space, since we’re working in two dimensions, however the insights could be generalized simply to 3 dimensions).
Defining the temperature is trickier — for the reason that temperature is proportional to the typical kinetic vitality of the molecules within the field, we take into account the sq. of the imply velocity of the preliminary distribution to be a proxy for temperature. To take away any stochasticity from this estimate, the preliminary velocities assigned to the molecules in these simulations are set to a single specified worth. For example, if the required worth is 1 m/s, the preliminary velocities of all molecules are both +1 m/s or -1 m/s. This ensures that the preliminary whole kinetic vitality has a well-defined worth that continues to be the identical throughout all simulations which have the identical temperature. Primarily, when the temperatures in two simulations are identical, their preliminary whole kinetic energies are identical, which ought to make sure that the typical kinetic vitality in the course of the simulation can be the identical.
The outcomes of the simulations are given in Determine 2. The typical stress on the partitions of the field will increase linearly with a rise within the inverse of the amount at a continuing temperature (see Determine 2a). The slope of every isotherm is proportional to the temperature, in keeping with the best fuel regulation. In a second set of simulations, it’s noticed that the stress will increase linearly with a rise within the temperature at a continuing quantity (see Determine 2b). On this case, the slope of every isochore is inversely proportional to the amount, additionally in keeping with the best fuel regulation. Due to this fact, these microscopic simulations are in a position to reproduce developments in thermodynamic variables which can be in keeping with macroscopic theories like the best fuel regulation.
The n-body simulation offered on this article is an instance of a easy molecular dynamics simulation with none interactive forces. That is, after all, a gross oversimplification of how molecules behave and work together with one another however, as we have now seen, it’s adequate to foretell properties of a great fuel. Nevertheless, the best fuel assumption is never used for calculating properties of gases for engineering functions, just like the growth of steam in a steam turbine. Extra complicated equation-of-state fashions that embody interactions between molecules are required to precisely mannequin such processes. Including short-range interactions between molecules on this code can result in higher replica of developments predicted by such fashions for actual gases. Additional, the utilization of potentials like Lennard-Jones and addition of thermostats can enable prediction of properties for liquids as nicely.
