Statistical mechanics
Develop a framework to simulate and visualize molecular collisions, and extract thermodynamic insights utilizing Python
Think about somebody continuously throwing a ball at your head. Relying upon the scale 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 instantaneous 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 strain on you (termed atmospheric strain). Since we’re “used to” the atmospheric strain, it doesn’t really feel like something out of the extraordinary, however if you’re ever in an airplane or a hyperbaric chamber, the place the strain deviates from the atmospheric strain, chances are you’ll discover completely different sensations, like ear-popping, which might be your physique’s responses to adjustments in strain.
A macroscopic understanding of strain, 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 supplied by Daniel Bernoulli to clarify strain, and in addition relate it to the temperature and kinetic power of molecules. This “kinetic concept,” nevertheless, 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 viewpoint, which was in line 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, strain, quantity, and temperature — by performing a collection of numerical simulations utilizing Python.
Kinetic concept of gases
The kinetic concept of gases is a mannequin that describes a fuel by way of molecules which might 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 power 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 amount 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’ll 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 common kinetic power 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; moderately, it’s an averaged property of a set of molecules.
Word that the kinetic concept of gases is a mannequin, it isn’t 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 super fuel very effectively, which is a really helpful macroscopic mannequin to foretell the conduct of gases below varied circumstances.
Ideally suited fuel regulation
A great fuel is a hypothetical fuel whose properties are associated to one another by a easy equation-of-state, referred to as the perfect fuel regulation.
Right here P refers back to the strain 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 effectively, and that is mirrored by R, which is the common fuel fixed. The best fuel regulation is an correct mannequin for a fuel that’s at a excessive temperature or low strain, since many of the assumptions of the kinetic concept of gases are true in these regimes.
We’ll check two relationships based mostly on the perfect fuel regulation with our mannequin.
- At a continuing temperature and variety of moles, the strain is inversely proportional to the amount.
- At a continuing quantity and variety of moles, the strain is straight 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 by experiments.
To simulate the movement of molecules in response to the kinetic concept of gases, we have to arrange a dynamic n-body simulation with collisions. Word that the target of this text is to provide 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 referred to 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 referred to as shade, which may be used to differentiate between completely different gases and set the colour of the factors within the animation.
class Molecule(object):
#Constructor
def __init__(self,molecule,place,velocity,mass,shade="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.shade=shade#Setters for place, velocity, mass and shade
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,shade):
self.shade=shade
def set_mass(self,mass):
self.mass=mass
#Getters for place, velocity, mass and shade
def get_position(self):
return self.place
def get_velocity(self):
return self.velocity
def get_color(self):
return self.shade
def get_mass(self):
return self.mass
Including molecules to a field
Subsequent, we create a Simulation class and outline key enter parameters, reminiscent of the size of the simulation field, that finally management the amount of the fuel. We additionally initialize variables for bookkeeping, reminiscent of for the wall momentum, that can enable us to calculate strain.
class Simulation(object):
#Constructor
def __init__(self,title,box_dim,t_step,particle_radius):
#Set simulation inputs
self.title=title #Title 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 record 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 have to 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 consumer. Word that it doesn’t matter what distributions the initialized positions and velocities observe (so long as there may be nothing unphysical, like a molecule exterior 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 short-term 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 certain v_mean and better certain 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 adverse with likelihood 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 short-term 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",shade="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,:],shade=shade,mass=1) for i in vary(n)]
self.molecules.prolong(_add_list)
self.n_molecules+=len(_add_list)
Lastly, we write a perform for common 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 space between each two molecules. This may turn out to be useful to detect collisions.
#Nonetheless in Simulation class
def make_matrices(self):
#Make empty matrices to retailer positions, velocities, colours, and lots more and plenty
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 aircraft, however extension to a 3D field is feasible with minor adjustments. We outline a perform that goes by a three-item guidelines:
- Verify if there are molecule-molecule collisions and if sure, replace velocities accordingly
- Replace positions of the molecules based mostly on their velocities
- Verify if there are molecule-wall collisions and if sure, replace velocities and wall momentum accordingly
To filter molecules which might be shut sufficient to collide, we scour the space matrix and return the pairs of indices of molecules which have distances lower than the sum of their radii. It’s doable 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 rate vector and the timestep to the earlier place to get the brand new place. To establish collisions with the partitions, we merely verify if the brand new place of a molecule exceeds both the decrease or higher bounds of the size of the field. If it does, we alter the signal of the rate 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: Verify molecule collisions#Discover molecule pairs that can 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()
#Verify 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 ultimate 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: Verify wall collisions
#Verify 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
#Verify 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:
#Replicate 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
#Verify collisions with the underside and left partitions
if _x_new[_x_new<=0].form[0] > 0:
#Increment collision counter
self.wall_collisions+=1
#Verify 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:
#Replicate 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 laborious half is finished! Lastly, we have to write a perform 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 info
_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 perform to deal with collisions and replace positions
self.update_positions()
#Caclulate ultimate strain
self.P=self.wall_momentum/(self.n_iters*self.t_step*np.sum(self.box_dim))
print("Common strain on wall: {0}".format(self.P))
return self.P
That concludes the code involving the physics of the simulation. Operating a simulation isn’t any enjoyable, nevertheless, for those who 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 perform to create a determine with a side ratio that’s in line with the supplied 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 perform that takes the iteration quantity (body) as enter and creates a plot. This perform is then supplied as an argument to the FuncAnimation perform within the Animation module.
#Nonetheless in Simulation class
def show_molecules(self,i):
#Clear axes
plt.cla()#Plot a line exhibiting the trajectory of a single molecule
plt.plot(self.x_dynamics[0,0,:i+1],self.x_dynamics[0,1,:i+1],shade="pink",linewidth=1.,linestyle="-")
#Plot a single molecule in pink that's being tracked
plt.scatter(self.x_dynamics[0,0,i],self.x_dynamics[0,1,i],shade="pink",s=20)
#Plot the remainder of the molecules
plt.scatter(self.x_dynamics[1:,0,i],self.x_dynamics[1:,1,i],shade=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 response to 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 perform 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’s one other animation that we are able to make, exhibiting histograms of velocities each iteration. This may enable us to look at the convergence of the rate 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,shade="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 laborious work. The trouble spent in writing the code in an object-oriented vogue will repay now, since we now have a generalized solver and might run simulations with completely different enter parameters with just some 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 strain 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 comparable to the above simulation are proven under. 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 pink, 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 assessments can be utilized to assist this quantitatively.
Extracting thermodynamic insights
We return to the perfect fuel regulation that relates varied thermodynamic variables to one another. As talked about beforehand, we check two relationships — strain in opposition to quantity and strain in opposition to temperature. We preserve the variety of molecules within the field fixed for all the following simulations. The three variables — strain, quantity, and temperature — are calculated as follows: the strain is the online momentum exchanged with the partitions throughout the complete simulation divided by the product of the whole simulation time and perimeter of the field. The amount 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 might be generalized simply to 3 dimensions).
Defining the temperature is trickier — for the reason that temperature is proportional to the common kinetic power 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 complete kinetic power has a well-defined worth that continues to be the identical throughout all simulations which have the identical temperature. Basically, when the temperatures in two simulations are identical, their preliminary complete kinetic energies are identical, which ought to be certain that the common kinetic power through the simulation can also be the identical.
The outcomes of the simulations are given in Determine 2. The typical strain 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 line with the perfect fuel regulation. In a second set of simulations, it’s noticed that the strain 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 line with the perfect fuel regulation. Due to this fact, these microscopic simulations are in a position to reproduce tendencies in thermodynamic variables which might be in line with macroscopic theories like the perfect fuel regulation.
The n-body simulation introduced 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 ample to foretell properties of a super fuel. Nevertheless, the perfect 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 embrace 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 tendencies 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 effectively.
