Native backends

Let us have a closer look at backends.

Interactions

There is a native backend in atooms.backends.f90 to compute interactions between particles. It is written in Fortran 90 and compiled at run-time using f2py-jit. For demonstration purposes, the core package provides a simple Lennard-Jones interaction potential.

As a simple example, we compute the interaction energy of three particles interacting with the Lennard-Jones potential.

from atooms.system import Particle, Cell, System
from atooms.backends.f90 import Interaction

model = {
    "potential": [{
        "type": "lennard_jones",
        "parameters": {"epsilon": [[1.0]], "sigma": [[1.0]]}
    }],
    "cutoff": [{
        "type": "cut_shift",
    "parameters": {"rcut": [[2.5]]}
    }]
}

particles = [Particle(position=[0.0, 0.0, 0.0], species=1),
             Particle(position=[1.0, 0.0, 0.0], species=1),
             Particle(position=[2.0, 0.0, 0.0], species=1)]
cell = Cell([10., 10., 10.])
system = System(particles, cell)
system.interaction = Interaction(model)
print(system.potential_energy())  # reference: -0.01257276409199999

-0.012572764091999991

The model dictionary specify the Lennard-Jones potential and its parameters, along with a cutoff to truncate and shift the potential at \(r_c=2.5\) in the units of \(\sigma\).

Additional interaction models are defined by the companion atooms.models package, which is just a pip install atooms.models away. They can be used to perform simulations with the dynamics package to which we turn now.

Dynamics

The atooms.dynamics package implements Newtonian and stochastic dynamics backends for atooms.

Currently, implemented integration algorithms

• Netwonian dynamics

  • velocity-Verlet

  • Nose-Poincaré

  • event-driven

• Stochastic dynamics

  • overdamped Langevin dynamics

To run a molecular dynamics simulation of a Lennard-Jones system from an existing xyz file

from atooms.trajectory import Trajectory
from atooms.simulation import Simulation
from atooms.dynamics.netwonian import VelocityVerlet

# Start from the last frame of input.xyz
trajectory = Trajectory('input.xyz')
system = trajectory[-1]
system.interaction = Interaction('lennard_jones')
backend = VelocityVerlet(system, timestep=0.002)
sim = Simulation(backend, steps=200)
sim.run()

We can do the same via the Python API, which wraps a typical template of a simulation. We store the configurations in output.xyz

from atooms.dynamics.api import md

md('input.xyz', 'output.xyz',
   method='velocity-verlet', model='lennard_jones',
   dt=0.002, nsteps=200, config_number=20)

The same can be achieved from the command line

md.py --method velocity-verlet -n 200 --dt 0.002 --config-number 20 input.xyz output.xyz

Active

There is a separate backend for simulation of active matter, developed by Iacopo Ricci.

This is an example of simulation of the Vicsek model in 2d. We first set the system up.

import numpy
import atooms
from atooms.system import System

dims=2
system = System(N=200, d=dims)
system.density = 0.5
for p in system.particle:
    p.orientation = numpy.random.uniform(-numpy.pi,numpy.pi)
    p.position = [L*numpy.random.uniform() for i in range(dims)]
    p.fold(system.cell)

Simulate the system with the Vicsek backend.

from atooms.active.vicsek import Vicsek
from atooms.active.neighbors import VicsekNeighbors
from atooms.simulation import Simulation

backend = Vicsek(system, eta=0.4, v0=0.5, noise='vectorial')
neighbors = VicsekNeighbors(system, method='kdtree')
bck.neighbors = neighbors

Simulation(backend).run(10)

This simulation employs the so-called ‘vectorial’ noise implementation (discussed in Grégoire and Chaté) and the scipy-based kD-tree neighbor search algorithm.