+from pylj import md, comp, sample
+
+[docs]def md_nvt(number_of_particles, temperature, box_length, number_of_steps, sample_frequency):
+ # Creates the visualisation environment
+ #%matplotlib notebook
+ # Initialise the system
+ system = md.initialise(number_of_particles, temperature, box_length, 'square')
+ # This sets the sampling class
+ sample_system = sample.Interactions(system)
+ # Start at time 0
+ system.time = 0
+ # Begin the molecular dynamics loop
+ for i in range(0, number_of_steps):
+ # At each step, calculate the forces on each particle
+ # and get acceleration
+ system.particles, system.distances, system.forces = comp.compute_forces(system.particles,
+ system.distances,
+ system.forces, system.box_length)
+ # Run the equations of motion integrator algorithm
+ system.particles = md.velocity_verlet(system.particles, system.timestep_length, system.box_length)
+ # Sample the thermodynamic and structural parameters of the system
+ system = md.sample(system.particles, system.box_length, system.initial_particles, system)
+ # Allow the system to interact with a heat bath
+ system.particles = comp.heat_bath(system.particles, system.temperature_sample, temperature)
+ # Iterate the time
+ system.time += system.timestep_length
+ system.step += 1
+ # At a given frequency sample the positions and plot the RDF
+ if system.step % sample_frequency == 0:
+ sample_system.update(system)
+ return system
+
+
+[docs]def md_nve(number_of_particles, temperature, box_length, number_of_steps, sample_frequency):
+ # Creates the visualisation environment
+ #%matplotlib notebook
+ # Initialise the system
+ system = md.initialise(number_of_particles, temperature, box_length, 'square')
+ # This sets the sampling class
+ sample_system = sample.Interactions(system)
+ # Start at time 0
+ system.time = 0
+ # Begin the molecular dynamics loop
+ for i in range(0, number_of_steps):
+ # At each step, calculate the forces on each particle
+ # and get acceleration
+ system.particles, system.distances, system.forces = comp.compute_forces(system.particles,
+ system.distances,
+ system.forces, system.box_length)
+ # Run the equations of motion integrator algorithm
+ system.particles = md.velocity_verlet(system.particles, system.timestep_length, system.box_length)
+ # Sample the thermodynamic and structural parameters of the system
+ system = md.sample(system.particles, system.box_length, system.initial_particles, system)
+ # Iterate the time
+ system.time += system.timestep_length
+ system.step += 1
+ # At a given frequency sample the positions and plot the RDF
+ if system.step % sample_frequency == 0:
+ sample_system.update(system)
+ return system
+
+
+import numpy as np
+from pylj import comp, util
+
+
+[docs]def initialise(number_of_particles, temperature, box_length, init_conf, timestep_length=5e-3):
+ """Initialise the particle positions (this can be either as a square or random arrangement), velocities (based on
+ the temperature defined, and calculate the initial forces/accelerations.
+
+ Parameters
+ ----------
+ number_of_particles: int
+ Number of particles to simulate.
+ temperature: float
+ Initial temperature of the particles, in Kelvin.
+ box_length: float
+ Length of a single dimension of the simulation square, in Angstrom.
+ init_conf: string, optional
+ The way that the particles are initially positioned. Should be one of:
+
+ - 'square'
+ - 'random'
+ timestep_length: float (optional)
+ Length for each Velocity-Verlet integration step, in seconds.
+
+ Returns
+ -------
+ System
+ System information.
+ """
+ system = util.System(number_of_particles, temperature, box_length, init_conf=init_conf,
+ timestep_length=timestep_length)
+ v = np.sqrt(2 * system.init_temp)
+ theta = 2 * np.pi * np.random.randn(system.particles.size)
+ system.particles['xvelocity'] = v * np.cos(theta)
+ system.particles['yvelocity'] = v * np.sin(theta)
+ system.particles, system.distances, system.forces = comp.compute_forces(system.particles, system.distances,
+ system.forces, system.box_length)
+ return system
+
+[docs]def velocity_verlet(particles, timestep_length, box_length):
+ """Uses the Velocity-Verlet integrator to move forward in time. The
+
+ Updates the particles positions and velocities in terms of the Velocity Verlet algorithm. Also calculates the
+ instanteous temperature, pressure, and force and appends these to the appropriate system array.
+
+ Parameters
+ ----------
+ particles: util.particle_dt, array_like
+ Information about the particles.
+ timestep_length: float
+ Length for each Velocity-Verlet integration step, in seconds.
+ box_length: float
+ Length of a single dimension of the simulation square, in Angstrom.
+
+ Returns
+ -------
+ util.particle_dt, array_like:
+ Information about the particles, with new positions and velocities.
+ """
+ xposition_store = particles['xposition']
+ yposition_store = particles['yposition']
+ [particles['xposition'], particles['yposition']] = update_positions([particles['xposition'],
+ particles['yposition']],
+ [particles['xvelocity'],
+ particles['yvelocity']],
+ [particles['xacceleration'],
+ particles['yacceleration']], timestep_length,
+ box_length)
+ [particles['xvelocity'], particles['yvelocity']] = update_velocities([particles['xvelocity'], particles['yvelocity']],
+ [particles['xacceleration'],
+ particles['yacceleration']], timestep_length)
+ particles['xprevious_position'] = xposition_store
+ particles['yprevious_position'] = yposition_store
+ return particles
+
+[docs]def sample(particles, box_length, initial_particles, system):
+ """Sample parameters of interest in the simulation.
+
+ Parameters
+ ----------
+ particles: util.particle_dt, array_like
+ Information about the particles.
+ box_length: float
+ Length of a single dimension of the simulation square, in Angstrom.
+ initial_particles: util.particle_dt, array-like
+ Information about the initial particle conformation.
+ system: System
+ Details about the whole system
+
+ Returns
+ -------
+ System:
+ Details about the whole system, with the new temperature, pressure, msd, and force appended to the appropriate
+ arrays.
+ """
+ temperature_new = util.calculate_temperature(particles)
+ pressure_new = comp.calculate_pressure(particles, box_length, temperature_new)
+ msd_new = util.calculate_msd(particles, initial_particles, box_length)
+ system.temperature_sample = np.append(system.temperature_sample, temperature_new)
+ system.pressure_sample = np.append(system.pressure_sample, pressure_new)
+ system.force_sample = np.append(system.force_sample, np.sum(system.forces))
+ system.msd_sample = np.append(system.msd_sample, msd_new)
+ return system
+
+[docs]def update_positions(positions, velocities, accelerations, timestep_length, box_length):
+ """Update the particle positions using the Velocity-Verlet integrator.
+
+ Parameters
+ ----------
+ positions: (2, N) array_like
+ Where N is the number of particles, and the first row are the x positions and the second row the y positions.
+ velocities: (2, N) array_like
+ Where N is the number of particles, and the first row are the x velocities and the second row the y velocities.
+ accelerations: (2, N) array_like
+ Where N is the number of particles, and the first row are the x accelerations and the second row the y
+ accelerations.
+ timestep_length: float
+ Length for each Velocity-Verlet integration step, in seconds.
+ box_length: float
+ Length of a single dimension of the simulation square, in Angstrom.
+
+ Returns
+ -------
+ (2, N) array_like:
+ Updated positions.
+ """
+ positions[0] += velocities[0] * timestep_length + 0.5 * accelerations[0] * timestep_length * timestep_length
+ positions[1] += velocities[1] * timestep_length + 0.5 * accelerations[1] * timestep_length * timestep_length
+ positions[0] = positions[0] % box_length
+ positions[1] = positions[1] % box_length
+ return [positions[0], positions[1]]
+
+[docs]def update_velocities(velocities, accelerations, timestep_length):
+ """Update the particle velocities using the Velocity-Verlet algoritm.
+
+ Parameters
+ ----------
+ velocities: (2, N) array_like
+ Where N is the number of particles, and the first row are the x velocities and the second row the y velocities.
+ accelerations: (2, N) array_like
+ Where N is the number of particles, and the first row are the x accelerations and the second row the y
+ accelerations.
+ timestep_length: float
+ Length for each Velocity-Verlet integration step, in seconds.
+
+ Returns
+ -------
+ (2, N) array_like:
+ Updated velocities.
+ """
+ velocities[0] += 0.5 * accelerations[0] * timestep_length
+ velocities[1] += 0.5 * accelerations[1] * timestep_length
+ return [velocities[0], velocities[1]]
+
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+[docs]class Scattering(object):
+ """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.
+ """
+ def __init__(self, system):
+ fig, ax = environment(4)
+ 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) + 0.5])
+ self.fig.canvas.draw()
+
+ 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([np.amin(iq) - np.amax(iq) * 0.05, np.amax(iq) + np.amax(iq) * 0.05])
+ self.ax[1, 1].set_xlim([np.amin(x), np.amax(x)])
+
+
+[docs]class Interactions(object):
+ """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.
+ """
+ def __init__(self, system):
+ fig, ax = environment(4)
+
+ setup_cellview(ax[0, 0], system)
+ setup_forceview(ax[0, 1])
+ setup_pressureview(ax[1, 0])
+ setup_tempview(ax[1, 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[0, 1], system)
+ update_tempview(self.ax[1, 1], system)
+ update_pressureview(self.ax[1, 0], system)
+
+ self.fig.canvas.draw()
+
+
+[docs]class JustCell(object):
+ """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.
+ """
+ def __init__(self, system):
+ fig, ax = environment(1)
+
+ setup_cellview(ax, 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, system)
+
+ self.fig.canvas.draw()
+
+
+[docs]class RDF(object):
+ """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.
+ """
+ def __init__(self, system):
+ fig, ax = environment(2)
+ 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) + 0.5])
+ self.fig.canvas.draw()
+
+
+[docs]def environment(panes):
+ """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.
+
+ 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.
+ """
+ fig, ax = plt.subplot()
+ if panes == 1:
+ fig, ax = plt.subplots(figsize=(4, 4))
+ elif panes == 2:
+ fig, ax = plt.subplots(1, 2, figsize=(8, 4))
+ elif panes == 4:
+ fig, ax = plt.subplots(2, 2, figsize=(8, 8))
+ else:
+ AttributeError('The only options for the number of panes are 1, 2, or 4')
+ return fig, ax
+
+
+[docs]def setup_cellview(ax, system):
+ """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.
+ """
+ xpos = system.particles['xposition']
+ ypos = system.particles['yposition']
+ mk = (1052.2 / (system.box_length - 0.78921) - 1.2174)
+ 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):
+ """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_rdfview(ax, system):
+ """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/Å', fontsize=16)
+
+
+[docs]def setup_diffview(ax):
+ """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('log(I[q])', fontsize=16)
+ ax.set_xlabel('log(q)/Å$^{-1}$', fontsize=16)
+
+
+[docs]def setup_pressureview(ax):
+ """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('Pressure/Pa', fontsize=16)
+ ax.set_xlabel('Time/s', fontsize=16)
+
+
+[docs]def setup_tempview(ax):
+ """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 update_cellview(ax, system):
+ """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):
+ """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.5, 100))
+ gr = hist / (system.number_of_particles * (system.number_of_particles / system.box_length ** 2) * np.pi *
+ (bin_edges[:-1] + 0.5 / 2.) * 0.5)
+ average_rdf.append(gr)
+ x = bin_edges[:-1] + 0.5 / 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.05])
+
+
+[docs]def update_diffview(ax, system, average_diff, q):
+ """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.
+ """
+ hist, bin_edges = np.histogram(system.distances, bins=np.linspace(0, system.box_length/2 + 0.5, 100))
+ gr = hist / (system.number_of_particles * (system.number_of_particles / system.box_length ** 2) * np.pi *
+ (bin_edges[:-1] + 0.5 / 2.) * 0.5)
+ x2 = np.log(np.fft.rfftfreq(len(gr))[5:])
+ y2 = np.log(np.fft.rfft(gr)[5:])
+ average_diff.append(y2)
+ q.append(x2)
+ line1 = ax.lines[0]
+ line1.set_xdata(x2)
+ line1.set_ydata(y2)
+ ax.set_ylim([np.amin(y2) - np.amax(y2) * 0.05, np.amax(y2) + np.amax(y2) * 0.05])
+ ax.set_xlim([np.amin(x2), np.amax(x2)])
+
+
+[docs]def update_forceview(ax, system):
+ """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 * 1e-10)
+ 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 * 1e-10)-np.amax(system.force_sample * 1e-10) * 0.05,
+ np.amax(system.force_sample * 1e-10)+np.amax(system.force_sample * 1e-10) * 0.05)
+
+
+[docs]def update_tempview(ax, system):
+ """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_pressureview(ax, system):
+ """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]
+ line.set_ydata(system.pressure_sample * 1e10)
+ 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.pressure_sample * 1e10) - np.amax(system.pressure_sample * 1e10) * 0.05,
+ np.amax(system.pressure_sample * 1e10) + np.amax(system.pressure_sample * 1e10) * 0.05)
+
+
+[docs]def setup_msdview(ax):
+ """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/Å$^2$', fontsize=16)
+ ax.set_xlabel('Time/s', fontsize=16)
+
+
+[docs]def update_msdview(ax, system):
+ """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)
+
+import numpy as np
+
+
+[docs]class System:
+ """Simulation system.
+
+ This class is designed to store all of the information about the job that is being run. This includes the particles
+ object, as will as sampling objects such as the temperature, pressure, etc. arrays.
+
+ Parameters
+ ----------
+ number_of_particles: int
+ Number of particles to simulate.
+ temperature: float
+ Initial temperature of the particles, in Kelvin.
+ box_length: float
+ Length of a single dimension of the simulation square, in Angstrom.
+ init_conf: string, optional
+ The way that the particles are initially positioned. Should be one of:
+
+ - 'square'
+ - 'random'
+ timestep_length: float (optional)
+ Length for each Velocity-Verlet integration step, in seconds.
+ """
+ def __init__(self, number_of_particles, temperature, box_length, init_conf='square', timestep_length=5e-3):
+ self.number_of_particles = number_of_particles
+ self.init_temp = temperature
+ if box_length <= 600:
+ self.box_length = box_length
+ else:
+ raise AttributeError('With a box length of {} the particles are probably too small to be seen in the '
+ 'viewer. Try something (much) less than 600.'.format(box_length))
+ self.timestep_length = timestep_length
+ self.particles = None
+ if init_conf == 'square':
+ self.square()
+ elif init_conf == 'random':
+ self.random()
+ else:
+ raise NotImplementedError('The initial configuration type {} is not recognised. '
+ 'Available options are: square or random'.format(init_conf))
+ self.step = 0
+ self.time = 0.
+ self.distances = np.zeros(self.number_of_pairs())
+ self.forces = np.zeros(self.number_of_pairs())
+ self.temperature_sample = np.array([])
+ self.pressure_sample = np.array([])
+ self.force_sample = np.array([])
+ self.msd_sample = np.array([])
+ self.initial_particles = np.array(self.particles)
+
+[docs] def number_of_pairs(self):
+ """Calculates the number of pairwise interactions in the simulation.
+
+ Returns
+ -------
+ int:
+ Number of pairwise interactions in the system.
+ """
+ return int((self.number_of_particles - 1) * self.number_of_particles / 2)
+
+[docs] def square(self):
+ """Sets the initial positions of the particles on a square lattice.
+ """
+ part_dt = particle_dt()
+ self.particles = np.zeros(self.number_of_particles, dtype=part_dt)
+ m = int(np.ceil(np.sqrt(self.number_of_particles)))
+ d = self.box_length / m
+ n = 0
+ for i in range(0, m):
+ for j in range(0, m):
+ if n < self.number_of_particles:
+ self.particles[n]['xposition'] = (i + 0.5) * d
+ self.particles[n]['yposition'] = (j + 0.5) * d
+ n += 1
+
+[docs] def random(self):
+ """Sets the initial positions of the particles in a random arrangement.
+ """
+ part_dt = particle_dt()
+ self.particles = np.zeros(self.number_of_particles, dtype=part_dt)
+ self.particles['xposition'] = np.random.uniform(0, self.box_length, self.number_of_particles)
+ self.particles['yposition'] = np.random.uniform(0, self.box_length, self.number_of_particles)
+
+[docs]def pbc_correction(position, cell):
+ """Correct for the periodic boundary condition.
+
+ Parameters
+ ----------
+ position: float
+ Particle position.
+ cell: float
+ Cell vector.
+
+ Returns
+ -------
+ float:
+ Corrected particle position."""
+ if np.abs(position) > 0.5 * cell:
+ position *= 1 - cell / np.abs(position)
+ return position
+
+[docs]def calculate_temperature(particles):
+ """Determine the instantaneous temperature of the system.
+
+ Parameters
+ ----------
+ particles: util.particle_dt, array_like
+ Information about the particles.
+
+ Returns
+ -------
+ float:
+ Calculated instantaneous simulation temperature.
+ """
+ k = 0
+ for i in range(0, particles['xposition'].size):
+ v = np.sqrt(particles['xvelocity'][i] * particles['xvelocity'][i] + particles['yvelocity'][i] *
+ particles['yvelocity'][i])
+ k += 0.5 * v * v
+ return k / particles['xposition'].size
+
+[docs]def calculate_msd(particles, initial_particles, box_length):
+ """Determines the mean squared displacement of the particles in the system.
+
+ Parameters
+ ----------
+ particles: util.particle_dt, array_like
+ Information about the particles.
+ initial_particles: util.particle_dt, array_like
+ Information about the initial state of the particles.
+ box_length: float
+ Size of the cell vector.
+
+ Returns
+ -------
+ float:
+ Mean squared deviation for the particles at the given timestep.
+ """
+ dx = particles['xposition'] - initial_particles['xposition']
+ dy = particles['yposition'] - initial_particles['yposition']
+ for i in range(0, particles['xposition'].size):
+ if (np.abs(dx[i]) > 0.5 * box_length):
+ dx[i] *= 1 - box_length / np.abs(dx[i])
+ if (np.abs(dy[i]) > 0.5 * box_length):
+ dy[i] *= 1 - box_length / np.abs(dy[i])
+ dr = np.sqrt(dx * dx + dy * dy)
+ return np.average(dr ** 2)
+
+[docs]def particle_dt():
+ """Builds the data type for the particles, this consists of:
+
+ - xposition and yposition
+ - xvelocity and yvelocity
+ - xacceleration and yacceleration
+ - xprevious_position and ypresvious_position
+ - xforce and yforce
+ - energy
+ """
+ return np.dtype([('xposition', np.float64), ('yposition', np.float64), ('xvelocity', np.float64),
+ ('yvelocity', np.float64), ('xacceleration', np.float64), ('yacceleration', np.float64),
+ ('xprevious_position', np.float64), ('yprevious_position', np.float64), ('xforce', np.float64),
+ ('yforce', np.float64), ('energy', np.float64)])
+