Crosslinking MPD-TMA using bond/react

Hello! I am trying to crosslink MPD and TMA using “fix bond/react”. In my pre-reacted and post-reacted files I have 37 and 34 atoms. During the crosslinking process, unreacted MPD amide groups within 4.5 ̊A of an unreacted TMA carboxyl group (amide nitrogen to carboxyl carbon distance) “react” to form a crosslink. Moreover, three atoms of HA (MPD amide), OH and HO (TMA) are deleted during crosslinking. When I run my input file, I get an error saying that “ERROR on proc 0: Bond/react: Number of equivalences in map file must equal number of atoms in reaction templates (…/fix_bond_react.cpp:3779)”. I would highly appreciate any suggestions.

Best,
Majid
mpd_tma_map.txt (2.0 KB)
input_mpd_tma_crossLinking.in (2.3 KB)
pre_reacted_mpd_tma.data (5.5 KB)
post_reacted_mpd_tma.data (5.0 KB)
MPD_TMA_crosslinked_FULL.data (140.4 KB)

For questions on fix bond/react you just mention @jrgissing in your post so he notices sooner.

Deleted atoms should be included in your post-reaction template, in which they cannot be connected to non-deleted atoms. So, if you delete three atoms and create three atoms, your post-reaction template should have three more atoms than your pre-reaction template. The equivalences section should have exactly as many atoms as the pre-reaction template.

Also, there are quite a few typos or instances of incorrect formatting in your map file. I would suggest looking at example map files here: examples/PACKAGES/reaction

Hi Dr. Gissinger @jrgissing

I am also trying to simulate the MPD–TMC crosslinking reaction using the fix bond/react protocol, but I am encountering “bond atoms missing” errors. I observe the errors in two different ways:

  1. When reaction site stabilization is enabled (xmax = 0.03) and the number of stabilization steps is less than the run length for fix bond/react, the simulation crashes immediately after the stabilization steps are exhausted.

  2. When the number of stabilization steps is greater than the run length for fix bond/react, the simulation completes successfully. However, the resulting system cannot be used for subsequent NVT or NPT simulations, as the same error reappears.

I believe the issue may be related to the reaction templates and atom mapping, but I am unsure what I am missing. I have attached my input files along with a brief document describing my atom-type labeling.

I would greatly appreciate any guidance on this issue.

map.txt (287 Bytes)

ready.data (1.3 MB)

react.lammps (1.4 KB)

pre.txt (2.5 KB)

post.txt (2.5 KB)

molecule mapping and typelabels notes.pdf (167.7 KB)

TMC.lt (8.9 KB)

NVT.lammps (572 Bytes)

MPD.lt (7.9 KB)

In the future, please create a new post for a new question/topic. A couple notes:

  1. The number of stabilization steps is always meant to be less than the number of run steps, because it a local, temporary stabilization.
  2. If your simulation is unstable, it suggests that your stabilization parameters are inappropriate for your system. For example, note that you have 0.003 instead of 0.03 for xmax in your input script, which could prevent atoms from moving enough to stabilize with the standard number of stabilization steps (typically 60-100).
  3. Interrupted ‘fix bond/react’ runs can be continued via restart files. If using data files, they may need to be minimized first, if a reaction was in the middle of stabilization when it was interrupted.

Thank you Dr. Gissinger @jrgissing

I changed xmax to 0.03 and number of stabilization steps to 100.

Initially, I was running into errors whilst the Fix bond/react is going on.

I decided to change the time step from 1fs to 0.01fs and the system is pretty stable.

It is even able to perform an NVT or NPT for more than 50,000 runs after fix bond/react is off.

I am not facing any critical problem right now, but I will let you know if I encounter one.

Thank you very much for your time and response.

Please note that 0.01 fs is a very small timestep for typical fixed-bond force fields. You are likely masking the issue rather than addressing it.

Thank you, I changed the timestep back to 1fs and kept number of stabilization steps to 100.

I am producing Bond atoms missing errors in these occassions

  1. when xmax = 0.03, bond between atom 3569 and 2216 which corresponds to atom types CST and NM (where the new bond is being formed between reactants) gets missing.
  2. When xmax = 0.02 (same order of magnitude as eg 1), bond between atom 4594 and 2231 which also corresponds to atom types CST and NM (where the new bond is being formed between reactants) gets missing.
  3. When xmax = 0.3 or 0.5 (10 order of magnitude bigger than eg 1) the missing bonds occur between atom pairs (2038, 2047) and (342, 351), respectively. These bonds correspond to atom types CM and HBM, which are not involved in the reaction.

A drawing of my atom type label map is attached below:

molecule mapping and typelabels notes.pdf (167.7 KB)

One immediate potential issue I notice is that your pairwise coefficients for type HNM are all zero. Were these settings on purpose?

No sir, please those pair parameters were obtained from OPLSAA force field, it was not on purpose.

Dr. Gissinger, I compared my parameters with those reported in one of your publications (https://doi.org/10.1021/acsami.4c16229) and found that the Lennard-Jones parameter for HNM (739 in the paper) is also set to zero in your work. However, I noticed a few additional differences that I am currently working to resolve.

  1. A new atom type was introduced for Cl after it became an ion.
  2. The partial charges in my system are slightly different from those reported in the paper
  3. It appears that the molecular templates used in the paper consist only of atoms labelled in the molecule structure, making them considerably smaller than the template I am using.

I will attempt to rebuild the system to match the paper exactly and will follow up if I am unable to resolve the issue. Thank you very much for your time and assistance.

Please also see a discussion of using REACTER with OPLS-AA in this recent paper: https://doi.org/10.1080/08927022.2025.2611025

Thank you Dr. Gissinger, I will read it.