Basics¶
Atooms provides a high-level interface to the main objects of particle-based simulations. It mostly focuses on classical molecular dynamics and Monte Carlo simulations, but it is not limited to that. It can be used to simulate and analyze lattice models such as TASEP or kinetically constrained models.
We will start by having a look at the basic objects of particle-based simulations and how to store them on a file.
Particles’ properties¶
Particles’ positions are stored as numpy arrays, but we can pass a simple list with x, y, z coordinates when we create them
from atooms.system.particle import Particle
particle = Particle(position=[1.0, 0.0, 0.0])
print(particle.position, type(particle.position))
[1. 0. 0.] <class 'numpy.ndarray'>
Particles can live in an arbitrary number of spatial dimensions
particle = Particle(position=[1.0, 0.0, 0.0, 0.0, 0.0])
print(len(particle.position))
5
By default, particles have a few more properties such as velocity, chemical species, mass and radius. They can all be altered at will or even set to None.
import numpy
particle = Particle(position=[1.0, 0.0, 0.0], velocity=[1.0, 0.0, 0.0])
particle.species = 'Na'
particle.position += numpy.array([0.0, 1.0, 1.0])
particle.velocity *= 2
particle.radius = None # point particles have no radius
print(particle)
Particle(species=Na, mass=1.0, position=[1. 1. 1.], velocity=[2. 0. 0.], radius=None)
You may want to add physical properties to particles, like charge or whatever. Of course, in python you can do it very easily
particle.charge = -1.0
This won’t break anything!
Dealing with velocities¶
You may not need velocities at all (for instance because you are working with Monte Carlo simulations) but if you do, atooms provides a few useful methods and functions. For instance, you can assign velocity from a Maxwell-Boltzmann distribution at a temperature T.
particle = [Particle() for i in range(1000)]
for p in particle:
p.maxwellian(T=1.0)
ekin = sum([p.kinetic_energy for p in particle])
ndim = 3
ndof = len(particle) * ndim
T = 2.0 / ndof * ekin
print(T)
0.9978434233794155
Doing so will leave a non-zero total momentum, but we can fix it (note that all masses are equal)
from atooms.system.particle import fix_total_momentum, cm_velocity
print(cm_velocity(particle))
fix_total_momentum(particle)
print(cm_velocity(particle))
[ 0.00927695 0.0473026 -0.02404831]
[3.39450690e-17 5.55111512e-17 3.35287353e-17]
Boundary conditions¶
To avoid major finite size effects, we enclose particles in a cell with periodic boundary conditions. By convention, the cell origin is at the origin of the reference frame.
from atooms.system.cell import Cell
L = 2.0
cell = Cell(side=[L, L, L])
print(cell.side, cell.volume)
[2. 2. 2.] 8.0
Atooms provides means to fold particles back in the “central” simulation cell, i.e. the one centered at the origin at the reference frame. For simplicity, let us work with particles in 1d.
cell = Cell(side=[1.0])
particle = Particle(position=[2.0]) # particle outside the central cell
particle.fold(cell)
print(particle.position)
[0.]
The particle is now folded back at the origin.
A related method returns the nearest periodic image of a given particle with respect to another particle
particle_1 = Particle(position=[-0.45])
particle_2 = Particle(position=[+0.45])
image = particle_1.nearest_image(particle_2, cell, copy=True)
print(image)
Particle(species=A, mass=1.0, position=[0.55], velocity=[0. 0. 0.], radius=0.5)
The System object¶
Objects like particles and the simulation cell can be gathered in an instance of a god-like class called System. The system contains all the relevant physical objects of your simulation. Reservoirs like thermostats, barostats and particle reservoirs can be added as well. These objects are placeholders for thermodynamic state variables like temperature, pressure or chemical potential. Any class meant to describe the interaction between particles also belongs to the system.
Let us build a system with a few particles in a cell and use the system methods to modify the system density and temperature. Note that density and temperature are python properties and thus modify the attributes of particles and cell under the hoods using the set_density and set_temperature methods respectively
from atooms.system import System
system = System(particle=[Particle() for i in range(100)],
cell=Cell([10.0, 10.0, 10.0]))
system.density = 1.2 # equivalent to system.set_density(1.2)
system.temperature = 1.5 # equivalent to system.set_temperature(1.2)
print(system.density, system.temperature)
1.1999999999999997 1.4999999999999998
Note that the system temperature is the kinetic one and need not coincide with the one of the thermostat.
from atooms.system import Thermostat
system.thermostat = Thermostat(temperature=1.0)
system.temperature = 1.5 # equivalent to system.set_temperature(1.2)
print(round(system.temperature, 2), system.thermostat.temperature)
1.5 1.0
Interaction and backends¶
Classical particles interact with each other via a potential \(u(\{r_i\})\), where \(\{r_i\}\) is the set of particles’ coordinates. Atooms relies on third-party efficient backends written in C, Fortran or CUDA to actually compute the interaction between the particles. Here we will use the LAMMPS backend, see Molecular dynamics ith LAMMPS for further details. It accepts a string variable that defines the interaction potential using the LAMMPS syntax, see https://lammps.sandia.gov/doc/pair_style.html, and stores a reference to the system object of which we want to compute the energy.
As proof of principle, we compute the interaction energy between two Lennard-Jones particles
from atooms.system import System, Particle, Cell
from atooms.backends import lammps
lammps.lammps_command = 'lmp'
x = 1.122 # Minimum of the potential
system = System(particle=[Particle(position=[0.0, 0.0, 0.0]),
Particle(position=[x, 0.0, 0.0])],
cell=Cell([10.0, 10.0, 10.0]))
cmd = """
pair_style lj/cut 2.5
pair_coeff 1 1 1.0 1.0 2.5
"""
# The backend will add an interaction to the system
backend = lammps.LAMMPS(system, cmd)
# Compute and get the potential energy
# The cache option allows to get the potential energy without recalculating it
print(system.potential_energy(), system.potential_energy(cache=True))
-0.99999388 -0.99999388
The energy and forces are stored in system.interaction.energy and system.interaction.forces.
Trajectory files¶
To write the state of the system to a file, we use a Trajectory class. Trajectories are composed of multiple frames, each one holding the state of the system at a given step during the simulation. We use a basic xyz format to write the state of the system and then parse the trajectory file we produced to see how it looks like.
from atooms.trajectory import TrajectoryXYZ
system = System(N=4)
system.cell = Cell([10.0, 10.0, 10.0])
# Open the trajectory in write mode and write the state of the system
# at step 0
with TrajectoryXYZ('test.xyz', 'w') as th:
th.write(system, step=0)
# Read the xyz file back as plain text
with open('test.xyz') as fh:
print(fh.read())
4
step:0 columns:species,position dt:1 cell:10.0,10.0,10.0
A -0.250000 0.250000 -0.250000
A -0.250000 0.250000 0.250000
A 0.250000 -0.250000 0.250000
A 0.250000 0.250000 0.250000
Note that trajectories are file-like objects: they must be opened and closed, preferably using the with syntax.
We can write multiple frames by calling write() repeatedly.
with TrajectoryXYZ('test.xyz', 'w') as th:
for i in range(3):
th.write(system, step=i*10)
To get the system back we read the trajectory. Trajectories support iteration and indexing, just like lists.
with TrajectoryXYZ('test.xyz') as th:
# First frame
system = th[0]
print(system.particle[0].position, system.cell.side)
# Last frame
system = th[-1]
print(system.particle[0].position, system.cell.side)
# Iterate over all frames
for i, system in enumerate(th):
print(th.steps[i], system.particle[0].position)
[-0.25 0.25 -0.25] [10. 10. 10.]
[-0.25 0.25 -0.25] [10. 10. 10.]
0 [-0.25 0.25 -0.25]
10 [-0.25 0.25 -0.25]
20 [-0.25 0.25 -0.25]
Particles on a lattice¶
Suppose we want to simulate a system where particles can only be located at discrete sites, say a one-dimensional lattice or perhaps a network with a complex topology. Particle positions can then be described as plain integers, holding the index of the site on which a particle is located. We create such a system and then write it to a file in xyz format
import numpy
from atooms.system import System, Particle
# Build model system with integer coordinates
particle = [Particle() for i in range(3)]
particle[0].position = 0
particle[1].position = 1
particle[2].position = 2
system = System(particle=particle)
# Write xyz trajectory
from atooms.trajectory import TrajectoryXYZ
with TrajectoryXYZ('test.xyz', 'w') as th:
th.write(system, 0)
# Read the xyz file back as plain text
with open('test.xyz') as fh:
print(fh.read())
3
step:0 columns:species,position dt:1
A 0
A 1
A 2
Everything went fine. However, we have to tweak things a bit when reading the particles back, to avoid positions being transformed to arrays of floats instead of integers. This can be done with the help of a callback that transforms the system accordingly as we read the trajectory.
# Read file as an xyz trajectory
with TrajectoryXYZ('test.xyz') as th:
# We add a callback to read positions as simple integers
# Otherwise they are read as numpy arrays of floats.
def modify(system):
for p in system.particle:
p.position = int(p.position[0])
p.velocity = None
p.radius = None
return system
th.add_callback(modify)
for p in th[0].particle:
print(p)
Particle(species=A, mass=1.0, position=0, velocity=None, radius=None)
Particle(species=A, mass=1.0, position=1, velocity=None, radius=None)
Particle(species=A, mass=1.0, position=2, velocity=None, radius=None)
Our particles have now integer coordinates. Note that, on passing, we have set to None velocities and radii as they are not relevant in this case.