## 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.