Re: [lammps-users] request: allow fix bond/break to change atom types
Andrew Jewett <jewett@...1937...>
Thu, 6 Jul 2017 18:42:47 -0700
On Thu, Jul 6, 2017 at 7:19 AM, Steve Plimpton <sjplimp@...24...> wrote:
> Bad idea to try to invoke an extra neighbor rebuild from inside your fix.
> Instead, set the next_reneighbor flag (see fix.h) to trigger a neighbor
> rebuild in the normal timestep. That comes immeditately after
> post_intergrate() anyway. That's what fix bond/create does.
Dear Steve, Jacob, and Tim
Thanks for making code suggestions and sending me example code.
I'll try Jacob's version get back to you if it works. (I tried
something similar to Jacob's solution without success, but my version
did other things was much more complicated. My fault. I should
always start with something simple.)
I've been struggling with perplexing crashes since I started toying
around with these kinds of fixes a week ago. I probably should not
have attempted it. Forgive my impatience. I've got a hard deadline
and I'm not supposed to be working on this stuff at all. So it's
fantastic to have examples of code that does something close to what
I'm trying to do. I really appreciate it.
P.S. If I use any of your code, I will contact you before any paper
submission that use it are made. Perhaps we can also work together to
make a version of fix bond create/break which can handle all of the
goals that we have, or at least share as much code as possible.
(Jacob, have you had a look at Tim's "topo.cpp",
"fix_bond_create_rxn.cpp", and "fix_bond_break_rxn.cpp" files?)
Here are some things to think about in the future:
1) There are some cases where it's very convenient for the changes
made by one fix need to be finalized before the next one begins,
(without any atom movement in between, for example to prevent an atom
from drifting too far after a bond break). When this happens, is
there a way to tell LAMMPS that you don't want to move any atom
positions for the next N verlet iterations? (Where "N" is the time
interval between when the fix is applied)
2) Minimization following topology changes:
It would be great if there was a way to define a dynamic group of
atoms which were "recently_influenced" by the formation or deletion of
nearby bonds, and send only that group of atoms to be minimized.
(Perhaps this could be done using a for-loop that run dynamics for N
steps, and then invokes "minimize" combined with "setforce" to
immobilize atoms which aren't in the "recently_influenced" group).
The sudden changes in energy caused by creating a new angle
interaction can be surprisingly large. An alternative approach might
be to slow them down temporarily by changing their mass or friction
coefficients for the next N verlet steps.
(I like to pontificate.)
> Re: limitations of fix bond create/delete. Again I suggest conferring
> with Tim and Jake - they have been working on a bigger picture
> solution to these kinds of issues.
> On Wed, Jul 5, 2017 at 8:58 PM, Andrew Jewett <jewett@...1937...> wrote:
>> Again, I'm hoping to bribe anybody who helps me solve this problem with
>> authorship on the "moltemplate" paper (the contents of which is also related
>> to this topic). I think this will be a well cited paper. More details
>> On Jul 5, 2017 7:03 AM, "Steve Plimpton <sjplimp@...24...> wrote:
>> Hi Andrew - you could look at fix bond/create, which
>> has some logic to change an atom type under certain
>> Thanks for the prompt reply. Detailed comments below...
>> In principle changing atom types for atoms in a bond
>> should be OK to do. You need to insure that all other
>> procs know about the change, e.g in their ghost or owned
>> atoms. The forward comm should accomplish that
>> for ghost atoms.
>> Thanks? That's very reassuring to hear.
>> 1) Is it possible to forgo all this simply by invoking "neighbor->build()"
>> at the end of the fix's "post_integrate()". (I can worry about optimizing
>> speed later.).
>> 2) If so, then, after modifying the contents of the local atom->type
>> array, would I need to do anything else before doing this?
>> Interestingly, un/pack_forward_comm() from "fix_bond_create.cpp" does not
>> pack or unpack atom type information (just partner lists, probabilities, and
>> special bond information).
>> 3) (I'm curious, how does it transmit atom type changes to other procs?)
>> If you had one proc change the
>> type of a ghost atom owned by another proc, then
>> you'd need a reverse comm as well, to propagate that.
>> Good point.
>> I will have to modify both atoms in the bond, and one of them could be
>> owned by another processor. Unfortunately neither fix_atom_swap.cpp, nor
>> fix_bond_create.cpp can do this.
>> 4) Actually, I was thinking of reporting this as a bug in
>> fix_bond_create.cpp: Users can specify either "iparam" OR "jparam", but not
>> both. (If they specify both, then only one of the atoms will be changed.)
>> This significantly limits what "fix bond/create" can do, and makes it makes
>> it practically impossible to use it for many polymerization problems (for
>> (I also have a selfish motive to complain about this: If someone solved
>> this issue, then I could probably steal that code and put it in the fix I am
>> working on. I might post a separate question about this on the list...)
>> Either way, I am interested to learn of any other fixes I can study which
>> can modify properties of atoms belonging to other processors.
>> Thanks for engaging me in a discussion on this topic. I think it is
>> really important (and not just for my own work). LAMMPS already has
>> fantastic capabilities and infrastructure. With some really trivial changes
>> like this, LAMMPS will be able to do exiting, amazing new things.
>> I'm CCing this to Tim Sirk who was also modifying trying to modify
>> fix_bond_create.cpp. (Hi Tim)
>> Re: bigger picture, you should talk to Jake Gissinger, CCd.
>> He's been working on a general "fix react" command to
>> allow more complex reactions than fix bond break/create allow,
>> and your model might fit into that framework.
>> On Wed, Jul 5, 2017 at 7:36 AM, Andrew Jewett <jewett@...1937...> wrote:
>>> Suppose I want to modify "fix_bond_break.cpp" so that whenever a bond
>>> is broken, both atoms from the bond are also changed to atom type 2
>>> (for example). How would I go about this?
>>> --- Goal ----
>>> Actually, I have already written a much more general fix (based on
>>> "fix_bond_break.cpp") which allows LAMMPS users to simulate
>>> "intelligent molecules". These molecules can make or break bonds and
>>> make arbitrarily complex decisions to alter their properties in
>>> response to contact with other molecules. (Essentially, this fix
>>> allows LAMMPS users to run arbitrary cellular automaton on networks of
>>> bonded atoms. There's a Conway's game of life example which uses this
>>> The reason we created this fix was to study coarse-grained machines in
>>> cells which are driven out of equilibrium and exhibit complex behavior
>>> (such as microtubules which are dynamically unstable, as well as
>>> proteins which read and walk along DNA controlling gene expression and
>>> cell fate). Because the fix is general, the hope is that anybody
>>> could use it to model new biomolecules we don't understand yet. We
>>> want to eventually get many biologists interested in using LAMMPS.
>>> Alas, I have not released this fix because when I run it in parallel,
>>> it crashes, (unsurprisingly). But it only seems to crash when I ask
>>> it to alter the atom types. Rather than post that code, just knowing
>>> the answer to the question above would be incredibly useful.
>>> Thanks to everyone who read this far.
>>> P.S. We are currently writing a paper on this, in case anyone wants be
>>> part of that. Help on this issue would definitely earn authorship.
>>> P.P.S. More specific questions (feel free to ignore these):
>>> 1) I'm curious: Is it a problem to change the atom types from within
>>> FixBondBreak::post_integrate()? (I ask because I noticed that
>>> "fix_atom_swap.cpp" places the code that changes the atom types inside
>>> "pre_exchange()" instead)
>>> 2) Is there anything else I have to do after changing the atom types
>>> besides invoke comm->forward_comm_fix() and update
>>> "un/pack_forward_comm()"? (To make it easier, let's assume, for the
>>> moment, that the distance c.utoff is the same for all atoms types.
>>> I'll worry about that issue later.)