Hi all,
I noticed that if I make a 2-atomic molecule with a long (big R_0) FENE bond then the bond breaks if I don’t use a long enough cutoff in my non-bonded potentials. The minimal example is on my github. Sample output is in bad.out and good.out. The bad.out has the fail message about missing atoms. It failed after a few minutes while the good.* was able to run for 1h and I stopped it after. The only difference (that I see) is the cutoff r_c of the gaussian potentials in good/bad.def - the good.* has a longer cutoff. I had many similar runs before making this post and this looks reproducible - not having long enough r_c (looks like compared to R_{0, FENE}) leads to bond errors. All non-bonded potentials give negligible force even at the smaller cutoff. So the only effect of a bigger r_c I can think of is the construction of neighboring lists. I thought (for no particular reason except my physical intuition on how it should work physically) that bonded interactions were always computed regardless of the r_c. However, this does not seem to be the case:
The bond has R_0 \sim 16, the smaller r_{c,s} \sim 5 (bad) and the bigger is r_{c,b} \sim 9 (good). So if the atoms beyond r_c don’t feel each other even if they are bonded then it makes sense that atoms farther than 5 don’t feel the bond and thus fly away freely until something bad happens.
The effect is approximately reproducible (not to a step but the approximate physical time of the fail stays given the initial conditions, the seed, and everything), so I tried simulating almost until the effect and then dumping at every frame and the explanation above looks plausible: I see the pair that causes the crash and it is very far apart before the crash and its flying looks unconstrained short before the crash.
You can ask - what is the system that has such long bonds compared to unbonded interactions? My particles represent globular protein parts that are connected by a linker made of disordered protein. The non-bonded interactions are an attraction due to R_g \sim a \sqrt{N} and a repulsion due to \sigma \sim a N^{1/3}. At the same time, the maximum linker length is R_0 = L \sim a N, so it is expected to be bigger than all bonded interactions.
My questions are:
- Do bonded interactions have an r_c?
- Are bonded interactions accounted for beyond “global” r_c (the one reported by default at the beginning of a run)?
2.1) If no, then is there a way to enable them without increasing global r_c?
2.2) If yes, then what could be the reason for the crash in my example on github?
Bonds do not have a cutoff, but any bond that crosses subdomain boundaries has to have its atoms within the communication cutoff, which - by default - is the same as the “master neighbor list cutoff”, i.e. the longest non-bonded cutoff plus the neighborlist skin. If the second atom is not within the ghost region it will be “lost”.
You can boost the communication cutoff, but that can come with a serious performance penalty. For a coarse grain system that was somewhat similar than what you describe, we implemented pair style list, which is not subject to those limitations but now makes the interaction subject to minimum image conventions, i.e. the bond in the list must not be longer than half the simulation cell.
I see, thank you for the suggestion. I looked at pair_style list but it is for non-bonded interactions, right? Do you mean I could add a dummy interaction between atoms with long bonds so they don’t get lost?
Please read the documentation of pair style list again. It computes forces only between the pairs of atoms in the list. It does this in a unique way. It is not like other pair styles.
Sorry, I don’t think I understood.
Yes, I see that pair_style list needs an explicit list of pairs and their interactions like
15 259 lj126 1.0 1.0 50.0
15 603 morse 10.0 1.2 2.0 10.0 # and another comment
18 470 harmonic 50.0 1.2 5.0
19 332 quartic 10.0 5.0 -1.2 1.2
but it does not have a FENE bond that I am using for physical reasons (it behaves ~like an entropic polymer of finite length). pair_style list has harmonic and I could try to choose harmonic parameters to reproduce fene behavior as much as possible but it will not be great.
So was I suggesting that if I had a long bond between say atoms 1 and 2 then I could still have a fene bond between 1 and 2 but also add something like 1 2 harmonic 1e-10 0 100 in the pair_style list so that 1 and 2 are never lost (within r_c = 100). But how I understood your last answer is that this will not work and the interactions from pair_style list will not make other pair_style-s to keep track of selected atoms - is that right?
And another question about a more complicated molecule. I found a related thread here and you said there that “Perhaps, with setting special_bonds lj/coul 0.0 1.0 1.0 you can avoid it, but I am not certain.” and the guy answered that you were “spot on” but did not give much detail.
So, my questions:
- Do you now think your suggestion about effectively enabling 1-3 and 1-4 relieves the requirement for them to be inside
comm_modify cutoff Rc? I see you were “not certain”, so - did anything change?
I tried (using special_bonds fene is the same as you suggested, right?) and it seems to be working but I am not sure if I can rely on this enough to launch many tasks and not check each one for crushes manually.
- If doing
special_bonds fene does not really relieve the worries about 1-3 and 1-4, then can I not use special_bonds explicitly at all - will it help? It again seems to be working on a sample run but I am not sure if that’s just luck.
If you want a different functional form for the pairwise interaction, you have to program it into the C++ code of pair style list.
Correct, pair style list avoids the whole issue of having both atoms of a bond as local or ghost atoms in a subdomain, by computing each “half” of the interaction separately.
This does not apply to you, since your problem is quite different from what is mentioned in the referenced post.
At this point, I don’t know how to further assist you since I don’t understand well enough what you are doing and at the same time, you don’t seem to understand what I am trying to explain. Thus we are in a bit of a deadlock situation.
Since you say you now has something that “works”, then my only recommendation is: make extensive tests to confirm the results you are getting from LAMMPS are what you should be getting. This is something were you don’t need much understanding of LAMMPS’ internals. If this turns out positive, you can go on with your project.