What goes around comes around

• Topic: periodic boundary conditions

• Task:

  1. Implement the function that calculates the squared distance between particles, taking into account the periodicity of the system.

  2. Run the simulation. The small particle should bounce off the larger ones. Does it?

  3. Draw the system so that a few periodic images are also visible.

• Template: pbc.py

• Data: supercell.txt

• Further reading: https://en.wikipedia.org/wiki/Periodic_boundary_conditions

pbc.py

class pbc.Atom(position, velocity, mass=1.0)

A point like object.

An atom has a position (a 3-vector), a velocity (3-vector) and a mass (a scalar).

Parameters:

• position (array) – coordinates \([x, y, z]\)

• velocity (array) – velocity components \([v_x, v_y, v_z]\)

• mass (float) – mass \(m\)

accelerate(force, dt)

Set a new velocity for the particle as

\[\vec{v}(t+\Delta t) = \vec{v}(t) + \frac{1}{2m}\vec{F} \Delta t\]

Parameters:

• force (array) – force acting on the planet \([F_x, F_y, F_z]\)

• dt (float) – time step \(\Delta t\)

move(force, dt)

Move the atom.

Parameters:

shift (array) – coordinate change \([\Delta x, \Delta y, \Delta z]\)

save_position()

Save the current position of the particle in the list ‘trajectory’.

Note: in a real large-scale simulation one would never save trajectories in memory. Instead, these would be written to a file for later analysis.

class pbc.PeriodicBox(lattice)

Class representing a simulation box with periodic boundaries.

The box is orthogonal, i.e., a rectangular volume. As such, it is specified by the lengths of its edges (lattice constants).

Parameters:

lattice (array) – lattice constants

distance_squared(particle1, particle2)

Calculates and returns the square of the distance between two particles.

\[r^2_{ij} = (x_i-x_j)^2 + (y_i-y_j)^2 + (z_i-z_j)^2.\]

In a periodic system, each particle has an infinite number of periodic copies. Therefore the distance between two particles is not unique. The function returns the shortest such distance, that is, the distance between the the periodic copies which are closest to each other.

Note

This function is incomplete!

Parameters:

• particle1 (pbc.Atom) – the first body

• particle2 (pbc.Atom) – the second body

Returns:

the squared distance \(r^2_{ij}\)

Return type:

float

shift_inside_box(position)

If the given position (3-vector) is outside the box, it is shifted by multiple of lattice vectors until the new position is inside the box. That is, the function transforms the position vector to an equivalent position inside the box.

Parameters:

position (array) – the position to be shifted

vector(particle1, particle2)

Returns the vector pointing from the position of particle1 to the position of particle2.

\[\vec{r}_{i \to j} = \vec{r}_j - \vec{r}_i\]

In a periodic system, each particle has an infinite number of periodic copies. Therefore the displacement between two particles is not unique. The function returns the shortest such displacement vector.

Parameters:

• particle1 (pbc.Atom) – the first body

• particle2 (pbc.Atom) – the second body

Returns:

components of \(\vec{r}_{i \to j}\), \([x_{i \to j}, y_{i \to j}, z_{i \to j}]\)

Return type:

array

pbc.animate(particles, box, multiply=[3, 3])

Animates the simulation.

Parameters:

• particles (list) – list of pbc.Atom objects

• box (pbc.PeriodicBox) – supercell

• multiply (array) – number of periodic images to draw in x and y directions

pbc.calculate_forces(particles, box, parameters=[1.0, 0.1])

Calculates the total force applied on each atom.

The forces are returned as a numpy array where each row contains the force on an atom and the columns contain the x, y and z components.

[ [ fx0, fy0, fz0 ],
  [ fx1, fy1, fz1 ],
  [ fx2, fy2, fz2 ],
  ...               ]

Parameters:

• atoms (list) – a list of pbc.Atom objects

• box (pbc.PeriodicBox) – supercell

• parameters (list) – a list of physical parameters

Returns:

forces

Return type:

array

pbc.calculate_kinetic_energy(particles)

Calculates the total kinetic energy of the system.

\[K_\text{total} = \sum_i \frac{1}{2} m_i v_i^2.\]

Parameters:

particles (list) – a list of pbc.Atom objects

Returns:

kinetic energy \(K\)

Return type:

float

pbc.calculate_momentum(particles)

Calculates the total momentum of the system.

\[\vec{p}_\text{total} = \sum_i \vec{p}_i = \sum_i m_i \vec{v}_i\]

Parameters:

particles (list) – a list of pbc.Atom objects

Returns:

momentum components \([p_x, p_y, p_z]\)

Return type:

array

pbc.calculate_potential_energy(particles, box, parameters=[1.0, 0.1])

Calculates the total potential energy of the system.

The potential energy is calculated using the Lennard-Jones model

\[U = \sum_{i \ne j} 4 \epsilon \left( \frac{ \sigma^{12} }{ r^{12}_{ij} } - \frac{ \sigma^6 }{ r^6_{ij} } \right).\]

Parameters:

• particles (list) – a list of pbc.Atom objects

• box (pbc.PeriodicBox) – supercell

• parameters (list) – list of parameters

Returns:

potential energy \(U\)

Return type:

float

pbc.draw(frame, xtraj, ytraj, ztraj, bounds)

Draws a representation of the particle system as a scatter plot.

Used for animation.

Parameters:

• frame (int) – index of the frame to be drawn

• xtraj (array) – x-coordinates of all particles at different animation frames

• ytraj (array) – y-coordinates at all particles at different animation frames

• ztraj (array) – z-coordinates at all particles at different animation frames

• bounds (array) – list of lower and upper bounds for the plot as [[xmin, xmax], [ymin, ymax]]

pbc.main()

The main program.

The program reads the system from a file, runs the simulation, and plots the trajectory.

pbc.print_progress(step, total)

Prints a progress bar.

Parameters:

• step (int) – progress counter

• total (int) – counter at completion

pbc.read_particles_from_file(filename)

Reads the properties of planets from a file.

Each line should define a single Planet listing its position in cartesian coordinates, velocity components and mass, separated by whitespace:

x0 y0 z0 vx0 vy0 vz0 m0
x1 y1 z1 vx1 vy1 vz1 m1
x2 y2 z2 vx2 vy2 vz2 m2
x3 y3 z3 vx3 vy3 vz3 m3
...

Parameters:

filename (str) – name of the file to read

Returns:

list of pbc.Atom objects

Return type:

list

pbc.show_trajectories(particles, box, tail=10, skip=10, multiply=[3, 3])

Plot a 2D-projection of the trajectory of the system.

The function creates a plot showing the current and past positions of particles.

Parameters:

• particles (list) – list of Planet objects

• box (pbc.PeriodicBox) – supercell

• tail (int) – the number of past positions to include in the plot

• skip (int) – only every nth past position is plotted - skip is the number n, specifying how far apart the plotted positions are in time

• multiply (array) – number of periodic images to draw in x and y directions

pbc.update_positions(particles, forces, dt)

Update the positions of all particles using pbc.Atom.move() according to

\[\vec{r}(t+\Delta t) = \vec{r}(t) + \vec{v} \Delta t + \frac{1}{2m}\vec{F} (\Delta t)^2\]

Parameters:

• particles (list) – a list of pbc.Atom objects

• force (array) – array of forces on all bodies

• dt (float) – time step \(\Delta t\)

pbc.update_positions_no_force(particles, dt)

Update the positions of all particles using pbc.Atom.move() according to

\[\vec{r}(t+\Delta t) = \vec{r}(t) + \vec{v} \Delta t\]

Parameters:

• particles (list) – a list of pbc.Atom objects

• dt (float) – time step \(\Delta t\)

pbc.update_velocities(particles, forces, dt)

Update the positions of all particles using pbc.Atom.accelerate() according to

\[\vec{v}(t+\Delta t) = \vec{v}(t) + \frac{1}{m}\vec{F} \Delta t\]

Parameters:

• particles (list) – a list of pbc.Atom objects

• force (array) – array of forces on all bodies

• dt (float) – time step \(\Delta t\)

pbc.velocity_verlet(particles, box, dt, time, trajectory_dt=1.0)

Verlet algorithm for integrating the equations of motion, i.e., advancing time.

There are a few ways to implement Verlet. The leapfrog version works as follows: First, forces are calculated for all particles and velocities are updated by half a time step, \(\vec{v}(t+\frac{1}{2}\Delta t) = \vec{v}(t) + \frac{1}{2m}\vec{F} \Delta t\). Then, these steps are repeated:

• Positions are updated by a full time step using velocities but not forces,

  \[\vec{r}(t+\Delta t) = \vec{r}(t) + \vec{v}(t+\frac{1}{2}\Delta t) \Delta t.\]

• Forces are calculated at the new positions, \(\vec{F}(t + \Delta t)\).

• Velocities are updated by a full time step using the forces

  \[\vec{v}(t+\frac{3}{2}\Delta t) = \vec{v}(t+\frac{1}{2}\Delta t) + \frac{1}{m}\vec{F}(t+\Delta t) \Delta t\]

These operations are done using the methods calculate_forces(), update_velocities() and update_positions_no_force().

Because velocities were updated by half a time step in the beginning of the simulation, positions and velocities are always offset by half a timestep. You always use the one that has advanced further to update the other and this results in a stable algorithm.

Parameters:

• particles (list) – a list of pbc.Atom objects

• box (pbc.PeriodicBox) – supercell

• dt (float) – time step \(\Delta t\)

• time (float) – the total system time to be simulated

• trajectory_dt (float) – the positions of particles are saved at these these time intervals - does not affect the dynamics in any way

pbc.write_particles_to_file(particles, filename)

Write the configuration of the system in a file.

The format is the same as that specified in read_particles_from_file().

Parameters:

• particles (list) – list of pbc.Atom objects

• filename (str) – name of the file to write

pbc.write_xyz_file(particles, filename)

Write the configuration of the system in a file.

The information is written in so called xyz format, which many programs can parse.

Parameters:

• particles (list) – list of pbc.Atom objects

• filename (str) – name of the file to write