ERROR: Non-numeric pressure - simulation unstable

LAMMPS
1

Hello all,

I’m using LAMMPS (2 Aug 2023 - Update 3) on my Macbook M2. I’m working on creating a crystalline ice at 1 K (H2O) in lammps using the TIP3P potential. I have downloaded and build a large cell of over 1500 water molecules with bond and angles included in the file description. I have modeled my code after the one at 8.4.3. TIP3P water model — LAMMPS documentation. However, I’m using a ready made data file which of a different format from that of a molecule. Initially, I had trouble with the Coulomb long interactions so I expanded my cell boundaries and changed my commands around to bond_style zero and angle_style zero, using the kspace_style ewald 1.0e-4. I also extended by neighbour command to accommodate my system. Here is the snippet of the code:

units real
dimension 3
processors * * *
boundary p p p
atom_style full
box tilt large
atom_modify map array
read_data large_box.data

variable H equal 1
variable O equal 2

neighbor 2 bin
neigh_modify delay 0 every 1 check yes page 100000 one 10000

timestep 1

variable albite_equil equal 5000

variable albite_temp equal 1 #kelvin

Pair Styling

pair_style lj/cut/coul/long 8.0
pair_coeff 1 1 0.0 1.0 # No LJ interaction between H atoms
pair_coeff 2 2 0.1521 3.1507 # LJ parameters for O atoms
pair_coeff 1 2 0.0 1.0 # No LJ interaction between H and O

Bonds

bond_style zero
bond_coeff 1 0.9572 # TIP3P bond length (H-O) and force constant

Angles

angle_style zero
angle_coeff 1 104.52 # TIP3P angle (H-O-H) and force constant

kspace_style ewald 1.0e-4

mass 1 1.008 # Mass of Hydrogen
mass 2 15.999 # Mass of Oxygen

set the charge for H and O manually (from lammps manual on TIP3P) as it isn’t contained in the data file

set type 1 charge 0.417
set type 2 charge -0.834

perform these steps to mitigate errors and get some configuration that’s stable

fix rigid all shake 0.001 10 10000 b 1 a 1
minimize 0.0 0.0 1000 10000

unfix rigid
reset_timestep 0

My error occurs because the dynamics of my system cause the pressure to reach -2.2260919e+16 during the minimisation phase. Also note that I have to assign charges (also taken from the manual for TIP3) to particles separately as they are 0.0 in the original data file. I would really appreciate insight on how to validate where my system is failing. I have read the previous forum answers and implemented the suggestions to my code. Tell me if I missed anything.

Thank you in advance,
Vikentiy

2

Please see Please Read This First: Guidelines and Suggestions for posting LAMMPS questions for information about how to correctly quote text in the forum.

 neigh_modify delay 0 every 1 check yes page 100000 one 10000

You should not need to change “page” and “one” for a system at a normal density.

Such high-pressures are a sign of a bad geometry or a bad topology.

3

I’m not sure what “trouble” you had and how expanding your cell would have helped. But a negative pressure indicates a system that “wants” to compress itself into a smaller volume. Coincidentally, you started your LAMMPS run from an “expanded cell”. You may just be seeing the effects of starting a system with physically unfeasible density that can’t converge in a stable manner back to its optimal density.