Core-shell using neigh_modify exclude

Dear Developers,

I am using Dec 2018 LAMMPS. My system is ionic solid involving one atom type governed by core-shell potentials, which is oxygen atom. As per ‘Adiabatic core/shell model’ document, I have followed all the steps and ran many NVT simulations at different temperatures with various random seed for velocity creation. So, my system is working fairly well for many seeds for simulations upto nano seconds; but for few seeds it is giving an error like “ERROR on proc 2: Bond atoms 1711 1712 missing on proc 2 at step 14372”. As per the information given in document and earlier threads, the error is because of bad dynamics. I would have understood the case of bad dynamics if most of my simulations would have given me this error, but the issue is only 2 out of 5 are showing this error. So I doubt, the potentials are fine but something has to do with processor communications or some other aspect of lammps.

Searching for alternative ways to troubleshoot this issue, I landed on many commands, like comm_modify, newton, neigh_modify. The most helpful command and the one which finished without any error, is “neigh_modify exclude molecule/intra”. Used this command for group of oxygen core-shell, which are described as molecules with bonds in structure file. But the issues with this option is the total energy calculated is different from the simulations which are successful without using ‘neigh_modify exclude’, and I understand, this is fairly consistent with what is mentioned in documents regarding use of ‘neigh_modify exclude’ in presence of long range coulomb forces. What I doubt here is, if the simulation is successful with ‘neigh_modify exclude’ specifically for bonded core-shell within molecule containing only 2 atoms (core and shell), why is it affecting total energy? And if the simulation is working fine with neigh_modify exclude, is there an issue with “special bond” command? I feel the issue is with neighbor list and even with using command “neigh_modify check no every 1 delay 0” the favorable results are not achieved. I just want to know, if there is an issue with numerical collapse or some other ways to achieve, what I want to achieve. Please help, thank you.

Thank you,

Saurabh

PS: Please find the attachment of my input code in.lfp and structure file 545cs.lfpfp

in.lfp (1.63 KB)

545cs.lfpfp (417 KB)