Jumos: customisable molecular simulation package¶

System setup¶

The system will contain all the physical informations about what we are simulating. In Jumos, a system is represented by the Universe type. It should at least contain data about the simulation cell, the system topology, the initial particles coordinates and the interactions between the particles.

Simulation setup¶

During the simulation, the system will be propagated by a propagator: MolecularDynamics, MonteCarlo, Minimization are examples of propagators. Other algorithms can take place here, either algorithms to compute thermodynamic properties or algorithms to modify the behaviour of the simulation.

Running the simulation¶

The simulation run is the main part of the simulation, and consist in three main steps:

• Use the propagator to genrerate an update in the particles positions;
• Compute thermodynamical properties of the system;
• Output data to file for later analysis.

Lennard-Jones fluid¶

Here is a simple simulation script for running a simulation of a Lennard-Jones fluid at \(300K\).

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

#!/usr/bin/env julia

# Loading the Jumos module before anything else
using Jumos

# Molecular Dynamics with 1.0fs timestep
sim = Simulation(:MD, 1.0)

# Create a cubic cell with a width of 10A
cell = UnitCell(unit_from(10.0, "A"))
# Read a topology from a file
topology = Topology("lennard-jones.xyz")
# Create an universe from the cell and the topology

universe = Universe(cell, topology)
positions_from_file!(universe, "lennard-jones.xyz")
# Initialize random velocities at 300K
create_velocities!(universe, 300)

# Add Lennard-Jones interactions between He atoms
add_interaction!(
    universe,
    LennardJones(unit_from(0.2, "kJ/mol"), unit_from(2.0, "A")),
    "He", "He"
)

# You can either bind algorithms to variables ...
out_trajectory = TrajectoryOutput("LJ-trajectory.xyz", 1)
push!(sim, out_trajectory)
# ... or create them directly in function call
push!(sim, EnergyOutput("LJ-energy.dat", 10))

propagate!(sim, universe, 5000)

# Simulation scripts are normal Julia scripts !
println("All done")

Each simulation script should start by the using Jumos directive. This imports the module and the exported names in the current scope.

Then, in this script, we create two main objects: a simulation, and an universe. The topology and the original positions for the universe are read from the same file, as an .xyz file contains topological information (mainly the atomics names) and coordinates. The velocities are then initialized from a Boltzmann distribution at 300K.

The only interaction — a Lennard-Jones interaction — is also added to the universe before the run. The next lines add some outputs to the simulation, namely a trajectory and an energy output. Finally, the simulation runs for 5000 steps.

Other example¶

Some other examples scripts can be found in the example folder in Jumos’ source tree. To go there, use the julia prompt:

julia> Pkg.dir("Jumos")
"~/.julia/v0.4/Jumos"

julia>; cd $ans/examples
~/.julia/v0.4/Jumos/examples

Universe¶

The Universe type contains data. It is built around four other basic types:

• The Topology contains data about the atomic organisation, i.e. the particles, bonds, angles and dihedral angles in the system;
• The UnitCell contains data about the bounding box of the simulation;
• The Interactions type contains data about the potentials to use for the atoms in the system;
• The Frame type contains raw data about the positions and maybe velocities of the particles in the system;

Simulation¶

The Simulation type contains algorithms. The progagator algorithm is the one responsible for propagating the Universe along the simulation. It also contains some analysis algorithms (called compute in Jumos); and some output algorithms, to save data during the simulation run.

Creating simulation cell¶

UnitCell(A[, B, C, alpha, beta, gamma, celltype])¶

Creates an unit cell. If no celltype parameter is given, this function tries to guess the cell type using the following behavior: if all the angles are equals to \(\pi/2\), then the cell is an OrthorombicCell; else, it is a TriclinicCell.

If no value is given for alpha, beta, gamma, they are set to \(\pi/2\). If no value is given for B, C, they are set to be equal to A. This creates a cubic cell. If no value is given for A, a cell with lenghts of 0 Angström and \(\pi/2\) angles is constructed.

julia> UnitCell() # Without parameters
OrthorombicCell
    Lenghts: 0.0, 0.0, 0.0

julia> UnitCell(10.) # With one lenght
OrthorombicCell
    Lenghts: 10.0, 10.0, 10.0

julia> UnitCell(10., 12, 15) # With three lenghts
OrthorombicCell
    Lenghts: 10.0, 12.0, 15.0

julia> UnitCell(10, 10, 10, pi/2, pi/3, pi/5) # With lenghts and angles
TriclinicCell
    Lenghts: 10.0, 10.0, 10.0
    Angles: 1.5707963267948966, 1.0471975511965976, 0.6283185307179586

julia> UnitCell(InfiniteCell) # With type
InfiniteCell

julia> UnitCell(10., 12, 15, TriclinicCell) # with lenghts and type
TriclinicCell
    Lenghts: 10.0, 12.0, 15.0
    Angles: 1.5707963267948966, 1.5707963267948966, 1.5707963267948966

UnitCell(u::Vector[, v::Vector, celltype])

If the size matches, this function expands the vectors and returns the corresponding cell.

julia> u = [10, 20, 30]
3-element Array{Int64,1}:
 10
 20
 30

julia> UnitCell(u)
OrthorombicCell
    Lenghts: 10.0, 20.0, 30.0

Indexing simulation cell¶

You can access the cell size and angles either directly, or by integer indexing.

getindex(b::UnitCell, i::Int)¶

Calling b[i] will return the corresponding length or angle : for i in [1:3], you get the i^th lenght, and for i in [4:6], you get [avoid get] the angles.

In case of intense use of such indexing, direct field access should be more efficient. The internal fields of a cell are : the three lenghts a, b, c, and the three angles alpha, beta, gamma.

Boundary conditions and cells¶

Only fully periodic boundary conditions are implemented for now. This means that if a particle crosses the boundary at some step, it will be wrapped up and will appear at the opposite boundary.

Distances and cells¶

The distance between two particles depends on the cell type. In all cases, the minimal image convention is used: the distance between two particles is the minimal distance between all the images of theses particles. This is explicited in the Periodic boundary conditions and distances computations part of this documentation.

Atom type¶

An Atom instance is a representation of a type of particle in the system. It is paramerized by an AtomType, which can be one of:

• Element : an element in the perdiodic classification;
• DummyAtom : a dummy atom, for particles without interactions;
• CorseGrain : corse-grain particle, for corse-grain simulations;
• UnknownAtom : All the other kind of particles.

You can access the following fields for all atoms:

• label::Symbol : the atom name;
• mass : the atom mass;
• charge : the atom charge;
• properties::Dict{String, Any} : any other property: dipolar moment, etc.;

Atom([label, type])¶

Creates an atom with the label label. If type is provided, it is used as the Atom type. Else, a type is guessed according to the following procedure: if the label is in the periodic classification, the the atom is an Element. Else, it is a CorseGrain atom.

julia> Atom()
"Atom{UnknownAtom}"

julia> Atom("He")
"Atom{Element} He"

julia> Atom("CH4")
"Atom{CorseGrain} CH4"

julia> Atom("Zn", DummyAtom)
"Atom{DummyAtom} Zn"

mass(label::Symbol)¶

Try to guess the mass associated with an element, from the periodic table data. If no value could be found, the 0.0 value is returned.

mass(atom::Atom)

Return the atomic mass if it was set, or call the function mass(atom.label).

Topology¶

A Topology instance stores all the information about the system : atomic types, atomic composition of the system, bonds, angles, dihedral agnles and molecules.

Topology([natoms = 0])¶

Create an empty topology with space for natoms atoms.

Atoms in the system can be accesed using integer indexing. The following example shows a few operations available on atoms:

# topology is a Topology with 10 atoms

atom = topology[3]  # Get a specific atom
println(atom.label) # Get the atom label

atom.mass = 42.9    # Set the atom mass
topology[5] = atom  # Set the 5th atom of the topology

Periodic table information¶

The Universes module also exports two dictonaries that store information about atoms:

• ATOMIC_MASSES is a Dict{String, Float64} associating atoms symbols and atomic masses, in internal units ;
• VDW_RADIUS is a Dict{String, Integer} associating atoms symbols and Van der Waals radii, in internal units.

Using non-default computation¶

By default, the computation algorithm is automatically determined by the potential function type. ShortRangePotential are computed with CutoffComputation, and all other potentials are computed by DirectComputation. If we want to use another computation algorithm, this can be done by providing a computation keyword to the add_interaction! function. The following values are allowed:

• :direct to use a DirectComputation;
• :cutoff to use a CutoffComputation. The cutoff can be specified with the cutoff keyword argument;
• :table` to use a TableComputation. The table size can be specified with the numpoints keyword argument, and the maximum distance with the rmax keyword argument.

Potential functions¶

Potential computation¶

As stated at the begiging of this section, we need two informations to compute interactions between particles: a potential function, and a potential computation. The potential compuatation algorithms come in four flavors:

• The DirectComputation is only a small wrapper on the top of the potential functions, and directly calls the potential function methods for energy and force evaluations.
• The CutoffComputation is used for short range potentials. All interactions at a longer distance than the cutoff distance are set to zero. The default cutoff is \(12\ A\), and this can be changed by passing a cutoff keyword argument to the add_interaction function. With this computation, the energy is shifted so that their is a continuity in the energy at the cutoff distance.
• The TableComputation use table lookup to extrapolate the potential energy and the forces at a given point. This saves computation time at the cost of accuracy. This algorithm is parametrized by an integer, the size of the undelying array. Increases in this size will result in more accuracy, at the cost of more memory usage. The default size is 2000 points — which corresponds to roughly 15kb. TableComputation has also a maximum distance for computations, rmax. For any bigger distances, the TableComputation will returns a null energy and null forces. So TableComputation can only be used if you are sure that the particles will never be at a greater distance than rmax.
• The LongRangeComputation is not implemented yet.

Which computation for which potential ?¶

Not all computation algorithms are suitable for all potential functions. The usable associations are in the table below.

Function	DirectComputation	CutoffComputation	TableComputation
ShortRangePotential
BondedPotential
AnglePotential
DihedralPotential

Time integration¶

An integrator is an algorithm responsible for updating the positions and the velocities of the current frame of the universe. It is represented by a subtype of the Integrator type. You can set the integrator to use with your simulation using the set_integrator!() function.

Force computation¶

The NaiveForces algorithm computes the forces by iterating over all the pairs of atoms, and calling the appropriate interaction potential. This algorithm is the default in Jumos.

Controlling the simulation¶

While running a simulation, we often want to have control over some simulation parameters: the temperature, the pressure, etc. This is the goal of the control algorithms, all subtypes of the Control type.

Checking the simulation consistency¶

Molecular dynamic is usually a garbage in, garbage out set of algorithms. The numeric and physical issues are not caught by the algorithm themselves, and the physical (and chemical) consistency of the simulation should be checked often.

In Jumos, this is achieved by the Check algorithms, which are presented in this section. Checking algorithms can be added to a simulation by using the push! function.

Default algorithms¶

Default algorithms for molecular dynamic are presented in the following table:

Functions for algorithms selection¶

The following functions are used to to select specific algorithms for the simulation. They allow to add and change all the algorithms, even in the middle of a simulation run.

set_integrator!(sim, integrator)¶

Sets the simulation integrator to integrator.

Usage example:

sim = Simulation(MolecularDynamic(0.5))

run!(sim, 300)   # Run with a 0.5 fs timestep

set_integrator!(sim, Verlet(1.5))
run!(sim, 3000)  # Run with a 1.5 fs timestep

push!(simulation, check)¶

Adds a check or control algorithm to the simulation list and issues a warning if the algorithm is already present.

Usage example:

# Note the parentheses, needed to instanciate the new check.
push!(sim, AllPositionsAreDefined())

push!(sim, RescaleVelocities(300, 10))

Energy related values¶

type TemperatureCompute¶

Computes the temperature of the universe.

Key: :temperature

Return value: The current frame temperature.

type EnergyCompute¶

Computes the potential, kinetic and total energy of the current simulation step.

Keys: :E_kinetic, :E_potential, :E_total

Return value: A tuple containing the kinetic, potential and total energy.

energy = EnergyCompute()

# unpacking the tuple
E_kinetic, E_potential, E_total = energy(universe)

# accessing the tuple values
E = energy(sim)

E_kinetic = E[1]
E_potential = E[2]
E_total = E[3]

Volume¶

type VolumeCompute¶

Computes the volume of the current unit cell.

Key: :volume

Return value: The current cell volume

Pressure¶

type PressureCompute¶

TODO

Key:

Return value:

Existing outputs¶

Defining a new output¶

Adding a new output with custom values, can be done either by using a custom output or by subtyping the Output type to define a new output. The the former wayis to be prefered when adding a one-shot output, and the latter when adding an output which will be re-used.

Minimal images¶

These functions take a vector and wrap it inside a cell, by finding the minimal vector fitting inside the cell.

minimal_image(vect, cell::UnitCell)¶

Wraps the vect vector inside the cell unit cell. vect can be a Vector (i.e. a 1D array), or a view from an Array3D. This function returns the wrapped vector.

minimal_image!(vect, cell::UnitCell)¶

Wraps vect inside of cell, and stores the result in vect.

minimal_image!(A::Matrix, cell::UnitCell)

If A is a 3xN Array, wraps each one of the columns in the cell. The result is stored in A. If A is not a 3xN array, this throws an error.

Distances¶

Distances are computed using periodic boundary conditions, by wrapping the \(\vec r_i - \vec r_j\) vector in the cell before computing its norm.

Reading and writing trajectories¶

One can read or write frames from a trajectory. In order to do so, more information is needed : namely an unit cell and a topology. Both are optional, but allow for better computations. Some file formats already contain this kind of informations so there is no need to provide it.

Trajectories can exist in differents formats: text formats like the XYZ format, or binary formats. In Jumos, the format of a trajectory file is automatically determined based on the file extension.

traj_reader = Reader("filename.xyz", cell=[10., 10., 10.])

for frame in eachframe(traj_reader)
    # Do stuff here
end

close(traj_reader)

traj_reader = Reader("filename-in.nc", topology="topology.xyz")
traj_writer = Writer("filename-out.xyz")

for frame in eachframe(traj_reader)
    write(traj_writer, frame)
end

close(traj_writer)
close(traj_reader)

Supported formats¶

The following table summarizes the formats supported by Jumos, giving the reading and writing capacities of Jumos, as well as the presence or absence of the unit cell and the topology information in the files. The last column indicates the accepted keywords.

Format	Extension	Read	Write	Cell	Topology	Keywords
XYZ	.xyz					cell
Amber NetCDF	.nc					topology

Readind and writing topologies¶

Topologies can also be represented and stored in files. Some functions allow to read directly these files, but there is usally no need to use them directely.

Supported formats for topology¶

Topology reading supports the formats in the following table.

Format	Reading ?	Writing ?
XYZ
LAMMPS data file

If you want to write a toplogy to a file, the best way for now is to create a frame with this topology, and write this frame to an XYZ file.

User potential¶

The easier way to define a new potential is to create UserPotential instances, providing potential and force functions. To add a potential, for example an harmonic potential, we have to define two functions, a potential function and a force function. These functions should take a Float64 value (the distance) and return a Float64 (the value of the potential or the force at this distance).

UserPotential(potential, force)¶

Creates an UserPotential instance, using the potential and force functions.

potential and force should take a Float64 parameter and return a Float64 value.

UserPotential(potential)

Creates an UserPotential instance by automatically computing the force function using a finite difference method, as provided by the Calculus package.

Here is an example of the user potential usage:

# potential function
f(x) = 6*(x-3.)^2 - .5
# force function
g(x) = -12.*x + 36.

# Create a potential instance
my_harmonic_potential = UserPotential(f, g)

# One can also create a potential whithout providing a funtion for the force,
# at the cost of a less effective computation.
my_harmonic_2 = UserPotential(f)

force(my_harmonic_2, 3.3) == force(my_harmonic_potential, 3.3)
# false

isapprox(force(my_harmonic_2, 3.3), force(my_harmonic_potential, 3.3))
# true

Subtyping PotentialFunction¶

A more efficient way to use custom potential is to subtype the either PairPotential, BondedPotential, AnglePotential or DihedralPotential, according to the new potential from.

For example, we are goning to define a Lennard-Jones potential using an other function:

This is obviously a ShortRangePotential, so we are going to subtype this potential function.

To define a new potential, we need to add methods to two functions: call and force. It is necessary to import these two functions in the current scope before extending them. Potentials should be declared as immutable, this allow for optimizations in the code generation.

# import the functions to extend
import Base: call
import Jumos: force

immutable LennardJones2 <: PairPotential
    A::Float64
    B::Float64
end

# potential function
function call(pot::LennardJones2, r::Real)
    return pot.A/(r^12) - pot.B/(r^6)
end

# force function
function force(pot::LennardJones2, r::Real)
    return 12 * pot.A/(r^13) - 6 * pot.B/(r^7)
end

The above example can the be used like this:

# Add a LennardJones2 interaction to an universe
universe = Universe(...)
add_interaction!(universe, LennardJones2(4.5, 5.3), "He", "He")

# Directly compute values
pot = LennardJones2(4.5, 5.3)
pot(3.3) # value of the potential at r=3.3
force(pot, 8.12) # value of the force at r=8.12

Computing values¶

Algorithms to compute properties from a simulation are called Compute in Jumos. They all act by taking an Universe as parameter, and then setting a value in the Universe.data dictionary.

To add a new compute algorithm (we will call it MyCompute), we have to subtype Jumos.Simulations.BaseCompute and provide specialised implementation for the Base.call function; with the following signature:

Base.call(::Compute, ::Universe)

This function can set a Universe.data entry with a Symbol key to store the computed value. This value can be anything, but basic types (scalars, arrays, etc.) are to be prefered.

Outputing values¶

An other way to create a custom output is to subtype Output. The subtyped type must have at least one field: frequency::OutputFrequency, which is responsible for selecting the write frequency of this output. Two functions can be overloaded to set the output behaviour.

Base.write(::Output, context::Dict{Symbol, Any})

Write the ouptut. This function will be called every n simulation steps according to the frequency field.

The context parameter contains all the values set up by the computation algorithms, and a special key :frame refering to the current simulation frame.

Jumos.setup(::Output, ::Simulation)¶

This function is called once, at the begining of a simulation run. It should do some setup job, like adding the needed computations algorithms to the simulation.

As an example, let’s build a custom output writing the x position of the first atom of the simulation at each step. This position will be taken from the frame, so no specific computation algorithm is needed here. But this position will be writen in bohr, so some conversion from Angstroms will be needed.

# File FirstX.jl

using Jumos

import Base.write
import Jumos.setup

type FirstX <: Output
    file::IOStream
    frequency::OutputFrequency
end

# Default values constructor
function FirstX(filename, frequency=1)
    file = open(filename, "w")
    return FirstX(file, OutputFrequency(frequency))
end

function write(out::FirstX, context::Dict)
    frame = context[:frame]
    x = frame.positions[1][1]
    x = x/0.529 # Converting to bohr
    write(out.file, "$x \n")
end

# Not needed here
# function setup(::FirstX, ::Simulation)

This type can be used like this:

using Jumos
require("FirstX.jl")

sim = Simulation(:md, 1.0)
# ...

push!(sim, FirstX("The-first-x-file.dat"))

Writing a new integrator¶

To create a new integrator, you have to subtype the Integrator type, and provide the call method, with the following signature: call(::Integrator, ::MolecularDynamic).

The integrator is responsible for calling the getforces!(::MolecularDynamics, ::Universe) function when the MolecularDynamics.forces field needs to be updated. It should update the two Array3D: Universe.frame.positions and Universe.frame.velocities with appropriate values.

The Universe.masses field is a Vector{Float64} containing the particles masses. Any other required information should be stored in the new Integrator subtype.

Computing the forces¶

To create a new force computation algorithm, you have to subtype ForcesComputer, and provide the method call(::ForcesComputer, ::Universe, forces::Array3D).

This method should fill the forces array with the forces acting on each particles: forces[i] should be the 3D vector of forces acting on the atom i. In order to do this, the algorithm can use the Universe.frame.posistions and Universe.frame.velocities arrays.

Potentials to use for the atoms can be optained throught the following functions:

pairs(::Interactions, i, j)¶

Get a Vector{PairPotential} of interactions betweens atoms i and j.

bonds(::Interactions, i, j)¶

Get a Vector{BondPotential} of interactions betweens atoms i and j.

angles(::Interactions, i, j, k)¶

Get a Vector{AnglePotential} of interactions betweens atoms i and j and k.

dihedrals(::Interactions, i, j, k, m)¶

Get a Vector{DihedralPotential} of interactions betweens atoms i, j, k and m.

Adding new checks¶

Adding a new check algorithm is as simple as subtyping Check and extending the call(::Check, ::Universe, ::MolecularDynamic) method. This method should throw an exception of type CheckError if the checked condition is not fullfiled.

type CheckError¶

Customs exception providing a specific error message for simulation checking.

julia> throw(CheckError("This is a message"))
ERROR: Error in simulation :
This is a message
in __text at no file (repeats 3 times)

Adding new controls¶

To add a new type of control to a simulation, the main way is to subtype Control, and provide two specialised methods: call(::Control, ::Universe) and the optional setup(::Control, ::Simulation). The call method should contain the algorithm inplementation, and the setup method is called once at each simulation start. It should be used to add add some computation algorithm to the simulation, as needed.