Dear All,
I want to do a simulation for many proliferating cells. So I need to
keep track of the cell parts(made up of atoms), and during the
I knew somebody would attempt lammps for something like this
eventually. Cool! Extra points!
If LAMMPS does not have what you need, you may have to write your own
code (LAMMPS fix) just to get the particles in your system to move
semi-realistically (to say nothing of making the cells divide).
The movement of cells can be pretty complex stuff. So when I saw your
post, I got worried you were just going to blindly use some kind of
simple integrator like Langevin (Brownian) dynamics (or SRD, or SPH,
or even USER-SELM).
I guess motility is actively driven in various ways by motors which
consume ATP. Movement is not random. It happens in response to
gradients in background nutrient concentration (for example, which
could depend on the presence of other cells) and perhaps other
material (created and left behind by other cells). Cells also have
internal states (eg "straight" or "tumbling" states in e-coli) which
change stochastically over time:
http://en.wikipedia.org/wiki/Motility
http://en.wikipedia.org/wiki/Biased_random_walk_(biochemistry)
(There are also hydrodynamic effects which people obsess about.)
Perhaps in your case, you don't need to worry about any of these
fine-grained details, but I would check what other people have
modelled cell movement before proceeding. Forgive me if you know all
this better than I do (probably) and have done this already.
simulation I have to create many new atoms and many new bonds in a way
that I am able to remove them any time I need.
My question is: as far as I could do with lammps, it is not so easy to
create bonds in the middle of the simulation by scripting.
Not to say anything against LAMMPS, but some software like ESPResSo++
and OpenMM have fancy python-based scripting interface which makes it
easy to do things like run for N steps, check something, and depending
on the outcome invoke a function which adds indivudual atoms and
bonds. Presumably you could define a molecule as a python class whose
constructor invokes the necessary ESPResSo++ or OpenMM commands for
creating the atoms and bonds you need each time a new molecule-object
(cell) is created.
If you wanted to do the same thing in LAMMPS, it is harder. You would
have to delve into the LAMMPS code and see how Paul Crozier
implemented fix gcmc (at least that's what comes to my mind,
fix_gcmc.cpp).
-- wish list --
It would be cool if the LAMMPS C++ code internally had a simple
function for creating a single new atom (or bonded interaction), or a
new molecule (a copy of an existing one) and moving it into position.
That would make it easier for users to make these kinds of fixes. (A
related problem is the automatic generation of angle interactions when
new bonds are added. This is a more difficult problem. I don't know
if ESPResSo++ or OpenMM can do that either. That is asking a lot.
But a molecule duplication function should not be too hard to
implement. It's not clear if we need this, so I hesitate to suggest
to Paul or Steve to implement any of these features.)
Do you
suggest me to continue using lammps with scripting. Or write my own
code beside LAMMPS to control the creation/removing of bonds/atoms
very well, or forget lammps for this use and wrrite my own code(which
I dont like to do)?
If you do not need to divide cells too frequently, and you can
tolerate a little bit of file IO, I think the easiest thing to do
(quick and dirty) is to write a loop which will periodically halt the
simulation and writes the coordinates (using the "write_data" or
"dump" commands, for example). In each iteration of the loop use the
"shell" command to invoke a program that will read the old data file
and build a new one (with additional molecules, atoms, bonds etc. as
needed).
http://lammps.sandia.gov/doc/jump.html
http://lammps.sandia.gov/doc/shell.html
http://lammps.sandia.gov/doc/read_data.html
http://lammps.sandia.gov/doc/write_data.html
I have used a combination of these LAMMPS input-script commands with
my own DATA-file builder to implement my own exotic Monte-Carlo steps.
(My monte-carlo lammps script is not very elegant or fast, though.
I'm happy to post it as a proof-of-principle.)
As for writing your own script which can read and write molecules,
it's probably best that you do this yourself. There's an example
"create_new_molecule.sh" script below, but I suggest you create your
own instead. (The moltemplate way is general, but it is probably very
slow.)
--- moltemplate ---
I certainly did not intend this email to be another advertisement for
moltemplate. moltemplate is a molecule-building tool used for setting
up a simulation. It is probably too slow invoke it frequently during
a simulation. Nevertheless, I'll show you the moltemplate way because
knowing one way of doing this might help you think of your own better
way.
If you are already have a "system.lt" file describing the molecules in
your system, then you could write a bash script which appends a line
to your "system.lt" file (telling it to create a new molecule) and
then run moltemplate again to create a new data file with that
molecule:
------- create_new_molecule.sh ----------
# extract out the coordinates of the data file you (presumably) just created
# in LAMMPS using the "write_data old.data" command
# Save these coords in the file "old_coords.raw" (simple 3-column format)
extract_lammps_data.py Atoms < old.data \
> awk '{print $1" "$5" "$6" "$7}' \
> sort -g | awk '{print $2" "$3" "$4}' \
> old_coords.raw
# (I assume you are using "atom_style full". If not, change
# the column numbers accordingly. "extract_lammps_data.py" is
# a script which comes with moltemplate, but you could use pizza.py)