ReaxFF and Hybrid Pair style

Dear all

I’m working on a system composed of one layer in the middle (ReaxFF) sandwiched between two layers that are modeled using a MEAM potential. Interactions between the layers are approximated with a Morse potential.

The simulation starts with minimizing the middle layer and this process completes successfully. The next step is adding the two layers on top and bottom using the lattice command with region to create atoms then minimizing the system. The upper and bottom layer are smaller than the middle and hence the volume doesnt change. However, I observed that when the atoms are added, the pressure increases considerably because the upper layer is tending to expand.

The issue is that the energy of the middle layer modeled by ReaxFF decreases by ~2eV. When i broke down the energy components the elp (Coulomb) energy dropped significantly just after adding the atoms. Looking at the potential energy of the atoms, the middle layer changes just when adding the other atoms. The charges of new atoms are zero but it seems like the change is from the QEQ. I tried turning off the interlayer interactions and the issue persists. Strangely when both bottom and top layers are far from the middle layer, the energy does not change. So it seems like the QEQ is somehow calculating energy for these atoms even the charge is zero. I attached a sample script to reproduce this issue using lammps 23Jun2022 compiled with oneapi.

I conducted several tests and I came to the conclusion that, just the addition of these atoms causes this issue. I tried adding using a read data command and yet the issue persists. you can confirm this by looking at the elp energy and by the pe of atoms in ovito. I also performed an NPT run to relax the pressure but the energy did not change. I’m not sure how the coulomb energy can be affected by the presence of other atoms not modeled by the same potential. Also please note that the issue persists regardless of the other layers size. Even one unit cell would change the energy.

I hope I can get some help regarding this. Thanks in advance.

I cant upload to the fourm since i registered recently. the attached files are here: WeTransfer - Send Large Files & Share Photos Online - Up to 2GB Free

The ReaxFF model by its very nature is not a good choice for use in a hybrid pair style.
As you have already identified, the problem originates with the required charge equilibration. You have an inconsistency here, since you need to provide parameters for the atoms that are not contributing to the ReaxFF forces and you are asking the ReaxFF pair style to provide them. Of course it cannot. Fix qeq/reaxff does not contain a check for it. If you would replace it with fix qeq/shielded (which is the same method but implemented separately), it will refuse to run due to the invalid parameters.

Also MEAM is not a good choice for use in a hybrid styles because its model requires to account for contributions from neighboring atoms (to the embedding term and the crystal orientation), but in a hybrid model those atoms are present but not visible to MEAM and thus have no effect even though they should.

My recommendation is to either find a different pair style than reaxff to model that part of the system (possibly with a lower accuracy if considered only by itself but likely with better results in the context of a hybrid pair style), or look for a reaxff parameterization that can handle the entire system.

I tried to make the qeq fix read the parameters from an external file (made from the ff file), and set the fix group to be only the atoms modeled by the ReaxFF. The fix does not seem to read any parameters for the other atoms when the group is set, as expected. If the problem is in qeq not reading the parameters for the extra atoms from the other potential, limiting the group should fix the problem ? But it does not.

I would like to add that after more testing. Strangely, only the presence of the bottom layer does this behavior. The structure with just one layer at the top is normal. The presence of only bottom layer or both layers does this. It seems to me not related to the qeq parameters since it worked with only a top layer present.

Perhaps @RayShan has some idea. If I remember correctly, it was him that implemented and tested hybrid pair style support in ReaxFF. My personal recommendation is to either use ReaxFF for everything or don’t use ReaxFF.

A more detailed description for @RayShan:

I took the minimized structure and performed a ‘run 0’ after adding the upper and lower layer separately.
The first row is for only bottom layer present:

Step Temp Lx Volume Pxx Pyy Pzz PotEng TotEng Fmax v_MXPE v_GPE v_CPE v_OPE v_TiPE v_COUL v_eb v_ea v_elp v_emol v_ev v_epen v_ecoa v_ehb v_et v_eco v_ew v_ep v_eqeq

    0   0              15.182497      19816.649     -113795.31     -108881.78     -21979.374     -102472.47     -102472.47      106.15653     -9.6751715     -6.9678003     -7.2252344     -9.3769367     -11.385205     -2.9994033     -7.3036732      0.589918       0.00026354655  0              0.93954934     4.510672e-27  -1.4205095e-07  0              0              0              1.8988818     -4.8181163      1.818713     

    0   0              15.182497      19816.649     -10201.695     -4343.5205     -0.046630596   -82897.274     -82897.274      76.324419     -7.2770705     -6.9678003     -6.7874644     -9.6356293     -5.9305912     -1.292217      -7.3036732      0.589918       0.00026354655  0              0.93954934     4.510672e-27  -1.4205095e-07  0              0              0              1.8956095     -3.1109301      1.818713     

As it seems, the QEQ is doing its job just fine, but the ep calculates different coulomb energy for the two cases while they should be identical.