Use the force

  • Topic: atomistic models, pair potentials, conservative forces

  • Task:

    1. Implement the method that calculates forces on all atoms.

    2. Run the program and compare your forces against values the program calculates by numerically differentiating the potential energy. If your function works, all values should be a close match. If there are any notable discrepancies, there is a bug.

    3. Change the \(\Delta\) parameter and study how it affects numeric accuracy. Does the numeric result approach your analytic solution?

    4. Make the program create random configurations and check that you always get matching results.

  • Template: forces.py

  • Data: atoms.txt

  • Further reading:

forces.py

class forces.Atom(position, velocity, mass=1.0)[source]

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\)

move(shift)[source]

Move the atom.

Parameters:

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

forces.calculate_forces(atoms, parameters)[source]

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 ],
  ...               ]

Note

This function is incomplete!

Parameters:

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

  • parameters (list) – a list of physical parameters

Returns:

forces

Return type:

array

forces.calculate_potential_energy(atoms, parameters)[source]

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:

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

  • parameters (list) – a list of physical parameters

Returns:

potential energy

Return type:

array

forces.create_random_configuration(n=10)[source]

Randomly creates a set of Atom objects.

The atoms are placed uniformly in a \([-3, 3] \times [-3, 3] \times [-3, 3]\) box. They are also given random velocities and masses.

Parameters:

n (int) – number of atoms

Returns:

the atoms

Return type:

list

forces.distance_squared(atom1, atom2)[source]

Calculates the squared distance between two atoms,

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

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

Return type:

float

forces.draw_particles_and_forces(atoms, forces_A, forces_B=None)[source]

Draws the particles and forces as dots and arrows in the xy plane.

The function accepts two sets of forces and draws them both so that the forces can be compared to each other.

Parameters:

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

  • forces_A (array) – forces on all atoms

  • forces_B (array) – forces on all atoms

forces.energy_gradient(atoms, alpha, delta, parameters)[source]

Calculates a numeric approximation for the gradient of the potential energy with respect to the position of one atom.

The potential energy of the system is defined via calculate_potential_energy(), and it is a function of the coordinates of all atoms,

\[U = U(x_0, y_0, z_0, x_1, y_1, z_1, \ldots).\]

The force of a conservative interaction is given by the gradient of the potential energy. More precisely, the force on atom \(\alpha\) is the gradient with respect to the coordinates of that particular atom,

\[\vec{F}_\alpha = - \nabla_\alpha U = \frac{ \partial U }{\partial x_{\alpha} } \hat{i} + \frac{ \partial U }{\partial y_{\alpha} } \hat{j} + \frac{ \partial U }{\partial z_{\alpha} } \hat{k}.\]

This function calculates this vector numerically. Each partial derivative is calculated using the symmetric finite difference result

\[\frac{ \partial U }{\partial x } = \frac{ U(x+\Delta) - U(x-\Delta) }{2 \Delta}.\]

In practice, the function shifts atom \(\alpha\) and calls calculate_potential_energy() to obtain the needed energies.

Parameters:

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

  • alpha (int) – index of the atom with respect to which the gradient is calculated

  • delta (float) – the shift \(\Delta\) used for calculating finite differences

Returns:

force components \([ F_{x,\alpha}, F_{y,\alpha}, F_{z,\alpha} ]\)

Return type:

array

forces.main(delta=0.001, random_atoms=0)[source]

The main program.

Reads the particle data from the file atoms.txt or creates a random configuration and calculates the forces acting on all particles both directly with calculate_forces() and by numerically differentiating the potential energy with energy_gradient(). The method reports these results as well as the error

\[\Delta \vec{F} = \vec{F}_\text{analytic} - \vec{F}_\text{numeric},\]

its absolute value \(|\Delta \vec{F}|\) and the error relative to the magnitude of the force, \(|\Delta \vec{F}| / |\vec{F}_\text{numeric}|\).

If the force calculation is correct, the results should be approximately the same and the error should be close to zero. Due to numeric inaccuracy, one expects a difference which approaches zero as \(\Delta\) goes to zero.

Parameters:

  • delta (float) – the shift \(\Delta\) used for calculating finite differences

  • random_atoms (int) – If zero, configuration is read form file. If positive, a random system with that many atoms is created.

forces.read_atoms_from_file(filename)[source]

Reads the properties of atoms from a file.

Each line should define a single forces.Atom 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 forces.Atom objects

Return type:

list

forces.vector(atom1, atom2)[source]

The vector pointing from one atom, \(i\), to another, \(j\),

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

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

Return type:

array