I wondered if performing a Reaxff simulation ignoring one of the atoms defined by the force field is possible. For example, I want to carry out this simulation, but I don’t want that the atoms in green react with my other molecules. However, I wish these green atoms experience the other modifications I’m doing in my simulation (like the increase in pressure to compress the system size). Is that possible? And if so, how can I do it?
Please check the neigh_modify exclude keyword (neigh_modify command — LAMMPS documentation), I think exclude type would fit your purpose. I tested on the AuO example and I think neigh_modify exclude type works as expected, though I don’t know much about ReaxFF so someone else may have better comments on this.
If you can represent the green-green atom interaction with a different potential, say lj/cut, then you could experiment with using pair_style hybrid and use either “none” or “zero” for the interactions of the green atom type with everything else. You could use lj/cut for that, too, but then you need to use the lj/cut/soft version, so there are no crashes with short distances between the two subsystems. Due to removing the interactions, the green atoms won’t “see” anything but other “green” atoms and thus can “move through” the other and thus triggering division by zero problems or similar with lj/cut.
neigh_modify exclude would reliably only work for pair-wise additive potentials. While the exclusion would certainly work for the pairwise non-bonded and bonded terms, I am not so certain about the others. Some very careful testing would be required. Since the picture shows a molecular system, I would assume that angle and dihedral interactions might be more important, and that is where inconsistencies are possible.
It would help if you’re a bit more precise about what exactly you’re modelling, and what exactly should the interactions be between green atoms and green atoms, as well as between green atoms and ReaxFF atoms.
Worth noting: ReaxFF usually comes with some kind of charge equilibration, so if the green atoms are charged, you do need to think very carefully about whether their charges are included in the charge equilibration of nearby ReaxFF atoms.
I’m modeling the formation of biochar. I did the same procedure without the green atoms to compress the other molecules until the density was close to the experimental one. But now I want to “constraint” the final molecule by putting something that represents the pore of the biochar so that, for example, the molecules are constrained to be close because of the pore (green atoms). But, not sure if this is okay and if it would work.
Okay. If this were my project, here’s how I’d start to approach it:
Search (and search again) for other simulations that have done similar things in the past (eg “void formation char ReaxFF”). There’s nothing worse than putting in the work only to find that you don’t know whether your results are plausible or not.
Calculate the box size needed for your target density, accounting for void volumes, and stick to using an NVT integrator.
Augment your ReaxFF potential file with parameters for a “void” particle type that results in zero forces on other atoms, and also doesn’t affect their bond coordination numbers. It will help if your voids are large enough (have a large enough interaction radius) that they always fall outside the ReaxFF bond cutoff. (I don’t know ReaxFF well enough to tell you exactly how to do this, but this is a matter of mathematics.)
Set up your simulation with pair_style hybrid/overlay, with ReaxFF applied to all particles, and lj/cut used to apply a WCA potential (Weeks-Chandler-Andersen – a purely repulsive cutoff LJ) between the void type and all particles: