'ERROR on proc 7: Non-numeric atom coords - simulation unstable' error

Hello,

Recently I am trying to reproduce a solid electrolyte interface between lithium slab and ethylene carbonate (EC), but I got the error ‘Non-numeric atom coords - simulation unstable’. I constructed the model following the literature:

  1. lithium slab and EC are simulated for 100 ps under NPT to reach equilibrium state separately using reaxff;
  2. then gathering them together to form the interface model;
  3. the model was conducted energy minimization at first, then run 100 ps under NVT with a timestep of 0.1 fs.

I checked the log file and found two warnings:
WARNING: Van der Waals parameters for element LI indicate inner wall+shielding, but earlier atoms indicate a different van der Waals method. This may cause division-by-zero errors. Keeping van der Waals setting for earlier atoms. (src/REAXFF/reaxff_ffield.cpp:251)
WARNING: Changed valency_val to valency_boc for F (src/REAXFF/reaxff_ffield.cpp:296)
I checked the related posts as well and it seems doesn’t matter as long as the simulation can run. And after step 1, the density of EC I got was the same as the value the literature reported. But unluckily, when I started the interface model simulation, the error was reported after 100,000 timesteps.

I have checked the related posts on the forum and try to figure it out. The first method I tried is restart the simulation, but I got the same error within 100,000 steps. The second method is revising the atomic parameter (in the 4th line, 6th parameter), 0 to 1 following the suggestion from a post. This is a way to fix the error ‘division-by-zero’ according to the post, but it is not working neither.
Then I simulated a small system containing one EC and lithium slab. Its dimension is the same as the one reporting the error. However, this model can run normally within the simulated 8 hours. This may show the force field is working. Additionally, I was thinking if the non-periodic in Z direction or the fixed region affected the simulation, even though these setup followed the literature. So I tested the system under the periodic situation. Then I got the error ‘Non-numeric pressure - simulation unstable’ or the same error as reported under non-periodic. Then I tried to adjust kinds of parameters in the input file (including timestep, neighbor, Tdamp, thermostat, this caused the different error), but none of them worked.

I really tried my best to figure it out but I am failed. Therefore I posted it here to look for some help. I have uploaded the input file mentioned above to help you check it. I would appreciate it so much if you could provide any help!

Best regards,
Linquan

Li-EC-SEI.data (334.1 KB)
SEI_nonperiodic.in (1.6 KB)
SEI_periodic_npt.in (1.5 KB)
ffield.reax.CHOSLiFN (31.1 KB)

Hi @Linquan, a few points:

  • You should provide the version of LAMMPS you’re working with and the version of the Reax you’re working with. There are different versions of the forcefiel and each of them has different known issues.
  • This potential is not part of the LAMMPS potential repository. You should provide the article you got the file from. The authors are also very likely to be the best persones to help you on which parameters and setup to use if needed.
  • You do not state at which point of your simulation the error happens. Neither do you provide a log file.
  • The step 2 of your procedure has to be more detailed. It is very easy for something to go wrong here and lead to unstable simulation, even when minimizing the energy after merging.
  • The input files you provide are bogus.
    • None of them is an NPT input. That is launching them will not lead to a constant pressure simulation. They are NVT as far as I can tell.
    • There are a lot of useless or commented commands making them globally unreadable for understanding what is happening.
    • In the case of the non-periodic z direction, using both boundary p p f and fix 1 all npt temp 300 300 100 iso 1.0 1.0 100.0 (which is actually commented) should result in an error because of the iso keyword. You should also adjust your Tdamp and Pdamp parameter to correspond to your change in timestep. See the NVT command documentation.
    • You have an anisotropic system anyway. Consider barostatting each direction independently. See the same documentation for how to do that.

This gives you some more leads to enhance your scripts, test more things and come back with more details to get meaningful help. Anyway since you are a “beginner at LAMMPS” as you said in your November 2024 post, know that ReaxFF and reactive forcefield in general are hard to use and not very transferable. Be sure to use a forcefield fitted to your needs. You might require to reach out to experts colleagues to guide you on proper usage or comment on results since it is very easy to produce garbage with those models when using them in a wrong way.

Hi Germain,

Thanks for your detailed reply. I would like to clarify the points you mentioned.

  1. The version of LAMMPS is 2 Aug 2023 kokkos. In terms of the version of Reax, It should be a package with C++ code integrated into LAMMPS.
  2. I am sorry that I didn’t upload the literature where I get the force field from. I will attach it at the end of this reply. In fact, I emailed to the corresponding authors of the related works, but I didn’t get any reply so far.
  3. Sorry for the unclear state. The error happens in the step 3, where the energy minimization is finished normally, and the error happens after 100,000 steps under NVT. The log file will be uploaded at the end of this reply as well. To be more specific for the step 2, the box dimensions of lithium slab and EC after step 1 are 34.13334.234143.964 (the real dimension in Z direction is about 69 angstrom, I put the vacuum zone in Z direction. Will it affect the equilibrium of lithium slab?) and 33.30733.16170.563 angstrom, respectively. Then I put the lithium slab into a box (0 to 34.993, 0 to 35.904, 0 to 70) and put the EC into the same box (0 to 34.993, 0 to 35.904, 71.5 to 145) by PACKMOL.
  4. In terms of input files.
    4.1. I am sorry making you confused about the input files. I should have stated this in my post. As I conducted a lot of adjustment of parameters, including changing the NPT into NVT, or NVT into NPT, this makes my input files look messy, and that is why my input files are unreadable. The input files you saw are the latest version of my adjustment. I am really sorry for the inconvenience.
    4.2. This has been explained in 4.1. I should delete commented commands before uploading them.
    4.3. Yes, I got an error as well when applying the NPT for the case of the non-periodic z direction. Additionally, I have tested the systems by adjusting the Tdamp to larger or smaller under NVT, but I even got an error within the shorter steps.
    4.4. Thanks for your suggestion. It is really inspiring. I didn’t realize barostatting each direction independently.

I really appreciate your help. Your words enlightened me a lot. About the force field, to be honest, I have read a lot of related papers and most say they apply the force field C/H/O/N/S/Li/F/N (probably the element order is different) for their systems. So I just use this force field to reproduce one of these results. Even though I do not know if it is suitable for the current system, but it should be. Anyway, I will follow your suggestion to revise my input files and try to get more results.

Sincerely appreciate your help!

Best regards,
Linquan

log-nonperiodic-NVT.out (7.3 KB)
log-periodic-non-numeric coords.out (5.5 KB)
log-periodic-non-numeric pressure.out (4.4 KB)
SEI_nonperiodic.in (999 Bytes)
ff-COHFLiSN-aem.pdf (2.1 MB)
forcefield work.pdf (3.2 MB)
The original forcefield paper but no forcefield files.pdf (1.7 MB)
reproduced work the cited forcefield work is also uploaded.pdf (5.2 MB)
SEI_periodic_npt.in (684 Bytes)

Non-numeric atom coords - simulation unstable error typically means you have bad dynamics, i.e. your simulation blew up and atoms flew apart. Your timestep could be too large, or you could have a problem with your force field definition. You could try visualizing the simulation and see what is going wrong.

Hi Stamoor,

Thanks for your reply. According to your suggestion, can I assume that the issue is caused by the forcefield definition? As I personally think that 0.1 fs is small enough. However, the forcefield parameters are the same as the one in the literature. I will restart the file and visualize it to see what happens.

Best regards,
Linquan

Your personal opinion on your simulation parameter is of little value if it is not backed up by the literature or solid experience. The simulation parameters should also be the same as the one presented for the forcefield for your systems to remain stable. I agree with @stamoor that visualizing the trajectory should almost be a reflex in this kind of situation.

Thanks for your comment. I will output the trajectory of each step to visualize it