Flex
Simulation of a flexible cell membrane
Flex is an application to simulate particle dynamics. Particle Dynamics refers to systems composed of objects with physical properties of mass, position, velocity, force, etc., where all of the mass is located at a single point. This simplifies the simulation of the system by eliminating torques and rotations. However, the particles can have a simple extended boundary to allow for collisions between particles as well as particle-boundary interactions.
Although the main goal of the project is to simulate a flexible cell membrane that can be ruptured with increasing pressure, the application has been designed so that more general spring-connected particle models can also be generated, loaded and simulated. To facilitate the use and investigation of such models, including simple test models consisting of only one or two particles, the project includes a command-line application, Flex, which can be used in conjunction with the graphical app FlexGui. The command-line application supports a number of options for generating various kinds of models and running them, as described below.
As in previous projects, the code is partitioned in modules: numerics, util, model and view. In this project the util model contains a number of utility classes that support XML file i/o, command-line option parsing and numeric checks.
The Flex project implements a number of numerical integration algorithms which provided varying degrees of accuracy and performance. These include Runge-Kutta 4 (rk4), Euler-Richardson and Verlet. The goal is to enable simulation of a simple cell membrane that is represented by spherical particles connected by springs. Consequently, the solver must properly handle spring forces in addition to collision forces.
Testing of the numerical algorithms is divided into two components: a command-line application and a GUI application. These two applications complement each other in the development process. The command-line application enables script-driven testing of models with generation of test results that can be easily compared against reference data. The GUI app provides a visual display of particle and boundary positions and interactions that facilitates determination of suitable initial conditions and simulation parameters.
The GUI app is called 'FlexGui'. It displays a graphical representation of particle movement, along with some input and output information, as seen in the screenshot above.
The view classes are similar to those in the previous demo, Leap:
The main difference is the replacement of the phase plot panel by the ParticleGraphics panel. This class is responsible for drawing the particles, springs and system boundary. As in the Leap project, the function worldToView() converts model position coordinates to panel window coordinates.
In addition to the model components, the ParticleGraphics class can optionally draw DebugVector objects. These are simply lines with a length and origin that represent the application of various kinds of forces. The storing of DebugVector objects must be enabled by a call to ParticleModel.storeDebugVectors().
The main class in the model package is ParticleModel. This class is similar to the model class in previous projects. Its main component is ParticleSystemState, which contains the position and velocity variables for the particles in the system:
Other important components in ParticleModel are two arrays, one for Particle objects and one for Spring objects. The Particle class contains all persistent data for a specific particle other than position and velocity (in this case, just mass and position). A Spring object defines a connection between two particles along with an equilibrium distance.
A useful feature in this project is the ability to interconvert the model data to and from a readable text format, a simple XML format. This allows for the definition and inspection of simple test models for debugging purposes. The conversion is achieved through an interface, XmlObject, that defines two methods: toXml() and fromXml(). The xml conversion makes use of standard classes in org.w3.dom and javax.xml. The actual file i/o methods are provided by an associated utility class XmlFile.
The significant numerics classes and related classes are shown below:
The numerical integration of motion is calculated by methods in the Numerics class in the numerics module. The methods and interface are similar to those used in previous projects such as bounce and leap, but with some important differences, as described in more detail below.
The new classes to note include the Vector class, which combines x and y values so that they can be treated as a single state variable in the SystemState. Another new class is VectorList. This can contain a list of Vector objects, just like the new SystemState, and is used by algorithms such as rk4, in which a set of intermediate values (k1, k2, k3, k4) need to be stored. In previous projects, another SystemState object was used to hold such values, but since that was a 'bogus' use of the SystemState class, a new, more-appropriately-named class was created for the purpose.
The same Function class is present, although with a modification to the function signature. The two model classes that define the derivative functions for position and velocity are also shown here as derived classes from this class.
One motivation for changing the numerical integration interface is for performance reasons. More specifically, a particle dynamics simulation with a large number of particles will correspondingly requires a large number of state variables, two for each particle (position and velocity). In addition, there are only two derivative functions, one for position and one for velocity, that are shared by all the particles.
Consequently, the previous way of passing derivatives to the numerical integration method in an array, one for each state variable, is redundant and wasteful. Instead, an array of just two derivative functions is passed in, and the appropriate function is obtained and applied to each state variable in succession:
// rk4, step 1:
curState.calcDerivatives( dt );
VectorList f1s = new VectorList(); // eval of deriv at time t, using state at time t
SystemState midState1 = curState.createNew();
for( int k = 0; k < curState.size(); ++k )
{
// calc F1
Function derivFn = (Function) stateDerivFns.get( k % numDerivFns );
Vector ydot = derivFn.eval( curState, k, t );
Vector f1 = ydot.mul( dt );
f1s.addVariable( f1 );
// add midstate val
Vector yc = curState.getVariable( k );
Vector ym = yc.add( f1.mul( 1.0/2.0 ));
midState1.addVariable( ym );
}
A couple of other changes to note are 1) the use of a Vector data type in place of a double, so that both x and y coordinates are contained within a single state variable. 2) the index of the state component is now passed as an additional parameter to the eval() function. This is needed to calculate the correct derivative value for each particle.
Another interface change was needed to calculate the velocity derivative function (acceleration). For this kind of application, the acceleration is determined by the sum of a set of forces, rather than a simple analytical equation. Since this force summation calculation needs to be done for every integrated state in the numerics method, a new method was added to the SystemState interface, 'calcDerivatives()'. This must be called for every state prior to making use of the derivative functions. In the implementation of this method, the sum of all forces on each particle are calculated and used to determine its acceleration.
Also, because the derivative functions need access to objects of type ParticleSystemState, which is derived from SystemState, the numerical integration method cannot simply create a new SystemState, but instead obtains a SystemState reference to a new ParticleSystemState object by calling curState.createNew().
A further optimization was to only return the final integrated state to the model instead of an array of all states, as these are normally just discarded by the model anyway.
| Name | Descrip |
|---|---|
| singleWall | single particle moving towards wall |
| pairWall | two particles moving towards wall |
| pairSpring | two particles connected by a spring |
| openMembrane | ring of particles connected by springs, with free end |
| closedMembrane | ring of particles connected by springs, closed |
The command-line app is called 'Flex'. It can run a simulation with a variety of options, specified by command-line arguments:
mhackslab:mdperf $ java Flex -help
usage: Flex [-option [optionVal]] [...]
options:
-genModel <S> generates model <S>
option value is one of the following:
singleWall
pairWall
pairSpring
openMembrane
closedMembrane
-integMethod <S> (=rk4) integration method to use for particle dynamics
option value is one of the following:
euler
euler-rich
rk4
-outFile <S> (=model-out.xml) model output file
-numBodies <N> (=10) number of bodies in the model
-simTime <N> (=10) run app until this simulation time limit (in seconds) is reached
-inFile <S> model input file
-help show usage
mhackslab:mdperf $
Below is a sample output that is produced by running the Flex command-line app on a 'singleWall' model for a sim time of 40 seconds:
mhackslab:mdperf $ java Flex -inFile doc/singleWall.xml -simTime 40 --- params --- numBodies 10 simTimeLimit 40.0000 integMethod euler-rich timestepsPerUpdate 100 timestepLen 0.00100000 --- error stats --- beg KE = 0.500000 end KE = 0.499827 err KE = -0.000172820 --- perf stats --- run time = 0.112000 sec sim time = 40.1000 sec perf = 358.036 bss (body*sim_sec/sec) mhackslab:mdperf $
In addition to the report on the command-line, the app also writes the state of the model to the file 'model-out.xml'. The contents of this output file show the resulting particle position and velocity:
<?xml version="1.0" encoding="UTF-8"?> <ParticleModel integMethod="euler-rich" timestepLen="0.0010" timestepsPerUpdate="100"> <boundary> <pos x="-25.0" y="-25.0"/> <dim x="50.0" y="50.0"/> </boundary> <particles> <Particle mass="1.0" radius="5.0"/> </particles> <ParticleSystemState> <integVars> <pos x="1.056446985940132" y="0.0"/> <vel x="-0.9998271652722177" y="0.0"/> </integVars> </ParticleSystemState> </ParticleModel>