Homework: It’s a gas at the bridge

• Topic: thermodynamics, ideal gas, structural design, constraints

• Template: A4_2dmd.py

• Data:

  • A4_gas_simulation.txt

  • A4_gas_system.txt

  • A4_bridge_simulation.txt

  • A4_bridge_system.txt

• Further reading:

  • https://en.wikipedia.org/wiki/Ideal_gas_law

  • https://en.wikipedia.org/wiki/Langevin_equation

  • https://en.wikipedia.org/wiki/Equipartition_theorem

This assignment consists of two parts. Both parts will be simulated using the same code, a 2-dimensional molecular dynamics simulator, but the systems to be simulated are very different.

In the first part, you simulate a gas consisting of atoms and molecules. It is a simple system and a good place to test your understanding of both molecular dynamics and thermodynamics.

In the second part, your task is to build a structure using a particle model. This is a design problem with no right or wrong answers. There are, however, better and worse solutions.

Part I: Theory

Ideal gas is often the first model students learn in introductory thermodynamics, since it is so simple yet useful. The ideal gas model describes gases as a collection of particles in a large volume.

The ideal gas law, derivable both from experiments and statistical physics, states that the pressure \(p\), volume \(V\), temperature \(T\) and number of molecules \(N\) are connected by the equation

\[p V = N k_B T,\]

where \(k_B\) is the Boltzmann constant.

The ideal gas law is not exact. Gases obey it approximately, when their density is small enough and temperature high enough. When gases are cooled down or compressed, they typically turn into a liquid or a solid, and the ideal gas law is only valid when the gas is not even near such a phase transformation.

Your first task is to verify that this law approximately holds in a simulation. Note that our simulation is run in 2D, which means that volume of the gas is actually the area of the simulation box and pressure is defined as force per unit length (not area). The simulation also uses scaled units where \(k_B = 1\). Do not expect to see realistic values for the physical quantities.

The gas consists of 20 atoms and 10 diatomic molecules. The atoms and molecules interact via a Lennard-Jones potential,

\[U(r) = 4 \epsilon \left( \frac{ \sigma^{12} }{ r^{12} } - \frac{ \sigma^6 }{ r^6 } \right),\]

and the two atoms in each molecule connect to each other via spring-like bonds with potential energy

\[U(r) = \frac{1}{2} k (r - r_0)^2.\]

For a classical gas, statistical physics states that all quadratic (i.e., second power, \(v^2\) or \(r^2\)) forms of energy should, on average, contain the same amount of energy in each degree of freedom (DOF),

\[\langle E_\text{DOF} \rangle = \frac{1}{2} k_B T.\]

This is also known as the energy equipartition principle.

Here “a degree of freedom” roughly means the same thing as “a distinct way to move”. An atom has two degrees of freedom in 2D: movement left-right and up-down. A diatomic molecule can also rotate around its center of mass and vibrate by having the bond connecting the two atoms stretch and contract. Translation and rotation are just movement so they only store kinetic energy, but vibration involves a spring that also stores potential energy. Therefore a diatomic molecule has four degrees of freedom with kinetic energy (left-right, up-down, rotation, vibration) and one degree of freedom with potential energy (vibration).

Your second task is to check if energy is distributed as theory suggests. Note that the Lennard-Jones potential is not quadratic and so the energy it stores behaves differently to what was described above. You do not need to pay attention to it in this task.

The simulation can be run in both microcanonical (constant energy) and canonical (constant average temperature) ensembles. The microcanonical simulation is done using normal Newtonian dynamics whereas the canonical simulation is done using Langevin dynamics. In Langevin dynamics, one adds a friction force and a random force to each particle, so that the equations of motion become

\[\vec{a} = \frac{1}{m} \vec{F} - \gamma \vec{v} + \frac{1}{m} \vec{F}_\text{random}.\]

Here \(\vec{a}\) is acceleration, \(m\) mass, \(\vec{F}\) real force, \(\gamma\) friction coefficient, \(\vec{v}\) velocity and \(\vec{F}_\text{random}\) a normally distributed random force with no correlations and standard deviation \(\sigma = \sqrt{ 2 \gamma m T }\).

The friction force opposes particle movement, so it tends to slow them down and lower the temperature. The random force pushes the particles around adding to their speed and increasing the temperature. The magnitudes of these forces are tuned so that the average temperature of the system will approach \(T\). Magnitude of the friction coefficient \(\gamma\) detemines how strongly dynamics are altered. For \(\gamma = 0\), you have normal Newtonian dynamics. If \(\gamma > 0\), the system is driven towards temperature \(T\) and this coupling grows stronger as \(\gamma\) is made larger. You could try a value of about \(\gamma = 0.1\) to simulate a canonical ensemble. Feel free to experiment with different values and see what happens.

Part II: Instructions

The second problem is more of an engineering challenge than a fundamental physics problem. The files with “bridge” in their name define a system where an elastic wire is taut horizontally between immovable supports. A single heavy particle is dropped on the wire, and depending on how hard the impact is, the particle either bounces back or the wire snaps.

This happens because the wire is made of a chain of particles where each pair of adjacent particles is connected by a spring-like interaction, just like the diatomic molecules in the gas simulation. However, the springs have a maximum length, and if the particles are separated farther than that length, no force should be applied between the particles. This essentially breaks the spring (although the bond will reform if the particles are brought back together).

The elastic band in the given template snaps quite easily. Your task is to design and build a more robust bridge between the two supports.

You must build your bridge using similar particles as those in the template (mass: 1, name: bridge), and you are not allowed to change the properties of the falling particle or the supports. Your final solution must not contain more than 50 particles and 100 bonds in total. Eight particles are already taken, so you have 42 to work with.

Your solution will be tested with an ever increasing pulling force applied on the heavy particle and scored according to how hard an impact the bridge can take. For reference, the elastic band holds, when the pulling force magnitude is 5.0, but snaps for 10.0. Also note that the program that checks the solutions will treat any simulation failed, if the heavy particle reaches the botom of the simulation box. So you cannot make your bridge too elastic either.

Input syntax

The simulator reads the system and simulation data from two input files.

The system file should contain

• the size of the simulation box (<box>)

• properties of the particles (<atom>).

The simulation file should contain

• potential energy parameters (<interactions>)

• temperature and thermostat parameters (<temperature>)

• timing parameters (<time>)

• data recording and visualization options (<recording>)

• optionally, constraints (<constraint>).

There is no need to list the full syntax of these files here. The given input files already contain suitable values for most parameters, and comments point out the values you may want to change. In the gas simulation, you only need to change the temperature and friction coefficient \(\gamma\). In the bridge simulation, however, you need to make a substantial addition to A4_bridge_system.txt, since you are required to design and add your own structure there. You may also want to change the force pulling the heavy particle down to test your solution.

You are supposed to build your structure using particles, and each particle is defined by wrapping its properties in the <atom> tag. The most important properties are coordinates x and y, and connect, which defines a spring-like connection between the atom in question and the atom whose index matches the number given in the connect order. For instance,

<atom>
x: 1
y: 1
name: bridge
index: 10
connect: 11
</atom>

<atom>
x: 2
y: 1
name: bridge
index: 11
</atom>

will create two atoms with indices 10 and 11, and a spring interaction between them.

You only need to declare a connection once. You could also write connect: 10 in the properties of the second atom, but you don’t have to. The program automatically makes the connections reciprocal.

Task

To complete this assignment, you need to complete the following tasks, submit your code as a plain .py file, submit your final bridge desing as the A4_bridge_system.txt file, and answer the questions. In the computing tasks, you must implement the incomplete functions so that they operate exactly as described in the documentation.

• Computing task I:

  1. Implement a function to calculate the force due to the spring interaction. Make sure to test your implementation properly!

  2. Run a microcanonical simulation. Check that the total energy and momentum is conserved and observe how the gas behaves.

  3. Run a canonical simulation in at least two different temperatures. Observe the gas again. Compare to the previous simulation.

  4. Check the ideal gas law and energy equipartition in all cases.

• Computing task II:

  1. Design a bridge structure using at most 50 particles and 100 bonds in total. Build it in A4_bridge_system.txt.

  2. Test the robustness of your bridge by simulating a heavy object falling on it. Increase the force pulling the heavy particle down and see how hard you can pull before your bridge breaks apart.

  3. If you are not content with the result, try making some adjustments in your design and try again.

• Questions to answer:

  1. Explain in detail how you implemented the spring force.

  2. Explain how you measure the necessary quantities in the gas simulation. Include plots in your explanation and consider the role of statistics.

  3. What happens in the gas simulation, when you change from Newtonian to Langevin dynamics? How do the particles move? Is there a difference in measured quantities?

  4. How well does the gas simulation follow ideal gas law and energy equipartition? Are these laws valid in all of the simulations you ran or only in some? If there are discrepancies, can you explain them? Include data in your explanation and be mindful of the accuracy and error estimates of your results.

  5. Explain, in your own words, your strategy in the bridge design problem and what lead you to the final design. Include a picture of your bridge in your answer. You may also include pictures of other designs you considered.

A4_2dmd.py

class A4_2dmd.Animation(root, particles, lattice_parameters, force_parameters)

A class for creating an animation based on simulation results.

The animation is drawn in a graphical window generated using tkinter.

Parameters:

• root (tkinter.canvas) – graphical engine object from tkinter

• particles (list) – list of Particle objects

• lattice_parameters (list) – lattice parameters [Lx, Ly]

• force_parameters (list) – interaction parameters

clock()

Updates the animation at a constant frequency.

The graphical representation is updated to the next frame using Animation.shift_particles(). This is repeated after a short delay to run the animation.

get_bond_width_and_colour(r)

Gives the width and colour of a bond.

To visualize bond straining, we draw bonds that are close to breaking as thin red lines. This function calculates the proper width and colour for a bond of given length.

Parameters:

r (float) – length of the bond

Returns:

line width, colour hex

Return type:

float, string

get_particle_size_and_colour(particle)

Gives the radius and colour of a particle.

To make different particles more easily distinguishable, we draw them with various sizes and colours.

Particles are distinguished according to their name. This function only knows the names used in this exercise.

Parameters:

particle (Particle) – the particle

Returns:

circle radius, colour string

Return type:

float, string

rgb_to_hex(rgb)

Transforms RGB colour to hexadecimal representation.

Parameters:

rgb (vector) – 3-component vector storing Red-Green-Blue values [0-255]

Returns:

colour hex

Return type:

string

shift_particles(frame)

Updates the graphical representation.

Shifts the graphical elements representing particles and bonds to the positions corresponding to the given frame.

If the given frame number is too high, the animation simply loops back to beginning.

Parameters:

frame (int) – number of the frame to be drawn

class A4_2dmd.Particle(position, velocity, mass, name, index, connections)

A class representing a point-like particle.

Parameters:

• position (list) – coordinates [x, y]

• velocity (list) – components [vx, vy]

• mass (float) – particle mass

• name (str) – a name for the particle, to help distinguish it

• index (int) – an identifying number

• connections (list) – indices of the atoms that connect to this atom by a spring

position

coordinates [x, y]

Type:

list

velocity

components [vx, vy]

Type:

list

force

components [Fx, Fy]

Type:

list

mass

particle mass

Type:

float

name

a name for the particle, to help distinguish it

Type:

str

index

an identifying number

Type:

int

connected_atoms

list of Particle objects connected to this atom

Type:

list

constraint_type

one of FREEZE_TAG, F_TAG or V_TAG

Type:

str

constraint_value

value of the force or velocity constraint

Type:

list

accelerate(dt, gamma=0)

Set a new velocity for the particle as

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

By default, the force \(F\) is the total force applied on the particle by all other particles. It should be precalculated and stored in the attributes of the particle.

If a non-zero gamma is given, a drag force \(\vec{F}_\text{drag} = - \gamma m \vec{v}\) is also applied.

If the particle is constrained, the constraints are also applied:

• A frozen particle always has zero velocity.

• A velocity constrained particle has constant velocity.

• If an external force is applied, it is added to the total force.

Parameters:

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

• gamma (float) – coefficient \(\gamma\) for the drag force

distance_squared_to(other_particle, lattice)

Calculates and returns the square of the distance between this and another particle using vector_to().

Parameters:

• other_particle (Particle) – the end point of the vector

• lattice (list) – lattice parameters [Lx, Ly]

Returns:

squared distance from this to the other particle

Return type:

float

distance_to(other_particle, lattice)

Calculates and returns the distance between this and another particle using vector_to().

Parameters:

• other_particle (Particle) – the end point of the vector

• lattice (list) – lattice parameters [Lx, Ly]

Returns:

distance from this to the other particle

Return type:

float

kinetic_energy()

Calculates the kinetic energy of the particle.

Returns:

kinetic energy

Return type:

float

move(dt)

Set a new position for the particle as

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

Parameters:

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

move_linearly(dt)

Set a new position for the particle as

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

Parameters:

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

save_position()

Save the current position of the particle.

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.

unit_vector_to(other_particle, lattice)

Returns the unit vector pointing from the position of this particle towards the position of other_particle using vector_to().

Parameters:

• other_particle (Particle) – the end point of the vector

• lattice (list) – lattice parameters [Lx, Ly]

Returns:

unit vector pointing from this to the other particle

Return type:

array

vector_to(other_particle, lattice)

Returns the vector pointing from the position of this particle to the position of other_particle.

Takes periodic boundary conditions into account.

Parameters:

• other_particle (Particle) – the end point of the vector

• lattice (list) – lattice parameters [Lx, Ly]

Returns:

vector pointing from this to the other particle

Return type:

array

wrap(lattice)

If the particle is outside of the simulation area, its position is shifted by a suitable multiple of lattice vectors so that the particle ends up back inside the simulation area.

Parameters:

lattice (list) – lattice parameters [Lx, Ly]

A4_2dmd.animate(particles, lattice_parameters, force_parameters)

Animates the simulation.

Parameters:

• particles (list) – list of Particle objects

• lattice_parameters (list) – list of lattice parameters

• margin (float) – each axis will be drawn in the range [-margin , lattice parameter + margin]

A4_2dmd.apply_constraints(particles, constraints)

Saves constraint information in Particle objects.

Parameters:

• particles (list) – list of Particle objects

• constraints (list) – list of constraints

A4_2dmd.calculate_average_and_error(values, start=0)

Calculates the average and standard error of mean of a sequence.

The values in the sequence are assumed to be uncorrelated.

If the beginning of the sequence cannot be used in the analysis (equilibrium has not yet been reached), one can ignore the early values by specifying a starting index.

Parameters:

• values (array) – values to analyse

• start (int) – index of the first value to be included in the analysis

A4_2dmd.calculate_forces(particles, force_parameters, lattice_parameters)

Calculate the total force acting on every atom.

Saves the result to each Particle object in particle.force.

Parameters:

• particles (list) – list of Particle objects

• force_parameters (list) – interaction parameters

• lattice_parameters (list) – lattice parameters [Lx, Ly]

Returns:

the virial

Return type:

float

A4_2dmd.calculate_kinetic_energy(particles)

Calculate the total kinetic energy of the system.

Parameters:

particles (list) – list of Particle objects

Returns:

total kinetic energy

Return type:

float

A4_2dmd.calculate_lj_potential_energy(particles, force_parameters, lattice_parameters)

Calculate the total potential energy of all Lennard-Jones interactions.

Parameters:

• particles (list) – list of Particle objects

• force_parameters (list) – interaction parameters

• lattice_parameters (list) – lattice parameters [Lx, Ly]

Returns:

total LJ potential energy

Return type:

float

A4_2dmd.calculate_momentum(particles)

Calculate the total momentum of the system.

Parameters:

particles (list) – list of Particle objects

Returns:

momentum [px, py]

Return type:

array

A4_2dmd.calculate_pressure(particles, lattice_parameters, virial)

Calculate the current pressure.

For a molecular simulation with constant pressure, volume and temperature, one can derive the relation

\[pV = Nk_B T + \frac{1}{d} \langle \sum_i \vec{F}_i \cdot \vec{r}_i \rangle,\]

where \(p, V, N, k_B, T, d, \vec{F}_i, \vec{r}_i\) are, respectively, pressure, volume, number of atoms, Boltzmann constant, temperature, number of dimensions, force acting on atom \(i\) and position of atom \(i\). The sum of the products of forces and positions is known as the virial and it should be calculated with calculate_forces().

The function uses this relation to solve the effective instantaneous pressure as

\[p = \frac{1}{V} Nk_B T + \frac{1}{dV} \sum_i \vec{F}_i \cdot \vec{r}_i.\]

This is not necessarily the true instantaneous pressure, but calculating the average of this quantity over an extended simulation should converge towards the true pressure.

Parameters:

• particles (list) – list of Particle objects

• lattice_parameters (list) – lattice parameters [Lx, Ly]

• virial (float) – the virial

Returns:

pressure

Return type:

float

A4_2dmd.calculate_spring_potential_energy(particles, force_parameters, lattice_parameters)

Calculate the total potential energy of all spring-like bonds.

Parameters:

• particles (list) – list of Particle objects

• force_parameters (list) – interaction parameters

• lattice_parameters (list) – lattice parameters [Lx, Ly]

Returns:

total spring potential energy

Return type:

float

A4_2dmd.calculate_temperature(particles)

Calculate the current instantaneous temperature.

The calculation is based on the kinetic energies of particles.

According to the equipartition principle, each quadratic degree of freedom (DOF) stores, on average, the energy

\[\langle E_\text{DOF} \rangle = \frac{1}{2} k_B T,\]

where \(k_B\) is the Boltzmann constant and \(T\) is the temperature.

In 2D, every unconstrained atom has 2 degrees of freedom for their linear movement, so the total kinetic energy of the system is, on average,

\[K = 2 N_\text{atoms} \langle E_\text{DOF} \rangle = N_\text{atoms} k_B T.\]

For simplicity, we set \(k_B = 1\), so the temperature is

\[T = \frac{1}{ N_\text{atoms}} K.\]

At the macroscopic scale, the kinetic energy of a microscopic system may be observed either as kinetic energy (movement) or internal energy (hotness). This function assumes there is no macroscopic kinetic energy so that all microscopic kinetic energy can be interpreted as internal energy. That is, this function assumes the system as a whole is at rest and the movement of particles is random.

If the system has a moving center of mass or rotation around a center, there is collective motion which is observed as macroscopic movement. In such a case this function systematically reports too high temperatures, because not all of the microscopic energy is internal energy at the macro scale.

Parameters:

particles (list) – list of Particle objects

Returns:

temperature

Return type:

float

A4_2dmd.count_connections(particles, printout=False)

Counts the number of atomic pairs connected via springs.

Parameters:

• particles (list) – list of Particle objects

• printout (bool) – if True, all connections are listed

A4_2dmd.find_connected_atoms(particles)

Builds connections between Particle objects.

Particles can be defined with a connection to another particle, which means that there is a spring-like bond connecting these two particles. As particle data is read, only the indices of connected particles are read and saved in the Particle objects. This function goes through all Particles and adds links to the Particle objects that are connected.

After this operation, each Particle has a list named connected_atoms containing the other Particle objects that are connected to it.

Connections are always reciprocal: If A is connected to B, then B is connected to A. The original particle data needs not fulfill this condition, but the function will always form connections to both connected Particles even if only one of them orginally declared the connection.

Parameters:

particles (list) – list of Particle objects

A4_2dmd.find_info(lines, tag)

Searches for the information wrapped in the given tag among the given lines of text.

If tag is, e.g., “foo”, the function searches for the start tag <foo> and the end tag </foo> and returns the lines of information between them.

The function only finds the first instance of the given tag. However, in order to catch multiple instances of the tag, the function also returns all the information in lines following the end tag.

For instance, if lines contains the strings:

aa

<foo>
bb
cc
</foo>

dd
ee
ff

the function will return two lists: [“bb”, “cc”], [“”, “dd”, “ee”, “ff”].

Parameters:

• lines (list) – the information as a list of strings

• tag (str) – the tag to search

Returns:

the lines between start and end tags, the lines following the end tag

Return type:

list, list

A4_2dmd.langevin_force(particles, dt, gamma, t_external)

Applies a random Gaussian force to all particles.

In Langevin dynamics, the Newtonian dynamics are extended by adding a drag force

\[\vec{F}_\text{drag} = - \gamma m \vec{v}\]

and a random Gaussian force \(\vec{F}_\text{random}\) with standard deviation \(\sigma = \sqrt{ 2 \gamma m T }\) and no correlation between forces at different times. This function adds such random force to all particles.

Physically, this could represent a system where the simulated particles are surrounded by an evironment of other particles that are not explicitly included in the simulation. The drag force represents flow resistance from moving through this environment (cf. air resistance). The random force represents random collisions between the simulated particles and the particles of the environment (cf. Brownian motion).

This approach also leads to correct sampling of the canonical ensemble at temperature T, so Langevin dynamics can also be used as a thermostat.

Parameters:

• particles (list) – list of Particle objects

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

• gamma (float) – coefficient \(\gamma\) for the drag force

• t_external (float) – external temperature

A4_2dmd.lj_energy(atom_A, atom_B, sigma_sixth, epsilon, lattice_parameters)

Calculate the Lennard-Jones potential energy of two atoms.

Denote the distance between the atoms as \(r\). The potential energy is calculated as:

\[U = 4 \epsilon \left( \frac{ \sigma^{12} }{ r^{12} } - \frac{ \sigma^6 }{ r^6 } \right),\]

where \(\sigma\) and \(\epsilon\) are parameters of the model.

Parameters:

• atom_A (Particle) – atom taking part in the interaction

• atom_B (Particle) – atom taking part in the interaction

• sigma_sixth (float) – parameter \(\sigma^6\)

• epsilon (float) – parameter \(\epsilon\)

• lattice_parameters (list) – lattice parameters [Lx, Ly]

Returns:

potential energy \(U\)

Return type:

float

A4_2dmd.lj_force(atom_A, atom_B, sigma_sixth, epsilon, lattice_parameters)

Calculate the Lennard-Jones force that atom B applies on atom A.

Returns the force associated with the potential energy U given by the function lj_energy():

\[\vec{F}_A = - \nabla_A U,\]

where the energy is differentiated with respect to the coordinates of atom A.

Parameters:

• atom_A (Particle) – atom taking part in the interaction

• atom_B (Particle) – atom taking part in the interaction

• sigma_sixth (float) – parameter \(\sigma^6\)

• epsilon (float) – parameter \(\epsilon\)

• lattice_parameters (list) – lattice parameters [Lx, Ly]

Returns:

force acting on atom A, [Fx, Fy]

Return type:

array

A4_2dmd.main(system_file='system.txt', simu_file='simulation.txt')

The main program.

Reads system and simulation information from files, runs a simulation, and calculates statistics.

Possibly also shows an animation of the simulation.

Parameters:

• system_file (str) – name of the file containing system info

• simu_file (str) – name of the file containing physical and simulation info

A4_2dmd.parse_line(line)

Separates tag and info on a line of text.

The function also removes extra whitespace and comments separated with #.

For instance if line is ” x : 1.23 # the x coordinate”, the function returns (“x”, “1.23”).

Parameters:

line (str) – a string of information

Returns:

tag, info

Return type:

str, str

A4_2dmd.plot_function_and_average(time, values, average, error, ylabel, show=True)

Plots a time series of a quantity as well as its average and error marginal.

Parameters:

• time (array) – time at each data point

• values (array) – recorded value at each data point

• average (float) – average of the recorded values

• error (float) – error estimate for the average

• ylabel (str) – name of the variable

• show (bool) – Of True, the plot is shown on screen, otherwise only created in memory.

A4_2dmd.print_progress(step, total)

Prints a progress bar.

Parameters:

• step (int) – progress counter

• total (int) – counter at completion

A4_2dmd.randomize_velocities(particles, temperature_parameters)

Replace the velocities of all particles with random velocities drawn from the Maxwell-Boltzmann velocity distribution.

The function makes sure the total momentum from the newly assigned velocities is exactly zero.

Parameters:

• particles (list) – list of Particle objects

• temperature_parameters (list) – temperature parameters

A4_2dmd.read_atom(lines)

Reads the properties of a single particle.

Parameters:

lines (list) – information as a list of strings

Returns:

a new Particle object created from the given information

Return type:

Particle

A4_2dmd.read_box(lines, default=10.0)

Reads lattice parameter info from given lines.

Parameters:

• lines (list) – information as a list of strings

• default (float) – the default lattice parameter in all directions

Returns:

lattice parameters [Lx, Ly]

Return type:

list

A4_2dmd.read_constraint(lines)

Reads one constraint.

Parameters:

lines (list) – information as a list of strings

Returns:

constraint info as (name, type, value)

Return type:

str, str, float

A4_2dmd.read_constraints(filename)

Read constraint settings from a file.

Parameters:

filename (str) – the file containing cosntraint information

Returns:

list of constraints

Return type:

list

A4_2dmd.read_interactions(lines)

Reads interaction parameters.

The information is returned as a list with these elements:

• params[CK_INDEX] : spring constant

• params[CR_INDEX] : spring equilibrium length

• params[CM_INDEX] : spring maximum length (no force beyond this separation)

• params[SIGMA_INDEX] : Lennard-Jones sigma parameter

• params[EPS_INDEX] : Lennard-Jones epsilon parameter

Parameters:

lines (list) – information as a list of strings

Returns:

interaction parameters

Return type:

list

A4_2dmd.read_physics(filename)

Read interaction and temperature data from a file.

Parameters:

filename (str) – the file containing physical information

Returns:

interaction parameters, temperature parameters

Return type:

list, list

A4_2dmd.read_recording(lines)

Reads data recording parameters.

The information is returned as a list with these elements:

• params[PLANE_INDEX] : animation switch: If “yes”, an animation will be drawn at the end

• params[PLOT_INDEX] : plot switch: If “yes”, temperature and pressure will be plotted at the end

• params[ANIDT_INDEX] : simulation time between animation frames

• params[RECDT_INDEX] : simulation time between recording of physical data

• params[THERMAL_INDEX] : simulation time before data collecting begins

Parameters:

lines (list) – information as a list of strings

Returns:

timing parameters

Return type:

list

A4_2dmd.read_system(filename)

Read particle and lattice data from a file.

Parameters:

filename (str) – the file containing system information

Returns:

lattice parameters, list of Particle objects

Return type:

list, list

A4_2dmd.read_temperature(lines)

Reads temperature parameters.

The information is returned as a list with these elements:

• params[T_INDEX] : external temperature

• params[TAU_INDEX] : thermostat strength (0 for no thermostat)

• params[RS_INDEX] : random start switch: If “yes”, all atoms will be given new random velocities following the Maxwell-Boltzmann distribution at the start of the simulation.

Parameters:

lines (list) – information as a list of strings

Returns:

temperature parameters

Return type:

list

A4_2dmd.read_timescale(filename)

Read simulation and recording timing data from a file.

Parameters:

filename (str) – the file containing simulation information

Returns:

simulation timescale parameters, data recording parameters

Return type:

list, list

A4_2dmd.read_timing(lines)

Reads timing parameters.

The information is returned as a list with these elements:

• params[DT_INDEX] : simulation time step

• params[RUNTIME_INDEX] : total simulation time

Parameters:

lines (list) – information as a list of strings

Returns:

timing parameters

Return type:

list

A4_2dmd.run_simulation(particles, force_parameters, lattice_parameters, time_parameters, rec_parameters, temperature_parameters)

Run a molecular dynamics simulation.

Parameters:

• particles (list) – list of Particle objects

• force_parameters (list) – interaction parameters

• lattice_parameters (list) – lattice parameters

• time_parameters (list) – timing parameters

• rec_parameters (list) – recording parameters

• temperature_parameters (list) – temperature parameters

Returns:

arrays containing measured physical quantities at different times,

(time, temperature, pressure, momentum, kinetic energy, spring energy, LJ energy)

Return type:

tuple

A4_2dmd.show_system(particles, lattice_parameters, force_parameters, margin=0.1, frame=- 1, filename='particles.png')

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

Parameters:

• particles (list) – list of Particle objects

• lattice_parameters (list) – list of lattice parameters

• force_parameters (list) – interaction parameters

• margin (float) – each axis will be drawn in the range [-margin , lattice parameter + margin]

• frame (int) – index of the frame to show - the last one by default

• filename (str) – name of the file where the system will be drawn

A4_2dmd.spring_energy(atom_A, atom_B, k, r0, rmax, lattice_parameters)

Calculate the spring potential energy of two connected atoms.

Denote the distance between the atoms as \(r\) and the maximum spring length as \(r_\max\).

If the atoms are close enough, \(r < r_\max\), the energy is

\[U = \frac{1}{2} k (r - r_0)^2,\]

where \(k\) is the spring constant and \(r_0\) is the equilibrium distance.

If the atoms are too far apart, \(r \ge r_\max\), the energy is

\[U = \frac{1}{2} k (r_\max - r_0)^2.\]

With these definitions, the potential energy has the minimum \(U(r_0) = 0\) and increases parabolically. At large separations, \(r > r_\max\), the energy is constant.

The function does not check if the particles should interact. It assumes they always do.

Parameters:

• atom_A (Particle) – atom taking part in the interaction

• atom_B (Particle) – atom taking part in the interaction

• k (float) – spring constant \(k\)

• r0 (float) – spring equilibrium length \(r_0\)

• rmax (float) – maximum spring length \(r_\max\)

• lattice_parameters (list) – lattice parameters [Lx, Ly]

Returns:

potential energy \(U\)

Return type:

float

A4_2dmd.spring_force(atom_A, atom_B, k, r0, rmax, lattice_parameters)

Calculate the spring force that atom B applies on atom A.

Returns the force associated with the potential energy given by the function spring_energy():

\[\vec{F}_A = - \nabla_A U,\]

where the energy is differentiated with respect to the coordinates of atom A.

The function does not check if the particles should interact. It assumes they always do.

Note

This function is incomplete!

Parameters:

• atom_A (Particle) – atom taking part in the interaction

• atom_B (Particle) – atom taking part in the interaction

• k (float) – spring constant \(k\)

• r0 (float) – spring equilibrium length \(r_0\)

• rmax (float) – maximum spring length \(r_\max\)

• lattice_parameters (list) – lattice parameters [Lx, Ly]

Returns:

force acting on atom A, [Fx, Fy]

Return type:

array

A4_2dmd.update_positions(particles, dt)

Update the positions of all particles according to

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

using Particle.move().

Parameters:

• particles (list) – list of Particle objects

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

A4_2dmd.update_positions_without_force(particles, dt)

Update the positions of all particles according to

\[\Delta \vec{r} = \vec{v} \Delta t\]

using Particle.move_linearly().

Parameters:

• particles (list) – list of Particle objects

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

A4_2dmd.update_velocities(particles, dt, gamma=0)

Update the velocities of all particles according to

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

If a non-zero gamma is given, a drag force \(\vec{F}_\text{drag} = - \gamma m \vec{v}\) is also applied.

using Particle.accelerate().

Parameters:

• particles (list) – list of Particle objects

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

• gamma (float) – coefficient \(\gamma\) for the drag force

A4_2dmd.velocity_verlet(particles, force_parameters, lattice_parameters, time_parameters, rec_parameters, temperature_parameters)

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

The algorithm works as follows:

• First, forces are calculated at current time, \(\vec{F}(t)\).

• Second, velocities are calculated half a time step in the future, \(\vec{v}(t + \frac{1}{2}\Delta t) = \vec{v}(t) + \frac{1}{m}\vec{F}(t) \frac{1}{2}\Delta t\).

• Then, the following steps are repeated for as long as the simulation runs,

  • Positions are updated by one time step using the velocities, \(\vec{r}(t + \Delta t) = \vec{r}(t) + \vec{v}(t + \frac{1}{2}\Delta t) \Delta t\).

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

  • Velocities are updates by one 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\)

Note that the algorithm uses velocities “half a time step from the future” to update positions and forces “half a time step from the future” to update the velocities. This approach effectively averages upcoming and previous values and leads to a stable algorithm that is symmetric with respect to time reversal.

Parameters:

• particles (list) – list of Particle objects

• force_parameters (list) – interaction parameters

• lattice_parameters (list) – lattice parameters

• time_parameters (list) – timing parameters

• rec_parameters (list) – recording parameters

• temperature_parameters (list) – temperature parameters

Returns:

average temperature, pressure

Return type:

float, float

A4_2dmd.write_system_file(particles, lattice_parameters, filename='new_system.txt')

Write system information in a file.

Parameters:

• particles (list) – list of Particle objects

• lattice_parameters (list) – lattice parameters

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