An error was reported in the gas hydrate simulation

LAMMPS LAMMPS General Discussion
,
1

Good evening everyone,

I’m currently using LAMMPS to simulate a polycrystalline CO₂ hydrate system. The workflow involves:

  1. Initial relaxation using NVE with a Langevin thermostat, which runs without issues.
  2. Then switching to NPT ensemble for further equilibration.

However, upon starting the NPT run, I encounter the following error:

Setting up Verlet run ...
Unit style    : real
Current step  : 200000
Time step     : 2
ERROR: Non-numeric pressure - simulation unstable (../fix_nh.cpp:1049)
Last command: run     1000000
WARNING: Shake determinant < 0.0 (../fix_shake.cpp:1825)

After reducing the timestep to 0.1 fs, a new error occurs:

ERROR on proc 20: TIP4P hydrogen is missing (../pair_lj_cut_tip4p_long.cpp:224)
  • The water model is TIP4P, and I’ve reviewed the molecular definitions and force field settings (including pair_style lj/cut/tip4p/long) for errors.
  • Atom ordering and group assignments appear correct.
  • fix shake is applied to water, but CO₂ molecules are treated as rigid or flexible in different trials (results similar).

Thank you very much for your valuable time and guidance!

1.data (2.0 MB)
in.poly (3.6 KB)
in.nvt (3.7 KB)
log.lammps (17.7 KB)

2

You need to re-read the forum guidelines and follow them, since your post is near unreadable.
It is missing crucial information, and I have told you before that partial input are next to useless.

Despite all these shortcomings there are three issues to note:

  • you say that you are simulating CO2 hydrate, but your comments and commands suggest that there are also Ca and He atoms.
  • you are using fix shake to have rigid water, but I see no equivalent for the CO2 molecules
  • you should not be using pair style hybrid, but use lj/cut/tip4p/long for all interactions. With your setup all atoms modeled with lj/cut/coul/long “see” the negative charge of the water molecule at the O atom and not the point M as they should.
3

Thanks for the feedback, and sorry for the confusion in my previous post. I’d like to clarify a few points:

  1. The system does include some methane. I used “Ca” and “He” as placeholder atom types by mistake — I’ll correct the labeling.
  2. I originally applied fix rigid to both CO₂ and CH₄, but fix npt gave errors when combined with rigid fixes, so I temporarily removed them. Any advice on how to stabilize fix rigid with npt would be appreciated.
  3. During minimize, the system energy increased sharply. I expected minimization to reduce energy even if atoms overlap. Could this be due to bad initial geometry or force field incompatibility?
  4. Thanks for pointing out lj/cut/tip4p/long. I didn’t realize it was necessary for the M site in TIP4P. I’ll revise my setup to avoid using pair_style hybrid.
  5. I got WARNING: Inconsistent image flags. I exported the data file from OVITO after wrapping. Does that break periodicity or cause image flag issues?
4

When looking at your log file, there are some indications of problems and warnings that you should not ignore.

LAMMPS (28 Mar 2023)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
  using 1 OpenMP thread(s) per MPI task
Loaded 0 plugins from D:\ LAMMPS-64bit-latest-MPI\plugins
Loaded 1 plugins from D:\LAMMPS 64-bit 28Mar2023-MPI\plugins

This indicates that a) you are using a more than 2 year old LAMMPS version and b) that you may have a mixed installation of two LAMMPS packages on top of each other. That is bad. Nobody will be looking into issues unless they are confirmed with the latest release (12 June 2025 at this time).

   Step          Temp          PotEng         KinEng         E_vdwl         E_coul         E_pair         E_bond        E_angle        E_dihed        E_impro         E_mol         Enthalpy        TotEng          Pxx            Pyy            Pzz           Volume           Lx             Ly             Lz          Density    
         0   0              81532.749      0              190783.05      640692.13      75590.893      2584.7338      22.361195      0              0              2607.095       881247.44      81532.749      137565.42      152644.81      150530.26      373248         72             72             72             0.90799595   
        10   0             -96027.892      0              25487.134      632214.32     -98365.714      360.26102      80.923894      0              0              441.18492     -4919.7028     -96027.892      16358.158      16703.113      17150.472      373248         72             72             72             0.90799595   
        20   0             -4.6344116e+10  0              25328.859     -4.6343387e+10 -4.6344118e+10  338.6464       81.549301      0              0              420.1957      -6.1792054e+10 -4.6344116e+10 -4.3850099e+09 -7.2964088e+08 -3.3990506e+09  373248         72             72             72             0.90799595   
        30   0             -3.2503001e+16  0              25328.859     -3.2503001e+16 -3.2503001e+16  338.6464       81.549301      0              0              420.1957      -3.2503001e+16 -3.2503001e+16  15681          15845.417      16503.91       373248         72             72             72             0.90799595   
        31   0             -3.2503001e+16  0              25328.859     -3.2503001e+16 -3.2503001e+16  338.6464       81.549301      0              0              420.1957      -4.45126e+16   -3.2503001e+16  15681          15845.417     -6.6187562e+15  373248         72             72             72             0.90799595

Your Coulomb energy becomes very large. That is usually an indication of a problem with the force field (LJ?) parameters that allow oppositely charged atoms come too close to each other. That in turn will lead to large forces and a lot of problems when starting MD.

WARNING: Should not use fix nve/limit with fix shake or fix rattle (src/fix_nve_limit.cpp:78)
WARNING: Shake determinant < 0.0 (src/RIGID/fix_shake.cpp:1796)

The first warning should not be ignored and the second is likely the consequence, which means that the SHAKE equations are degenerate. This is likely a consequence of the previous issue.

BTW: there also is a warning about mixing tip4p/long with coul/long:

WARNING: Using a tip4p sub-style with other sub-styles that include coulomb interactions can result in inconsistent calculation of the coulomb interactions (src/pair_hybrid.cpp:375)

In general, if you ignore warnings, you better know very well, what you are doing.

Impossible to say without having a closer look (which makes no sense with all the other issues that need resolving). Some general explanation is here: 11.2. Errors and warnings details — LAMMPS documentation

You should probably read the entire page top to bottom to be better prepared for errors and warnings that LAMMPS will issue.

5

Thank you very much for your kind guidance. I have identified and resolved the previous issues, and I have also downloaded and compiled the latest version of LAMMPS. The image flags problem has now been fixed.
However, I am still encountering an abnormal energy issue at the beginning of the simulation. At this point, I’m a bit unsure how to proceed with further debugging. I would greatly appreciate it if you could offer some advice or suggestions.
in.test (2.3 KB)
1.data (2.0 MB)
log.lammps (5.5 KB)

6

I already mentioned this issue and also speculated about its origin. That is all I can do. This is not my research and thus not my job to figure out all the details. All we can do here is point you into a direction. The rest is up to you.

7

I really appreciate your help. I’ll keep working on it. Thanks again!