dihedral ownership

I am trying to modify fixes bond/break and bond/create for a special case, but I have trouble understanding the ‘ownership’ for dihedrals. I have a game plan, any insight is appreciated.

Consider the case ‘newton on’ for simplicity.

I understand every dihedral involves 4 atoms but is only owned by 1 atom.

I am inferring that the owner atom is the ‘central’ atom. Is that so? (There is no real documentation on this, which part of the code would I have to read? )

If so, which one is the central (owner) atom? For an improper it is the 2nd atom, but how about a proper dihedral? Is it ALWAYS the 2nd atom?

I have system with a periodic topology. I can infer all connectivity for an atom from only the atom tag. I would like to do the following.

When I break a bond, I know this is the central bond of 2 dihedrals and it is involved in no other dihedrals.

I would like to ‘disable’ these dihedrals by zeroing atom->num_dihedral[i] for the dihedral-owning atom, but wihout touching atom->dihedral_atomX[i] values.

I will decrement atom->ndihedrals accordingly.

There is a good chance I will re-create the same bond and dihedrals later on in the simulation.

After creating the bond,

I will only revert atom->num_dihedral[i] to the original value.
I will increment atom->ndihedrals accordingly.
I am hoping atom->dihedral_atomX[i] values will then already be there because I never removed them in the first place.

This way I am hoping to flip the dihedrals on/off without actually adding/removing them which is completely unnecessary in this case.

I have no use for special neighbours, so I will not be updating them.

Is this a bad way of implementing this, and if so, why?

Thank you very much.

Murat

I am trying to modify fixes bond/break and bond/create for a special case,
but I have trouble understanding the 'ownership' for dihedrals. I have a
game plan, any insight is appreciated.

Consider the case 'newton on' for simplicity.

I understand every dihedral involves 4 atoms but is only owned by 1 atom.

I am inferring that the owner atom is the 'central' atom. Is that so? (There
is no real documentation on this, which part of the code would I have to
read? )

If so, which one is the central (owner) atom? For an improper it is the 2nd
atom, but how about a proper dihedral? Is it ALWAYS the 2nd atom?

yes. see Atom::data_dihedrals() in atom.cpp

I have system with a periodic topology. I can infer all connectivity for an
atom from only the atom tag. I would like to do the following.

connectivity follows the bond topology. everything else is doomed to be wrong.

When I break a bond, I know this is the central bond of 2 dihedrals and it
is involved in no other dihedrals.

I would like to 'disable' these dihedrals by zeroing atom->num_dihedral[i]
for the dihedral-owning atom, but wihout touching atom->dihedral_atomX[i]
values.

*very* bad idea.

I will decrement atom->ndihedrals accordingly.

There is a good chance I will re-create the same bond and dihedrals later on
in the simulation.

After creating the bond,
I will only revert atom->num_dihedral[i] to the original value.
I will increment atom->ndihedrals accordingly.
I am hoping atom->dihedral_atomX[i] values will then already be there
because I never removed them in the first place.

nope. not going to work. totally incompatible with parallel runs and
domain decomposition anyway.

This way I am hoping to flip the dihedrals on/off without actually
adding/removing them which is completely unnecessary in this case.

I have no use for special neighbours, so I will not be updating them.

Is this a bad way of implementing this, and if so, why?

why can't you use the existing feature in fix bond/break and fix bond/create?

axel.

Thanks for the response. I cannot use the current fix bond/create because I need to create 2 dihedrals of distinct type and specific order for each bond that is created.

I don’t understand what you mean by “connectivity follows the bond topology. everything else is doomed to be wrong.”. Let me briefly explain my system.

I have n total atoms and you ask me about the i’th atom up until n/2, I know exactly that this atom is bound to atoms i-1, i+1 by harmonics and atom n-i by a breakable bond.

I also know that two dihedrals exist : (i-1, i, n-i+1, n-i) of type 1 and (i+1, i , n-i+1, n-i+2) of type 2. Dihedrals go away if the bond (i, n-i) is broken, they come back if it is created. There is nothing else in the system, it is that simple. So when I have i and n, I have complete knowledge of possible connectivity for i’th atom.

I don’t think I should be parsing the topology for such a system. I am not using any special neighbours or at least I can do away with them (take weights 1 1 1). In the worst case I can hard code the bonds and dihedrals to be created. I was hoping that I could do even better by flipping them on/off rather than destroying/creating them.

I can make sure all related atoms have ghosts in the related processor. This is straightforward by the way bond/create works if I set a large enough skin distance.

This is a stiff polymer (dna) with a very high relaxation time so I am striving for performance. I don’t want to do things if I don’t absolutely have to.

I don’t see why this is “totally incompatible with parallel runs and domain decomposition”. What breaks it exactly? Are dihedral_atomX arrays migrated based on if num_dihedral values? If so, I see your point, otherwise, I don’t. (I can’t find it in the code)

You talk as if this is non-sense. I believe either you are correct or I am not communicating my problem properly. My breakable bond/dihedral implementation was badly designed before and I understood when you warned me. This time however I dont’ understand what I am missing. I am new to parallel programming but I am enjoying studying lammps, so forgive me if this is actually non-sense.

Thanks for the response. I cannot use the current fix bond/create because
I need to create 2 dihedrals of distinct type and specific order for each
bond that is created.

I don't understand what you mean by "connectivity follows the bond
topology. everything else is doomed to be wrong.". Let me briefly explain
my system.

I have n total atoms and you ask me about the i'th atom up until n/2, I
know exactly that this atom is bound to atoms i-1, i+1 by harmonics and
atom n-i by a breakable bond.

​please keep in mind that in LAMMPS "i" is not guaranteed to remain in the
same order. atoms migrate between domains and the i index only represents
the number of local (and ghost) atoms.

I also know that two dihedrals exist : (i-1, i, n-i+1, n-i) of type 1 and

(i+1, i , n-i+1, n-i+2) of type 2. Dihedrals go away if the bond (i, n-i)
is broken, they come back if it is created. There is nothing else in the
system, it is that simple. So when I have i and n, I have complete
knowledge of possible connectivity for i'th atom.

​well, it is not exactly that simple. you only get the atoms you want
through using a map. i.e. you need to look at atom->map(atom->tag[i]-1) and
so on.
you also have to make certain, that your communication cutoff i large
enough to have sufficient suitable ghost atoms included.

I don't think I should be parsing the topology for such a system. I am not

using any special neighbours or at least I can do away with them (take
weights 1 1 1). In the worst case I can hard code the bonds and dihedrals
to be created. I was hoping that I could do even better by flipping them
on/off rather than destroying/creating them.

if the bond topology will always remain the same with only bonds turned on
or off, then this process is best handled via setting the involved bond and
dihedral types to negative or positive, respectively. this way, there is no
problem with looking up old data. if the connectivity changes, then it
would be cleaner to do things properly. i have learned in long years that
making some optimistic assumptions to save time and effort rarely pays off.
you are more likely to end up with a piece of code that may be difficult to
debug and create problems at a time when you least expect it and will
create the most damage.

I can make sure all related atoms have ghosts in the related processor.
This is straightforward by the way bond/create works if I set a large
enough skin distance.

This is a stiff polymer (dna) with a very high relaxation time so I am
striving for performance. I don't want to do things if I don't absolutely
have to.

I don't see why this is "totally incompatible with parallel runs and
domain decomposition". What breaks it exactly? Are dihedral_atomX arrays
migrated based on if num_dihedral values? If so, I see your point,
otherwise, I don't. (I can't find it in the code)

yes, all "owned" data is migrated between domains, please see the
correponding AtomVecXXX:pack_exchange() methods.​

but my objection goes beyond that: *any* code that relies on the fact that
data when flagged as unused retains its value is fundamentally broken.​ it
is a very bad idea and bad programming style to boot.

You talk as if this is non-sense. I believe either you are correct or I am
not communicating my problem properly. My breakable bond/dihedral
implementation was badly designed before and I understood when you warned
me. This time however I dont' understand what I am missing. I am new to
parallel programming but I am enjoying studying lammps, so forgive me if
this is actually non-sense.

​this is not so much an issue of parallel programming, but understanding
the fundamental data model in LAMMPS, which is based on domain
decomposition and atoms taking "their" data with them (which is then used
to recreate derived information like the list of bonds etc).

​the more i think about how you describe your model, the more ​it seems to
me, it might also be implemented through a pair style used with
hybrid/overlay. check out pair_style list as an example. you would just
pass through a list of all possible bonds to flag all of the bonds that are
currently active and then in a second pass process this list to compute
force contributions to from the active bonds and the resulting active
dihedral contributions. of course, you'd need to have a way to restart the
bond/dihedral status, but that could be added rather easily.

​axel.​

I think setting bond and dihedral types to negative will actually do the trick for me. That looks exactly like what I needed.

Regarding the data model of lammps, I have read the developer's manual and 2 very helpful blog posts which I will include here as reference for others. Are there any other documents that describe this model? I will take a deeper look at the pack_exchange() methods.

Thank you very much for the input.

Murat

Blog posts :
https://sites.google.com/site/scienceuprising/code-packages/lammps/a-dissection-of-lammps-classes
http://kirilllykov.github.io/blog/2012/10/13/writing-fixes-for-lammps/

I think setting bond and dihedral types to negative will actually do the
trick for me. That looks exactly like what I needed.

Regarding the data model of lammps, I have read the developer's manual
and 2 very helpful blog posts which I will include here as reference for
others. Are there any other documents that describe this model? I will
take a deeper look at the pack_exchange() methods.

having a suitable developers guide is one of the TODO items that are
increasingly important as LAMMPS is growing and particularly the
number of "downstream developers". the pdf on the LAMMPS home page is
at best a starting point, and it has to be fairly vague so that it
doesn't mislead people when specific details of the implementation
change. a better place for a suitable documentation would be to embed
this into the source code and then extract it with a tool like
doxygen.
however, that would be a *MASSIVE* undertaking that would require a
larger number of people.

we have had discussions in the past to try and organize "LAMMPS
developer meetings" that try to address such "housekeeping" tasks, but
would also offer the opportunity for people actively developing with
LAMMPS to discuss their project with other (more or less experienced)
developers. we had something along the lines of the latter set up this
spring as part of a scientific software development workshop. it was
quite interesting, stimulating and productive (at least as far as i
could tell).

i am curious whether there would be a critical mass of people
interested in setting up such an activity again, perhaps for a longer
period (1 or even 2 weeks) where then the deal would be that
participants get training in LAMMPS implementation specifics and in
exchange the entire group would work on some large refactoring effort
or some other related project, e.g. embedding a suitable
implementation guide (quite a bit of the information is already
present as comments).

axel.