Mechanica real-time interactive biological and active-matter simulation¶
Contents:
Introduction¶
Mechanica is an interactive particle based physics, chemistry and biology simulation environment, with a heavy emphasis towards enabling users to model and simulate complex sub-cellular and cellular biological physics problems. Mechanica is part of the Tellurium http://tellurium.analogmachine.org project.
Mechanica is designed first and foremost to enable users to work interactively with simulations – so they can build, and run a simulation in real-time, and interact with that simulation whilst it’s running. The goal is to create an SolidWorks type environment where users can create and explore virtual models of soft condensed matter physics, with a emphasis towards biological physics.
Mechanica is a native compiled C++ shared library with a native and extensive Python API, that’s designed to used from an ipython console (or via scripts of course).
Biological cells and biochemical molecules and particles are the prototypical example of active matter. These are all actively driven agents that transduct free energy from their environment. These agents can sense and respond to mechanical, chemical and electrical environmental stimuli with a range of behaviors, including dynamic changes in morphology and mechanical properties, chemical uptake and secretion, cell differentiation, proliferation, death, and migration.
One of the greatest challenges at these medium length scales is that these dynamics and behaviors typically cannot be derived from first principles. Rather, the observed behaviors are described phenomenologically or empirically. Thus, the scientist exploring these kinds of phenomena needs a great deal of flexibility to propose and experiment with different kinds of interactions.
This presents a significant challenge, as simulation environments that make it simple for the end user to write models without resorting to hard-coding C++ or FORTRAN usually are very limited in the level of flexibility they provide the end user. For example, if users want to write a standard molecular dynamics model, there are many different, really good choices of simulation engines and these kinds of models can easily be specified by human readable configuration files. However, as the kinds of interactions are not well standardized or formalized at medium length scales, users almost always are forced to resort to hard-coding FORTRAN or C++.
Particle Dynamics enables modeling a wide range of length scales
Our goal here is to deliver a modeling and simulation framework that lets users INTERACTIVELY create, simulate and explore models at biologically relevant length scales. We believe that interactive simulation is key to increasing scientific productivity, much like interactive modeling environments such as SolidWorks has revolutionized engineering practice.
We thus present Mechanica, an interactive modeling and simulation environment based on off-lattice formalism that lets users create models for a wide range of biologically relevant problems, and we enable users to write models using any combination of the following modeling methodologies:
- Coarse Grained Molecular Dynamics
- Discrete Element Method (DEM). DEM particles add rotational degrees-of-freedom as well as stateful contact and often complicated geometries (including polyhedra).
- Dissipative Particle Dynamics (DPD) is a particle-based method, where particles represent whole molecules or fluid regions rather than single atoms, and atomistic details are not considered relevant to the processes addressed. The particles’ internal degrees of freedom are averaged out and represented by simplified pairwise dissipative and random forces, so as to conserve momentum locally and ensure correct hydrodynamic behavior. DPD allows much longer time and length scales than are possible using conventional MD simulations.
- Sub-Cellular Element (SCM). Frequently used to model complex sub-cellular active mechanics. SCM are similar to DPD, where each particle represents a region of space and is governed by empirically derived potentials, but adds active response.
- Smoothed particle hydrodynamics (SPH)is a particle method very similar to DPD and is frequently used to model complex fluid flows, especially large fluid deformations, fluid-solid interactions, and multi-scale physics.
- Reactive Molecular Dynamics. In RMD, particles react with other particles and form new molecules, and can absorb or emit energy into their environment. Mechanica is designed to support reactive particles, as one of our main goals is very efficient particle creation and deletion. Very few classical molecular dynamics packages support reactive MD, as they are almost all highly optimized towards conserved number of particles.
- Perhaps most uniquely, Mechanica allows users to attach a chemical cargo to each particle, and host a chemical reaction network at each element. Furthermore, we allow users to write fluxes between particles. A flux defines a movement of material from one site to another. Furthermore, we also allow users to attach their own handlers to a variety of different events that particles (or other objects) can emit. Therefore, we also support developing full Transport Dissipative Particle Dynamics simulations.
- Flux Networks. The concept of a flux is extremly general, and this lets us define a connector type that lets users connect different model elements. Flux networks allow us to define a wide range of problems, from biological fluid flow in areas like the liver and the eye, to physiologically based pharmacokinetic (PBPK) modeling, and even to electric circuits and pipe flow networks.
Warning
Only a subset of these features are presently available, and we encourage users to look at the Status page, and PLEASE LET US KNOW WHAT FEATURES YOU WANT. We can only deliver the kind of software users want if you let us know what features you want to see. Please contact us at <somogyie@indiana.edu> or on Twitter at @AndySomogyi
Once we have a well-defined, and user tested API for generalized particle dynamics, we will integrate our existing Vertex Model code into Mechanica. Vertex Model is another specialized form of classical Molecular Dynamics, but with instead of the traditional bonded relationships of bonds, angles, dihedrals, impropers, Vertex Models add some new kinds of bonded relationships such as polygons and volumes to represent surface and volume forces.
The kinds of problems Mechanica is designed to enable users to model.
Getting Mechanica¶
The easiest way to install Mechanica for most users is via PIP. We presently have PIP binaries for Windows and MacOS. We do have a Ubuntu 18.04 pip wheels, available on GitHub, but these are not on PyPi because they are not “manylinux” compatible.
Mechanica requires at least Python version 3.7.0, on Mac you can install Python a variety of ways, but we typicaly use Brew, https://brew.sh.
Note
Spyder has a pretty complex architecture, and Spyder support is new. We’ve tried to test as best as we can, but this is alpha code, you may encounter issues.
Installing via pip for Mac and Windows¶
Note, we presently only have Windows and Mac PyPi packages.
Python comes with an inbuilt package management system, pip. Pip can install, update, or delete any official package. The PyPi home page for mechancia is
https://pypi.org/project/mechanica/.
You can install packages via the command line by entering:
python -m pip install --user mechanica
We recommend using an user install, sending the --user flag to pip.
pip installs packages for the local user and does not write to the system
directories.
Installing via pip for Linux¶
We presently only support Ubuntu 18.04, pip wheel binaries are here:
Download and run pip directly on this file.
Preferably, do not use sudo pip, as this combination can cause problems.
Quickstart¶
Argon¶
This example will create a complete simulation of a set of argon atoms.
First we simply import the packages, we will use Numpy to create initial conditions:
import mechanica as m
import numpy as np
Define some variables that define the potential cutoff distance and size of the simulation domain. Potential cutoff is the maximum distance that for which the runtime evaluates potentials, for performance reasons, it should be kept reasonably small
# potential cutoff distance
cutoff = 1
# dimensions of universe
dim=[10., 10., 10.]
The first thing we have to do, before we create any simulation objects is initialize the simulator. This essentially sets up the simulation environment, and gives us a place to create our model.
m.Simulator(dim=dim)
Particles interact via Potentials, Potential. We supply a variety of
pre-built potentials, and users can create their own, but for now, the
Lennard-Jones 12-6 is one of the most commonly used potentials. All we have to
do is create an instance of one, and bind it to an an object type. To create a
LJ potential
# create a potential representing a 12-6 Lennard-Jones potential
# A The first parameter of the Lennard-Jones potential.
# B The second parameter of the Lennard-Jones potential.
# cutoff
pot = m.Potential.lennard_jones_12_6(0.275 , cutoff, 9.5075e-06 , 6.1545e-03 , 1.0e-3 )
It’s most common to sub-class the base particle type. By subclassing it, users can specify a variety of attributes that define that particular type, such as mass.:
# create a particle type
# all new Particle derived types are automatically registered with the universe
class Argon(m.Particle):
mass = 39.4
Whenever a new Particle derived type is created, it is automatically registerd with the Universe.
The total force on any phyiscal object such as a particle is simply the sum of all potentials that act on that object. To make a potential act on a particular kind of type, we bind the potential to the type of object we want it to act on:
# bind the potential with the *TYPES* of the particles
m.Universe.bind(pot, Argon, Argon)
For our simulation, we want to fill the simulation domain with a uniformly random distrubted particles. The simplest way to do this is the use the numpy random function to generate an array of uniform random distibuted positions:
# uniform random cube
positions = np.random.uniform(low=0, high=10, size=(10000, 3))
Then with the array of positions, we simply itterate over the postions, and create a new particle; When we create a new particle derived object, it gets automatically added to the universe, we don’t need to do anything else with it:
for pos in positions:
# calling the particle constructor implicitly adds
# the particle to the universe
Argon(pos)
Now all that’s left is to run the simulation. The Simulator object has
two methods on it, run and irun. The run method runs the simulator, and
continues until the window is closed, or some stop condition. The irun method
starts the simulation, but leaves the console open for further input. The irun
method is only used when running from iPython.:
# run the simulator interactive
m.Simulator.run()
The complete simulation script is here, and can be downloaded here:
Download: this example script:
import mechanica as m
import numpy as np
# potential cutoff distance
cutoff = 1
# dimensions of universe
dim=[10., 10., 10.]
m.Simulator(dim=dim)
# create a potential representing a 12-6 Lennard-Jones potential
# A The first parameter of the Lennard-Jones potential.
# B The second parameter of the Lennard-Jones potential.
# cutoff
pot = m.Potential.lennard_jones_12_6(0.275 , cutoff, 9.5075e-06 , 6.1545e-03 , 1.0e-3 )
# create a particle type
# all new Particle derived types are automatically
# registered with the universe
class Argon(m.Particle):
mass = 39.4
# bind the potential with the *TYPES* of the particles
m.Universe.bind(pot, Argon, Argon)
# uniform random cube
positions = np.random.uniform(low=0, high=10, size=(10000, 3))
for pos in positions:
# calling the particle constructor implicitly adds
# the particle to the universe
Argon(pos)
# run the simulator interactive
m.Simulator.run()
A basic argon simulation
Examples¶
Mitosis and Events¶
In cell biology, mitosis is what the cell division process is called. This is where a cell undergoes fission and splits into two. The mitosis process splits the contents of the parent cell into two daughter cells. Usually the two daughter cells are the same type, but sometimes differen, and the contents usually get split semi-uniformly, but not always.
The mitosis process is part of the cell cycle, and usually gets triggered when the parent cell reaches some specified growth stage. Mitosis is extremly complicated when we look at every single detail, but in this example, we model an extremly simplified mitosis, which is triggered at regular time intervals. This is toy model simply to illustrate how to use time based events.
A basic mitosis simulation, about 15 lines of Python.
We start with the normal imports, and setup the domain size, cutoff distance, and make a simulator
import mechanica as m
import numpy as np
# potential cutoff distance
cutoff = 1
# new simulator, don't load any example
m.Simulator(dim=[20., 20., 20.])
Make a soft sphere interaction potential. The kappa here is the strength of the repussion (peak repussion energy), and the epsilon is the depth of the attractive well, r0 is the rest-length of this potential, i.e. if we set a r0 to 5, two particles will rest at center separation distance of 5. min and max are the cutoff ranges, and tol is the tolerance, or how accuratly this potential matches the exact potential, a value of around 0.1 is fine:
pot = m.Potential.soft_sphere(kappa=400, epsilon=5,
r0=0.7, eta=3, tol = 0.1, min=0.05, max=1.5)
Make a Cell class that extends particle, set the usuall class
properties. Here we add an on_time event, and bind it to the base
Particle.fission method.
What happens here, is that for every instance of the Cell class, the on_time event will get triggered at a period of 1 time unit, and call the fission method on that object. Basically, we’re addign an event handler to every single instance of the object when we define an event in the class definition:
class Cell(m.Particle):
mass = 20
target_temperature = 0
radius=0.5
events = [m.on_time(m.Particle.fission, period=1, distribution='exponential')]
dynamics=m.Overdamped
Here, we also added the Here, we also added the dynamics=m.Overdamped option to the cell definition. This essentially makes the object move with overdamped motion, in that the only time the object moves is when there are forces acting on it. Overdamped is generally how all biological materials behave.
Now we go ahead and proceed as usual, and simply bind the potential to the Cell objects. Again, this will create a small attractive adhesion between the objects, and and a spring like repulsion to enable the objects to push on each other:
m.Universe.bind(pot, Cell, Cell)
Make a random force. This essentailly gives each cell a certian motility, in that it will move around small random motions. The first arg is the mean, and the second is the strength of the random force. Bind the random force just like any other force:
rforce = m.forces.random(0, 0.5)
m.Universe.bind(rforce, Cell)
Create a single cell instance at the center of our domain, and start the simulation:
Cell([10., 10., 10.])
# run the simulator interactive
m.Simulator.run()
The complete simulation script is here, and can be downloaded here:
Download the overdamped example:
mitosis.od.py
We also included a version that does not use overdamped dynamics. Run this one, and see how the individaul objects move around and never really settle down. In conventional dynamics, each object has an inertia and will continue to move unless it is stopped by some other force. In overdamped dyanamics, object move only where there is a force on them.
A mitosis simulation usign conventional inertial dymamcis.
You can download the complete dynamics example here:
mitosis.dyn.py
Two Types¶
We can make a model with more than one type simply by making different
Particle derived classes:
import mechanica as m
import numpy as np
# potential cutoff distance
cutoff = 8
count = 1000
# dimensions of universe
dim=np.array([20., 20., 20.])
center = dim / 2
# new simulator, don't load any example
m.Simulator(dim=dim, cutoff=cutoff)
Make a Big and Small types:
class Big(m.Particle):
mass = 500000
radius = 3
class Small(m.Particle):
mass = 0.1
radius = 0.2
target_temperature=0
Make a pair on interaction potentials, one between the Small particle types, and one betwen the Big and Small particles. Here, we choose the r0 paramter according to the particle types radius. TODO: next version will automatically pick up the radisu from the type definitions. We bind these two potentials between the different particle types.:
pot_bs = m.Potential.soft_sphere(kappa=10000, epsilon=100, r0=3.2, \
eta=3, tol = 0.1, min=0.1, max=8)
pot_ss = m.Potential.soft_sphere(kappa=50, epsilon=10, r0=0.2, \
eta=2, tol = 0.05, min=0.01, max=4)
# bind the potential with the *TYPES* of the particles
m.bind(pot_bs, Big, Small)
m.bind(pot_ss, Small, Small)
Make a single Big particle in the middle of our domain:
Big(position=center, velocity=[0., 0., 0.])
Make a disk of small particles above the big particle. We use the built-in :function::random_point function to fill regions of space with randomly distributed points:
for p in m.random_point(m.Disk, count) * \
1.5 * Big.radius + center + [0, 0, Big.radius + 1]:
Small(p)
And finally run the simulation:
# run the simulator interactive
m.Simulator.run()
A basic two-type simulation, about 15 lines of Python. We are using Newtonian instead of overdamped dynamics, so each object has a good deal of inertia. We notice that single there is an attractive potential between the big and small objcts, most of the small objects are clusterd aournd the single big particle in the middle. Also notice that the small particles tend to produce clusters, as they too are attractive to each other.
The complete simulation script is here, and can be downloaded here:
Download: this example script:
Cell Sorting¶
In this section, we use cell sorting as a motivating example to develop a mechanica model. We first introduce the biological phenomena, then we start our analysis by first identifying what are the key objects and processes in the biological descrioption, and how we represent these physical concepts with Mechanica objects and processes.
Cell sorting between biological cells of different types is one of the basic mechanisms creating tissue domains during develoment and wound healing and in maintaining domains in homeostasis. Cells of two different types, when dissociated and randomly mixed can spontaneously sort to reform ordered tissues. Cell sorting can generate regular patterns such as checkerboards, engulfment of one cell type by another, or other more complex patterns [Ste70].
A combination of cells from 3-germ layers are dissociiated and re-organize, each cell type sorts out into its own region [TH55].
Both complete and partial cell sorting (where clusters of one cell type are trapped or contained inside a closed envelope of another cell type) have been observed experimentally in vitro in embryonic cells. Cell sorting does not involve cell division nor differentiation but only spatial rearrangement of cell positions.
A basic two-type simulation, about 25 lines of Python, as a simple extension of the two type example.
In a classic in vitro cell sorting experiment to determine relative cell adhesivities in embryonic tissues, mesenchymal cells of different types are dissociated, then randomly mixed and reaggregated. Their motility and differential adhesivities then lead them to rearrange to reestablish coherent homogenous domains with the most cohesive cell type surrounded by the less-cohesive cell types [AP72] [AA84].
Cell-sorting behavior of cell aggregates is similar to liquid surface tension, in the spontaneous separation of immiscible liquids (water vs. oil). Adhesive forces between mixed cells play a similar role in cell sorting that intermolecular attractive (cohesive) forces play in liquid surface tension. In cell sorting, the cells with the strongest highest adhesivities will be sorted to the center, while the less cohesive ones will remain outside.
Modeling and Abstraction¶
To develop a computational model of our biological system, we must first identify the key objects and processes of our physical system. If we look at the left side of the following figure, we can see a sheet of biogical cells. From the previous description of the cell sorting, we also know that cells move about. We know that epethelial sheets are essentially a sheet of cells that form a surface. Here we can identify our first biogical object, a cell. From the previous discussion, we know that there are more than one type of cell, so lets call our two cell types, cell type A and cell type B.
On the left is a confocal microscope image of a developing Drosophila wing, where the outlines of the cells have been highlighte with a florescent protein, which binds to E-cadherin, a surface proteion involved in cell adhesion. We can represent this sheet of biological cells with a set of polygons constrained to a two dimensional surface. Taken from [FOBS14]
Example Model¶
We start the example just like other Mechania models:
import mechanica as m
import numpy as np
# total number of cells
A_count = 5000
B_count = 5000
# potential cutoff distance
cutoff = 3
# dimensions of universe
dim=np.array([20., 20., 20.])
center = dim / 2
# new simulator, don't load any example
m.Simulator(dim=dim, cutoff=cutoff)
We make two different cell types, A and B:
class A(m.Particle):
mass = 1
radius = 0.5
dynamics = m.Overdamped
class B(m.Particle):
mass = 1
radius = 0.5
dynamics = m.Overdamped
To represent the cell interactions, we create three different interacation potentials, pot_aa, pot_bb, and pot_ab. These represent the strength of interaction between cell types:
# create three potentials, for each kind of particle interaction
pot_aa = m.Potential.soft_sphere(kappa=400, epsilon=40, r0=1.5, \
eta=2, tol = 0.05, min=0.01, max=3)
pot_bb = m.Potential.soft_sphere(kappa=400, epsilon=40, r0=1.5, \
eta=2, tol = 0.05, min=0.01, max=3)
pot_ab = m.Potential.soft_sphere(kappa=400, epsilon=5, r0=1.5, \
eta=2, tol = 0.05, min=0.01, max=3)
And bind those potentials to our cell types:
# bind the potential with the *TYPES* of the particles
m.Universe.bind(pot_aa, A, A)
m.Universe.bind(pot_bb, B, B)
m.Universe.bind(pot_ab, A, B)
We need a random force to keep the cells moving around, this represents basically a cell motility:
# create a random force. In overdamped dynamcis, we neeed a random force to
# enable the objects to move around, otherwise they tend to get trapped
# in a potential
rforce = m.forces.random(0, 50)
# bind it just like any other force
m.bind(rforce, A)
m.bind(rforce, B)
Create the particle instances:
# create particle instances, for a total A_count + B_count cells
for p in np.random.random((A_count,3)) * 15 + 2.5:
A(p)
for p in np.random.random((B_count,3)) * 15 + 2.5:
B(p)
And finally run the simulation:
# run the simulator
m.Simulator.run()
The complete simulation script is here, and can be downloaded here:
Download: this example script:
Writing Output Data With Events¶
Writing output data is a frequently done task with simulations. The
built-in on_time event provides a very convenient system to register an
output function that you want to get called at regular intervals.
During event execution, the simulation is halted, so it’s important to do whatever processing quickly as no to slow down the simulation.
An example of an output function would be to iterate over the particles, and write their positions to a file. Starting with out previous mitosis events example, Mitosis and Events, we can add a function to write output data, and add that to a time event. We use the Python csv writer to write this to an output file:
def write_data(time):
print("time is now: ", time)
positions = [list(p.position) for p in m.Universe.particles]
with open(fname, "a") as f:
writer = csv.writer(f)
writer.writerow([time, positions])
Here we iterate over the particles, grab the particle positions, convert it to a list, and concatenate that together with the time, and write that out as a row to the output file.
And we attach that function to the on_time event. We only specify a
period, and do not give any distribution as we want this to be called at
regular intervals:
m.on_time(write_data, period=0.05)
The complete simulation script is here, and can be downloaded here:
Download: this example script:
Worm Like Chain¶
Cyanobacteria often form worm-like filaments and may grow in large masses or “tufts” one meter or more in length. Some are unicellular, a few form branched filaments, and a few form irregular plates or irregular colonies. Cyanobacterial cells usually divide by binary fission, and the resulting cells may separate to form new colonies. Some of the filamentous cyanobacteria are motile by means gliding or rotating around a longitudinal axis,
|
|
This model starts a worm-like chain and introduces the Mechanica Bond and Angle Bonded Interactions to create a very simply bacterial filament model. This model does not include self-propelled particles, we will introduce these later.
In many bacterial filaments, each bacterial cell sticks to at most one or two
neighbors. We can not use a non-bonded type potential for this, as a non-bonded
interaction would stick to all bacteria with a range. Instead, we use a
Bond to explicitly bond pairs of bacterial cells. Bacterial filaments
tend to have a certain stiffness, we use a Angle to model this resistance
to bending.
We can create a filament starting with the standard simulation setup:
# potential cutoff distance
cutoff = 8
# dimensions of universe
dim=np.array([20., 20., 20.])
center = dim / 2
# new simulator, don't load any example
m.Simulator(dim=dim, cutoff=cutoff)
We approximate a bacteria with a simple spherical Bead. For the time being, we ignore the fact that bacteria in fact tend to be lozenge shaped. We’ll cover non-spherical particles later:
class Bead(m.Particle):
mass = 0.4
radius = 0.2
dynamics = m.Overdamped
Use over-damped dynamics to model Brownian motion. To represent steric occlusion, add a soft-sphere interaction potential to the bacterial particles. This form of the potential creates a slightly compressible bacteria, and with a non-zero epsilon term, is slightly attractive to other bacteria:
pot_bb = m.Potential.soft_sphere(kappa=20, epsilon=5.0, \
r0=0.2, eta=4, tol=0.01, min=0.01, max=0.5)
Use a harmonic potential to model the attachment of neighboring cells:
# harmonic bond between particles
pot_bond = m.Potential.harmonic(k=40., r0=0.2, max = 2)
And a harmonic_angle potential to represent this stiffness of the filament:
# angle bond potential
pot_ang = m.Potential.harmonic_angle(k=20, theta0 = 0.85*np.pi, tol=0.1)
Attach the non-bonded potential to the Bead type, this will cause all instances of the bead to be able to push on each other, and give them a definite shape:
# bind the potential with the *TYPES* of the particles
m.bind(pot_bb, Bead, Bead)
Make a random force to apply to the beads, In over-damped dynamics, we need a random force to enable the objects to move around, otherwise they tend to get trapped in a potential:
rforce = m.forces.random(0, 1)
Bind it just like any other force:
m.bind(rforce, Bead)
We want to start our filament basically oriented in a line, so make a numpy array of contiguous positions for the x axis, and hole the y and z fixed. We will iterate over this array of positions, and create a new Bead particle for each one. Recall that a Bond connects adjacent pairs of bacteria, but a Angle is a 3-body term, so we need to attach it to each bacterial cell, but we need to keep track of the previous and next cells:
# make a array of positions
xx = np.arange(4., 16, 0.15)
p = None # previous bead
bead = Bead([xx[0], 10., 10.0]) # current bead
for i in range(1, xx.size):
n = Bead([xx[i], 10.0, 10.0]) # create a new bead particle
m.Bond(pot_bond, bead, n) # create a bond between prev and current
if(i > 1):
m.Angle(pot_ang, p, bead, n) # make an angle bond between prev, cur, next
p = bead
bead = n
And finally run the simulator:
m.Simulator.run()
The simulation output shows filament, you can play around with the initial positions, natural angle, bond and sphere stiffness.
The complete simulation script is here, and can be downloaded here:
Download: this example script:
Epiboly¶
Epiboly is one of the five major types of cell movement processes that occur in early embryonic development in many organisms, such as sea urchins, tunicates, amphibians. In epiboly, a layer of epithelial cells divides and spreads around a yolk sac. The exact mechanicsm of epiboly are an active reasarch area, so an interactive simulation enviorment such as Mechanica is an ideal testbed for exploring different hypotheses.
We can create a very simple epiboly type model starting with our two type model, the main change is use overdamped instead of Newtonian dynamics. Start with the normal model:
import mechanica as m
import numpy as np
# potential cutoff distance
cutoff = 8
# number of small cells
count = 3000
# dimensions of universe
dim=np.array([20., 20., 20.])
center = dim / 2
# new simulator, don't load any example
m.Simulator(dim=dim, cutoff=cutoff)
Make big and small particles just like before, execept we now set the small particle dynamics to Overdamped. We use a Big cell to represent a yolk sac, and small particles to represent epithelial cells:
class Big(m.Particle):
mass = 500000
radius = 3
class Small(m.Particle):
mass = 0.1
radius = 0.2
target_temperature=0
dynamics = m.Overdamped
Create a couple interaction potentials that act between the agent types:
pot_bs = m.Potential.soft_sphere(kappa=10, epsilon=50, r0=2.9, \
eta=3, tol = 0.1, min=0.1, max=9)
pot_ss = m.Potential.soft_sphere(kappa=20, epsilon=0.0001, r0=0.2, \
eta=2, tol = 0.05, min=0.01, max=3)
# bind the potential with the *TYPES* of the particles
m.bind(pot_bs, Big, Small)
m.bind(pot_ss, Small, Small)
# create a random force. In overdamped dynamcis, we neeed a random force to
# enable the objects to move around, otherwise they tend to get trapped
# in a potential
rforce = m.forces.random(0, 0.00000001)
# bind it just like any other force
m.bind(rforce, Small)
Make a single big particle (yolk sac) in the middle:
Big(position=center, velocity=[0., 0., 0.])
Make a disk of particles above the big one:
for p in m.random_point(m.Disk, count) * \
1.5 * Big.radius + center + [0, 0, Big.radius + 1]:
Small(p)
Run the simulation:
m.Simulator.run()
A basic epiboly simulation. The epithelial cells start in a disc above the yolk, and are attracted to the yolk. Since the epethialial cells both adhear to each other, and the yolk, and they are soft spheres, so the push other cells out of the way. Thus, they start as a blob, and relax to a sheet coverign the complete surface of the yolk sac.
The complete simulation script is here, and can be downloaded here:
Download: this example script:
Advanced Epiboly With Clusters and Events¶
We can write a much more detailed model of epiboly using clusters and events:
import mechanica as m
import numpy as np
# dimensions of universe
dim=np.array([30., 30., 30.])
center = dim / 2
m.Simulator(
dim=dim,
cutoff=12,
integrator=m.FORWARD_EULER,
dt=0.0005)
class C(m.Cluster):
radius=2.3
class B(m.Particle):
radius=0.25
dynamics = m.Overdamped
mass=15
style={"color":"skyblue"}
def split(self, event):
print("C.split(" + str(self) + ", event: " + str(event) + ")")
axis = self.position - C.yolk_pos
print("axis: " + str(axis))
m.Cluster.split(self, axis=axis)
events = [m.on_time(split,
period=0.2,
predicate="largest")]
class Yolk(m.Particle):
radius = 10
mass = 1000000
dynamics=m.Overdamped
style={"color":"gold"}
total_height = 2 * Yolk.radius + 2 * C.radius
yshift = total_height/2 - Yolk.radius
cshift = total_height/2 - C.radius - 1
yolk = Yolk(position=center-[0., 0., yshift])
c = C(position=center+[0., 0., cshift])
C.yolk_pos = yolk.position
c.B(4000)
pb = m.Potential.soft_sphere(kappa=300, epsilon=6, r0=0.5, \
eta=2, tol = 0.05, min=0, max=3)
pub = m.Potential.soft_sphere(kappa=400, epsilon=0, r0=0.5, \
eta=2, tol = 0.05, min=0.001, max=1.5)
py = m.Potential.soft_sphere(kappa=300, epsilon=25, r0=1, \
eta=2, tol = 0.04, min=0, \
max=10, shift=True)
rforce = m.forces.random(0, 1)
m.bind(rforce, C.B)
m.bind(pb, C.B, C.B, bound=True)
m.bind(pub, C.B, C.B, bound=False)
m.bind(py, Yolk, C.B)
m.Simulator.irun()
The complete simulation script is here, and can be downloaded here:
Download: this example script:
Cytoplasmic Particles¶
The motion of particles inside a cell is a very interesting scientific
problem. We have a number of ways of representing particles bound to a region,
we introduce the simplest method here. Mechanica provides a well
potential that creates a hard potential wall, and constrains particles inside
it. We can use the square well potential to model receptors in a cell’s
cytoplasm.
Start with the usual model setup:
import mechanica as m
import numpy as np
# potential cutoff distance
cutoff = 9
receptor_count = 10000
# dimensions of universe
dim=np.array([20., 20., 20.])
center = dim / 2
# new simulator, don't load any example
m.Simulator(dim=dim, cutoff=cutoff, cells=[4, 4, 4], threads=8)
Make a couple particle types, one to represent a nucleus, and another for receptors. We want the nucleus fixed in space, so make it very heavy:
class Nucleus(m.Particle):
mass = 500000
radius = 1
class Receptor(m.Particle):
mass = 0.2
radius = 0.05
target_temperature=1
#dynamics = m.Overdamped
Create a array of random positions within a solid sphere for the initial receptor positions:
# locations of initial receptor positions
receptor_pts = m.random_point(m.SolidSphere, receptor_count) * 5 + center
Make a well potential, and a soft_sphere potential. The well
confines the receptors to within a fixed distance of the nucleus, and the soft
sphere potential gives the receptors a definite shape, and prevents overlapping:
pot_nr = m.Potential.well(k=15, n=3, r0=7)
pot_rr = m.Potential.soft_sphere(kappa=15, epsilon=0, r0=0.3, eta=2, tol=0.05, min=0.01, max=1)
# bind the potential with the *TYPES* of the particles
m.bind(pot_rr, Receptor, Receptor)
m.bind(pot_nr, Nucleus, Receptor)
Create a random force (Brownian motion), zero mean of given amplitude, this gives the receptors some random motion.
tstat = m.forces.random(0, 3)
# bind it just like any other force
m.bind(tstat, Receptor)
Position the nucleus at the middle, and loop over the receptor points, and make a receptor, and finally run the simulation:
n=Nucleus(position=center, velocity=[0., 0., 0.])
for p in receptor_pts:
Receptor(p)
# run the simulator interactive
m.Simulator.run()
The square well binds particles to a region of space
The complete simulation script is here, and can be downloaded here:
Download: this example script
Plotting Simulation Metrics¶
stuff:
import mechanica as m
import numpy as np
import math
import matplotlib
import matplotlib.pyplot as plt
from time import time
from random import random
print ( matplotlib.__version__ )
# potential cutoff distance
cutoff = 10
# number of particles
count = 6000
# number of time points we avg things
avg_pts = 3
# dimensions of universe
dim=np.array([30., 30., 30.])
center = dim / 2
# new simulator, don't load any example
m.Simulator(dim=dim, cutoff=cutoff, integrator=m.FORWARD_EULER,
dt=0.002)
clump_radius = 2
# number of bins we use for averaging
avg_bins = 3
prev_pos = np.zeros((count, 3))
avg_vel = np.zeros((count, 3))
avg_pos = np.zeros((count, 3))
avg_index = 0;
# the output display array.
# matplotlib uses (Y:X) axis arrays instead of normal (X:Y)
# seriously, WTF matplotlib???
display_velocity = np.zeros((100, 200))
display_velocity_count = np.zeros((100, 200))
# keep a handle on all the cells we've made.
cells = count * [None]
def cartesian_to_spherical(pt, origin):
"""
Convert a given point in cartesian to a point in spherical (r, theta, phi)
relative to the origin.
"""
pt = pt - origin
r = np.linalg.norm(pt)
theta = math.atan2(pt[1], pt[0])
if theta < 0:
theta = theta + 2 * math.pi
phi = math.acos(pt[2] / r)
return np.array((r, theta, phi))
def spherical_index_from_cartesian(pos, origin, theta_bins, phi_bins):
"""
calculates the i,j indexex in a 2D array of size irange,jrange
from a given cartesian coordinate and origin.
theta_bins is number of bins in the theta direction
phi_bins is number of bins in the phi direction
"""
# get spherical coordinates in (r, theta, phi)
sph = cartesian_to_spherical(pos, origin)
# i is the horizontal axis, theta axis, and j is vertica, the phi
i = math.floor(sph[1] / (2. * np.pi) * theta_bins)
j = math.floor(sph[2] / (np.pi + np.finfo(np.float32).eps) * phi_bins)
return (i, j)
class Yolk(m.Particle):
mass = 500000
radius = 6
class Cell(m.Particle):
mass = 5
radius = 0.25
target_temperature=0
dynamics = m.Overdamped
total_height = 2 * Yolk.radius + 2 * clump_radius
yshift = total_height/2 - Yolk.radius
cshift = total_height/2 - 1.3 * clump_radius
pot_yc = m.Potential.soft_sphere(kappa=300, epsilon=50, r0=1, \
eta=2, tol = 0.03, min=0.1, max=8, shift=True)
pot_cc = m.Potential.soft_sphere(kappa=600, epsilon=0.5, r0=1, \
eta=3, tol = 0.05, min=0, max=2.5, shift=True)
# bind the potential with the *TYPES* of the particles
m.Universe.bind(pot_yc, Yolk, Cell)
m.Universe.bind(pot_cc, Cell, Cell)
# create a random force. In overdamped dynamcis, we neeed a random force to
# enable the objects to move around, otherwise they tend to get trapped
# in a potential
rforce = m.forces.random(0, 1)
# bind it just like any other force
m.bind(rforce, Cell)
yolk = Yolk(position=center-[0., 0., yshift])
for i, p in enumerate(m.random_point(m.SolidSphere, count)):
pos = p * clump_radius + center+[0., 0., cshift]
cells[i] = Cell(position=pos)
def calc_avg_pos(e):
global avg_index
print("calc_avg_pos, index: ", avg_index)
for i, p in enumerate(cells):
avg_vel[i] += p.position - prev_pos[i]
prev_pos[i] = p.position
avg_pos[i] += p.position
if avg_index == (avg_bins - 1):
avg_pos[:] = avg_pos / avg_bins
for i in range(count):
# get the theta / phi index from the cartesian coordinate
# remeber, matplotlib is backwards and wants matricies in
# transposed order.
ii, jj = spherical_index_from_cartesian(avg_pos[i], \
yolk.position, \
display_velocity.shape[1], \
display_velocity.shape[0])
# counts of samples we have for this spherical coordinate
display_velocity_count[jj, ii] += 1
# velocity of the vertical (y) direction
display_velocity[jj, ii] += avg_vel[i][1] / avg_bins
display_velocity[:] = display_velocity / avg_bins
Z = display_velocity
plt.pause(0.01)
plt.clf()
#plt.contour(Z)
yy = np.linspace(0, np.pi, num=100)
xx = np.linspace(0, 2 * np.pi, num=200)
c=plt.pcolormesh(xx, yy, Z, cmap ='jet')
plt.colorbar(c)
plt.show(block=False)
avg_pos[:] = 0
avg_vel[:] = 0
display_velocity_count[:] = 0
display_velocity[:] = 0
# bump counter where we store velocity info to be averaged
avg_index = (avg_index + 1) % avg_bins
m.on_time(calc_avg_pos, period=0.01)
# run the simulator interactive
m.Simulator.run()
The complete simulation script is here, and can be downloaded here:
Download: this example script:
Concepts¶
Running a Simulation¶
In most cases, running simulation is as simple as initializing the simulator
with the init function, building the physical model, and calling either the
run or irun methods to run the
simulation. However sometimes users might want more control.
A running simulation has two key components: interacting with the operating
system, and listening for user input, and time stepping the simulation (model)
itself. Whenever the run() or irun() are
invoked, they automatically start time-stepping the simulation. Users can
however explicitly control the time-stepping of the model directly. To display
the window, and start the operating system message loop, you can call the
show() method. This works just like the MatPlotLib show, in that
it displays the windows, but does not time step the simulation. The
Universe.start(), Universe.step(), Universe.stop() methods
start the universe time evolution, perform a single time step, and stop the time
evolution. If the universe is stopped, you can simply call the
Universe.start() method to continue where it was stopped. All methods to build
and manipulate the universe are available either with the universe stopped or
running.
With run(), the function will open the main window, run the simulation, and
will return when the window closes. With irun(), it will open the window,
run the simulation, and return right away, and users can then interact with
simulation objects whilst the simulation is running. Note, only works in the
command line ipython.
For convince, all the the running / showing / closing methods are aliased as top-level methods, i.e.:
>>> import mechanica as m # import the package
>>> m.init() # initialize the simulator
>>> # create the model here
>>> ...
>>> m.irun() # run in interactive mode (only for ipython console)
>>> m.run() # display the window and run
>>> m.close() # close the main window
>>> m.show() # display the window
Building A Model¶
The first step in formalizing knowledge is writing it down in such a way that it has semantic meaning for both humans and computers. This section will cover the key concepts Mechanica provides that enable users to build models of physical things.
The two key concepts we cover here are objects and processes. Objects are logical representations of physical matter. We use a particle to represent either individual things such as cells, molecules, etc, or particles can represent clumps of matter, such a volume of fluid. Processes are the ways objects interact with each other. Here we will cover the basic interaction potentials that we provide, and will also cover reactions, fluxes and events.
Making Things Move¶
Newton’s first law states that an object either remains at rest or continues to move at a constant velocity, unless acted upon by a force. That is true regardless if we are considering atoms or galaxies. In order make any object in Mechanica move, we must apply a force to it. To make objects move, Mechanica sums up all of the forces that act on an object, and uses that to calculte the object’s velocity and position.
The nature of forces in Mechanica are incredibly flexible, but we provide a variety of built-in forces to enable common behaviors.
Conservative forces are usually a kind of Potential object, where the
force is described in terms of it’s potential energy function. Long-range,
fluid, and most bonded interactions are examples of conservative potential
energy function based forces. All potential based forces contribute to the total
potential energy of the system, and we can read the total potential energy
either via the Universe.potential_energy attribute, or we can also read
the potential energy of all objects of a type, via the type’s
potential_energy attribute.
We make it easy to create forces, and apply them to objects:
# create a potential, for a simple lennard-jones fluid:
fluid_potential = Potential.lennard-jones-12-6(…)
# bind it to ALL types
m.bind(fluid_potential, Particle, Particle)
This example creates a simple potential, and binds it to ALL objects. As all
objects in our modeling world are either an instance of the base Particle
type, or a instance of a subclass of it.
Coordinate Systems¶
We can access the cartesian coordinates of any particle via the
Particle.position() method:
>>> p = MyParticleType.items()[0]
>>> print(p.position)
Out[5]: array([8.147237 , 1.3547701, 9.0579195], dtype=float32)
For spherical coordinates, we use the spherical coordinates \((r, \theta, \phi)\) as often used in mathematics: radial distance \(r\), azimuthal angle \(\theta\), and polar angle \(\phi\).
We can get the spherical coordinates of any particle via the
Particle.spherical_position() method, or for a list of partices, we can
call the spherical_positions() method on the list:
>>> In [3]: particles = Cell.items()
>>> In [4]: print(particles.spherical_positions())
[[11.88918877 1.38879788 0.87157094]
[13.48884487 1.18363881 1.71954644]
[13.44023037 -0.2029787 1.40948272]
...
[10.63551331 -1.95453715 1.43269634]
[13.24325752 -1.98586023 1.21553445]
[12.08293533 -0.55892485 1.20240927]]
Controlling Temperature¶
Right now, I have the concept of a ‘Potential’, these are objects that are specified in terms of potential function, and internally, the integrator does a bit of magic with them, and uses them calculate the conservative force that gets added to the total force. Things like bonds, angles, long-range non-bonded forces are all specified in terms of potentials. This works great for conservative forces, and is numerically actually faster then specifying a force function directly. Also, but specifying conservative forces as a potential, that lets me have both a ‘potential_energy’ and a ‘kinetic_energy’ attributes on the universe (and also the object type, i.e. if a user creates a ‘MyParticleType’, they can call MyParticleType.kinetic_energy and this returns the total kinetic energy of all objects of this type).
However, for non-conservative forces, like temperature, friction, etc, these are almost always defined as forces. We can associate a potential energy with a conservative force, but not a non-conservative (or random) force.
That would imply that we need have to allow the user to represent both potentials and forces. I would have preferred to just work in potential or forces, as this simplifies the things for the users, but I don’t really see a way around it.
import mechanica as m
import numpy as np
# potential cutoff distance
cutoff = 1
# dimensions of universe
dim=[10., 10., 10.]
# new simulator
m.init(dim=dim)
# create a potential representing a 12-6 Lennard-Jones potential
# A The first parameter of the Lennard-Jones potential.
# B The second parameter of the Lennard-Jones potential.
# cutoff
pot = m.Potential.lennard_jones_12_6(0.275 , cutoff, 9.5075e-06 , 6.1545e-03 , 1.0e-3 )
# create a particle type
# all new Particle derived types are automatically
# registered with the universe
class A(m.Particle):
mass = 39.4
target_temperature = 10000
radius=0.25
class B(m.Particle):
mass = 39.4
target_temperature = 0
radius=0.25
# bind the potential with the *TYPES* of the particles
m.Universe.bind(pot, A, A)
m.Universe.bind(pot, A, B)
m.Universe.bind(pot, B, B)
# create a thermostat, coupling time constant determines how rapidly the
# thermostat operates, smaller numbers mean thermostat acts more rapidly
tstat = m.forces.berenderson_tstat(10)
# bind it just like any other force
m.Universe.bind(tstat, A)
m.Universe.bind(tstat, B)
size = 10000
# uniform random cube
positions = np.random.uniform(low=0, high=10, size=(size, 6))
velocities = np.random.normal(0, 0.1, size=(size,6))
for (pos,vel) in zip(positions, velocities):
# calling the particle constructor implicitly adds
# the particle to the universe
A(pos[:3], vel[:3])
A(pos[3:], vel[3:])
# run the simulator interactive
m.Simulator.run()
The complete simulation script is here, and can be downloaded here:
Download: this example script:
So, user experience would be like this:
# create a thermostat force, effectively maintains the temperate of a set of things thermostat = Force.langevin_thermostat(298)
# bind it to all objects of type MyParticle m.bind(thermostat, MyParticle)
# create a friction force friction = Force.friction(…)
# bind it to all objects of type SomeOtherParticle m.bind(friction, SomeOtherParticle)
Binding¶
Binding objects and processes together is one of the key ways to create a Mechanica simulation. Binding is a very generic concept, but essentially it serves to connect a process (such a potential, flux, reaction, etc..) with one or more objects that that process should act on.
Binding Families of Objects¶
Binding Individual Objects¶
Particles¶
Universe¶
Collisions / Reactions¶
Adjective Materials / Diffusive / Dissolved Chemicals¶
Interacting With The Operating System¶
The Simulator is manages all of the interaction between the operating
system, end user input, external messaging and the physical model (which resides
in the Universe object.
In order to display a window (s), receive user input, and listen for external
messages, the simulator needs to run an event loop. These are handled by the
Simulator.run() and Simulator.irun() methods.
Types and Subtypes¶
Finding and Accessing Objects¶
One of the more comming user activities is finding, accessing and interacting with system objects after we’ve created them.
Most mechanica methods that return a list of particles actually return a
specialized list called a ParticleList. This is a special list that can
only contain particle derived types and has number of convenience methods for
dealign with spatial information.
Suppose we want the average position, or average velocity for a list of particles, then we can simply:
>>> parts = A.items()
>>> type(parts)
ParticleList
>>> parts.positions()
Out[5]:
array([[10.97836208, 7.98962736, 16.90347672],
[ 9.46043396, 4.44753504, 17.40228081],
[13.12018967, 11.84001255, 6.71417236],
...,
[13.93455601, 4.16581154, 4.48115969]])
>>> parts.positions().mean(axis=0)
Out[6]: array([ 9.9923632 , 10.01337742, 9.92124116])
Each Particle derived type has an Particle.items() method on the
type that returns all of the objects of that type:
class MyType(m.Particle):
...
# make a few instances...
MyType()
# get all the instances of that type:
parts = MyType.items()
We can acess ALL of the particles in the entire simulation via the
universe.particles() function.
Frequently we might want to grab all the objects in a grid, for example, if we
want to display some quantity as a function of spatial position. We can use the
universe.grid() function to get all the particles binned in on a regular
grid. Each element in the returned 3D array is a particle list of the particles
at that grid site. For example, if we wanted a 10x10x10 grid, of particles, we
would:
parts = universe.grid([10, 10, 10])
Finding Neighbors¶
Each object spatial object is aware if it’s neighbors, and we can get a list of
all the neghbors of an object by calling the Particle.neighbors method.
The neighbors method accepts two optional arguments, distance and types. Distance is the distance away from the surface of the present particle to search for neibhbors, and types is a tuple of particle types to restrict the search to. For example, to search for all objects of type A and B a distance of 1 unit away from a particle p, we would:
>>> nbrs = p.neighbors(distance=1, types=(A, B))
>>> print(len(nbrs))
Boundary Conditions¶
We can specify a variety of different boundary conditions via the bc argument to the Simulator constructor. We offer a range of boundary condition options in the Boundary Condition Constants section.
We specify boundary conditions via the bc argument to the main
init() function call. Boundary conditions can be one of the simple kinds
if we use the numeric argument from the Boundary Condition Constants, or can be a dictionary to use flexible boundary conditions.
For flexible boundary conditions, we pass a dictionary to the bc
init() argument like:
m.init(bc={'x':'periodic', 'z':'no_slip', 'y' : 'periodic'})
The top-level keys in the bc dictionary can be either "x", "y", "z", or
"left", "right", "top", "bottom", "front", or "back". If
we choose one of the axis directions, "x", "y", "z", then
boundary condition is symmetric, i.e. if we set "x" to some value, then both
"left" and "right" get set to that value. Valid options for axis
symmetric boundaries are "periodic", "freeslip", "noslip" or
"potential".
Note, as a user convenience, we include some synonym, i.e. "no_slip" and
"noslip" are equivlent, and "freeslip" and "free_slip" are
equivalent.
Periodic¶
The periodic is one of the most widely used boundary condition types in general. This effectively simulates an infinite domain where any agents that leave one side automatically go back in the other. Also, any agents near a boundary will interact with the agents of the opposite boundary, i.e. if there is a repulsive interaction, the agent on the left boundary will be repulsed by an agent on the right boundary. Periodic boundary conditions also determine how chemical fluxes in the Spatial Transport and Flux section operate. With periodic boundaries will, flux will interact with agents on the opposing boundary, but with non-periodic, fluxes will bounce back, or not see any agents on the opposing boundary.
We can specify periodic boundary conditions using either the 'periodic'
dictionary argument, or via numeric constants from the Boundary Condition Constants section, i.e.
m.init(bc={'x':'periodic', 'y':'no_slip', 'z' : 'periodic'})
Freeslip and NoSlip¶
Free-slip and no-slip boundary conditions reflect impacting particles back into
the simulation domain. Free-slip boundaries are essentially equivalent to a
boundary moving at the same tangential velocity with the simulation
objects. We can think it as each impacting agent collides with an equivlent
ghost agent with the same tangent velocity. No-slip boundaries are equivalent to
a stationary wall, in that impacting partricles bounce straight back, inverting
their velocity. We can specify these using the 'freeslip' or 'noslip'
options to the boundary condition dictionary.
Velocity¶
No-slip boundary conditions are just a specialization of velocity boundaries with zero velocity. Velocity boundaries model a simulation domain with a moving wall boundary. We define a velocity boundary with a nested dictionary to the individual boundary argument, i.e.:
m.init(bc={'top':{'velocity':[-1, 0, 0]},
'bottom':'noslip',
'x':'noslip',
'y':'noslip'})
Potential¶
We can attach a potential to most boundary conditions, and we define a
potential only boundary condition using "potential" option, and then bind a
potential to the boundary just like we bind potentials between agent types:
m.init(bc={'left':'potential', 'y':'no_slip', 'z' : 'periodic'})
pot = m.Potential.coulomb(q=1, tol=0.0001, min=0.05)
class MyAgent(m.Particle): pass
m.bind(pot, MyAgent, m.universe.boundary_conditions.left)
In this example, we define the left boundary to be a potential type, then add a coulomb potential between that boundary and our agent type, this any instance of that agent will be repulsed by the boundary whenever it is near it.
Potential boundary condtions default to 'free_slip', but we can also create
other kinds of potentials, such as say a 'velocity' and bind a potential to
that.
Reset¶
The reset boundary condition is a special form of a periodic boundary condition. The basic idea is we have two separate movment processes going on: we have the physical parcels of space (agents) moving around, but they also carry with them a chemical cargo. We have advection (movment of the physical parcels), and diffusion (movment of chemical cargo between agents).
The reset boundary conditions enables us to model driven, advective flow through a domain, where we reset the chemical cargo at the domain boundaries. This is usefull if we model flow where we have chemical sources or sinks in the simulation domain, and do not want any secreted material to re-enter the domain. In effect, this enables us to model an section of domain with driven flow, where the chemical state of the entering flow is always the same.
We enable reset peridic boundaries by setting the boundary condition value to a list of the keywords “periodic” and “reset”, i.e.:
m.init(dt=0.1, dim=[15, 5, 5], cutoff = 3,
bc={'x':('periodic','reset')})
The reset style boundary condition enables the fluid parcels to pass through periodic boundary conditions normally, howerver it re-sets any attached chemical cargo to the their initially specified values.
A complete example of reset boundary conditons is here with two particles. Here we have two particles, and we can see that the chemical cargo of the moving particle gets reset every time it passes through the periodic domain boundary:
import mechanica as m
m.init(dt=0.1, dim=[15, 5, 5], cutoff = 3,
bc={'x':('periodic','reset')})
class A(m.Particle):
species = ['S1', 'S2', 'S3']
style = {"colormap" : {"species" : "S1",
"map" : "rainbow",
"range" : "auto"}}
m.flux(A, A, "S1", 2)
a1 = A(m.universe.center - [0, 1, 0])
a2 = A(m.universe.center + [-5, 1, 0], velocity=[0.5, 0, 0])
a1.species.S1 = 3
a2.species.S1 = 0
m.run()
Clusters¶
Clusters and Hierarchical Orgianization¶
A Cluster is a special kind of particle that contains other particles, including more clusters.
Defining Clusters¶
We can create an instance of a cluster simply via…
Clusters are most usefull when they contain nested particle types, We define cluster with embedded types with standard python syntax:
import mechanica as m
class MyCell(m.Cluster):
class A(m.Particle)
class B(m.Particle)
This defines a new cluster type that contains two nested particle types, A and B.
Creating Clusters¶
Simply:
c = MyCell()
# creates an instance of the MyCell.A particle type in the cluster.
c.A()
# same for B
c.B()
Defining Interactions¶
Create a potential just like for free particles:
pot = m.Potential.soft_sphere(kappa=20, epsilon=5.0, \
r0=0.2, eta=4, tol=0.01, min=0.01, max=0.5)
We allow interactions between particles that belong to a cluster, and between particles belonging to different clusters, use the bound argument to the bind fuction:
m.Universe.bind(pot, MyCell.A, MyCell.A, bound=True)
Or to create an interaction potential that is between ALL instances of the particle type, reguardless of wether or not it is in cluster, simply set bound to false:
m.Universe.bind(pot, MyCell.A, MyCell.A, bound=False)
Stuff¶
Splitting and Cleavage¶
With the Cluster concept previously introduced, we can do two very interesting things with them: Splitting and Cleavage. We define a split, we start with a cluster, and split it into two clusters, the original cluster, and a new daughter cluster, and we proportion the contents of the original cluster between itself, and the new dauger cluster. Cleavage on the other hand occurs in-place in a cluster, and divides the contents of the cluster, and creates new clusters contained within the orignal one.
The split operation is typically used to apply a cell division process, where single cell splits into itself, and a new daughter. The cleavalge operation models the early embryo cleavage processes.
Splitting¶
The Cluster.split() method splits a given cluster into itself, and a new
daughter cluster. Split accepts optional normal and point arguments to define a
cleavage plane. If only a normal vector is given, and no point, the split
uses the center of mass of the cluster as the point. For example, to split a
cluster along the x axis, we can:
c = MyCell(...)
d = c.split(normal=[1., 0., 0.])
where d is the new daughter cell. This form uses the center of mass of the cluster as the cleavage plane point. We can also specify the full normal/point form as:
d = c.split(normal=[x, y, z], point=[px, py, pz])
If no named arguments are given, split interprets the first argument as a cleavage plane normal:
c.split([x, y, z])
We frequently want to split a cell along an axis, i.e. generate a cleavage plane coincident with an axis. In this case, we use the optional axis argument. Here, the split will generate a cleavage plane that contains the given axis:
c.split(axis=[x, y, z])
The default version of split uses a random cleavage plane that intersects the cell center:
d = c.split()
Split can also randomly pick half the particles in a cluster, and assign them to a new cluster with the random agument as in:
c.split(random=True)
Cleavage¶
The Cluster.cleave() method is the single entry point for all cleavage
operations, and it supports a lot of different options.
Cleavage is built on the concept of embryo cleavage. Embryology introduced some rather interesting terminolgy that we adopt here. First off, embryologists have a concept called an animal and vegatble poles. We represent these concepts using the built-in base Particle.orientation vector, and here we use this vector to represent a direction from the vegatable pole to the animal pole. All of the forms of the Cluster.cleave method will use the cluster’s orientation vector in this way. So, setting a cluster’s orientation vector, if the object is used to represent an embryo sets the direction of the pole. All of the contined cells can then refer to this property to get their orientation, i.e.:
embryo = MyEmbryo(orientation=[1.0, 0.0, 0.0], ...)
creates an new embryo with it’s pole aligned along the positive x axis.
The Cluster.cleave method is designed to perform the kinds of operations in:
Radial Cleavage¶
We perform a radial cleavage using the kind=RADIAL argument, and a number that specifies what level of cleavage to perform. As in the above figure, level 1 splits a single cluster into two, along the orientation axis, level 2 splits it into 4 new clusters, orientated along the orietation axis, level 3 splits it into 8 new clusters, usign a cleavage plane as in
For example, if we create an embryo, and call:
embryo = MyEmbryo(orientation=[1.0, 0.0, 0.0], ...)
embryo.cleave(kind=RADIAL, level=3, cell_type=MyEmbryo.Basic)
this creates eight new clsuters inside the parent embryo cluster, and assigns all of them the MyEmbryo.Basic type. Or optionallly, we may use any Cluster derived type. If the cell_type option is left blank, the new cluster default to the top level Cluster type.
Subsequent Radial Cleavage¶
We don’t have to perform the cleavage in a single step, rather we can call:
embryo.cleave(kind=RADIAL, cell_type=MyEmbryo.Basic)
without specifing a level, in this case, ‘cleave` looks at the current number of contained clusters, and perfoerms the next pre-defined cleavage operations. For example, if the embryo only contained a single cell, then it pefrorms the first cleavage operation. If there are two cells, it splits them along the orientation axis, if there are 4, then it splits them perpendicular to the orientation axis, as in the previous figure.
Potentials and Effect On Mechanics¶
One of the main goals of Mechanica is to enable users to rapidly develop and explore empirical or phenomenological models of active and biological matter in the 100nm to multiple cm range. There are no pre-defined force fields in this length scale, as such users need a good deal of flexibility to create and calibrate potential functions to model material rheology and particle interactions.
Mechanica provides a wide range of potentials in the Potential class, we can create any of these potential objects using the static methods in this class.
Creating, Plotting and Exploring Potentials¶
We create a potential object simply by calling one of the static methods on the Potential class. Potential objects conveniently have a plot method on them that displays a ‘matplotlib’ plot, for example:
p = m.Potential.glj(10)
p.plot()
results in
Out:
<Potential at 0x7fb7f6110210>
The Generalized Lennard-Jones Potential, Potential.glj is one of the
most potential we offer, and supports a wide range of parameters. The glj
potential has the form
where \(r_e\) is the effective radius, which is automatically computed as the sum of the interacting particle radii.
In most traditional molecular dynamics packages, it is important to specify the
r0 parameter in a potential, however, Mechanica automatically picks up the rest
length from the particle radius, For example, say we have two particles, one
with radius 1 and the other with 2, then, if these particles interact with a
potential like glj, they will automatically have a rest separation
distance of 3. We can make a simple model like to illustrate this:
class B(m.Particle):
mass = 1
dynamics = m.Overdamped
# make a glj potential, this automatically reads the
# particle radius to determine rest distance.
pot = m.Potential.glj(e=1)
m.bind(pot, B, B)
p1 = B(center + (-2, 0, 0))
p2 = B(center + (2, 0, 0))
p1.radius = 1
p2.radius = 2
m.Simulator.run()
The two particles start out a with a little separation, but the potential pulls them in such that they just touch at the specified radii, and we get an output of:
This complete model is can be downloaded here:
glj-1.py
The Mechanica runtime automatically shifts the potential, such that the minimum of the potential will always be when the separation distance between particles is \(r_1 + r_2\) where \(r_1\) and \(r_2\) are the radii of the two interacting particles.
The most important parameters to glj are the e, r0, m, and min
parameters. The e represents the strength of the potential, it is the depth of
the potential well, and ultimately determines how strongly objects stick
together. r0 is the the lowest point of the potential, and m and n are the
exponents which determine how sharp or shallow the potential is.
e : potential energy¶
Appropriate values of e should be close to the mass of an object, if e is much larger than about two to three times the mass at most, as this will result in unrealistic behavior and numeric instability. e however can be small as needed to create weakly bound or interacting objects.
r0 : effective radius¶
The r0 parameter is the effective radius, and determines the equilibrium distance of the potential.
Mechanica however automatically shifts the potential relative to particle radius, so r0 doesn’t determine separation distance as with a traditional molecular dynamics package. Rather r0 influences how soft the potential is. Smaller values of r0 lead to a stiff potential, with a harder shoulder, i.e. it makes the particles harder. Whereas larger values of r0 spread the potential well out over a larger distance, and also allows the potential to be longer ranging, and attract objects from longer distances, as in (A) of Fig. 12.
m : potential exponent¶
The m and n parameter are the exponents of the potential, and these are principally responsible for how soft or hard the potential is. Small values of m results in a softer, potential, and also allow a longer ranging potential, much like the r0 parameter. It is normally not required to specify the n parameter, as the potential constructor automatically sets this to \(2 * m\), because this format is widely used in many meso-scale particle dynamics based simulations, and results in well-behaved numerics and realistic physical behavior.
(A) Varying the \(r_0\) parameter in a glj potential, while holding all other parameters constant, with e=1, m=3, n=6. (B) Varying the m parameter in a glj potential, holding e=1, \(r_0=1\) and with m=2n.
max : potential cutoff distance¶
The potential cutoff distance, max has a very strong influence on material properties and dynamics. Longer cutoff distances enable us to model things such as an membrane pulling longer distance particles inwards, or long-range coulomb or magnetic potentials, or possibly some more phenomenological interactions such as boids, however longer ranged potentials are typically not very realistic in mesoscopic materials.
In mesoscopic and biological materials, we typically have an interaction that has a short range repulsion – that keeps things from being pushed together, and also a short-range attraction that allows object near each other to stick together, and also lets us model things like fluids, where they like to stick to each other. If we allow longer-range potentials, such as in Fig. 13, we end up with a situation where the particles that make up an object pull so hard on other particles that it causes the object to collapse down to a point. This is especially true when we use a softer potential that does not provide as much short-range repulsion.
For example, we can construct the model that’s shown in Fig. 13
by creating a Cluster type, adding a few particle types to it, and
creating two different potentials to act on them as:
class C(m.Cluster):
radius=3
class A(m.Particle):
radius=0.5
dynamics = m.Overdamped
mass=10
style={"color":"MediumSeaGreen"}
class B(m.Particle):
radius=0.5
dynamics = m.Overdamped
mass=10
style={"color":"skyblue"}
c1 = C(position=center - (3, 0, 0))
c2 = C(position=center + (7, 0, 0))
c1.A(2000)
c2.B(2000)
p1 = m.Potential.glj(e=7, m=1, max=1)
p2 = m.Potential.glj(e=7, m=1, max=2)
m.bind(p1, C.A, C.A, bound=True)
m.bind(p2, C.B, C.B, bound=True)
rforce = m.forces.random(0, 10)
m.bind(rforce, C.A)
m.bind(rforce, C.B)
Here we created two very soft glj potentials, we made them soft because
we used m=1 to set the exponent, as discussed earlier. We set one potential to
a cutoff of 1 (twice the particle radius) and the other to a cutoff of 2 (four
times the particle radius). The longer range potential pulls on many more
particles, and the net attraction of these overcomes the relatively soft short
range repulsion. The complete code to create Fig. 13 can be
downloaded here:
glj-1.py
To overcome the this issue of particle clumps collapsing, or overlapping particles, we recommend cutoff distances of around 2-3 times the smallest particle radius. The max cutoff distance defaults to three times the effective radius, r0, we found that this distance generally works well, so users typically would not need to explicitly set this parameter.
Two clusters, both with identical particles, and identical initial conditions. The left cluster (green) has a shorter cutoff (twice particle radius), on the potential, ad the right one (blue) has a longer cutoff (four times particle radius). (A) Initializing two clusters of identical particles, both start out the same size. (B) After a few seconds of simulation time, the left
Metrics and Derived Quantities¶
Mechanica provides numerous methods to compute a range of derived quantities. Some of these are top-level metrics, and depend on the entire simulation volume, and others are localized to individual objects. Simulator
Pressure and Virial Tensors¶
For a system of N particles in a volume V, we can compute the surface tension from the diagonal components of the pressure tensor \(P_{\alpha,\alpha}(\alpha=x,y,z)\). The \(P_{xx}\) components are:
where \(N\) is the number of particles, \(\rho\) is the particle density density, \(k\) Boltzmann constant, \(T\) is the temperature, \(\mathbf{r}_{ij}\) is the vector between the particles \(i\) and \(j\), and \(\mathbf{f}_{ij}\) is the force between them. Some important concepts here is that the force here is only between the particles used in the pressure tensor calculation, it specifically excludes any external force. The pressure tensor here is a measure of how much internal force exists in the specified set of particles.
We comonly refer to the right hand side above as the virial, it represents half of the the product of the stress due to the net force between pairs of particles and the distance between them. We formally define the virial tensor components as
The volume of a group of particles is not well defined, as such we separate out computing the virial component, and the volume, and give users the flexiblity of using different volume metrics.
We provide a number of different options for calculating the virial
tensor. You can compute the pressure tensor for the entire simulation domain, or
a specific region using the Universe.virial() method. Can compute the
pressure tensor for a specific cluster using the Cluster.virial() method,
or can compute the tensor at a specific particle location using
Particle.virial() method.
Radius of Gyration¶
In the radius of gyration is measure of the dimensions of a group
(Cluster) of particles such as a polymer chain, macro-molecule or some
larger object. The radius of gyration of group of particles at a given time is
defined as:
We can compute the radius of gyration for a cluster of particles using the
Cluster.radius_of_gyration() method.
Center of Mass¶
The center of mass of a system of particles, \(P_i, i-1, \ldits, n\), each with mass \(m_i\), at locations \(\mathbf{r}_i, i=1, \ldots, n\), with the center of mass \(\mathbf{R}\) satisfy the condition
with \(\mathbf{R}\) defined as:
where \(M\) is the sum of the masses of all of the particles.
We can compute the center of mass of a cluster particles with the
Cluster.center_of_mass() method.
Center of Geometry¶
Computes the geometric center of a group of particles with the
Cluster.center_of_geometry() method, or equivalently, with the
Cluster.centroid() method.
Moment of Inertia¶
For a system of \(N\) particles, the moment of inertia tensor is a symmetric tensor, and is defined as:
Its diagonal elements are defined as
and the the off-diagonal elements, also called the are:
We can compute the inertia tensor for a group of particles using the
Cluster.moment_of_inertia() or Cluster.inertia() methods.
Changing Type¶
One of the hallmarks of biological cells, and biological objects in general is that they change type over time. For example a XXX cell change change type to become a YYY cell. One way to think about what’s going on here is that you have a physical object of a certain kind, which has a certain set of physical substances in it. Then their contents generally stays the same, but the set of processes that act on the object change, the dynamics of the object changes, although it’s physical state may not be altered significantly.
In nature, An object’s genotype is the set of genetic makeup of that object, whereas an object’s phenotype is the set of observable characteristics, which are influenced both by its genotype it’s environment.
In Mechanica, we use the concept of a Python language type analogously to a phenotype. For a particular phenotype, we have an instance of an object, and the type defines the set of interactions and processes that that object participates in.
Type change is very common in nature, however is exceedingly rarely used in compute languages. One example of a computer language with type change is Smalltalk. Smalltalk has a method called become, which is one if it’s most powerful and unique features, and is also quite obscure outside of the Smalltalk community,
Smalltalk’s become method swaps the identities of a pair of object instances:
a become: b
In Smalltalk, the the object a assumes the identity or type of b, and visa
versa. The contents of the of all object references to a and b get swapped,
so whenever any reference called a, it’s now referring to the contents of
b. This is not simply swapping a pair of variables, like we would in python as
(b,a) = (a,b), but rather every reference that pointed to the object a
now points to b. Thus in Smalltalk, both the state and type of the objects
change.
Mechanica’s become is slightly different in that we only change the type of the object, and generally leave the state of the object alone.
This of course comes with a caveat, in that say the variables defined in the respective type definitions is different, such as if one type defined a different set of Species than the the other. for example, if we have:
class A(m.Particle):
species = ['S1', 'S2', 'S3']
class B(m.Particle):
species = ['S2', 'S3', 'S4']
a = A(); b = B();
a.become(B)
We can see that types ‘A’ and ‘B’ only share species ‘S2’ and ‘S2’, but species ‘S1’ only exists in ‘A’ and ‘S4’ only exists in ‘B’. In such cases, we keep the values of the common state variables alone, destroy the species ‘S1’ (because it does not exist in ‘B’, and create a new ‘S4’ and initialize it to it’s default value (if specified in an initial assignment).
Changing the type of all the neighbors of an random object. We can do this
with a single line of python, [p.become(B) for p in x.neighbors(3)]
To create the model (Fig. 14), we first create a standard simulation, create a couple particle types, and define their size and color:
import mechanica as m
import numpy as np
m.Simulator()
class A(m.Particle):
radius=0.1
dynamics = m.Overdamped
mass=5
style={"color":"MediumSeaGreen"}
class B(m.Particle):
radius=0.1
dynamics = m.Overdamped
mass=10
style={"color":"skyblue"}
Make a potential, here we choose the coulomb potential as we’re using
periodic boundary conditions, and simply want the objects to uniformly push away
from each other, and uniformly fill space. We call the random_points
function to create a set of random points that fills a cube
p = m.Potential.coulomb(q=2, min=0.01, max=3)
m.bind(p, A, A)
m.bind(p, B, B)
m.bind(p, A, B)
r = m.forces.random(0, 1)
m.bind(r, A)
pos = m.random_points(m.SolidCube, 50000) * 10 + m.Universe.center
Create a new object of type A for every position in the list, and show the simulator:
[A(p) for p in pos]
m.show()
Now, to choose the neighbors of an object, and change their type, we simply call
the neighbors method on a randomly chosen particle. This returns a list
of objects near a, and we simply use the Python list syntax to call the
Particle.become() method on every object in the neighbor list. Simply to make the
object more visible, we set it’s radius to 2. Changing the radius of
particles that interact with a coulomb potential has no effect on the
force.
a = A.items()[0]
a.radius = 2
[p.become(B) for p in a.neighbors(3)]
You can download and run the complete model here:
Download: this example script:
Style¶
All renderable objects in Mechanica have a style attribute. This style is essentially like a CSS .style attribute in Javascript / HTML. The style object behaves like a container for a variety of style descriptors. Our intent is to eventually read and generate all of the style nodes from a CSS file, but for the time being, we have to manually set style information.
Like CSS, each instance of a object automatically inherits the style of it’s type, but users may override values.
Currently, the only supported style descriptor is a color value. This is an instance of the
Chemical Reactions, Flux and Transport¶
Unlike a traditional micro-scale molecular dynamics approach, where each computational particle represents an individual physical atom, a DPD is mesoscopic approach, where each computational particle represents a ‘parcel’ of a real fluid. A single DPD particle typically represents anywhere from a cubic micron to a cubic mm, or about \(3.3 \times 10^{10}\) to \(3.3 \times 10^{19}\) water molecules.
Transport dissipative particle dynamics (tDPD) adds diffusing chemical solutes to each classical DPD particle. Thus, each tDPD particle represents a parcel of bulk fluid (solvent) with a set of chemical solutes at each particle. In tDPD, the main particles represent the bulk medium, or the ‘solvent’, and these carry along, or advect attached solutes. We introduce the term ‘cargo’ to refer to the localized chemical solutes at each particle.
In general, the time evolution of the chemical species at each spatial object is given as:
where the rate of change of the vector of chemical species at an object is equal to the flux vector, \(Q_i\). This is the sum of the transport and reactive fluxes. \(Q^T\), is the transport flux, and \(Q^R_i\) is a local reactive flux, in the form of a local reaction network. We will cover local reaction fluxes in later paper, for now we restrict this discussion to the passive or ‘Fickian’ flux, secretion flux and uptake flux.
Before we cover spatial transport, we first cover adding chemical reaction networks to objects.
To attach chemical cargo to a particle, we simply add a species specifier to
the particle type definition as:
class A(m.Particle):
species = ['S1', 'S2', S3']
This defines the three chemical species, S1, S2, and S3 in the type
definition. Thus, when we create an instance of the object, that instance will
have a vector of chemical species attached to it, and is accessible via the
Particle.species attribute. Internally, we allocate a memory block for
each object instance, and users can attach a set of reactions to define the time
evolution of these attached chemical species.
If we access this species list from the type, we get:
>>> print(A.species)
SpeciesList(['S1', 'S2', 'S3'])
This is a special list of SBML species definitions. It’s important to note that once we’ve defined the list of species in each time, that list is immutable. Creating a list of species with just their names is the simplest example, if we need more control, we can create a list from more complex species definition strings in Species.
If a type is defined with a species definition, every instance of that type will get a StateVector, of these substances. Internally, a state vector is really just a contiguous block of numbers, and we can attach a reaction network or rate rules to define their time evolution.
Each instance of a type with a species identifier gets a species attribute, and we can access the values here. In the instance, the species attribute acts a bit like an array, in that we can get it’s length, and use numeric indices to read values:
>>> a = A()
>>> print(a.species)
StateVector([S1:0, S2:0, S3:0])
As we can see, the state vector is array like, but in addition to the numerical values of the species, it contains a lot of meta-data of the species definitions. We can access individual values using array indexing as:
>>> print(a.species[0])
0.0
>>> a.species[0] = 1
>>> print(a.species[0])
1.0
The state vector also automatically gets accessors for each of the species names, and we can access them just like standard Python attributes:
>>> print(a.species.S1)
1.0
>>> a.species.S1 = 5
>>> print(a.species.S1)
5.0
We can even get all of the original species attributes directly from the instance state vector like:
>>> print(a.species[1].id)
'S2'
>>> print(a.species.S2.id)
'S2'
In most cases, when we access the species values, we are accessing the concentration value. See the SBML documentation, but the basic idea is that we internally work with amounts, but each of these species exists in a physical region of space (remember, a particle defines a region of space), so the value we return in the amount divided by the volume of the object that the species is in. Sometimes we want to work with amounts, or we explicitly want to work with concentrations. As such, we can access these with the amount or conc attributes on the state vector objects as such:
>>> print(a.species.amount)
1.0
>>> print(a.species.conc)
0.5
Species¶
This simple version of the species definition defaults to create a set of floating species, or species who’s value varies in time, and they participate in reaction and flux processes. We also allow other kinds species such as boundary, or have initial values.
The Mechanica Species object is essentially a Python binding around the libSBML Species class, but provides some Pythonic conveniences. For example, in our binding, we use convential Python snake_case sytax, and all of the sbml properties are avialable via simple properties on the objects. Many SBML concepts such as initial_amount, constant, etc. are optional values that may or may not be set. In the standard libSBML binding, users need to use a variety of isBoundaryConditionSet(), unsetBoundaryCondition(), etc… methods that are a direct python binding to the native C++ API. As a convience to our users, our methods simply return a Python None if the field is not set, otherwise returns the value, i.e. to get an initial amount:
>>> print(a.initial_amount)
None
>>> a.initial_amount = 5.0
This internally updates the libSBML Species object that we use. As such, if the user wants to save this sbml, all of these values are set accordingly.
The simplest species object simply takes the name of the species as the only argument:
>>> s = Species("S1")
We can make a boundary species, that is, one that acts like a boundary condition with a “$” in the argument as:
>>> bs = Species("$S1")
>>> print(bs.id)
'S1'
>>> print(bs.boundary)
True
The Species constructor also supports initial values, we specify these by adding a “= value” right hand side expression to the species string:
>>> ia = Species("S1 = 1.2345")
>>> print(ia.id)
'S1'
>>> print(ia.initial_amount)
1.2345
Spatial Transport and Flux¶
Recall that the bulk or solvent particles don’t represent a single molecule, but rather a parcel of fluid. As such, dissolved chemical solutes (cargo) in each parcel of fluid have natural tendency to diffuse to nearby locations.
Dissolved solutes have a natural tendency to diffuse to nearby locations.
This micro-scale diffusion of solutes results in mixing or mass transport without directed bulk motion of the solvent. We refer to the bulk motion, or bulk flow of the solvent as advection, and use convection to describe the combination of both transport phenomena. Diffusion processes are typically either normal or anomalous. Normal (Fickian) diffusion obeys Fick’s laws, and anomalous (non-Fickian) does not.
We introduce the concept of flux to describe this transport of material (chemical solutes) between particles. Fluxes are similar similar to conventional pair-wise forces between particles, in that a flux is between all particles that match a specific type and are within a certain distance from each other. The only differences between a flux and a force, is that a flux is between the chemical cargo on particles, and modifies (transports) chemical cargo between particles, whereas a force modifies the net force acting on each particle.
We attach a flux between chemical cargo as:
class A(m.Particle)
species = ['S1', 'S2', 'S3']
class B(m.Particle)
species = ['S1, 'Foo', 'Bar']
q = m.fluxes.fickian(k = 0.5)
m.Universe.bind(q, A.S1, B.S)
m.Universe.bind(q, A.S2, B.Bar)
This creates a Fickian diffusive flux object q, and binds it between species
on two different particle types. Thus, whenever any pair of particles instances
belonging to these types are near each other, the runtime will apply a Fickian
diffusive flux between the species attached to these two particle instances.
Passive Flux: Diffusion¶
We implement a diffusion process of chemical species located at object instances
using the basic passive (Fickian) flux type, with the flux(). This flux
implements a passive transport between a species located on a pair of nearby objects of type a
and b. A Fick flux of the species \(S\) attached to object types
\(A\) and \(B\) implements the reaction:
\(B\) respectivly. \(S\) is a chemical species located at each object instances. \(k\) is the flux constant, \(r\) is the distance between the two objects, \(r_{cutoff}\) is the global cutoff distance, and \(d\) is the optional decay term.
Active Fluxes: Production and Consumption¶
For active pumping, to implement such processes like membrane ion pumps, or
other forms of active transport, we provide the produce_flux() and
consume_flux() objects.
The produce flux implements the reaction:
The consumer flux implements the reaction:
Here, the \(\left(1 - \frac{b.S}{b.S_{target}} \right)\) influences the
forward rate, where \([b.S]\) is the concentration of the substance S, and
\(b.S_{target}\) is the target concentration. The flux will continue forward
so long as there is both concentration of the reactant, \(a.S\), and
the product \(b.S\) remains below its target value. Notice that if the
present concentation of \(b.S\) is above its target, the reaction will
proceed in the reverse direction. Thus, the pumping_flux() can be used to
implement both secretion and uptake reactions.
Bonded Interactions¶
Bonded interactions are explicit bonds between particle instances.
We support the standard molecular dynamics bonds, angles, and torsions.
Bonds¶
bonds do stuff
some more text akdjklasjdk akjdalkjdsl k
aDJalsdj
LKADJalj
asdjfasdjfjkasdhdfj
fsfjhjk
Continuum Reactions or Reaction Kinetics¶
The Carbon Reaction Kinetics (CRK) is a set of data structures (types) and APIs for constructing, interacting, attachign to physical objects and visualizing reaction kinetics models.
- The CRK objects form a network and provide a way for users and developers to create and interact with RK models.
- CRK objects are similar to, and generally have SBML analogs.
- Internally implemnted as C++, have C and Python APIs.
- Represent both the API for creating networks, and for simulation and accessing values.
- Are a tree structure that the nework visulazion and visual editing tools operate on.
The libSBML API is only designed for manipulating the SBML DOM, for creating and manipulating SBML documents. It does not have any way of accessing values in a simulaiton.
Our objects here are a way for both creating and manipuating an RK network, and also for accessing the values in a simulation.
The Carbon reaction-kinetics DOM is the core of our model construction, simulation and visulazion workbench
Python Interaction¶
The CRK provides a clean, concise and very native way for users to attach RK networks to objects. Users can add ‘species’ and optionally, ‘reactions’ definitions to a Python type, to automatically generate a rection kinetics network on every instance of that type:
# create one type of particle, species default to floating.
class A(m.Particle):
species = ['S1', 'S2', 'S3']
mass = 1
# another type, with some reactions, and more complex species
class B(m.Particle):
species = {'S1' :Species(boundary=True, 1.0),
'Foo':5.0,
'Bar':10.0}
reactions = ['S1 + S1 -> S3; k*S1*S2', 'S3->0; k2*S3']
# create a flux between species at particles
f = m.fluxes.fick(k=0.5, omega='gaussian')
m.Universe.bind(f, A.S1, B.Foo)
m.Universe.bind(f, A.S3, B.Foo)
The real power here with the CRK objects is that they also provide a way for users to directly access and manipulate simulation values even whilst the simulation is running. Adding a species (or reaction, or other CRK defintion) to a Python class automatically generates accessor symbols on that python type so that the actual value if the species is direclty availabe. For example, say a user created a bunch of the aformention particles. Both of these defined a ‘S1’ species, so a user could access these values directly as:
for p in m.Universe.particles:
print(p.S1)
This works because internally, CRK objects implement the Python Descriptor Object protocal, and provide a way to add symbols to type definitions that access the underlying value of that species or object name. This very different from the libSBML API, which only provides a way to define a network but has no concept of simulation or accessing values.
Network Visialization¶
The Network Visialization Component (NetViz) and Canvas shown in (Fig. 16) interacts with both the C API of the CRK ojbcts, and the C API of the canvas component. Becaue the CRK C API supports events, when end users click or interact with any object in the canvas, the coresponding CRK object will emit an event. These events are things such as mouse over, mouse enter, etc.. The NetViz componenent is essentially the glue between the network of CRK objects and the tree of Canvas objects.
- Canvas API provides a series of methods to create visual, on-screen objects. Objects such as rectangles, text, shapes, etc..
- On-screen Canvas objects listen for user mouse events, key press events, etc and provide an API for attaching event listners to these evetns,
- The NetViz component has a reference to a set of CRK objects. It iterates over this collection, and generates a corresponding visual representation using the Canvas API.
- NetViz listens for user events from the visual Canvas objects. NetViz can then signal a corresponding event in a CRK object. This enables end-user Python code to recieve events from end-user mouse input.
Model Loading¶
End users typically create CRK models using the Python API, or the interactive NetViz component. We also support readign and writing both SBML and Antimony documents. CRK also supports it’s own human-readable and writable (unlike XML 😄) file format.
Network Object Model¶
Reaction kinetics, flux and ODE netorks are the (??? really important thing ???). Something about chemical processes…
The Network Object Model (NOM) is a core component of the Carbon library that
- Enables the visual network renderer to represent and edit reaction kinetics models
- Users to attach networks to spatial objects
- Have the visual network user interface display and interact with a live, running simulation.
- is a way to represent, store and edit both the underlying structure of a reaction network, as well as interact with the values from a running simulation
Formalizing Physical Knowledge¶
The first step in formalizing knowledge is writing it down in such a way that it has semantic meaning for both humans and computers. In order to specify biological phenomena, we need a modeling description that corresponding closely to natural phenomena.
Let’s start by looking at the folliwing diagram of a single biological cell. This cell could be moving on a surface, or inside some tissue.
A biological cell is a complex dynamic entity that participates in many active processes.
We can imedietly idenfity a number of things. Things such as cell nuclie, cell membrane, cytoplasm, actin fibers, lamellpoidium, filopodium, etc.. We call these things objects. Objects represent the “things” such as molecules, proteins, cells, fluids, or materials. We define an object as any instantiable physical or logical entity that has certain state-full or structural properties.
We also understand that for any living or active things, these objects typically not static. Cells can move around, active fibers contract, cells can excert forces on it’s neighboring cells and enviormnet, cells can sense and respond with chemical signaling, and so forth as in the following figure
Biological cells can sense a variaety of signals in their envionrment and respond with a variety of behaviors.
We call anything that causes something to change a process. Processes represent anything that causes a change in state. A process may, for example, represent a transformation in which sub-strates are consumed and products are produced, or it may alter the state of an object, either continuously or discretely. Unlike traditional programming languages, processes in Mechanica can operate concurrently and continuously.
Everything we’ve discussed also has a type or in other words, what kind of thing it is. When we look at a say a particular individual molecule, we say that this molecue is a kind or type of molecue, say water molecule, calcium ion, or a protein, etc. A type serves to identify the category that an instance of a thing belongs. to. Processes too have a type, we have have an adhesion process, a surface tension process, etc..
These notions of objects, process and types are discussed in more detail in the following sections, but for now, we have a basic understanding that we may proceeed with developing a model of this cell.
Mechanisms¶
Now that we have established a vocabulary and a formalism, lets again look at Fig. 17, and start formalizing what kinds of objects make up the cell, and developing models or abstractions to represent these concepts.
The first thing we likely notice is that the cell exists in space, it occoupies a specifiec region of space, it has a well defined (though iregulary shapped) boundary. At the molecular level of course, the cell boundary (membrane) is composed of lipid bylayer, and these lipids in turn are usually composed of a hydrocarbon chain and a terminal carboxyl group and phosphate group. The questions we are interested in are at the cell, rather than the molecular length scale, so we we develop a suitable cellular scale description of geometry, rather than a molecular description.
We could of course explicity model every single atom, as is normally done in molecular dynamics simulations, however this is simply not tractable from a computational perspespective. Rather, we will represent the surface of a cell as closed surface composed of polygons. Each polygon is connected to it’s neighboring polygons, and the boundary of a polygon is defined by a set of edges and vertices. Vertices define a single point in space, and an edge connects two vertices. We can see an example of a cell boundary defined as set of connected polygons in Fig. 19.
We can represent the surface of a cell as a set of interconnected polygons.
Mechanica represents physical concepts using two key concepts: objects and processes. This paradigm of representing the world in terms of objects and processes forms the basis of the Cellular Behavior Ontology (CBO) [SSS+14], and also the Object Process Methodology (OPM) [Dor02].
Objects¶
Objects are the nouns, the actual things being described. Objects are such things as molecules, cells, membranes, ion channels, the extra-cellular matrix, fluids, etc. Objects have quantifiable characteristics such as location, amount, concentration, mass and volume. Objects define a state; they are comparable to standard data structures in conventional programming languages. Objects may inherit and extend other objects, or may contain other objects. Objects are grouped into two categories: continuous and discrete. Continous objects describe things such as continuously valued chemical or molecular fields which have position-dependent concentrations. Chemical fields may be 3D in that they are bounded by a surface or 2D in that they exist on a surface. For reasons of numerical efficiency, users may specify fields as spatially uniform.
Processes¶
Processes are the verbs. Processes may create and destroy objects, alter the state of objects, transform objects, or consume and produce sets of objects. As in nature, multiple Mechannica processes act concurrently on objects and may act at different rates, but may only be active under certain circumstances. Processes may be continuously active, or may be explicitly triggered by specific conditions or invoked directly by some other process. Processes may also be combined or aggregated, such that a process may be hierarchically composed of child processes, and child processes may be ordered either concurrently or sequentially. Processes fall into two categories: continuous and discrete.
Types¶
Every thing that we have previosly discussed, objects and process has a well-defined type. Types are a basic concept common to most programming languages, serving to classify variable instances into categories. The type of a variable determines the kind of data which may be stored in that variable. Essentially, a type in programming languages is roughly analogous to genotype in biology. The type of an object defines what operations are valid on that instance, i.e. we can sum two numeric types, but adding a numeric type to, say, a string is ill defined. Most programming languages do not have a concept related to the biological notion of a phenotype. A phenotype in biology is a metric, an observable categorization that determines the cell type of a particular cell instance. A phenotype is defined by a set of rules, or conditions; when these conditions are met, we say that a cell is of such type.
The Mechanica extends the basic concept of dynamic or static types with a rule-based type. Here, the a type may be defined via a set of rules, and when all of these rules are met, we say that a variable instance is a certain type. This notion is important because biological cells frequently undergo phenotypic change, that is, they change type over time. Again, processes are defined to operate on a specific type. Here we can create a type definition based on a set of conditions; when all of these conditions are met, the type expression becomes true, and the processes corresponding to that type definition now automatically become active on all object instances for which the type evaluates to true.
The type of an object in a programming language must be unambiguous, since the type assigns meaning to a block of memory. A type formally defines the layout of a memory block and the types of data it stores. Strongly-typed programming languages require that we associate a type to each object or process when we define its identifier (symbol). Weakly-typed programming languages determine object and process types according to a set of rules and do not require that we declare the type of the object when we create it. In some weakly-typed languages, an object can even change type. A type system is a collection of rules that define the set of object types which can participate in a given process.
Most programming languages lack a concept corresponding to phenotype. A programming language cannot, in general, look at a memory block and determine the type of object it contains by analyzing its contents. A phenotype, on the other hand, is a list of conditions that determine how to categorize an object from the object’s properties. However both programming languages and biology have type systems. Indeed, biological modeling is the creation of a type system for a specific situation –collecting and defining a set of rules which specify which types of objects participate in which types of processes.
Mechanica extends the programming language concept of type with a rule-based definition of type. The type definition needs to be fuzzy so that we can ask how much an object instance participates in a type. Current programming languages support only Boolean type inquiries, in the Mechanica inquiring if an object is of a specified type returns a real value between 0 and 1.
Types serve to classify variable instances into categories. The type of a variable determines the kind of data that may be stored in that variable. The type of an object defines what operations are valid on instances of that object type, i.e., we can sum two numeric types, but adding a numeric type to a string is ill-defined. Most programming languages do not have a concept related to the biological notion of a phenotype. A phenotype in biology is a emph{metric}, an observable categorization that determines a cell’s type. A phenotype is defined by a set of rules or conditions such that when these conditions are met, we say that a cell is of such type.
Mechanica extends the basic concept of dynamic or static types with a rule-based type, which is related to the concept of typestate oriented programming [SY86]. Here, the type may be defined via a set of rules, and when all of those rules are met, we say that a variable instance is a certain type. This notion is important because biological cells frequently undergo phenotypic change, that is, they change type over time. Processes are defined to operate on a specific type, and the runtime automatically applies these processes to all instances of the specified type. Here we can create a type definition based on a set of conditions; when all of these conditions are met, the type expression becomes true, and the processes corresponding to that type definition now automatically become active on all object instances for which the type evaluates to true.
Forces¶
Forces are kind of process that applies a force to a physical object.
Force processes provide a way to describe spatial change, such as motion, deformation, adhesion or response to external forces. forces are similar to force functions in molecular dynamics. A force can be defined to act on one or between two spatial objects, hence a force may have one or two arguments, and both of them must be spatial object subtypes. forces return the force that acts on its arguments. Any motion processes (adhesion, taxis, deformation) can be specified via a suitable force process. For example, when an adhesion process is active between a surfaces of a pair of cells, the adhesion process applies a force between the cell surfaces at the locations where the surfaces are in contact. This adhesive force acts to keep the cells in contact and resists surface separation.
The language runtime automatically applies the force functions to spatial objects and calculates the net force acting on each spatial object. The runtime then calculates the time evolution of each spatial object, typically as $mathbf{v} propto mathbf{F}/m$, where velocity is proportional to the net force acting on each spatial object.
We provide a suite of pre-built forces, and users can easily develop their own custom forces.
Mechanica API Reference¶
This is the API Reference page for the module: mechanica
-
mechanica.init(...)¶ Initializes a simulation. All of the keyword arguments are the same as on the Config object. You can initialize the simulator via a config
Configobject , or via keyword arguments. The keywords have the same name as fields on the config.
-
mechanica.run()¶ Starts the operating system messaging event loop for application. This is typically used from scripts, as console input will no longer work. The run method will continue to run until all of the windows are closed, or the
quitmethod is called. By default,runwill automatically start the universe time propagation.
-
mechanica.irun()¶ Runs the simulator in interactive mode, in that in interactive mode, the console is still active and users can continue to issue commands via the ipython console. By default,
irunwill automatically start the universe time propagation.
-
mechanica.close()¶ Closes the main window, but the application / simulation will continue to run.
-
mechanica.show()¶ Shows any windows that were specified in the config. This works just like MatPlotLib’s
showmethod. Theshowmethod does not start the universe time propagation unlikerunandirun.Keyword Arguments: - dim (
vector3 like) – [x,y,z] dimensions of the universe, defaults to [10., 10., 10.] - bc (
number or dictionary) – Boundary conditions, use one of the constants in the Boundary Condition Constants section. Defaults toPERIODIC_FULL. For a detailed discussion of boundary condtion options, see Boundary Conditions - cutoff (
number) – Cutoff distance for long range forces, try to keep small for best performance, defaults to 1 - cells (
vector3 like) – [x,y,z] dimensions of spatial cells, how we partition space. Defaults to [4,4,4] - threads (
number) – number of compute threads, defaults to 4 - integrator (
number) – kind of integrator to use, defaults toFORWARD_EULER - dt (
number) – time step (\(\delta t\)). If you encounter numerical instability, reduce time step, defaults to 0.01
- dim (
Simulator¶
-
mechanica.version(str)¶ The current Mechanica version, as a string.
-
class
mechanica.Simulator(object)¶ The Simulator is the entry point to simulation, this is the very first object that needs to be initialized before any other method can be called. All the methods of the Simulator are static, but the constructor needs to be called first to initialize everything.
The Simulator manages all of the operating system interface, it manages window creation, end user input events, GPU access, threading, inter-process messaging and so forth. All ‘physical’ modeling concepts to in the
Universeobject.-
__init__(self, conf=None, **kwargs)¶ The
Simulator(...)construtor call is really just a synonum for the maininit()function. We keep this call here as earlier versions used this as the init.
-
static
run()¶ alias to top level
mechanica.run()function.
-
static
irun()¶ alias to top level
mechanica.irun()function.
-
static
close()¶ alias to top level
mechanica.close()function.
-
static
show()¶ alias to top level
mechanica.show()function.
-
-
class
Simulator.Config¶ An object that has all the arguments to the simulator,
Event / Message Processing¶
void glfwPollEvents ( void ) This function processes only those events that are already in the event queue and then returns immediately. Processing events will cause the window and input callbacks associated with those events to be called.
On some platforms, a window move, resize or menu operation will cause event processing to block. This is due to how event processing is designed on those platforms. You can use the window refresh callback to redraw the contents of your window when necessary during such operations.
Universe¶
-
class
mechanica.universe¶ The universe is a top level singleton object, and is automatically initialized when the simulator loads. The universe is a representation of the physical universe that we are simulating, and is the repository for all physical object representations.
All properties and methods on the universe are static, and you never actually instantiate a universe.
Universe has a variety of properties such as boundary conditions, and stores all the physical objects such as particles, bonds, potentials, etc..
-
static
bind(thing, a, b)¶
-
kinetic_energy¶ A read-only attribute that returns the total kinetic energy of the system
Type: double
-
dt¶ Get the main simulation time step.
Type: double
-
time¶ Get the current simulation time
-
temperature¶ Get / set the universe temperature.
The universe can be run with, or without a thermostat. With a thermostat, getting / setting the temperature changes the temperature that the thermostat will try to keep the universe at. When the universe is run without a thermostat, reading the temperature returns the computed universe temp, but attempting to set the temperature yields an error.
-
boltzmann_constant¶ Get / set the Boltzmann constant, used to convert average kinetic energy to temperature
-
static
particles()¶ List of all the particles in the universe. A particle can be removed from the universe using the standard python
delsyntaxdel universe.particles()[23]
Type: list
-
dim¶ Get / set the size of the universe, this is a length 3 list of floats. Currently we can only read the size, but want to enable changing universe size.
Type: Vector3
-
static
start()¶ Starts the universe time evolution, and advanced the universe forward by timesteps in
dt. All methods to build and manipulate universe objects are valid whether the universe time evolution is running or stopped.
-
static
stop()¶ Stops the universe time evolution. This essentially freezes the universe, everything remains the same, except time no longer moves forward.
-
static
step(until=None, dt=None)¶ Performs a single time step
dtof the universe if no arguments are given. Optionally runs untiluntil, and can use a different timestep ofdt.Parameters: - until – runs the timestep for this length of time, optional.
- dt – overrides the existing time step, and uses this value for time stepping, optional.
-
static
grid(shape)¶ Gets a three-dimesional array of particle lists, of all the particles in the system. Each
Parameters: shape – (length 3 array), list of how many grids we want in the x, y, z directions. Minimum is [1, 1, 1], which will return a array with a single list of all particles.
-
static
virial([origin][, radius][, types])¶ Computes the Virial Tensor for the either the entire simulation domain, or a specific local virial tensor at a location and radius. Optionally can accept a list of particle types to restrict the virial calculation for specify types.
Parameters: - origin – An optional length-3 array for the origin. Defaults to the center of the simulation domain if not given.
- radius – An optional number specifying the size of the region to compute the virial tensor for. Defaults to the entire simulation domain.
- types – An optional list of
Particletypes to include in the calculation. Defaults to every particle type.
-
static
Boundary Conditions¶
-
class
mechanica.BoundaryCondition¶ -
name¶ Gets the name of the boundary condition, i.e. “top”, “bottom”, etc..
-
kind¶ The ‘kind’ of the boundary condition, i.e. “PERIODIC”, “FREESLIP”, or “VELOCITY”
-
velocity¶ gets / sets the velocity of the boundary, a length-3 vector. Note a “no-slip” boundary is just a “velocity” boundary with a zero velocity vector.
-
normal¶ get the normal vector of the boundary, points into the simulation domain.
-
restore¶ gets / sets the restoring percent. When objects hit this boundary, they get reflected back at restore percent, so if restore is 0.5, and object hitting the boundary at 3 length / time recoils with a velocity of 1.5 lengths / time.
-
-
class
mechanica.Boundaryconditions¶ The BoundaryConditions class really just serves as a contianer for the six instances of the BoundaryCondition object:
-
left¶
-
right¶
-
front¶
-
back¶
-
bottom¶
-
top¶
-
Boundary Condition Constants¶
-
mechanica.BOUNDARY_NONE¶ no boundary conditions
-
mechanica.PERIODIC_X¶ periodic in the x direction
-
mechanica.PERIODIC_Y¶ periodic in the y direction
-
mechanica.PERIODIC_Z¶ periodic in the z direction
-
mechanica.PERIODIC_FULL¶ periodic in all directions
-
mechanica.PERIODIC_GHOST_X¶
-
mechanica.PERIODIC_GHOST_Y¶
-
mechanica.PERIODIC_GHOST_Z¶
-
mechanica.PERIODIC_GHOST_FULL¶
-
mechanica.FREESLIP_X¶ free slip in the x direction
-
mechanica.FREESLIP_Y¶ free slip in the y direction
-
mechanica.FREESLIP_Z¶ free slip in the z direction
-
mechanica.FREESLIP_FULL¶ free slip in all directions
Constants¶
Particles and Clusters¶
-
class
mechanica.Particle(object)¶ The particle is the most basic physical object we provide. All other physical objects either extend the base Particle, or are collections of particles such as the
Cluster-
static
items()¶ Returns a list of all of the instances of this type.
-
position¶ vector3– gets / sets the global position of the particle
-
velocity¶ vector3– returns the velocity of the particle, read-only
-
force¶ vector3– returns the net force that acts on this particle
-
charge¶ number– gets the total charge of the particle.
-
mass¶
-
radius¶
-
name¶
-
name2¶
-
dynamics¶ number– one of the Integrator Constants that specifies the time evolution of this particle. Particles can be either intertial particles that obey Newtonian dynamics, \(F=ma\) or overdamped dynamics, \(F \propto mv\).
-
age¶
-
style¶ Style– gets / sets the style of an object. When we create a new particle instance, it’s style points to the style attribute of the particle’s type, so that if we change something in the particle type, this changes every instance of that type. For more details, see the Style section.
-
frozen¶ Get / sets the frozen attribute. Frozen particles are fixed in place, and will not move if any force acts on them.
-
id¶
-
type_id¶
-
flags¶
-
become(type)¶ Dynamically changes the type of an object. We can change the type of a
Particlederived object to anyther pre-existingParticlederived type. What this means is that if we have an object of say type A, we can change it to another type, say B, and and all of the forces and processes that acted on objects of type A stip and the forces and processes defined for type B now take over. See section Changing Type for more details.Parameters: type – (Type)
-
split()¶ Splits a single particle into two, for more details, see section Splitting and Cleavage. The particle version of split is fairly simple, however the
Cluster.split()offers many more options.
-
destroy()¶ Destroys the particle, and removes it form inventory. The present object is handle that now references an empty particle. Calling any method after destroy will result in an error.
-
spherical([origin])¶ Calculates the particle’s coordinates in spherical coordinates (\([\rho, \theta, \phi]\)), where \(\rho\) is the distance from the origin, \(\theta\) is the azimuthal polar angle ranging from \([0,2 \pi]\), and \(phi\) is the declination from vertical, ranging from \([0,\pi]\)
Parameters: origin ([x,y,z]) – a vector of the origin to use for spherical coordinate calculations, optional, if not given, uses the center of the simulation domain as the origin.
-
virial([distance])¶ Computes the virial tensor, see Pressure and Virial Tensors.
Parameters: distance – (number (,optional)) distance from the center of this particle to include the other particles to use for the virial calculation. Return type: 3x3 matrix
-
neighbors([distance][, types])¶ Gets a list of all the other particles that are near the current one. By default, we list all the nearest particles that interact with the current one via forces.
Parameters: - distance – (number (,optional)) - An optional search distance, if specified will get all objects within the given distance. Defaults to the global simulation cutoff distance.
- types – (tuple, (,optional)) – If specified, can provide a tuple of types to include in the neighbor search. If types are provides, this method will return all non-cluster particles within a certain distance. Defaults to all types.
For example, to search for all objects of type A and B a distance of 1 unit away from a particle p, we would:
>>> nbrs = p.neighbors(distance=1, types=(A, B)) >>> print(len(nbrs))
-
static
-
class
mechanica.Cluster(Particle)¶ A Cluster is a collection of particles.
-
split([axis][, random][, normal][, point])¶ Splits the cluster into two clusters, where the first one is the original cluster and the new one is a new ‘daughter’ cluster.
split is discussed in detail in Splitting and Cleavage
Parameters: - axis – (length 3 vector (,optional)) - orientation axis for a split. If the ‘axis’ argument is given, the ‘split’ method chooses a random plane co-linear with this vector and uses this as the cleavage plane.
- random – (Boolean (,optional)) - ‘split’ chooses a random cleavage plane coincident with the center of mass of the cluster.
- normal – (length 3 vector (,optional)) - a normal direction for the cleavage plane.
- point – (length 3 vector (,optional)) - if given, uses this point to determine the point-normal form for the clevage plane.
-
virial()¶ Computes the Virial Tensor for the particles in this cluster.
-
radius_of_gyration()¶ Computes the Radius of Gyration for the particles in this cluster.
-
center_of_mass()¶ Computes the Center of Mass for the particles in this cluster.
-
center_of_geometry()¶ Computes the Center of Geometry for the particles in this cluster.
-
moment_of_inertia()¶ Computes the Moment of Inertia for the particles in this cluster.
-
centroid()¶ Convenience synonym for
center_of_geometry
-
inertia()¶ Convenience synonym for
moment_of_inertia
-
ParticleList¶
-
class
mechanica.ParticleList¶ Most Mechanica functions that return a list of particle return a ParticleList. This is a special list type that adds a lot of convenience methods for dealing with spatial information.
-
virial()¶ Returns the virial, or pressure tensor for all objects in the list.
-
radius_of_gyration()¶ Computes the radius of gyration
-
center_of_mass()¶ Computes the center of mass.
-
center_of_geometry()¶ Computes the center of geometry
-
centroid()¶ Computes the centroid
-
moment_of_inertia()¶ Computes the moment of inertia tensor
-
inertia()¶ Computes the moment of inertia tensor, a synonum for
moment_of_intertia()
-
copy()¶ Creates a deep copy of this list.
-
positions()¶ Returns a numpy 3xN array of all the particle positions.
-
spherical_positions([origin])¶ Computes the positions in spherical coordinates relative to an origin, returns a 3xN array.
-
velocities()¶ Returns a 3xN array of particle velocities
-
forces()¶ Returns a 3xN array of particle forces.
-
Potentials¶
A Potential object is a compiled interpolation of a given function. The Universe applies potentials to particles to calculate the net force on them.
For performance reasons, we found that implementing potentials as interpolations can be much faster than evaluating the function directly.
A potential can be treated just like any python callable object to evaluate it:
>>> pot = m.Potential.lennard_jones_12_6(0.1, 5, 9.5075e-06 , 6.1545e-03, 1.0e-3 )
>>> x = np.linspace(0.1,1,100)
>>> y=[pot(j) for j in x]
>>> plt.plot(x,y, 'r')
Potential is a callable object, we can invoke it like any Python function.
-
class
mechanica.Potential¶ -
static
power([k], [r0], [min], max=[10], [tol=0.001])¶ Parameters: - k – interaction strength, represents the potential energy peak value. Defaults to 1
- r0 – potential rest length, zero of the potential, defaults to 0.
- alpha – Exponent, defaults to 1.
- min – minimal value potential is computed for, defaults to r0 / 2.
- max – cutoff distance, defaults to 3 * r0.
- tol – Tolerance, defaults to 0.01.
some stuff, \(r_0 - r_0 / 2\).
The power potential the general form of many of the potential functions, such as
linear(), etc… power has the form:\[U(r,k) = k (r-r_0)^{\alpha}\]Some examples of power potentials are:
>>> import mechanica as m >>> p = m.Potential.power(k=1, r0=1, alpha=1.1) >>> p.plot(potential=True)
power potential function of \(U(r) = 1.0 (r - 1)^{1.1}\)
>>> import mechanica as m >>> p = m.Potential.power(k=1,alpha=0.5, r0=1) >>> p.plot(potential=True)
power potential function of \(U(r) = 1.0 (r - 1)^{0.5}\)
>>> import mechanica as m >>> p = m.Potential.power(k=1,alpha=2.5, r0=1) >>> p.plot(potential=True)
power potential function of \(U(r) = 1.0 (r - 1)^{2.5}\)
-
static
lennard_jones_12_6(min, max, a, b)¶ Creates a Potential representing a 12-6 Lennard-Jones potential
Parameters: - min – The smallest radius for which the potential will be constructed.
- max – The largest radius for which the potential will be constructed.
- A – The first parameter of the Lennard-Jones potential.
- B – The second parameter of the Lennard-Jones potential.
- tol – The tolerance to which the interpolation should match the exact potential., optional
The Lennard Jones potential has the form:
\[\left( \frac{A}{r^{12}} - \frac{B}{r^6} \right)\]
-
static
lennard_jones_12_6_coulomb(min, max, a, b, tol)¶ - Creates a Potential representing the sum of a
- 12-6 Lennard-Jones potential and a shifted Coulomb potential.
Parameters: - min – The smallest radius for which the potential will be constructed.
- max – The largest radius for which the potential will be constructed.
- A – The first parameter of the Lennard-Jones potential.
- B – The second parameter of the Lennard-Jones potential.
- q – The charge scaling of the potential.
- tol – The tolerance to which the interpolation should match the exact potential. (optional)
The 12-6 Lennard Jones - Coulomb potential has the form:
\[\left( \frac{A}{r^{12}} - \frac{B}{r^6} \right) + q \left(\frac{1}{r} - \frac{1}{max} \right)\]
-
static
ewald(min, max, q, kappa, tol)¶ Creates a potential representing the real-space part of an Ewald potential.
Parameters: - min – The smallest radius for which the potential will be constructed.
- max – The largest radius for which the potential will be constructed.
- q – The charge scaling of the potential.
- kappa – The screening distance of the Ewald potential.
- tol – The tolerance to which the interpolation should match the exact potential.
The Ewald potential has the form:
\[q \frac{\mbox{erfc}( \kappa r)}{r}\]-
static
well(k, n, r0[, min][, max][, tol])¶ Creates a continuous square well potential. Usefull for binding a particle to a region.
Parameters: - k (float) – potential prefactor constant, should be decreased for larger n.
- n (float) – exponent of the potential, larger n makes a sharper potential.
- r0 (float) – The extents of the potential, length units. Represents the maximum extents that a two objects connected with this potential should come appart.
- min (float) – [optional] The smallest radius for which the potential will be constructed, defaults to zero.
- max (float) – [optional] The largest radius for which the potential will be constructed, defaults to r0.
- tol (float) – [optional[ The tolerance to which the interpolation should match the exact potential, defaults to 0.01 * abs(min-max).
\[\frac{k}{\left(r_0 - r\right)^{n}}\]As with all potentials, we can create one, and plot it like so:
>>> p = m.Potential.well(0.01, 2, 1) >>> x=n.arange(0, 1, 0.0001) >>> y = [p(xx) for xx in x] >>> plt.plot(x, y) >>> plt.title(r"Continuous Square Well Potential $\frac{0.01}{(1 - r)^{2}}$ \n", ... fontsize=16, color='black')
A continuous square well potential.
-
static
harmonic_angle(k, theta0, [min, ]max[, tol])¶ Creates a harmonic angle potential
Parameters: - k – The energy of the angle.
- theta0 – The minimum energy angle.
- min – The smallest angle for which the potential will be constructed.
- max – The largest angle for which the potential will be constructed.
- tol – The tolerance to which the interpolation should match the exact potential.
returns A newly-allocated potential representing the potential
\[k(\theta-\theta_0)^2\]Note, for computational effeciency, this actually generates a function of r, where r is the cosine of the angle (calculated from the dot product of the two vectors. So, this actually evaluates internally,
\[k(\arccos(r)-\theta_0)^2\]
-
static
harmonic(k, r0[, min][, max][, tol])¶ Creates a harmonic bond potential
Parameters: - k – The energy of the bond.
- r0 – The bond rest length
- min – [optional] The smallest radius for which the potential will be constructed. Defaults to \(r_0 - r_0 / 2\).
- max – [optional] The largest radius for which the potential will be constructed. Defaults to \(r_0 + r_0 /2\).
- tol – [optional] The tolerance to which the interpolation should match the exact potential. Defaults to \(0.01 \abs(max-min)\)
return A newly-allocated potential
\[U(r) = k (r-r_0)^2\]We can plot the harmonic function to get an idea of it’s behavior:
>>> import mechanica as m >>> p = p = m.Potential.harmonic(k=10, r0=1, min=0.001, max=10) >>> p.plot(s=2, potential=True)
-
static
soft_sphere(kappa, epsilon, r0, eta, min, max, tol)¶ Creates a soft sphere interaction potential. The soft sphere is a generalized Lennard Jones type potentail, but you can varry the exponents to create a softer interaction.
Parameters: - kappa –
- epsilon –
- r0 –
- eta –
- min –
- max –
- tol –
-
static
glj(e[, m][, n][, r0][, min][, max][, tol][, shifted])¶ Parameters: - e – effective energy of the potential
- m – order of potential, defaults to 3
- n – order of potential, defaults to 2*m
- r0 – mimumum of the potential, defaults to 1
- min – minimum distance, defaults to 0.05 * r0
- max – max distance, defaults to 5 * r0
- tol – tolerance, defaults to 0.01
- shifted – is this shifted potential, defaults to true
Generalized Lennard-Jones potential.
\[V^{GLJ}_{m,n}(r) = \frac{\epsilon}{n-m} \left[ m \left( \frac{r_0}{r} \right)^n - n \left( \frac{r_0}{r} \right) ^ m \right]\]where \(r_e\) is the effective radius, which is automatically computed as the sum of the interacting particle radii.
The Generalized Lennard-Jones potential for different exponents \((m, n)\) with fixed \(n = 2m\). As the exponents grow smaller, the potential flattens out and becomes softer, but as the exponents grow larger the potential becomes narrower and sharper, and approaches the hard sphere potential.
-
static
overlapping_sphere(mu=1, [kc=1], [kh=0], [r0=0], [min=0.001], max=[10], [tol=0.001])¶ Parameters: - mu – interaction strength, represents the potential energy peak value.
- kc – decay strength of long range attraction. Larger values make a shorter ranged function.
- kh – Optionally add a harmonic long-range attraction, same as
glj()function. - r0 – Optional harmonic rest length, only used if kh is non-zero.
- min – Minimum value potential is computed for.
- max – Potential cutoff values.
- tol – Tolerance, defaults to 0.001.
The overlapping_sphere function implements the Overlapping Sphere, from [OFPF+17]. This is a soft sphere, from our testing, it appears too soft, probably better suited for 2D models. This potential appears to allow particles to collapse too closely, probably needs more paramater fiddling.
Note
From the equation below, we can see that there is a \(\log\) term as the short range repulsion term. The logarithm is the radial Green’s function for cylindrical (2D) geometry, however the Green’s function for 3D is the \(1/r\) function. This is possibly why Osborne has success in 2D but it’s unclear if this was used in 3D geometries.
\[\begin{split}\mathbf{F}_{ij}= \begin{cases} \mu_{ij} s_{ij}(t) \hat{\mathbf{r}}_{ij} \log \left( 1 + \frac{||\mathbf{r}_{ij}|| - s_{ij}(t)}{s_{ij}(t)} \right) ,& \text{if } ||\mathbf{r}_{ij}|| < s_{ij}(t) \\ \mu_{ij}\left(||\mathbf{r}_{ij}|| - s_{ij}(t)\right) \hat{\mathbf{r}}_{ij} \exp \left( -k_c \frac{||\mathbf{r}_{ij}|| - s_{ij}(t)}{s_{ij}(t)} \right) ,& \text{if } s_{ij}(t) \leq ||\mathbf{r}_{ij}|| \leq r_{max} \\ 0, & \text{otherwise} \\ \end{cases}\end{split}\]Osborne refers to \(\mu_{ij}\) as a “spring constant”, this controls the size of the force, and is the potential energy peak value. \(\hat{\mathbf{r}}_{ij}\) is the unit vector from particle \(i\) center to particle \(j\) center, \(k_C\) is a parameter that defines decay of the attractive force. Larger values of \(k_C\) result in a shaper peaked attraction, and thus a shorter ranged force. \(s_{ij}(t)\) is the is the sum of the radii of the two particles.
We can plot the overlapping sphere function to get an idea of it’s behavior:
>>> import mechanica as m >>> p = m.Potential.overlapping_sphere(mu=10, max=20) >>> p.plot(s=2, force=True, potential=True, ymin=-10, ymax=8)
We can plot the overlapping sphere function like any other function.
-
static
linear(k, [min], max=[10], [tol=0.001])¶ Parameters: - k – interaction strength, represents the potential energy peak value.
- min – Minimum value potential is computed for.
- max – Potential cutoff values.
- tol – Tolerance, defaults to 0.001.
The linear potential is the simplest one, it’s simply a potential of the form
\[U(r) = k r\]>>> import mechanica as m >>> p = p = m.Potential.linear(k=5, max=10) >>> p.plot(potential=True)
Linear potential function
-
static
Dissapative Particle Dynamics¶
-
Potential.dpd(alpha=1, gamma=1, sigma=1, cutoff=1)¶ Parameters: - alpha – interaction strength of the conservative force
- gamma – interaction strength of dissapative force
- sigma – strength of random force.
The dissapative particle dynamics force is the sum of the conservative, dissapative and random forces:
\[\mathbf{F}_{ij} = \mathbf{F}^C_{ij} + \mathbf{F}^D_{ij} + \mathbf{F}^R_{ij}\]The conservative force is:
\[\mathbf{F}^C_{ij} = \alpha \left(1 - \frac{r_{ij}}{r_c}\right) \mathbf{e}_{ij}\]The dissapative force is:
\[\mathbf{F}^D_{ij} = -\gamma \left(1 - \frac{r_{ij}}{r_c}\right)^{2}(\mathbf{e}_{ij} \cdot \mathbf{v}_{ij}) \mathbf{e}_{ij}\]And the random force is:
\[\mathbf{F}^R_{ij} = \sigma \left(1 - \frac{r_{ij}}{r_c}\right) \xi_{ij}\Delta t^{-1/2}\mathbf{e}_{ij}\]
-
class
mechanica.DPDPotential¶ A subtype of the Potential class that adds features specific to Dissapative Particle Dynamics.
Forces¶
Forces are one of the fundamental processes in Mechanica that cause objects to move. We provide a suite of pre-built forces, and users can create their own.
The forces module collects all of the built-in forces that Mechanica provides. This module contains a variety of functions that all generate a force object.
-
mechanica.forces.berendsen_thermostat(tau)¶ Creates a Berendsen thermostat
Param: tau: time constant that determines how rapidly the thermostat effects the system. The thermostat picks up the target temperature \(T_0\) from the object that it gets bound to. For example, if we bind a temperature to a particle type, then it uses the
The Berendsen thermostat effectively re-scales the velocities of an object in order to make the temperature of that family of objects match a specified temperature.
The Berendsen thermostat force \(\mathbf{F}_{temp}\) has a function form of:
\[\mathbf{F}_{temp} = \frac{\mathbf{p}_i}{\tau_T} \left(\frac{T_0}{T} - 1 \right),\]where \(T\) is the measured temperature of a family of particles, \(T_0\) is the control temperature, and \(\tau_T\) is the coupling constant. The coupling constant is a measure of the time scale on which the thermostat operates, and has units of time. Smaller values of \(\tau_T\) result in a faster acting thermostat, and larger values result in a slower acting thermostat.
..class:: ConstantForce(value, period=None)
A Constant Force that acts on objects, is a way to specify things like gravity, external pressure, etc…
param value: value can be either a callable that returns a length-3 vector, or a length-3 vector. If value is a vector, than we simply use that vector value as the constant force. If value is callable, then the runtime calls that callable every period to get a new value. If it’s a callable, then it must return somethign convertible to a length-3 array. param period: If value is a callable, then period defines how frequenly the runtime should call that to obtain a new force value.
Reactions and Species¶
Events are the principle way users can attach connect thier own functions, and built in methods event triggers.
-
mechanica.on_time(method, period[, start][, end][, distribution])¶ Binds a function or method to a elapsed time interval.
on_time must be called with all keyword arguments, since there are a large number of options this method accepts.
Parameters: - method – A Python function or method that will get called. Signature
of the method should be
method(time), in that it gets a single current time argument. - period (float) – Time interval that the event should get fired.
- start (float) – [optional] Start time for the event.
- end (float) – [optional] end time for the event.
- distribution (str) – [optional] String that identifies the statistical distribution for event times. Only supported distibution currently is “exponential”. If there is no distibution argument, the event is called deterministically.
- method – A Python function or method that will get called. Signature
of the method should be
Fluxes¶
-
mechanica.flux(a, b, species_name, k[, decay=0])¶ The basic passive (Fickian) flux type. This flux implements a passive transport between a species located on a pair of nearby objects of type a and b. See Spatial Transport and Flux for a detailed discussion.
Parameters: - a (type) – object type a
- b (type) – object type b
- species_name (string) – string, textual name of the species to perform a flux with
- k (number) – flux rate constant
- decay (number) – decay rate, optional, defaults to 0.
A Fick flux of the species \(S\) attached to object types \(A\) and \(B\) implements the reaction:
\[\begin{split}\begin{eqnarray} a.S & \leftrightarrow a.S \; &; \; k \left(1 - \frac{r}{r_{cutoff}} \right)\left(a.S - b.S\right) \\ a.S & \rightarrow 0 \; &; \; \frac{d}{2} a.S \\ b.S & \rightarrow 0 \; &; \; \frac{d}{2} b.S, \end{eqnarray}\end{split}\]\(B\) respectivly. \(S\) is a chemical species located at each object instances. \(k\) is the flux constant, \(r\) is the distance between the two objects, \(r_{cutoff}\) is the global cutoff distance, and \(d\) is the optional decay term.
-
mechanica.produce_flux(a, b, species_name, k, target[, decay=0])¶ An active flux that represents active transport, and is can be used to model such processes like membrane ion pumps. See Spatial Transport and Flux for a detailed discussion. Unlike the consume_flux, the produce flux uses the concentation of only the source to determine the rate.
Parameters: - a (type) – object type a
- b (type) – object type b
- species_name (string) – string, textual name of the species to perform a flux with
- k (number) – flux rate constant
- target (number) – target concentation of the \(b.S\) species.
- decay (number) – decay rate, optional, defaults to 0.
The produce flux implements the reaction:
\[\begin{split}\begin{eqnarray} a.S & \rightarrow b.S \; &; \; k \left(r - \frac{r}{r_{cutoff}} \right)\left(a.S - a.S_{target} \right) \\ a.S & \rightarrow 0 \; &; \; \frac{d}{2} a.S \\ b.S & \rightarrow 0 \; &; \; \frac{d}{2} b.S \end{eqnarray}\end{split}\]
-
mechanica.consume_flux(a, b, species_name, k, target_concentation[, decay=0])¶ An active flux that represents active transport, and is can be used to model such processes like membrane ion pumps. See Spatial Transport and Flux for a detailed discussion.
Parameters: - a (type) – object type a
- b (type) – object type b
- species_name (string) – string, textual name of the species to perform a flux with
- k (number) – flux rate constant
- target (number) – target concentation of the \(b.S\) species.
- decay (number) – decay rate, optional, defaults to 0.
The consume flux implements the reaction:
\[\begin{split}\begin{eqnarray} a.S & \rightarrow b.S \; &; \; k \left(1 - \frac{r}{r_{cutoff}}\right)\left(b.S - b.S_{target} \right)\left(a.S\right) \\ a.S & \rightarrow 0 \; &; \; \frac{d}{2} a.S \\ b.S & \rightarrow 0 \; &; \; \frac{d}{2} b.S \end{eqnarray}\end{split}\]
Bonded Potentials¶
Events¶
Events are the principle way users can attach connect thier own functions, and built in methods event triggers.
-
mechanica.on_time(method, period[, start][, end][, distribution]) Binds a function or method to a elapsed time interval.
on_time must be called with all keyword arguments, since there are a large number of options this method accepts.
Parameters: - method – A Python function or method that will get called. Signature
of the method should be
method(time), in that it gets a single current time argument. - period (float) – Time interval that the event should get fired.
- start (float) – [optional] Start time for the event.
- end (float) – [optional] end time for the event.
- distribution (str) – [optional] String that identifies the statistical distribution for event times. Only supported distibution currently is “exponential”. If there is no distibution argument, the event is called deterministically.
- method – A Python function or method that will get called. Signature
of the method should be
-
mechanica.on_keypress(callback)¶ Register for a key press event. The callback function will get called with a key press object every time the user presses a key. The Keypress object has a
key_nameproperty to get the key value:def my_keypress_handler(e): print("key pressed: ", e.key_name) m.on_keypress(my_keypress_handler)
Style¶
-
class
mechanica.Style(object)¶ The Style object contains all of the rendering visual style attribute of an object. Most simulation objects in Mechanica have a style attribute.
-
color¶ Color4– gets / sets the color of the style. This is actually aColor4object, but we can set it’s color value either as a Color object, or via one of the color strings from Colors, i.e.>>> c = m.Style() >>> c.color = "DeepSkyBlue"
The color string is not case sensitive.
-
visible¶ boolean– gets / sets the visible property. This toggles the visibility of any rendered object.
-
Colors¶
You can construct a color object, (either a Color3 or Color4) using one of the following web color names:
- “AliceBlue”,
- “AntiqueWhite”,
- “Aqua”,
- “Aquamarine”,
- “Azure”,
- “Beige”,
- “Bisque”,
- “Black”,
- “BlanchedAlmond”,
- “Blue”,
- “BlueViolet”,
- “Brown”,
- “BurlyWood”,
- “CadetBlue”,
- “Chartreuse”,
- “Chocolate”,
- “Coral”,
- “CornflowerBlue”,
- “Cornsilk”,
- “Crimson”,
- “Cyan”,
- “DarkBlue”,
- “DarkCyan”,
- “DarkGoldenRod”,
- “DarkGray”,
- “DarkGreen”,
- “DarkKhaki”,
- “DarkMagenta”,
- “DarkOliveGreen”,
- “Darkorange”,
- “DarkOrchid”,
- “DarkRed”,
- “DarkSalmon”,
- “DarkSeaGreen”,
- “DarkSlateBlue”,
- “DarkSlateGray”,
- “DarkTurquoise”,
- “DarkViolet”,
- “DeepPink”,
- “DeepSkyBlue”,
- “DimGray”,
- “DodgerBlue”,
- “FireBrick”,
- “FloralWhite”,
- “ForestGreen”,
- “Fuchsia”,
- “Gainsboro”,
- “GhostWhite”,
- “Gold”,
- “GoldenRod”,
- “Gray”,
- “Green”,
- “GreenYellow”,
- “HoneyDew”,
- “HotPink”,
- “IndianRed”,
- “Indigo”,
- “Ivory”,
- “Khaki”,
- “Lavender”,
- “LavenderBlush”,
- “LawnGreen”,
- “LemonChiffon”,
- “LightBlue”,
- “LightCoral”,
- “LightCyan”,
- “LightGoldenRodYellow”,
- “LightGrey”,
- “LightGreen”,
- “LightPink”,
- “LightSalmon”,
- “LightSeaGreen”,
- “LightSkyBlue”,
- “LightSlateGray”,
- “LightSteelBlue”,
- “LightYellow”,
- “Lime”,
- “LimeGreen”,
- “Linen”,
- “Magenta”,
- “Maroon”,
- “MediumAquaMarine”,
- “MediumBlue”,
- “MediumOrchid”,
- “MediumPurple”,
- “MediumSeaGreen”,
- “MediumSlateBlue”,
- “MediumSpringGreen”,
- “MediumTurquoise”,
- “MediumVioletRed”,
- “MidnightBlue”,
- “MintCream”,
- “MistyRose”,
- “Moccasin”,
- “NavajoWhite”,
- “Navy”,
- “OldLace”,
- “Olive”,
- “OliveDrab”,
- “Orange”,
- “OrangeRed”,
- “Orchid”,
- “PaleGoldenRod”,
- “PaleGreen”,
- “PaleTurquoise”,
- “PaleVioletRed”,
- “PapayaWhip”,
- “PeachPuff”,
- “Peru”,
- “Pink”,
- “Plum”,
- “PowderBlue”,
- “Purple”,
- “Red”,
- “RosyBrown”,
- “RoyalBlue”,
- “SaddleBrown”,
- “Salmon”,
- “SandyBrown”,
- “SeaGreen”,
- “SeaShell”,
- “Sienna”,
- “Silver”,
- “SkyBlue”,
- “SlateBlue”,
- “SlateGray”,
- “Snow”,
- “SpringGreen”,
- “SteelBlue”,
- “Tan”,
- “Teal”,
- “Thistle”,
- “Tomato”,
- “Turquoise”,
- “Violet”,
- “Wheat”,
- “White”,
- “WhiteSmoke”,
- “Yellow”,
- “YellowGreen”,
- For example, to make some colors::
>>> m.Color3("red") Vector(1, 0, 0)
>>> m.Color3("MediumSeaGreen") Vector(0.0451862, 0.450786, 0.165132)
>>> m.Color3("CornflowerBlue") Vector(0.127438, 0.300544, 0.846873)
>>> m.Color3("this is total garbage") /usr/local/bin/ipython3:1: Warning: Warning, "this is total garbage" is not a valid color name. #!/usr/local/opt/python/bin/python3.7 Vector(0, 0, 0)
As it’s easy to make a mistake with color names, we simply issue a warnign, instead of an error if the color name can’t be found.
-
class
mechanica.Color3¶ The class can store either a floating-point or integer representation of a linear RGB color. Colors in sRGB color space should not beused directly in calculations — they should be converted to linear RGB using fromSrgb(), calculation done on the linear representation and then converted back to sRGB using toSrgb().
You can construct a Color object
-
static
from_srgb(srgb(int))¶ constructs a color from a packed integer, i.e.:
>>> c = m.Color3.from_srgb(0xffffff) >>> print(c) Vector(1, 1, 1)
-
static
Logging¶
Mechanica has extensive logging system. Many internal methods will log extensive details (in full color) to either the clog (typically stderr) or a user specified file path. Internally, the Mechanica logging system is currently implemented by the Poco (http://pocoproject) logging system.
Future versions will include a Logging Listener, so that the Mechanica log will log messages to a user specified function.
The logging system is highly configurable, users have complete control over the color and format of the logging messages.
All methods of the Logger are static, they are available immediately upon loading the Mechanica package.
If one wants to display logging at the lowest level (LOG_TRACE), where every logging message is displayed, one would run:
m.Logging.set_level(m.Logging.LOG_TRACE)
Logging the following messages:
m.Logger.log(m.Logger.LOG_FATAL, "A fatal message")
m.Logger.log(m.Logger.LOG_CRITICAL, "A critical message")
m.Logger.log(m.Logger.LOG_ERROR, "An error message")
m.Logger.log(m.Logger.LOG_WARNING, "A warning message")
m.Logger.log(m.Logger.LOG_NOTICE, "A notice message")
m.Logger.log(m.Logger.LOG_INFORMATION, "An informational message")
m.Logger.log(m.Logger.LOG_DEBUG, "A debugging message.")
m.Logger.log(m.Logger.LOG_TRACE, "A tracing message. This is the lowest priority.")
will produce the following output:
If one wants different colors on the log, these can be set via:
rr.Logger.set_property("traceColor", "red")
rr.Logger.set_property("debugColor", "green")
The available color property names and values are listed below at the Logger.setProperty method.
Logging Levels¶
The Logger has the following logging levels:
-
Logger.LOG_CURRENT¶ Use the current level – don’t change the level from what it is.
-
Logger.LOG_FATAL¶ A fatal error. The application will most likely terminate. This is the highest priority.
-
Logger.LOG_CRITICAL¶ A critical error. The application might not be able to continue running successfully.
-
Logger.LOG_ERROR¶ An error. An operation did not complete successfully, but the application as a whole is not affected.
-
Logger.LOG_WARNING¶ A warning. An operation completed with an unexpected result.
-
Logger.LOG_NOTICE¶ A notice, which is an information with just a higher priority.
-
Logger.LOG_INFORMATION¶ An informational message, usually denoting the successful completion of an operation.
-
Logger.LOG_DEBUG¶ A debugging message.
-
Logger.LOG_TRACE¶ A tracing message. This is the lowest priority.
Logging Methods¶
-
static
Logger.set_level([level])¶ sets the logging level to one a value from Logger::Level
Parameters: level (int) – the level to set, defaults to LOG_CURRENT if none is specified.
-
static
Logger.get_level()¶ get the current logging level.
-
static
Logger.disable_logging()¶ Suppresses all logging output
-
static
Logger.disable_console_logging()¶ stops logging to the console, but file logging may continue.
-
static
Logger.enable_console_logging(level)¶ turns on console logging (stderr) at the given level.
Parameters: level – A logging level, one of the above listed LOG_* levels.
-
static
Logger.enable_file_logging(fileName[, level])¶ turns on file logging to the given file as the given level.
Parameters: - fileName (str) – the path of a file to log to.
- level – (optional) the logging level, defaults to LOG_CURRENT.
-
static
Logger.disable_file_logging()¶ turns off file logging, but has no effect on console logging.
-
static
Logger.get_current_level_as_string()¶ get the textural form of the current logging level.
-
static
Logger.get_filename()¶ get the name of the currently used log file.
-
static
Logger.set_formatting_pattern(format)¶ Internally, Mechanica uses the Poco logging framework, so we can custom format logging output based on a formatting pattern string.
The format pattern is used as a template to format the message and is copied character by character except for the following special characters, which are replaced by the corresponding value.
An example pattern of “%Y-%m-%d %H:%M:%S %p: %t” set via:
m.Logger.set_formatting_pattern("%Y-%m-%d %H:%M:%S %p: %t")
would produce the following output:
Mechanica supports the following format specifiers. These were copied from the Poco documentation:
- %s - message source
- %t - message text
- %l - message priority level (1 .. 7)
- %p - message priority (Fatal, Critical, Error, Warning, Notice, Information, Debug, Trace)
- %q - abbreviated message priority (F, C, E, W, N, I, D, T)
- %P - message process identifier
- %T - message thread name
- %I - message thread identifier (numeric)
- %N - node or host name
- %U - message source file path (empty string if not set)
- %u - message source line number (0 if not set)
- %w - message date/time abbreviated weekday (Mon, Tue, …)
- %W - message date/time full weekday (Monday, Tuesday, …)
- %b - message date/time abbreviated month (Jan, Feb, …)
- %B - message date/time full month (January, February, …)
- %d - message date/time zero-padded day of month (01 .. 31)
- %e - message date/time day of month (1 .. 31)
- %f - message date/time space-padded day of month ( 1 .. 31)
- %m - message date/time zero-padded month (01 .. 12)
- %n - message date/time month (1 .. 12)
- %o - message date/time space-padded month ( 1 .. 12)
- %y - message date/time year without century (70)
- %Y - message date/time year with century (1970)
- %H - message date/time hour (00 .. 23)
- %h - message date/time hour (00 .. 12)
- %a - message date/time am/pm
- %A - message date/time AM/PM
- %M - message date/time minute (00 .. 59)
- %S - message date/time second (00 .. 59)
- %i - message date/time millisecond (000 .. 999)
- %c - message date/time centisecond (0 .. 9)
- %F - message date/time fractional seconds/microseconds (000000 - 999999)
- %z - time zone differential in ISO 8601 format (Z or +NN.NN)
- %Z - time zone differential in RFC format (GMT or +NNNN)
- %E - epoch time (UTC, seconds since midnight, January 1, 1970)
- %[name] - the value of the message parameter with the given name
- %% - percent sign
Parameters: format (str) – the logging format string. Must be formatted using the above specifiers.
-
static
Logger.get_formatting_pattern()¶ get the currently set formatting pattern.
-
static
Logger.level_to_string(level)¶ gets the textual form of a logging level Enum for a given value.
Parameters: level (int) – One of the above listed logging levels.
-
static
Logger.string_to_level(s)¶ parses a string and returns a Logger::Level
Parameters: s (str) – the string to parse.
-
static
Logger.get_colored_output()¶ check if we have colored logging enabled.
-
static
Logger.set_colored_output(b)¶ enable / disable colored output
Parameters: b (boolean) – turn colored logging on or off
-
static
Logger.set_property(name, value)¶ Set the color of the output logging messages.
In the future, we may add additional properties here.
The following properties are supported:
- enableColors: Enable or disable colors.
- traceColor: Specify color for trace messages.
- debugColor: Specify color for debug messages.
- informationColor: Specify color for information messages.
- noticeColor: Specify color for notice messages.
- warningColor: Specify color for warning messages.
- errorColor: Specify color for error messages.
- criticalColor: Specify color for critical messages.
- fatalColor: Specify color for fatal messages.
The following color values are supported:
- default
- black
- red
- green
- brown
- blue
- magenta
- cyan
- gray
- darkgray
- lightRed
- lightGreen
- yellow
- lightBlue
- lightMagenta
- lightCyan
- white
Parameters: - name (str) – the name of the value to set.
- value (str) – the value to set.
-
static
Logger.log(level, msg)¶ logs a message to the log.
Parameters: - level (int) – the level to log at.
- msg (str) – the message to log.
Rendering and System Interaction¶
The system module contains various function to interact with rendering engine and host cpu.
-
mechanica.system.cpuinfo()¶ Returns a dictionary that contains a variety of information about the current procssor, such as processor vendor and supported instruction set features.
-
mechanica.system.gl_info()¶ Returns a dictionary that contains OpenGL capabilities and supported extensions.
-
mechanica.system.egl_info()¶ Gets a string of EGL info, only valid on Linux, we use EGL for Linux headless rendering.
-
mechanica.system.image_data()¶ Gets the contents of the rendering frame buffer back as a JPEG image stream, a byte array of the packed JPEG image.
-
mechanica.system.context_has_current()¶ Checks of the currently executing thread has a rendering context.
-
mechanica.system.context_make_current()¶ Makes the single Mechanica rendering context current on the current thread. Note, for multi-threaded rendering, the initialization thread needs to release the context via
context_release(), then aquire it on the rendering thread via this function.
-
mechanica.system.context_release()¶ Release the rendering context on the currently executing thread.
-
mechanica.system.camera_move_to([eye][, center][, up])¶ Moves the camera to the specified ‘eye’ location, where it looks at the ‘center’ with an up vector of ‘up’. All of the parameters are optional, if they are not give, we use the initial default values.
Parameters: - eye – ([x, y, z]) New camera location, where we re-position the eye of the viewer.
- center – ([x, y, z]) Location that the camera will look at, the what we want as the center of the view.
- up – ([x, y, z]) Unit up vector, defines the up direction of the camera.
-
mechanica.system.camera_reset()¶ Resets the camera to the initial position.
-
mechanica.system.camera_rotate_mouse(mouse_pos)¶ Rotates the camera according to the current mouse position. Need to call
camera_init_mouse()to initialize a mouse movement.Parameters: mouse_pos – ([x, y]) current mouse position on the view window or image.
-
mechanica.system.camera_translate_mouse(mouse_pos)¶ Translates the camera according to the current mouse position. You need to call
camera_init_mouse()to initialize a mouse movment, i.e. set the starting mouse position.Parameters: mouse_pos – ([x, y]) current mouse position on the view window or image.
-
mechanica.system.camera_init_mouse(mouse_pos)¶ Initialize a mouse movment operation, this tells the simulator that a mouse click was performed at the given coordinates, and subsequent mouse motion will refer to this starting position.
Parameters: mouse_pos – ([x, y]) current mouse position on the view window or image.
-
mechanica.system.camera_translate_by(delta)¶ Translates the camera in the plane perpendicular to the view orientation, moves the camera a given delta x, delta y distance in that plane.
Parameters: delta – ([delta_x, delta_y]) a vector that indicates how much to translate the camera.
-
mechanica.system.camera_zoom_by(delta)¶ Zooms the camera in and out by a specified amount.
Parameters: delta – number that indicates how much to increment zoom distance.
-
mechanica.system.camera_zoom_to(distance)¶ Zooms the camera to the given distance.
Parameters: distance – distance to the universe center for the camera.
-
mechanica.system.camera_rotate_to_axis(axis, distance)¶ Rotates the camera to one of the principal axis, at a given zoom distance.
Parameters: - axis – ([x, y, z]) unit vector that defines the axis to move to.
- distance – how far away the camera will be.
-
mechanica.system.camera_rotate_to_euler_angle(angles)¶ Rotate the camera to the given orientiation defined by three Euler angles.
Parameters: angles – ([alpha, beta, gamma]) Euler angles of rotation about the X, Y, and Z axis.
-
mechanica.system.camera_rotate_by_euler_angle(angles)¶ Incremetns the camera rotation by given orientiation defined by three Euler angles.
Parameters: angles – ([alpha, beta, gamma]) Euler angles of rotation about the X, Y, and Z axis.
-
mechanica.system.view_reshape(window_size)¶ Notify the simulator that the window or image size was changed.
Parameters: window_size – ([x, y]) new window size.
Setting Up a Development Enviorment¶
Git¶
Git submodule updates
To update submodule, need to do both the followign:
git pull –recurse-submodules
git submodule update –recursive –remote
Windows¶
Triplet: x64-windows-static
Cmake defaults to : Multi-threaded Debug DLL (/MDd)
but assimp built statically seems to be Multi-threaded Debug (/MTd)
vcpkg static linking
Python Debug Builds¶
It is very informative to run the Python interpreter in debug mode, especially when developing advanced language extensions.
Numpy¶
The default version of numpy, from PyPy does not work with debug builds of Python, at least it does not on MacOS with Python 3.7 and Numpy version 1.18. The standard source build of Numpy will not work with debug Python builds, confirmed with the above versions.
At least will debug builds, Numpy also will not work with the built in MacOS ‘Accelerate’ framework which provides a version of BLAS. Don’t know exaclty why, but it compiles and installs fine, but fails the sanity test:
def _sanity_check():
x = ones(2, dtype=float32)
if not abs(x.dot(x) - 2.0) < 1e-5:
raise AssertionError()
We don’t know why the dot product fails with the Accelerate framework lapack. The fix is to build Numpy with Open BLAS. This is simple, however poorly documented. Steps:
1: Grab the numpy source code from git, and go into that directory:
git clone https://github.com/your-user-name/numpy.git
cd numpy
git remote add upstream https://github.com/numpy/numpy.git
2: Install Open BLAS from Brew
localhost:$ brew reinstall openblas
It will say something like:
==> Downloading ...
...
==> Pouring openblas-0.3.9.high_sierra.bottle.tar.gz
==> Caveats
openblas is keg-only, which means it was not symlinked into /usr/local,
because macOS provides BLAS and LAPACK in the Accelerate framework.
For compilers to find openblas you may need to set:
export LDFLAGS="-L/usr/local/opt/openblas/lib"
export CPPFLAGS="-I/usr/local/opt/openblas/include"
For pkg-config to find openblas you may need to set:
export PKG_CONFIG_PATH="/usr/local/opt/openblas/lib/pkgconfig"
==> Summary
🍺 /usr/local/Cellar/openblas/0.3.9: 23 files, 120MB
Make note of the LDFLAGS and CPPFLAGS values, you will need to paste these
into the Numpy config file.
3: In the numpy directory, there is a file called example.site.cfg, rename
this to site.cfg, and open it up with your favorite editor. Look for the line
that says openblas, will look like this:
# [openblas]
# libraries = openblas
# library_dirs = /opt/OpenBLAS/lib
# include_dirs = /opt/OpenBLAS/include
# runtime_library_dirs = /opt/OpenBLAS/lib
Uncomment these lines, and replace them with the values of LDFLAGS and
CPPFLAGS from above:
[openblas]
libraries = openblas
library_dirs = /usr/local/opt/openblas/lib
include_dirs = /usr/local/opt/openblas/include
runtime_library_dirs = /usr/local/opt/openblas/lib
4: The really important part is to tell the Numpy buuild to actually use Open
BLAS. There does not seem to be any way in the site.cfg file to tell numpy
this, rather it apears you have to use enviornment variables (ugh!!!). So, set
the NPY_BLAS_ORDER enviorment variable to use Open BLAS:
export NPY_BLAS_ORDER=openblas
5: Then perform a configuration, build and install of Numpy with the following commands (make sure to use your correct debug build Python here):
python3 setup.py config
python3 setup.py build
python3 setup.py install
Pay particular attention to the output of the config step:
python3 setup.py config
It will look something like this:
gcc: /var/folders/p7/vcwd91x50j96v_yh75xtl4p40000gn/T/tmpr2lbwhqr/source.c
gcc /var/folders/p7/vcwd91x50j96v_yh75xtl4p40000gn/T/tmpr2lbwhqr/var/ ... 40000gn/T/tmpr2lbwhqr/a.out
FOUND:
libraries = ['openblas', 'openblas']
library_dirs = ['/usr/local/opt/openblas/lib']
language = c
define_macros = [('HAVE_CBLAS', None)]
runtime_library_dirs = ['/usr/local/opt/openblas/lib']
FOUND:
libraries = ['openblas', 'openblas']
library_dirs = ['/usr/local/opt/openblas/lib']
language = c
define_macros = [('HAVE_CBLAS', None)]
runtime_library_dirs = ['/usr/local/opt/openblas/lib']
lapack_opt_info:
lapack_mkl_info:
FOUND:
libraries = ['mkl_rt', 'pthread']
library_dirs = ['/usr/local/lib']
define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
include_dirs = ['/usr/local/include', '/usr/include', '/Users/andy/local/include']
FOUND:
libraries = ['mkl_rt', 'pthread']
library_dirs = ['/usr/local/lib']
define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
include_dirs = ['/usr/local/include', '/usr/include',
'/Users/andy/local/include']
The first sections must say that it found Open BLAS first.
Then finish with the build and install steps.
6: Test it BOTH from Python and IPython, very important. They seem to set up different library paths, and frequently if numpy works in one, it won’t work in the other, if Open BLAS is not setup correctly. In IPython, simply:
>>> import numpy
And test it from python via the command line as:
python3 -c "import numpy; print(numpy.get_include())"
Debugging Mechanica From Python¶
We provide the mx-pyrun app / projct. This is a trivial program that simply calls the main Python main routine, but the purpose of this probram is to serve as a target app, as an entry point in your IDE so that the mechanica python library can be loaded and stepped through.
Linux¶
Prerequisites:
sudo apt-get install libjpeg-dev
sudo apt-get install python-dev
sudo apt-get install python3-pip pip3 install ipython –user
pip3 install numpy –user
libgl1-mesa-dev
libassimp-dev libgl1-mesa-dev libxrandr-dev libxinerama-dev libxcursor-dev libgegl-dev libglfw3-dev
libxi-dev
Preparing a PyPi Release¶
Currently, we copy dist files to a dropbox folder. Need to implement cleaner soltuion.
From the command line:
>>> python3 setup.py sdist bdist_wheel
>>> cp dist/* ~/Dropbox/dist/
Then on upload computer:
>>> twine upload dist/*
Status¶
Mechancia is very early in the development cycle, as such, it’s EXTREMLY feature-incomplete, unstable, and will almost certainly crash.
However, as such, this is YOUR chance to try it out, and let us know what kind of features you’d like to see.
History¶
Version Alpha 1.0.18.1¶
- generalized passive, consumer and producer fluxes
- better OpenGL info reporting, gl_info(), egl_info()
- enable boundary conditions on chemical speices, bug fix parsing init conditions
- use species boundary value to enable source / sinks
- source / sinks in example
Version Alpha 1.0.17.1¶
- multi-threaded rendering fixes
Version Alpha 1.0.16.2¶
- Logging, standardized all logging output, python api for setting log level.
- fix kinetic energy reporting
- synchronize gl contexts between GLFW and Magnum for multi-thread rendering
Version Alpha 1.0.16.2¶
- initialize Mechanica either via m.init, m.Simulator, or m.simulator.init
Version Alpha 1.0.16.1¶
- finally, completly expunged pybind11! pybind11 is finally GONE!
- context managment methods for multi-threaded headless rendering.
- universe.reset() method, clears objects
- set window title to script name
- add ‘positions()’, ‘velocities()’ and ‘forces()’ methods to particle list.
- universe.particles() is now a method, and returns a proper list
Version Alpha 1.0.15.5¶
- bug fix with boundary condition constants
Version Alpha 1.0.15.5¶
- bug fix with force calculation when distance too short: pic random separation vector of with minimal distance. Seems to work…
- better diagnostic messages
- added normal to boundary vectors
Version Alpha 1.0.15.4¶
- generalized boundary conditions
- add potentials to boundary conditions
- velocity, free-slip, no-slip and periodic boundary conditions
- render updates, back face culling
- headless rendering, rendering without X11 using GLES on Linux
- generalized power potential
- much improved error handling, much more consistency
- particle list fixes
- Rigid Body Dynamics ! (only cuboids currently supported, but still rigid bodies)
- add potentials to rigid bodies
- python api rigid body updates
- rendering updates, more consistency, simplify
- rigid body particle interactions
- friction force
- more expunging pybind, soon, soon we will be rid of pybind.
- bond dissociation_energy (break strength)
- lattice initializer
- add bonds to lattice initliazer
- performance logging
- updates to dissapative particle dynamics forces
- enable adding DPD force to boundaries.
- generlized single body force (external force)
- fluid dynamics examples
- visco-elastic materials, with bond breaking
- single-body time-dependent force definitions in python
Version Alpha 1.0.15.2¶
- initial dissapative particle dynamics
- doc constant force, dpd
Version Alpha 1.0.15.1¶
Version Alpha 0.0.14.1¶
- added convenience methods to get spherical and cartesian coords from lists
- updated example models
- update docs
- added plot function in examples to plot polar angle velocity.
- code cleanup
Version Alpha 0.0.14¶
- All new FLUX / DIFFUSION / TRANSPORT, We’ve not got Transport-Dissipative-Dynamics working!!!
- secrete methods on particle to perform atomic secrete
- bug fixes in neighbor list, make sure neighbor don’t contain the particle
- bug fixes in harmonic potential
- new overlapped sphere potential
- new potential plotting method, lots of nice improvements
- new examples
- update become to copy over species values
- lattice initializers
- add decay to flux
- detect hardware concurrency
- bug fix in Windows release-mode CPUID crash
- multi-threaded integration
- all new C++ thread pool, working on getting rid of OpenMP / pthreads
- event system bug fixes
- documentation updates
Version Alpha 0.0.13¶
- preliminary SBML species per object support
- SBML parsing, create state vector per object
- cpuinfo to determine instruction set support
- neighbor list bug fixes
- improve and simplify events
- on_keypress event
- colormap support per SBML species
Version Alpha 0.0.12¶
- free-slip boundary conditions
- rendering updates
- energy minimizer in initial condition generator
- updates to init condition code
- initial vertex model support
Version Alpha 0.0.11¶
- new linear potential
- triagulated surface mesh generation for spheres, triangulate sphere surfaces with particles and bonds, returns the set.
- banded spherical mesh generation
- bug fixes in making particle list from python list
- points works with spherical geometry
- internal refactoring and updates
- Dynamic Bonds! can dynamically create and destory bonds
- lots of changes to deal with variable bond numbers
- rendering updates for dyanmic bonds
- particle init refactor
- added metrics (pressure, center of mass, etc…) to particle lists
- add properties and methods to Python bond API
- bond energy calcs avail in python
- bond_str and repr
- automatically delete delete bond if particle is deleted
Version Alpha 0.0.10-dev1¶
- bug fixes in bond pairwise search
- improved particle __repr__, __str__
- new style visible attribute to style to toggle visibility on any rendered object
- make show() work in command line mode
- internal changes for more consistent use of handles vs direct pointers
- bind_pairwise to search a particle list for pairs, and bind them with a bond.
- new points and random_points to generate position distributions
- spherical plot updates
- new distance method on particles
- implmement become – now allow dynamic type change
- big fixes in simulation start right away instead of wait for event
- basic bond rendering (still lines, will upgrade to cylinders in future
- render large particles with higher resolution
- new particle list composite structure, all particles returned to python in this new list type. fast low overhead list.
- major performance improvment, large object cutoff optimization
- numpy array conversion bug fix
- neighbor list for particles in range
- enumerate all particles of type with ‘items()’
- new c++ <-> python type conversions, getting rid of pybind.
- better error handling, check space cells are compatible with periodic boundary conditions.
- add start, stop, show, etc. methods to top-level as convenience.
- fix ipython interaction with show, default is universe not running when showing
- enable single stepping and visualization with ipython
- enable start and stop with keyboard space bar.
- pressure tensor calculations, add to different objects.
- new Universe.center property
- better error handling in Universe.bind
- clean up of importing numpy
- expose periodic boundary conditions to python.
- periodic on individual axis.
- new metrics calculations, including center of mass, radius of gyration, centroid, moment of inertia
- new spherical coords method
- frozen particles
- add harmonic term to generalized Lennard-Jones ‘glj’ potential
Version Alpha 0.0.9-dev4¶
- tweaks in example models
- more options (periodic, max distance) in simulator ctor
- add flags to potentials
- persistence time in random force
- frozen option for particles
- make glj also have harmonic potential
- in force eval, if distance is less than min, set eval force to value at min position.
- accept bound python methods for events
Version Alpha 0.0.9¶
- all new cluster dynamics to create sub-cellular element models
- cluster splitting
- splitting via cleavage plane
- splitting via cleavage axis
- other splitting options
- new potential system to deal with cluster and non-cluster interactions
- revamped generalized Lennard-Jones (glj) potential
- new ‘shifted’ potential takes into account particle radius
- updated potential plotting
- more examples
- fixed major integrator bug
Version Alpha 0.0.8¶
- explicit Bond and Angle objects
- new example apps
- new square well potential to model constrained particles
- bug fixes in potential
- thread count in Simulator init
Version Alpha 0.0.7¶
- lots of changes related to running in Spyder.
- force windows of background process to forground
- detect if running in IPython connsole – use different message loop
- fix re-entrancy bugs in ipython message loop.
- Spyder on Windows tested.
Version Alpha 0.0.6¶
- lots of changes to simulation running / showing windows / closing windows, etc..
- documentation updates
Version Alpha 0.0.5 Dev 1¶
- Add documentation to event handlers, and example programs
- fix bugs in creating event events
- add version info to build system and make available as API.
Version Alpha 0.0.4 Dev 1¶
- All new particle rendering based on instanced meshes. Rendering quality is dramatically improved. Now in a position to do all sorts of discrete elements like ellipsoids, bonds, rigid particles, etc…
- Implement NOMStyle objects. This is essentially the CSS model, but for 3D applications. Each object has a ‘style’ property that’s a collection of all sorts of style attributes. The renderer looks at the current object, and chain of parent objects to find style attributes. Basically the CSS approach.
- More demo applications.
- Memory bugs resolved.
Version Alpha 0.0.3 Dev 1¶
- Windows Build!
- lots of portability updates
- some memleak fixes
Version Alpha 0.0.2 Dev 5¶
- lots of new documentation
- reorganize utility stuff to utily file
- add performance timing info to particle engine
- add examples (multi-size particles, random force, epiboly, events with creation, destruction, mitosis, …)
- new dynamics options, include both Newtonian (Velocity-Verlet) and over-damped.
- new defaults to set space cell size, better threading
- New explicit bond object
- add creation time / age to particle
- particle fission (mitosis) method (simple)
- clean up potential flags
- harmonic potential
- new reactive potential to trigger (partial implementation)
- random points function to create points for geometric regions
- prime number generator
- Fixed major bug in cell pair force calculation (was in wrong direction)
- major bug fix in not making sure potential distance does not go past end of interpolation segments.
- new random force
- new soft-sphere interaction potential
- add radius to particle type def
- update renderer to draw different sized particles
- add number of space cells to simulator constructor
- configurable dynamics (Newtonian, Over-damped), more to come particle delete functionality, and fix particle events
- examples bind events to destroy, creation and mitosis methods
- new event model
Version Alpha 0.0.1 Dev 3¶
- Refactoring of Particle python meta-types, simpler and cleaner
- Upgrade to GLFW 3.3
- New single body generalized force system
- Berendsen thermostat as first example single body generalized forces
- Per-type thermostat
- Arc-ball user interaction
- Simplify and eliminate redundancy between C++ and Python apps.
Version Alpha 0.0.1 Dev 2¶
- First public release
Features to be implemented¶
Linux binaries
mouse interction – rotate, zoom simulation
Documentation
Event system to hook up simulation events to user objects
strike: User definable visualization style Nosé–Hoover thermostat
Destroying particles
Collision reactions (when particles collide, they react, and can create and destroy particles)
Particle mitois
strike: attach chemical cargo to particles strike: inter-particle flux of chemical cargo reaction-kinetics network at each particle
strike: Windows binaries Movable boundary conditions
strike: reflective boundary conditions (only have periodic now) mouse object picking
strike: Python API for bonded interactions (bonds, angles, dihedrals, impropers) strike: pre-made DPD potentials (conservative, friction, thermostat) strike: `With addition of particle chemical cargo, fluxes and above potentials, we will have complete transport-dissapative-particle-dynamics simulation. And reactions gives us reactive TDPD.`
strike: `Visualization: We will attach a ‘style’ attribute to the particle type that will let users define how they’re presented in the renderer. This will have attributes such as color, size, etc… We want to let users attach transfer functions here, that will read particle attributes, such as local chemical concentration and map this to a color. To get decent performance, we’ll have to compile user specified functions into pixel shader and run them on the GPU.`
strike: `Multi-process rendering. Jupyter notebooks now uses a two process model, so it’s probematic to create a window in a background process. We want to enable the simulator render to an off-screen buffer, and render to this buffer. Then copy the buffer to the foreground process and display it here.`
Known Bugs¶
strike: In ipython, closing the window does not work correctly strike: energy is not conserved – bug in integrator. - only a subset of features are implmented
Features for next release¶
- Python API for vertex model (polygons, cells)
References¶
bibtex_bibfiles = [‘refs.bib’]
| [AA84] | P B Armstrong and M T Armstrong. A role for fibronectin in cell sorting. J Cell Sci, 69(1):179–197, July 1984. |
| [AP72] | Peter B Armstrong and David Parenti. Cell sorting in the presence of cytochalasin B. The Journal of Cell Biology, 55(3):542–553, December 1972. |
| [Dor02] | Dov Dori. Object-Process Methodology: A Holistic Systems Paradigm. Springer, Berlin, 2002. |
| [FOBS14] | Alexander G Fletcher, Miriam Osterfield, Ruth E Baker, and Stanislav Y Shvartsman. Vertex Models of Epithelial Morphogenesis. Biophysical Journal, 106(11):2291–2304, 2014. |
| [OFPF+17] | James M Osborne, Alexander G Fletcher, Joe Pitt-Francis, Philip K Maini, and David J Gavaghan. Comparing individual-based approaches to modelling the self-organization of multicellular tissues. PLoS Computational Biology, 13(2):e1005387, 2017. |
| [SSS+14] | James P Sluka, Abbas Shirinifard, Maciej Swat, Alin Cosmanescu, Randy W Heiland, and James A Glazier. The cell behavior ontology: describing the intrinsic biological behaviors of real and model cells seen as active agents. Bioinformatics, 30:2367–2374, 2014. |
| [Ste70] | Malcolm S Steinberg. Does differential adhesion govern self-assembly processes in histogenesis? Equilibrium configurations and the emergence of a hierarchy among populations of embryonic cells. Journal of Experimental Zoology, 173(4):395–433, April 1970. |
| [SY86] | Robert E Strom and Shaula Yemini. Typestate: A Programming Language Concept for Enhancing Software Reliability. In IEEE Transactions on Software Engineering, 157–171. 1986. |
| [TH55] | Philip L Townes and Johannes Holtfreter. Directed movements and selective adhesion of embryonic amphibian cells. Journal of Experimental Zoology, 128(1):53–120, February 1955. |
Acknowledgements¶
We acknowledge generous funding from National Institutes of Health Grant 5U24EB028887-02