I am running a polymer simulation where, in addition to the fixed topology of the polymer, bonds can be dynamically created during the simulation. Once created these bonds can slide along the polymer chain and are eventually destroyed at certain points.
The simulation is implemented using a C++ code that integrates the LAMMPS library. Using the lmp object, I invoke “fixes” to create and delete bonds, then simulate for some time before applying the next update.
The simulation runs smoothly for a while, but eventually (typically after around 2e8 steps, although this number varies), I encounter the following error: New bond exceeded bonds per atom.
I have already increased the maximum bonds per atom to 82. Interestingly, the error occurs even when there are fewer than 80 extra bonds created in the polymer. To diagnose this issue, I modified the create_bonds.cpp function to print the actual number of bonds per atom when the function is called.
Here’s the output I obtained when running the simulation on 4 cores just before the error when the code execute the fix “create bonds single/bond 2 2000 2021” :
Atom1=120 Atom2=912 MaxBondperAtom=82 NbondAtom1=1 NbondAtom2=1
Atom1=1748 Atom2=1862 MaxBondperAtom=82 NbondAtom1=82 NbondAtom2=37
Atom1=1284 Atom2=1432 MaxBondperAtom=82 NbondAtom1=37 NbondAtom2=46
Atom1=1923 Atom2=1945 MaxBondperAtom=82 NbondAtom1=51 NbondAtom2=52
As I have understood the atom numbers are the representation of these two atoms in the different cores and this should be fine, what I don’t understand how these such high numbers of bonds on the atom in the different processors pops up. Do you have any idea? Is it possible that the list of the bonds of the atoms are not updated correctly for some reason?
My possible guess is that the first processor, the one giving the output NbondAtom1=1 NbondAtom2=1, is the one that is “simulating” the two atoms while for the others these two atoms are ghost and maybe for some reason the bond list is not updated correctly in this case.
Do you have any idea how can I proceed further to debug my code?
It is impossible to debug code from remote without actually seeing the source.
But if you keep repeatedly calling the create_bonds command with the same atoms, you will be adding bonds to them, so it is up to your code to decide when and for which atoms it is adding bonds. LAMMPS permits to add bonds (or angles or dihedrals) to the same pair (or triplet or quartet) of atoms multiple times. So if you are overflowing the list of bonds per atom then that must be due to a bug in your code and not a LAMMPS issue.
I am calling repeatedly yes but at the same time I delete bonds when they slide over the polymer so in theory the topology is updated contimuosly (at least I meant it to be like this). But how do you see the fact that there are such high number of bonds in some cores and in one core just 1 bond for those atoms?
Well, based on the information you provide, this is the part of your code that you need to debug and make certain that you create the bonds and also delete the bonds correctly. A common problem of people doing this is that they confuse the local indexing with the global atom ID information.
Bond information is stored with atoms (by default with only one of the two, only with the command newton off with both) and the per-MPI process bond lists are constructed from them. Since atoms are distributed across processors through domain decomposition and since you say you create bonds in a specific region in space, then not having many bonds per atom outside that region makes perfect sense.
This is not easy to read and I don’t have the time to dig into it deeply, but there is an issue with how you delete bonds. You keep defining a group, but I don’t see how you check, if the group exists and if yes, first delete it. The group command in LAMMPS does not reset a group, but add atoms to a group. The fact that you add 1 to the atoms indices where you desire to delete the bond is suspicious, too, as this hints you are passing the local atom index instead of the atom IDs (atom->tag[i]) to the group command and that may mean that you attempt to delete a bond from a different pair of atoms than where you added it.
I add the +1 because atom indeces in my c++ code go from 0 to N-1 while lammps go to 1 to N as far as I have understood. Thus if I have to delete bond between atom 10 and atom 12 (in my indexing) I will call the command delete bonds for atom 11 and 13.
The atom indeces 11 and 13 are the global indeces as far as I have understood so it should be okay? No?!?
As for the group command I alway delete the group at the end of the function of delete_bonds in my code.
I know it is hard to debug, thank you for your hints!
This still doesn’t sound right to me. You have to store and use the atom tag property and if you do so then you should not need to one those. The atom IDs for arbitrary sets of atoms are not likely to be consecutive or start from 1. LAMMPS thus uses the atom->tag property to get the global atom ID for the local array indices of atoms (which range from 0 to nlocal-1) and then the atom->map() function to reverse this process, which will return -1 if an atom ID is not present on a specific MPI process and values >= atom->nlocal for ghost atoms.
Thanks and I understand your point, but then I don’t understand how the command create bond would work as a fix inside a normal input file lammps. As I have understood, you define the indeces of the atoms of the simulation in the file you read with the command read_data. There, the atoms are defined by indeces which are the ones to be called if you want to create a bond between them. Am I wrong? Otherwise how could I extract the information you say it is needed in the lammps input file?
You seem to be confusing fix bond/create and the create_bonds command here.
We are going in circles now. I cannot give you more detailed advice unless I more carefully study your code, for which I have no time. And you don’t seem to fully understand what I am talking about and probably need to spend some more time reading the developer information parts of the LAMMPS manual and/or the recent LAMMPS paper.
Yes probably I don’t understand what you said but at least as written in the lammps page of create_bonds:
The single/bond style creates a single bond of type btype between two atoms with IDs batom1 and batom2. Btype must be a value between 1 and the number of bond types defined.
The single/angle style creates a single angle of type atype between three atoms with IDs aatom1, aatom2, and aatom3. The ordering of the atoms is the same as in the Angles section of a data file read by the read_data command (i.e., the three atoms are ordered linearly within the angle; the central atom is aatom2). Atype must be a value between 1 and the number of angle types defined.
At least for single/angle it is explicitly written that the ordering of the atoms is the same as in the Angles defined in the data file. For bond it is not written but I guess it works the same.
The order of atoms in a bond doesn’t matter, you just look at the distance between two atoms.
What I suspect that your problem is, is not what you issue as command to LAMMPS, but rather how you determine which atoms need to form a bond and that bonds between which atoms need to be removed. As I already wrote it is highly suspicious that this list starts with 0.
Atom IDs need not start from 1 and need not be consecutive (when you lose or delete atoms).
So you can have a list of N atoms that you inspect, but you cannot derive the atom ID from a number of 0 to N-1. Instead to need to store the value of atom->tag for each of those and then pass this information to either the bond_create or group command for delete_bonds.
As I already wrote, I have given the advice that I am willing to give and I don’t want to keep repeating the same things over and over again. The problem is definitely a problem in your code and the symptoms suggest that you create bonds and do not delete them. How exactly this happens is for you to debug and sort out.
Just to close the thread I want to add that I checked at run time the number of bonds present in the systems and these are exactly the number I would expect, even before the simulation would crash the number of extra bonds are smaller than the number of extra bond per atom allowed. I also printed the pair of bonded monomers bonded with the new extra bonds and these match the underlying mapping I have in my c++ code.
If you would manage to create a LAMMPS input file that does what your program does (e.g. by tracking the commands you issue and writing them to a file. while that does not provide the functionality you are after, it will be a sufficient trace of what happens to LAMMPS), so that I can test what you do without having to go through your custom code. I would be willing to look through it and see if there is an oversight somewhere in LAMMPS (it is not likely, but we find (sometimes obscure) bugs all the time, so it would be foolish of me to assume that there are no bugs in LAMMPS).
Hi! Sorry for the delay. I had to zip the input file in lammps here as it is too big with all the commands (create bonds and delete bonds ecc for 1billion steps) and the initial condition chain4000.dat. I run this input file with lammps2021 and got the error as well. polymer.lam.zip (8.1 MB) chain4000.dat (405.9 KB)
Why didn’t you reduce the extra/bond/per/atom to trigger the crash faster?
I have set it to 10 for that reason yet my run has reached step 10875103 without a problem so far.
Update:
It is hours later and the simulation has reached step 165996641, i.e. over 15x more than before and despite lowering the maximum number of bonds per atoms, there has been no error.
It usally happened to me around 2e8 steps but yes for 80 bonds. It might be the version of lammps indeed. I will try to launch it again with the new version, I haven’t tried it yet. Thanks!