import sys
import copy
from math import *
import numpy as np
import numpy.random as random
import matplotlib.pyplot as plt
[docs]class Atom:
"""
A point like object.
An atom has a position (a 3-vector), a velocity (3-vector)
and a mass (a scalar).
Args:
position (array): coordinates :math:`[x, y, z]`
velocity (array): velocity components :math:`[v_x, v_y, v_z]`
mass (float): mass :math:`m`
"""
def __init__(self, position, velocity, mass = 1.0):
self.position = np.array(position)
self.velocity = np.array(velocity)
self.mass = mass
self.trajectory = []
[docs] def move(self, shift):
"""
Move the atom.
Args:
shift (array): coordinate change :math:`[\Delta x, \Delta y, \Delta z]`
"""
self.position += shift
[docs]def read_atoms_from_file(filename):
"""
Reads the properties of atoms from a file.
Each line should define a single :class:`forces.Atom` listing its
position in cartesian coordinates, velocity components and mass,
separated by whitespace:
.. code-block::
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
...
Args:
filename (str): name of the file to read
Returns:
list: list of :class:`forces.Atom` objects
"""
f = open(filename)
lines = f.readlines()
f.close()
atoms = []
for line in lines:
parts = line.split()
if len(parts) > 0:
position = np.array( [0.0]*3 )
velocity = np.array( [0.0]*3 )
position[0] = float(parts[0])
position[1] = float(parts[1])
position[2] = float(parts[2])
velocity[0] = float(parts[3])
velocity[1] = float(parts[4])
velocity[2] = float(parts[5])
mass = float(parts[6])
atoms.append( Atom( position, velocity, mass ) )
return atoms
[docs]def distance_squared(atom1, atom2):
"""
Calculates the squared distance between two atoms,
.. math ::
r^2_{ij} = (x_i-x_j)^2 + (y_i-y_j)^2 + (z_i-z_j)^2.
Args:
atom1 (forces.Atom): the first atom
atom2 (forces.Atom): the second atom
Returns:
float: the squared distance :math:`r^2_{ij}`
"""
pos1 = atom1.position
pos2 = atom2.position
dx = pos1[0] - pos2[0]
dy = pos1[1] - pos2[1]
dz = pos1[2] - pos2[2]
return dx*dx + dy*dy + dz*dz
[docs]def vector(atom1, atom2):
"""
The vector pointing from one atom, :math:`i`, to another, :math:`j`,
.. math ::
\\vec{r}_{i \\to j} = \\vec{r}_j - \\vec{r}_i
Args:
atom1 (forces.Atom): the first atom
atom2 (forces.Atom): the second atom
Returns:
array: components of :math:`\\vec{r}_{i \\to j}`, :math:`[x_{i \\to j}, y_{i \\to j}, z_{i \\to j}]`
"""
pos1 = atom1.position
pos2 = atom2.position
return pos2 - pos1
[docs]def calculate_forces(atoms, parameters):
"""
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.
.. code-block::
[ [ fx0, fy0, fz0 ],
[ fx1, fy1, fz1 ],
[ fx2, fy2, fz2 ],
... ]
.. note ::
This function is incomplete!
Args:
atoms (list): a list of :class:`forces.Atom` objects
parameters (list): a list of physical parameters
Returns:
array: forces
"""
# parameters of the LJ potential
sigma = parameters[0]
epsilon = parameters[1]
forces = np.array( [[0.0, 0.0, 0.0]]*len(atoms) )
# Loop over all pairs of atoms:
# - i gets all values between 0 and N-1
# - j gets the values between 0 and i-1 for each i
# This way we get e.g. (i,j) = (1,0) but not (0,1),
# so the pair of atoms 0 and 1 is found only once.
for i in range(len(atoms)):
for j in range(0,i):
part_i = atoms[i]
part_j = atoms[j]
dist_sq = distance_squared(part_i, part_j)
dist = sqrt(dist_sq)
# todo
force_j_to_i = np.zeros(3)
# pairwise forces obey the law of action and reaction
forces[i, :] += force_j_to_i[:] # add force on atom i
forces[j, :] -= force_j_to_i[:] # add inverse force on atom j
return forces
[docs]def calculate_potential_energy(atoms, parameters):
"""
Calculates the total potential energy of the system.
The potential energy is calculated using the Lennard-Jones model
.. math ::
U = \\sum_{i \\ne j} 4 \\epsilon \\left( \\frac{ \\sigma^{12} }{ r^{12}_{ij} }
- \\frac{ \\sigma^6 }{ r^6_{ij} } \\right).
Args:
atoms (list): a list of :class:`forces.Atom` objects
parameters (list): a list of physical parameters
Returns:
array: potential energy
"""
sigma = parameters[0]
epsilon = parameters[1]
sigma_six = sigma*sigma*sigma*sigma*sigma*sigma
u = 0.0
for i in range(len(atoms)):
for j in range(0,i):
part_i = atoms[i]
part_j = atoms[j]
dist_sq = distance_squared(part_i, part_j)
dist_six = dist_sq * dist_sq * dist_sq
dist_12 = dist_six * dist_six
dist = sqrt(dist_sq)
u += 4.0 * epsilon * \
( sigma_six*sigma_six/dist_12 - sigma_six/dist_six)
return u
[docs]def energy_gradient(atoms, alpha, delta, parameters):
"""
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
:meth:`calculate_potential_energy`, and it is a function
of the coordinates of all atoms,
.. math ::
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
:math:`\\alpha` is the gradient *with respect to the coordinates
of that particular atom*,
.. math ::
\\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
.. math ::
\\frac{ \\partial U }{\\partial x } = \\frac{ U(x+\\Delta) - U(x-\\Delta) }{2 \\Delta}.
In practice, the function shifts atom :math:`\\alpha`
and calls :meth:`calculate_potential_energy` to obtain the
needed energies.
Args:
atoms (list): a list of :class:`forces.Atom` objects
alpha (int): index of the atom with respect to which the gradient is calculated
delta (float): the shift :math:`\\Delta` used for calculating finite differences
Returns:
array: force components :math:`[ F_{x,\\alpha}, F_{y,\\alpha}, F_{z,\\alpha} ]`
"""
deltax = np.array([delta,0,0])
deltay = np.array([0,delta,0])
deltaz = np.array([0,0,delta])
atoms[alpha].move(deltax)
e_xplus = calculate_potential_energy(atoms, parameters)
atoms[alpha].move(-2*deltax)
e_xminus = calculate_potential_energy(atoms, parameters)
atoms[alpha].move(deltax)
atoms[alpha].move(deltay)
e_yplus = calculate_potential_energy(atoms, parameters)
atoms[alpha].move(-2*deltay)
e_yminus = calculate_potential_energy(atoms, parameters)
atoms[alpha].move(deltay)
atoms[alpha].move(deltaz)
e_zplus = calculate_potential_energy(atoms, parameters)
atoms[alpha].move(-2*deltaz)
e_zminus = calculate_potential_energy(atoms, parameters)
atoms[alpha].move(deltaz)
return np.array( [ (e_xplus - e_xminus)/(2*delta),
(e_yplus - e_yminus)/(2*delta),
(e_zplus - e_zminus)/(2*delta) ] )
[docs]def draw_particles_and_forces(atoms, forces_A, forces_B = None):
"""
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.
Args:
atoms (list): a list of :class:`forces.Atom` objects
forces_A (array): forces on all atoms
forces_B (array): forces on all atoms
"""
plt.clf()
ax = plt.axes()
ax.set_xlim([-3,3])
ax.set_ylim([-3,3])
ax.set_aspect('equal')
xs = []
ys = []
ms = []
# plot the atoms
for atom in atoms:
xs.append(atom.position[0])
ys.append(atom.position[1])
ms.append(atom.mass*40)
plt.scatter( xs, ys, marker='o', s=ms )
# plot the forces
# scale so that the longest arrow has
# approximately unit length
fscale = 1 / np.max(forces_A)
for i in range(len(atoms)):
x0 = xs[i]
y0 = ys[i]
fxA = fscale*forces_A[i,0]
fyA = fscale*forces_A[i,1]
plt.arrow( x0, y0, fxA, fyA, color='r', lw = 2 )
if forces_B is not None:
fxB = fscale*forces_B[i,0]
fyB = fscale*forces_B[i,1]
plt.arrow( x0, y0, fxB, fyB, color='b', lw = 2 )
plt.show()
[docs]def create_random_configuration(n = 10):
"""
Randomly creates a set of :class:`Atom` objects.
The atoms are placed uniformly in a :math:`[-3, 3] \\times [-3, 3] \\times [-3, 3]`
box. They are also given random velocities and masses.
Args:
n (int): number of atoms
Returns:
list: the atoms
"""
rng = random.default_rng()
atoms = []
for i in range(n):
new_atom = Atom( 3.0*(2*rng.random(size=3)-1), rng.standard_normal(size=3), rng.integers(1,11) )
atoms.append( new_atom )
return atoms
[docs]def main(delta = 0.001, random_atoms = 0):
"""
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
:meth:`calculate_forces` and
by numerically differentiating the potential energy
with :meth:`energy_gradient`.
The method reports these results as well as the error
.. math::
\\Delta \\vec{F} = \\vec{F}_\\text{analytic} - \\vec{F}_\\text{numeric},
its absolute value :math:`|\\Delta \\vec{F}|` and the error relative to the
magnitude of the force, :math:`|\\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 :math:`\\Delta` goes to zero.
Args:
delta (float): the shift :math:`\\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.
"""
if random_atoms > 0:
atoms = create_random_configuration(random_atoms)
else:
atoms = read_atoms_from_file('atoms.txt')
n_atoms = len(atoms)
parameters = [1.0, 0.1] # sigma, epsilon
calculated_forces = calculate_forces(atoms, parameters)
numeric_forces = []
for i in range(n_atoms):
finite_difference = -energy_gradient(atoms, i, delta, parameters)
numeric_forces.append( finite_difference )
print( "force on atom ", i )
print( " direct calculation ", calculated_forces[i] )
print( " numeric differentiation ", finite_difference )
error = calculated_forces[i]-finite_difference
print( " difference ", error )
print( " difference, abs value ", np.sqrt( error @ error ) )
print( " difference, relative ", np.sqrt( (error @ error) / (finite_difference @ finite_difference) ) )
numeric_forces = np.array(numeric_forces)
draw_particles_and_forces(atoms, numeric_forces, calculated_forces)
if __name__ == "__main__":
# This command suppresses scientific output of numbers
# so that you get 0.000001 instead of 1.0E-06 for small errors.
# You can comment it out if you prefer scientific notation.
np.set_printoptions(suppress=True)
# Change delta to tune numeric differentiation.
# Set random_atoms > 0 to create a random system.
main(delta = 0.0001, random_atoms = 0)