I met a problem for the usage of bond_harmonic. I wrote my own fix command to add specific atoms during simulation, by calling create_atom() in atom_vec.cpp. Additionally, my model also contains a pre-defined bond lists. I found that after adding several tens of atoms, the bond list seems to have problem. In the visualization, some bonds appeared on the newly-added atoms and some old bonds just disappear, which is not correct.
I guess this behavior is because my code make the bond list (might be static in memory) broken? But after I read the source code I still cannot locate the reason. I was wondering if someone with experience can let me know how the bond list works (what is the data structure and are they dynamic during simulation). Thank you so much!
There is more to adding atoms to a simulations than calling AtomVec::create_atom()
You can see this from other fixes that are adding atoms to the system, e.g. fix deposit, fix pour, fix gcmc etc.
As a minimum, you need to add the added atom to Atoms::natoms, assign it a new atom tag (one that has not yet been used), rebuild atom maps by calling Atom::map_init() and Atom::map_set(), and make sure the neighbor list rebuild is triggered on the next time step you add atoms.
Thank you for the very quick response! Now I called Atom::map_init() and Atom::map_set() after adding a atom, but timestep did not go forward after rebuilding atom maps. If I delete this two lines of code, the timestep works okay. The neighbor list rebuild is triggered on by “neigh_modify every 1 delay 0 check yes”. Could you give me a clue on how to solve this? Thank you!
Sorry, but without knowing more details there is nothing else that I can tell you outside of what I already told you. The function calls I mentioned are the minimum required. That doesn’t mean that they are sufficient. You have to look at the other fixes I mentioned to gain further insight. As I mentioned adding atoms to a parallel application using domain decomposition for parallelization and thus distributed data structures is not trivial.
Thank you again! I think this error is far beyond my coding knowledge, maybe related to deadlock or something. I am thinking a way to bypass this complexity, since I just want the bond working correctly. I have a question on atom map. Can I avoid rebuilding atom map by setting the atoms IDs like this, in type 1 (with bonds), set the beginning of IDs to a large number like 10000, then in group 2 (adding this type of atom during simulation), set the IDs from 0 to 10000. Is this a way to avoid rebuilding atom maps?
Sorry, but there are no shortcuts and I most certainly am not interested in promoting and helping people with ugly hacks that are bound to create problems sooner or later.
The fact that LAMMPS gets “stuck” means that you are not updating internal data correctly.
If you add atoms with bonds then you must understand the data model in LAMMPS, i.e. that bond information is stored with atoms and differently so depending on whether the Force::newton_bond flag is on or off. The list of bonds is constructed at each re-neighboring step from that information.
You also seem to be unaware of the fact that there are two types of indexing atoms: one for the local list of atoms (including ghost atoms) and then the global index. The atom map (in combination with the Atom::tag property) is required translate one to other. So if you add a new atom with a new tag, you have to rebuilt the atom map to be able to tell from the tag which local (or ghost) atom it is pointing to or whether that atom is not present in the current sub-domain.
Bottom line, you either have to teach yourself these details and improve your programming skills accordingly or give up on implementing your own fix.
Very appreciate for your hints and suggestions. I do have limited knowledge on how local atom index are mapped to global. I think the problem arises from I am not writing the adding atom function correctly. I checked the fix_deposit.cpp and have several questions:
I found they use pre_exchange() to add atoms. It it necessary? My code uses final_integrate() to add atoms. It seems works correctly but calling atom map rebuild would lead to stuck.
What should I do if the atom has bonus data? For example, AtomVecBody::bonus stores data about rigid body geometry. Currently I wrote a function to directly append data to the AtomVecBody::bonus array then link it back to atom->body[] array, but am not sure this is all I need.
Should I worry about ghost atoms when adding atoms? During a single timestep, when it is safe to clear ghost atoms and ghost bonus?
Fix::final_integrate() is definitely the wrong place
Fix::pre_exchange() is a suitable place, since at that point LAMMPS is ready to redistribute atoms and update the ghost atoms. When you add atoms you obviously also need to add the corresponding ghost atoms. That is most easily done during reneighboring. Please see the pseudo code for the flow of control during a timestep in 4.3. How a timestep works — LAMMPS documentation
Are you using atom style body? Please note that this is not used with fix rigid. Again, you are asking for advice without providing sufficient information. It is not possible to give specific advice. To add atoms or molecules or rigid objects with fix rigid requires additional special care.
Got it. I will move the function of adding atoms to pre_exchange(). Another question is, can I manually code in my fix to tell LAMMPS to update neighbor lists then call pre_exchange()? I checked the documentation for neigh_modify but did not find suitable keywords.
More specifically, I am using the Body package, modeling a rigid body system with rod shape (bacteria), whose length is growing. I hope the neighbor list can rebuild once a particle length exceeds a critical value, then trigger pre_exchange() to divide this particle into two daughter cells.
The situation is the inverse. Fix::pre_exchange() will only be called when the neighbor lists will be updated later. Please check out the pseudo code I pointed out already.