Bond atoms missing occurs when the model is doubled using the same force field and modeling method

Hello all.
I know that bond atoms missing occurs when the initial configuration or the force field parameters are unreasonable. I first simulated a smaller system and successfully ran 2ns of NVE and 5ns of NVT, but when I used the same modeling method and force field and only doubled the number of molecules in the system, a bond atoms missing error occurred when 2ns of NVE was completed and NVT reached about 2ns. I don’t know what to do. Hope someone can help me, thanks a lot. (6.8 MB)
PARM.lammps (9.7 KB) (722 Bytes)

From what I see, the initial geometry is reasonable. Probably it’s a problem with your force field: the lost atoms occur in molecule 4087 and involve a carbon and hydrogen atom. Maybe one of the angle, dihedral, or improper terms affecting the pair has a small force constant and that produces instability at some point. The fact that you observe the error on a larger system is only due to statistics. If the simulation is long enough, it may happen in the small system too.

Do you mean that I need to replace the force field of the molecule that produces bond atoms missing or replace the force field of entire system?

It is not a matter of replacing bits of the force field. Ideally, the whole system should be described with the same approximations. Since the piperazine molecule produces the error, I would check which terms are used to describe it and whether they were correctly assigned and specified, in particular, the angular and dihedral terms.

I added the force fields of these organic compounds by LigParGen, so I don’t know how to check them. And I did not change the force field, only changed the method of calculating long-range forces from pppm to ewald, and the simulation was successful finally. But I found that ewald’s calculation time is basically three times than that of pppm’s calculation.

Have you tried to increase the accuracy of the PPPM forces to at least 1e-5 or smaller values?

No, I haven’t tried this yet. I wonder the purpose of reducing the calculation accuracy.

As stated in the documentation, the accuracy parameter controls the desired relative error in forces. Therefore, a smaller value lead to smaller error in the computed forces, and thus a higher accuracy.
A value of 1e-5 should not increase the computational cost too much, and leads to a much better force and pressure control.

I understand now that you explained it this way. Thank you very much for your help.