Bond atoms missing during bond/react fix

I would like to cross-link an amorphous polyethylene (PE) system with bond/react fix. My PE system contains 100 carbon atoms, which has been minimized and relaxed sufficiently before cross-linking reaction.
I defined dihedrals (cvff forcefield) and need edge atoms to represent the rest of the long PE chain, so I specified 13 carbon atoms for each chain to updata the topology of atoms within three bonds of the reacting atom (Is this necessary?)
The most important problem is that when I run the simulation, it prompted a warning of “system is not charge neutral” first, and then stopped with an error of “bond atoms 992 994 missing”. I thought that the waring of non-equilibrium of charge is due to the deletion of the by-product, i.e., the two hydrogen atoms. Does it lead to the missing atoms error? If not, why does the missing atom error occurs and how can I address this problem. Could you please give me some advice.
The attachments are the diagram figure of my reaction template and the scripts. Thank you very much! (781.7 KB)
2-pre.temp (6.3 KB)
2-post.temp (6.5 KB) (672 Bytes) (784 Bytes)
log.lammps (4.2 KB)

@jrgissing is the author of fix bond/react. He can hopefully provide some insight.

The non-charge neutral warning would make sense, if the two eliminated atoms (green) have a partial charge. Then their removal would need to require a change of the partial charges of the two atoms forming a link where those atoms were (yellow) and they charge has to be adjusted by the negative or the green atom partial charge. Also their atom type would likely have to change.

There is not much help that I can provide here. The data file reader in the VMD TopoTools package currently does not support information stored by fixes. You would have to make a copy of the data file and remove the offending line(s).

… or you would have to use the nofix option to the write_data command to prevent writing the information (which will prohibit restarting the fix, but that may even be desired in your case).

Thank you all the same!
As for the advice you gave on the non-charge neutral warning, i.e., “their charge has to be adjusted by the negative or the green atom partial charge”, I am still confused about how to implement it. Directly specify the charge distribution of all the atoms after reaction in the post-reaction template? Or are there any specific commands for me to accomplish that? Thank you very much!

This should be part of the processing that fix bond/react needs to do.
I am not familiar with its details. You can study the documentation or wait for @jrgissing to notice that I mentioned him and that he has time to read about your post and provide additional insight.

by default, fix bond/react updates all atomic charges to match those in the post-reaction template. so, yes, you can ‘directly specify the charge distribution of all the atoms after reaction in the post-reaction template’

1 Like

Thank you for your kind response!
I’ve tried to specify the charge distribution in the post-reaction template to solve the non-charge neutral problem. However, the problem is that the charge on the deleted by-product atoms is changing, but it seems that I can not adjust the charge distribution dynamically and accordingly in the post-reaction template over time. I also tried to set the charge on all the atoms to zero. But I can only specify the charge on the atoms at the reaction site, and the rest of the molecule represented by the edge atoms are not affected. Thus, the whole system is still non-charge neutral, and the warning remains.

I do have an experimental version of fix bond/react that updates charges according to the post-reaction template, then rescales the total charge of the reaction site to equal that of the site before the reaction. Would that type of rescaling make sense for your model?

otherwise, this is a nontrivial problem. A potentially elegant, but more advanced, solution is to use charge equilibration methods (such as fix qeq), but I have not heard of anyone trying that yet with fix bond/react

That’s exactly what I need! Could you please send it to me? Or when will it be released?

Besides, I found that the data file of the system after reaction (exported by write_data command) contained a strange section named “bond_react_props_internal # i_limit_tags i_react_tags” at the end of the data file, as shown in the attachment. What’s its meaning and functionality? Can I just delete it because the LAMMPS seems to be unable to read the data file containing the section? I did not find the description about the “bond_react_props_internal” section in the manual. (919.3 KB)

Please find a version of fix bond/react that rescales charges, attached. Each charge in the reaction site is incremented by the same amount, such that the total charge of the reaction site is equal to that of the site before the reaction. I didn’t add a keyword, so the charge rescaling will always happen. This code has been tested on the latest version of lammps (23 Jun 2022)

EDIT: the wrong files (with a different type of rescaling) were attached at first. they have been reuploaded

fix_bond_react.cpp (162.8 KB)
fix_bond_react.h (9.2 KB)

this section of the data file records custom groups that fix bond/react creates internally. the groups are used to relax a given reaction site, and subsequently release it back to the system-wide thermostat (this is described in the docs, just not using the variable names printed in the data file). LAMMPS should be able to read those sections, but they can be deleted or never even printed using the ‘nofix’ option of write_data. just note that if reactions are being relaxed when the simulation run ends, you may need to use the restart file to get a stable restart of the simulation

Thank you very much! I do appreciate your kind help!