Non-numeric pressure and unstable minimization in TIP4P water and OPLS-UA ethanol mixture

water_ethanol_box_3000_0.2.txt (1.5 KB)
water_ethanol_box_3000_0.15.txt (1.5 KB)
water_ethanol_box_3000_0.2.log (11.8 KB)
water_ethanol_box_3000_0.15.log (17.7 KB)

Dear LAMMPS users,

I am attempting to prepare a simulation box containing a mixture of ethanol and water. Ethanol is modeled using the OPLS-UA force field, and water is modeled using the TIP4P model.

As defined in the input script water_ethanol_box_3000_0.2.txt, I generated a system containing 600 ethanol molecules and 2400 water molecules. After energy minimization, I attempted a short npt simulation. However, the run failed with the following error:

ERROR: Non-numeric pressure - simulation unstable (src/fix_nh.cpp)

As defined in the input script water_ethanol_box_3000_0.2.log , this instability appears to originate from the high energy and NaN pressure reported at the end of the minimization step.

For comparison, I changed the composition to 450 ethanol and 2550 water molecules (keeping the total number constant). The simulation with this new configuration (water_ethanol_box_3000_0.15.txt) ran without errors, but the total energy after minimization remained high, and the temperature spiked immediately after entering the npt phase. If I run the npt simulation for a sufficiently long time, the temperature stabilizes at around 315 K.

Here are additional details:

  • Atomic overlaps were checked using delete_atoms overlap 0.1 all all, and no overlaps were found.
  • The minimization settings used were the LAMMPS default (minimize 1.0e-4 1.0e-6 100 1000).
  • The simulation box size was set to match the target density (0.97 g/cm³).

I would like to ask:

  • Is the high energy state after minimization expected behavior in these mixed systems?
  • Is it acceptable that temperature spikes and then settles during npt?
  • How can I improve the minimization so that the system becomes more stable prior to npt?

Any advice or suggestions would be greatly appreciated.

Best regards,
asas

Minimization will only find a (local) minimum that is close to the initial state. That will be a global minimum only in the most trivial of cases. For random initial placement this will definitely not be the case and there are no minimization settings that can change that.

Kinetic energy may raise, but you should take measures that it doesn’t rise too much:

  • use a short(er) timestep at the beginning
  • don’t run fix npt but fix nvt or fix langevin + fix nve until temperature and pressure have stabilized somewhat and potential energy has gone down significantly
  • initially run a shorter MD runs interleaved with minimizations to remove potential energy and get you closer to equilibrium

You cannot, see my answers to the previous questions.