import numpy as np
import matplotlib.pyplot as plt
[docs]class Scattering(object): # pragma: no cover
"""The Scattering class will plot the particle positions, radial
distribution function, mean squared deviation and scattering profile (as a
fft of the rdf). This sampling class is ideal for observing the phase
transitions between solid, liquid, gas.
Parameters
----------
system: System
The whole system information.
size: string
The size of the visualisation:
- 'small'
- 'medium' (default)
- 'large'
"""
def __init__(self, system, size='medium'):
fig, ax = environment(4, size)
self.average_rdf = []
self.r = []
self.average_diff = []
self.q = []
self.initial_pos = [
system.particles["xposition"],
system.particles["yposition"],
]
setup_cellview(ax[0, 0], system)
setup_rdfview(ax[0, 1], system)
setup_diffview(ax[1, 1])
setup_msdview(ax[1, 0])
plt.tight_layout()
self.ax = ax
self.fig = fig
[docs] def update(self, system):
"""This updates the visualisation environment. Often this can be slower
than the cythonised force calculation so used is wisely.
Parameters
----------
system: System
The whole system information.
"""
update_cellview(self.ax[0, 0], system)
update_rdfview(self.ax[0, 1], system, self.average_rdf, self.r)
update_diffview(self.ax[1, 1], system, self.average_diff, self.q)
update_msdview(self.ax[1, 0], system)
self.fig.canvas.draw()
[docs] def average(self):
gr = np.average(self.average_rdf, axis=0)
x = np.average(self.r, axis=0)
line = self.ax[0, 1].lines[0]
line.set_xdata(x)
line.set_ydata(gr)
self.ax[0, 1].set_ylim([0, np.amax(gr) + np.amax(gr) * 0.05])
iq = np.average(self.average_diff, axis=0)
x = np.average(self.q, axis=0)
line = self.ax[1, 1].lines[0]
line.set_ydata(iq)
line.set_xdata(x)
self.ax[1, 1].set_ylim([0, np.amax(iq) + np.amax(iq) * 0.05])
self.ax[1, 1].set_xlim([np.amin(x), np.amax(x)])
self.fig.canvas.draw()
[docs]class Phase(object): # pragma: no cover
"""The Phase class will plot the particle positions, radial
distribution function, mean squared deviation and total energy of the
simulation. This sampling class is ideal for observing the phase
transitions between solid, liquid, gas.
Parameters
----------
system: System
The whole system information.
size: string
The size of the visualisation:
- 'small'
- 'medium' (default)
- 'large'
"""
def __init__(self, system, size='medium'):
fig, ax = environment(4, size)
self.average_rdf = []
self.r = []
self.average_diff = []
self.q = []
self.initial_pos = [
system.particles["xposition"],
system.particles["yposition"],
]
setup_cellview(ax[0, 0], system)
setup_rdfview(ax[1, 1], system)
setup_energyview(ax[0, 1])
setup_msdview(ax[1, 0])
plt.tight_layout()
self.ax = ax
self.fig = fig
[docs] def update(self, system):
"""This updates the visualisation environment. Often this can be slower
than the cythonised force calculation so used is wisely.
Parameters
----------
system: System
The whole system information.
"""
update_cellview(self.ax[0, 0], system)
update_rdfview(self.ax[1, 1], system, self.average_rdf, self.r)
update_energyview(self.ax[0, 1], system)
update_msdview(self.ax[1, 0], system)
self.fig.canvas.draw()
[docs] def average(self):
gr = np.average(self.average_rdf, axis=0)
x = np.average(self.r, axis=0)
line = self.ax[1, 1].lines[0]
line.set_xdata(x)
line.set_ydata(gr)
self.ax[1, 1].set_ylim([0, np.amax(gr) + 0.01 * np.max(gr)])
self.fig.canvas.draw()
[docs]class Interactions(object): # pragma: no cover
"""The Interactions class will plot the particle positions, total force,
simulation pressure and temperature. This class is perfect for showing the
interactions between the particles and therefore the behaviour of ideal
gases and deviation when the conditions of an ideal gas are not met.
Parameters
----------
system: System
The whole system information.
size: string
The size of the visualisation:
- 'small'
- 'medium' (default)
- 'large'
"""
def __init__(self, system, size='medium'):
fig, ax = environment(4, size)
setup_cellview(ax[0, 0], system)
setup_forceview(ax[1, 1])
setup_pressureview(ax[1, 0])
setup_tempview(ax[0, 1])
plt.tight_layout()
self.ax = ax
self.fig = fig
[docs] def update(self, system):
"""This updates the visualisation environment. Often this can be slower
than the cythonised force calculation so used is wisely.
Parameters
----------
system: System
The whole system information.
"""
update_cellview(self.ax[0, 0], system)
update_forceview(self.ax[1, 1], system)
update_tempview(self.ax[0, 1], system)
update_pressureview(self.ax[1, 0], system)
self.fig.canvas.draw()
[docs]class JustCell(object): # pragma: no cover
"""The JustCell class will plot just the particles positions. This is a
simplistic sampling class for quick visualisation.
Parameters
----------
system: System
The whole system information.
size: string
The size of the visualisation:
- 'small'
- 'medium' (default)
- 'large'
scale: float (optional)
The amount by which to scale down the size of the particles
"""
def __init__(self, system, size='medium', scale=1):
fig, ax = environment(1, size)
setup_cellview(ax, system, scale)
plt.tight_layout()
self.ax = ax
self.fig = fig
[docs] def update(self, system):
"""This updates the visualisation environment. Often this can be slower
than the cythonised force calculation so use this wisely.
Parameters
----------
system: System
The whole system information.
"""
update_cellview(self.ax, system)
self.fig.canvas.draw()
[docs]class CellPlus(object): # pragma: no cover
"""The CellPlus class will plot the particles positions in addition to one
user defined one-dimensional dataset. This is designed to allow user
interaction with the plotting data.
Parameters
----------
system: System
The whole system information.
xlabel: string
The label for the x-axis of the custom plot.
ylabel: string
The label for the y-axis of the custom plot.
size: string
The size of the visualisation:
- 'small'
- 'medium' (default)
- 'large'
"""
def __init__(self, system, xlabel, ylabel, size='medium'):
fig, ax = environment(2, size)
setup_cellview(ax[0], system)
ax[1].plot([0], color="#34a5daff")
ax[1].set_ylabel(ylabel, fontsize=16)
ax[1].set_xlabel(xlabel, fontsize=16)
plt.tight_layout()
self.ax = ax
self.fig = fig
[docs] def update(self, system, xdata, ydata):
"""This updates the visualisation environment. Often this can be slower
than the cythonised force calculation so use this wisely.
Parameters
----------
system: System
The whole system information.
xdata: float, array-like
x-data that should be plotted in the custom plot.
ydata: float, array-like
y-data that should be plotted in the custom plot.
"""
update_cellview(self.ax[0], system)
line1 = self.ax[1].lines[0]
line1.set_xdata(xdata)
line1.set_ydata(ydata)
self.ax[1].set_ylim([np.amin(ydata), np.amax(ydata)])
self.ax[1].set_xlim(np.amin(xdata), np.amax(xdata))
self.fig.canvas.draw()
[docs]class Energy(object): # pragma: no cover
"""The energy class will plot the particle positions and potential energy
of the system.
Parameters
----------
system: System
The whole system information.
size: string
The size of the visualisation:
- 'small'
- 'medium' (default)
- 'large'
"""
def __init__(self, system, size='medium'):
fig, ax = environment(2, size)
setup_cellview(ax[0], system)
setup_energyview(ax[1])
plt.tight_layout()
self.ax = ax
self.fig = fig
[docs] def update(self, system):
"""This updates the visualisation environment. Often this can be slower
than the cythonised force calculation so used is wisely.
Parameters
----------
system: System
The whole system information.
"""
update_cellview(self.ax[0], system)
update_energyview(self.ax[1], system)
self.fig.canvas.draw()
[docs]class MaxBolt(object): # pragma: no cover
"""The MaxBolt class will plot the particle positions and a histogram
of the particle velocities.
Parameters
----------
system: System
The whole system information.
size: string
The size of the visualisation:
- 'small'
- 'medium' (default)
- 'large'
"""
def __init__(self, system, size='medium'):
fig, ax = environment(2, size)
self.velocities = np.array([])
setup_cellview(ax[0], system)
setup_maxbolthist(ax[1])
plt.tight_layout()
self.ax = ax
self.fig = fig
[docs] def update(self, system):
"""This updates the visualisation environment. Often this can be slower
than the cythonised force calculation so used is wisely.
Parameters
----------
system: System
The whole system information.
"""
update_cellview(self.ax[0], system)
self.velocities = update_maxbolthist(self.ax[1], system, self.velocities)
self.fig.canvas.draw()
[docs]class RDF(object): # pragma: no cover
"""The RDF class will plot the particle positions and radial distribution
function. This sampling class is can be used to show the relative RDFs for
solid, liquid, gas.
Parameters
----------
system: System
The whole system information.
size: string
The size of the visualisation:
- 'small'
- 'medium' (default)
- 'large'
"""
def __init__(self, system, size='medium'):
fig, ax = environment(2, size)
self.average_rdf = []
self.r = []
self.average_diff = []
self.q = []
self.initial_pos = [
system.particles["xposition"],
system.particles["yposition"],
]
setup_cellview(ax[0], system)
setup_rdfview(ax[1], system)
plt.tight_layout()
self.ax = ax
self.fig = fig
[docs] def update(self, system):
"""This updates the visualisation environment. Often this can be
slower than the cythonised force calculation so used is wisely.
Parameters
----------
system: System
The whole system information.
"""
update_cellview(self.ax[0], system)
update_rdfview(self.ax[1], system, self.average_rdf, self.r)
self.fig.canvas.draw()
[docs] def average(self):
gr = np.average(self.average_rdf, axis=0)
x = np.average(self.r, axis=0)
line = self.ax[1].lines[0]
line.set_xdata(x)
line.set_ydata(gr)
self.ax[1].set_ylim([0, np.amax(gr) + np.amax(gr) * 0.05])
self.fig.canvas.draw()
[docs]def environment(panes, size='medium'): # pragma: no cover
"""The visualisation environment consists of a series of panes (1, 2, or 4
are allowed). This function allows the number of panes in the
visualisation to be defined.
Parameters
----------
panes: int
Number of visualisation panes.
size: string
The size of the visualisation:
- 'small'
- 'medium' (default)
- 'large'
Returns
-------
Matplotlib.figure.Figure object:
The relevant Matplotlib figure.
Axes object or array of axes objects:
The axes related to each of the panes. For panes=1 this is a single
object, for panes=2 it is a 1-D array and
for panes=4 it is a 2-D array.
"""
scale = 1
if size == 'small':
scale = 2
elif size == 'large':
scale = 0.5
if panes == 1:
fig, ax = plt.subplots(figsize=(4/scale, 4/scale))
elif panes == 2:
fig, ax = plt.subplots(1, 2, figsize=(8/scale, 4/scale))
elif panes == 4:
fig, ax = plt.subplots(2, 2, figsize=(8/scale, 8/scale))
else:
AttributeError("The only options for the number of panes are 1, 2, or " "4")
return fig, ax
[docs]def setup_cellview(ax, system, scale=1): # pragma: no cover
"""Builds the particle position visualisation pane.
Parameters
----------
ax: Axes object
The axes position that the pane should be placed in.
system: System
The whole system information.
scale: float (optional)
The amount by which the particle size should be scaled down.
"""
xpos = system.particles["xposition"]
ypos = system.particles["yposition"]
# These numbers are chosen to scale the size of the particles nicely to
# allow the particles to appear to interact appropriately
mk = 6.00555e-8 / (system.box_length - 2.2727e-10) - 1e-10 / scale
ax.plot(xpos, ypos, "o", markersize=mk, markeredgecolor="black", color="#34a5daff")
ax.set_xlim([0, system.box_length])
ax.set_ylim([0, system.box_length])
ax.set_xticks([])
ax.set_yticks([])
[docs]def setup_forceview(ax): # pragma: no cover
"""Builds the total force visualisation pane.
Parameters
----------
ax: Axes object
The axes position that the pane should be placed in.
"""
ax.plot([0], color="#34a5daff")
ax.set_ylabel("Force/N", fontsize=16)
ax.set_xlabel("Time/s", fontsize=16)
[docs]def setup_energyview(ax): # pragma: no cover
"""Builds the total force visualisation pane.
Parameters
----------
ax: Axes object
The axes position that the pane should be placed in.
"""
ax.plot([0], color="#34a5daff")
ax.set_ylabel("Energy/J", fontsize=16)
ax.set_xlabel("Step", fontsize=16)
[docs]def setup_rdfview(ax, system): # pragma: no cover
"""Builds the radial distribution function visualisation pane.
Parameters
----------
ax: Axes object
The axes position that the pane should be placed in.
system: System
The whole system information.
"""
ax.plot([0], color="#34a5daff")
ax.set_xlim([0, system.box_length / 2])
ax.set_yticks([])
ax.set_ylabel("RDF", fontsize=16)
ax.set_xlabel("r/m", fontsize=16)
[docs]def setup_diffview(ax): # pragma: no cover
"""Builds the scattering profile visualisation pane.
Parameters
----------
ax: Axes object
The axes position that the pane should be placed in.
"""
ax.plot([0], color="#34a5daff")
ax.set_yticks([])
ax.set_ylabel("I(q)", fontsize=16)
ax.set_xlabel("q/m$^{-1}$", fontsize=16)
[docs]def setup_pressureview(ax): # pragma: no cover
"""Builds the simulation instantaneous pressure visualisation pane.
Parameters
----------
ax: Axes object
The axes position that the pane should be placed in.
"""
ax.plot([0], color="#34a5daff")
ax.set_ylabel(r"Pressure/$\times10^6$Pa m$^{-1}$", fontsize=16)
ax.set_xlabel("Time/s", fontsize=16)
[docs]def setup_tempview(ax): # pragma: no cover
"""Builds the simulation instantaneous temperature visualisation pane.
Parameters
----------
ax: Axes object
The axes position that the pane should be placed in.
"""
ax.plot([0], color="#34a5daff")
ax.set_ylabel("Temperature/K", fontsize=16)
ax.set_xlabel("Time/s", fontsize=16)
[docs]def setup_maxbolthist(ax): # pragma: no cover
"""Builds the simulation velocity histogram
visualisation pane.
Parameters
----------
ax: Axes object
The axes position that the pane should be placed in.
"""
ax.step([0], [0], color="#34a5daff")
ax.set_ylabel("PDF", fontsize=16)
ax.set_xlabel("Velocity/ms$^{-1}$", fontsize=16)
[docs]def update_cellview(ax, system): # pragma: no cover
"""Updates the particle positions visualisation pane.
Parameters
----------
ax: Axes object
The axes position that the pane should be placed in.
system: System
The whole system information.
"""
x3 = system.particles["xposition"]
y3 = system.particles["yposition"]
line = ax.lines[0]
line.set_ydata(y3)
line.set_xdata(x3)
[docs]def update_rdfview(ax, system, average_rdf, r): # pragma: no cover
"""Updates the radial distribution function visualisation pane.
Parameters
----------
ax: Axes object
The axes position that the pane should be placed in.
system: System
The whole system information.
average_rdf: array_like
The radial distribution functions g(r) for each timestep, to later be
averaged.
r: array_like
The radial distribution functions r for each timestep, to later be
averaged.
"""
hist, bin_edges = np.histogram(
system.distances, bins=np.linspace(0, system.box_length / 2 + 0.5e-10, 100)
)
gr = hist / (
system.number_of_particles
* (system.number_of_particles / system.box_length ** 2)
* np.pi
* (bin_edges[:-1] + 0.5e-10 / 2.0)
* 0.5
)
average_rdf.append(gr)
x = bin_edges[:-1] + 0.5e-10 / 2
r.append(x)
line = ax.lines[0]
line.set_xdata(x)
line.set_ydata(gr)
ax.set_ylim([0, np.amax(gr) + np.amax(gr) * 0.01])
[docs]def update_diffview(ax, system, average_diff, q): # pragma: no cover
"""Updates the scattering profile visualisation pane.
Parameters
----------
ax: Axes object
The axes position that the pane should be placed in.
system: System
The whole system information.
average_diff: array_like
The scattering profile's i(q) for each timestep, to later be averaged.
q: array_like
The scattering profile's q for each timestep, to later be averaged.
"""
# the range of q is chosen to give a representive range for the
# interactions minimum is the reciprocal of the box length and the maximum
# is the reciprocal of the van der Waals diameter of the argon atom
qw = np.linspace(2 * np.pi / system.box_length, 10e10, 1000)[20:]
i = np.zeros_like(qw)
for j in range(0, len(qw)):
i[j] = np.sum(
3.644 * (np.sin(qw[j] * system.distances)) / (qw[j] * system.distances)
)
if i[j] < 0:
i[j] = 0
x2 = qw
y2 = i
average_diff.append(y2)
q.append(x2)
line1 = ax.lines[0]
line1.set_xdata(x2)
line1.set_ydata(y2)
ax.set_ylim([0, np.amax(y2) + np.amax(y2) * 0.05])
ax.set_xlim(np.amin(x2), np.amax(x2))
[docs]def update_forceview(ax, system): # pragma: no cover
"""Updates the total force visualisation pane.
Parameters
----------
ax: Axes object
The axes position that the pane should be placed in.
system: System
The whole system information.
"""
line = ax.lines[0]
line.set_ydata(system.force_sample)
line.set_xdata(np.arange(0, system.step) * system.timestep_length)
ax.set_xlim(0, system.step * system.timestep_length)
ax.set_ylim(
np.amin(system.force_sample) - np.amax(system.force_sample) * 0.05,
np.amax(system.force_sample) + np.amax(system.force_sample) * 0.05,
)
[docs]def update_energyview(ax, system): # pragma: no cover
"""Updates the total force visualisation pane.
Parameters
----------
ax: Axes object
The axes position that the pane should be placed in.
system: System
The whole system information.
"""
line = ax.lines[0]
if system.force_sample != []:
y = system.energy_sample + 1.3806e-23 * system.temperature_sample
line.set_xdata(np.arange(0, system.step) * system.timestep_length)
ax.set_xlim(0, system.step * system.timestep_length)
ax.set_xlabel("Time/s", fontsize=16)
else:
y = system.energy_sample
line.set_xdata(np.arange(0, system.step + 1))
ax.set_xlim(0, system.step)
line.set_ydata(y)
ax.set_ylim(
np.amin(y) - np.abs(np.amax(y)) * 0.05, np.amax(y) + np.abs(np.amax(y)) * 0.05
)
[docs]def update_tempview(ax, system): # pragma: no cover
"""Updates the simulation instantaneous temperature visualisation pane.
Parameters
----------
ax: Axes object
The axes position that the pane should be placed in.
system: System
The whole system information.
"""
line = ax.lines[0]
line.set_ydata(system.temperature_sample)
line.set_xdata(np.arange(0, system.step) * system.timestep_length)
ax.set_xlim(0, system.step * system.timestep_length)
ax.set_ylim(
np.amin(system.temperature_sample) - np.amax(system.temperature_sample) * 0.05,
np.amax(system.temperature_sample) + np.amax(system.temperature_sample) * 0.05,
)
[docs]def update_maxbolthist(ax, system, velocities): # pragma: no cover
"""Updates the simulation velocity histogram visualisation pane.
Parameters
----------
ax: Axes object
The axes position that the pane should be placed in.
system: System
The whole system information.
"""
velocities = np.append(velocities, np.sqrt(np.square(system.particles['xvelocity']) + np.square(system.particles['yvelocity'])))
line = ax.lines[0]
hist, bin_edges = np.histogram(velocities, bins=25, density=True)
line.set_ydata(hist)
line.set_xdata(bin_edges[:-1])
ax.set_xlim(0, bin_edges[-2])
ax.set_ylim(
0,
np.amax(hist) + np.amax(hist) * 0.05,
)
return velocities
[docs]def update_pressureview(ax, system): # pragma: no cover
"""Updates the simulation instantaneous pressure visualisation pane.
Parameters
----------
ax: Axes object
The axes position that the pane should be placed in.
system: System
The whole system information.
"""
line = ax.lines[0]
# Scaling the pressure to give the numbers with a nice number of
# significant figures
data = system.pressure_sample * 1e6
line.set_ydata(data)
line.set_xdata(np.arange(0, system.step) * system.timestep_length)
ax.set_xlim(0, system.step * system.timestep_length)
ax.set_ylim(
np.amin(data) - np.amax(data) * 0.05, np.amax(data) + np.amax(data) * 0.05
)
[docs]def setup_msdview(ax): # pragma: no cover
"""Builds the simulation mean squared deviation visualisation pane.
Parameters
----------
ax: Axes object
The axes position that the pane should be placed in.
"""
ax.plot([0], color="#34a5daff")
ax.set_ylabel("MSD/m$^2$", fontsize=16)
ax.set_xlabel("Time/s", fontsize=16)
[docs]def update_msdview(ax, system): # pragma: no cover
"""Updates the simulation mean squared deviation visualisation pane.
Parameters
----------
ax: Axes object
The axes position that the pane should be placed in.
system: System
The whole system information.
"""
line = ax.lines[0]
line.set_ydata(system.msd_sample)
line.set_xdata(np.arange(0, system.step) * system.timestep_length)
ax.set_xlim(0, system.step * system.timestep_length)
ax.set_ylim(0, np.amax(system.msd_sample) + np.amax(system.msd_sample) * 0.05)