Extreme pressure drop

ms_research · February 4, 2026, 3:22am

Hi everyone.

I was working on a LAMMPS project where the goal is to simulate the diffusion coefficients of oxygen and carbon dioxide in polyurethane membranes.

While running my project, the density dropped from a good 1.22g/cm^3 to a 0.4g/cm^3. Alongside it, there was an extreme pressure drop and overall, it is very concerning.

Below is the LAMMPS input script. Hopefully someone can help.

# ============================================================
# INITIALIZATION
# ============================================================
units           real
atom_style      full
boundary        p p p

# Initial force field settings
pair_style      lj/cut/coul/long 12.0
bond_style      harmonic
angle_style     harmonic
dihedral_style  opls
kspace_style    pppm 1.0e-4
special_bonds   lj/coul 0.0 0.0 0.5
comm_modify     cutoff 6.0

# ============================================================
# READ DATA & FORCE FIELDS
# ============================================================
read_data       polyurethane_realistic.data

# Pair coefficients
include         file.txt

# Bond coefficients
bond_coeff  1   320.0000 1.4100  # C_alk-O_eth
bond_coeff  2   340.0000 1.0900  # C_alk-H_alk
bond_coeff  3   268.0000 1.5290  # C_alk-C_alk
bond_coeff  4   214.0000 1.3270  # C_carb-O_eth
bond_coeff  5   570.0000 1.2290  # C_carb-O_carb
bond_coeff  6   490.0000 1.3350  # C_carb-N_amide
bond_coeff  7   300.0000 1.5000  # H_amide-N_amide
bond_coeff  8   300.0000 1.5000  # C_arom-N_amide
bond_coeff  9   469.0000 1.4000  # C_arom-C_arom
bond_coeff  10  367.0000 1.0800  # C_arom-H_arom
bond_coeff  11  300.0000 1.5000  # C_alk-C_arom
bond_coeff  12  1000.0000 1.2100  # O_gas-O_gas (stronger for O2)
bond_coeff  13  1000.0000 1.1600  # C_gas-O_gas (stronger for CO2)

# Angle coefficients
angle_coeff 1   50.0000 109.5000  # H_alk-C_alk-O_eth
angle_coeff 2   50.0000 109.5000  # C_alk-C_alk-O_eth
angle_coeff 3   50.0000 109.5000  # H_alk-C_alk-H_alk
angle_coeff 4   50.0000 109.5000  # C_alk-C_alk-H_alk
angle_coeff 5   50.0000 109.5000  # C_alk-C_alk-C_alk
angle_coeff 6   60.0000 109.5000  # C_alk-O_eth-C_alk
angle_coeff 7   60.0000 109.5000  # C_alk-O_eth-C_carb
angle_coeff 8   80.0000 120.0000  # O_carb-C_carb-O_eth
angle_coeff 9   80.0000 120.0000  # N_amide-C_carb-O_eth
angle_coeff 10  80.0000 120.5000  # N_amide-C_carb-O_carb
angle_coeff 11  50.0000 121.9000  # C_carb-N_amide-H_amide
angle_coeff 12  50.0000 121.9000  # C_arom-N_amide-C_carb
angle_coeff 13  50.0000 121.9000  # C_arom-N_amide-H_amide
angle_coeff 14  70.0000 120.0000  # C_arom-C_arom-N_amide
angle_coeff 15  70.0000 120.0000  # C_arom-C_arom-H_arom
angle_coeff 16  70.0000 120.0000  # C_arom-C_arom-C_arom
angle_coeff 17  70.0000 120.0000  # C_alk-C_arom-C_arom
angle_coeff 18  50.0000 109.5000  # C_arom-C_alk-H_alk
angle_coeff 19  50.0000 109.5000  # C_arom-C_alk-C_arom

# Dihedral coefficients
dihedral_coeff  * 1.0 1.0 1.0 0.0

# ============================================================
# GROUPS
# ============================================================
group           oxygen type 10
group           co2    type 11
group           polymer type 1 2 3 4 5 6 7 8 9

# ============================================================
# SETTINGS
# ============================================================
neighbor        2.0 bin
neigh_modify    delay 0 every 1 check yes
timestep        0.1

# ============================================================
# STAGE 1: OVERLAP REMOVAL (Soft Potential Push-off)
# ============================================================
print "========================================="
print "STAGE 1: Removing overlaps with soft potential"
print "========================================="

# Temporarily disable long-range electrostatics and switch to soft potential
kspace_style    none
pair_style      soft 2.5
pair_coeff      * * 0.0

variable        pre equal ramp(0,100)
fix             SOFT all adapt 1 pair soft a * * v_pre
fix             LIMIT all nve/limit 0.05

thermo          500
thermo_style    custom step temp density press pe ke
run             10000

unfix           SOFT
unfix           LIMIT

# ============================================================
# CRITICAL: RESTORE ORIGINAL FORCE FIELD
# ============================================================
print "========================================="
print "RESTORING FULL FORCE FIELD"
print "========================================="

pair_style      lj/cut/coul/long 12.0
kspace_style    pppm 1.0e-4
include         file.txt  # Restore pair coefficients

# ============================================================
# STAGE 2: ENERGY MINIMIZATION
# ============================================================
print "========================================="
print "STAGE 2: Energy minimization"
print "========================================="

minimize        1.0e-4 1.0e-6 1000 10000

# ============================================================
# STAGE 3: NVT EQUILIBRATION (Temperature)
# ============================================================
print "========================================="
print "STAGE 3: NVT equilibration at 310.15 K"
print "========================================="

reset_timestep  0
timestep        0.5

fix             NVT all nvt temp 310.15 310.15 50.0

thermo          1000
thermo_style    custom step temp density press pe ke etotal vol

run             100000

unfix           NVT

# ============================================================
# STAGE 4: NPT EQUILIBRATION (Pressure & Density)
# ============================================================
print "========================================="
print "STAGE 4: NPT equilibration at 1 atm"
print "========================================="

reset_timestep  0
timestep        1.0

# Use moderate damping parameters for stable equilibration
fix             NPT all npt temp 310.15 310.15 100.0 iso 1.0 1.0 500.0 drag 1.0

thermo          2000
thermo_style    custom step temp density press pe ke etotal vol

run             500000

# Continue NPT with finer control
unfix           NPT
fix             NPT2 all npt temp 310.15 310.15 100.0 iso 1.0 1.0 1000.0 drag 2.0

run             500000

unfix           NPT2

# ============================================================
# STAGE 5: PRODUCTION RUN (Diffusion Calculation)
# ============================================================
print "========================================="
print "STAGE 5: Production run - Calculating diffusion"
print "========================================="

reset_timestep  0

# Compute MSD with center-of-mass correction
compute         mO2 oxygen msd com yes
compute         mCO2 co2 msd com yes

# Calculate diffusion coefficients on-the-fly
# D = MSD / (6*t), units: A^2/fs → multiply by 1e-1 for cm^2/s
variable        diff_O2  equal c_mO2[4]/6/(step*dt+1.0e-6)
variable        diff_CO2 equal c_mCO2[4]/6/(step*dt+1.0e-6)

# Production NVT ensemble
fix             PROD all nvt temp 310.15 310.15 100.0

# Output MSD data
fix             WRITE all ave/time 1 1 5000 c_mO2[4] c_mCO2[4] v_diff_O2 v_diff_CO2 &
                file msd_results.txt

# Trajectory output
dump            1 all custom 10000 trajectory.lammpstrj id type x y z

# Thermodynamic output
thermo          5000
thermo_style    custom step temp density press c_mO2[4] c_mCO2[4] v_diff_O2 v_diff_CO2

run             5000000

print "========================================="
print "SIMULATION COMPLETE"
print "========================================="

Thank you very much in advance.

Micah

Germain · February 4, 2026, 2:01pm

Hi @ms_research,

You’ll find general advice on what to investigate and how in the “read this first” post.

You should also indicate your LAMMPS version and investigate your system through visual analysis. At the moment without a more restrained problem and all the input and output files (logs, plot etc.), it is nearly impossible to provide meaningful help on what could have gone wrong in such a complicated system.

Have you checked that your oxygen, carbon dioxide and polyurethane membranes are independently stable with the same simulation procedure? Are all your atoms correctly connected to one another? Have you plotted all the different energy components before your box explosion?

This is the best I can personally come up with at the moment with the few elements you have here. Note that the best you cand do is to investigate locally with someone from your group that knows you procedure or how to set up such simulations and will be better at telling you where to look.

akohlmey · February 5, 2026, 11:17am

There is very little that can be learned from an input file alone. Below are some comments.
But without being able to (easily!!) reproduce what you see or having the output and logs to look at, there is not much substantial that can be said. As @Germain already mentioned, the most likely sources of problems are the force field and the initial geometry.

This is pointless. LAMMPS will not permit the communication cutoff to be smaller than the largest pairwise cutoff plus neighbor list skin. This is required for the pair-wise interactions to be computed correctly.

This looks very suspicious. Where does this come from?

Is this long enough to switch to NPT? This seems rather short. Polymers tend to be very slow to equilibrate and often need help from some Monte Carlo steps to avoid so-called “jamming” issues. What is the pressure like during NVT?

Using the COM correction only makes sense for modeling flows. Your system should not have a COM drift. If it has, then you have a bigger problem and it cannot be solved by just ignoring the drift: you must find the origin and remove it.

Have you checked if your system has a stable COM?