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
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
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
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.