Developping new algorithms¶
Writing a new integrator¶
To create a new integrator, you have to subtype the BaseIntegrator
type, and
provide the call
method, with the following signature:
call(::BaseIntegrator, ::MolecularDynamic)
.
The integrator is responsible for calling the get_forces!(::MolecularDynamic)
function if it need the MolecularDynamic.forces
field to be updated. It should
update the two Array3D: MolecularDynamic.frame.positions
and MolecularDynamic.frame.velocities
with appropriate values.
The MolecularDynamic.masses
field is a Vector{Float64}
containing the particles
masses, in the same order than in the current simulation frame.
Any other required information should be stored in the new BaseIntegrator
subtype.
Computing the forces¶
Naive force computation¶
The NaiveForceComputer
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.
New algorithm for forces computations¶
To create a new force computation algorithm, you have to subtype BaseForcesComputer
,
and provide the method call(::BaseForcesComputer, forces::Array3D, frame::Frame,
interactions)
.
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 frame.posistions
and frame.velocities
.
interactions
is a dictionary associating tuples of integers (the atoms types)
to potential. The get_potential
function can be useful
to get a the potential function acting on two particles.
-
get_potential
(interactions, topology, i, j)¶ Returns the potential between the atom i and the atom j in the topology.
Due to the internal unit system, forces returned by the potentials are in
\(kJ/(mol \cdot A)\), and should be in \(uma \cdot A / fs^2\) for being
used with the newton equations. The conversion can be handled by the unexported
Simulations.force_array_to_internal!
function, converting the values of an
Array3D from \(kJ/(mol \cdot A)\) to \(uma \cdot A / fs^2\).
Adding a new check¶
Adding a new check algorithm is as simple as subtyping BaseCheck
and extending
the call(::BaseCheck, ::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
BaseControl
, and provide two specialised methods: call(::BaseControl,
::MolecularDynamic)
and the optional setup(::BaseControl, ::MolecularDynamic)
.
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.
Computing values¶
To add a new compute algorithm (MyCompute
), we have to subtype BaseCompute
and provide specialised implementation for the call
function; with the
following signature:
-
call
(::MyCompute, ::MolecularDynamic)¶ This function can set a
MolecularDynamic.data
entry with any kind of key to store the computed value.
Outputing values¶
An other way to create a custom output is to subtype BaseOutput
. The subtyped
type must have two integer fields: current
and frequency
, and the constructor
should initialize current
to 0. The write
function should also be overloaded
for the signature write(::BaseOutput, ::Dict)
. The dictionairy parameter
contains all the values set up by the computation algorithms,
and a special key :frame
refering to the current simulation frame.
BaseOutput
subtypes can also define a setup(::BaseOutput, ::MolecularDynamic)
function to do some setup job, like adding the needed computations 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 <: BaseOutput
file::IOStream
current::Integer
frequency::Integer
end
# Default values constructor
function FirstX(filename, frequency=1)
file = open(filename, "w")
return FirstX(file, 0, 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, ::MolecularDynamic)
This type can be used like this:
using Jumos
require("FirstX.jl")
sim = MolecularDynamic(1.0)
# ...
add_output(sim, FirstX("The-first-x-file.dat"))