# The relaxation process using NPT results in density changes that fall short of requirements

I want to simulate HDPE, it has a density of 0.95g/cm3, I first use Materials Studio to build the initial model and use PCFF force field. After that, I use the NPT ensemble of LAMMPS to relax the initial structure, cooling down from 600K to 300K. The results show that his density will suddenly decrease, then increase with the decrease of temperature, and then stabilize at about 0.77g/cm3.
This density is too low for my requirements, has anyone encountered the same problem?

There is likely something wrong with your input or data file.

Hello，My settings are as follows, can you help me see where the problem is?
606060.data (7.4 MB)

# ***** Define a few parameters *******

variable T equal 600
variable DT equal 0.4 # 0.4fs

# ******* Potential Function **********

velocity all create \${T} 12345 mom yes rot yes dist gaussian

#能量最小化
timestep \${DT}
min_style cg
minimize 1.0e-10 1.0e-10 10000 100000

# equilibration and thermalization

fix 1 all nvt temp \$T \$T 5
thermo 4000
thermo_style custom step temp
run 1000000
unfix 1
write_restart restart.nvt.first

fix 1 all npt temp \$T 450 5 iso 1 1 50
#dump 2 all atom 4000 dump.atom
thermo 4000
thermo_style custom step temp vol density
run 3000000
unfix 1
#undump 2

This is incomplete.

My complete relaxation process is as follows, thank you！

units real
dimension 3
boundary p p p
atom_style full
bond_style class2
angle_style class2
dihedral_style class2
improper_style class2
kspace_style ewald 1.0e-4
pair_style lj/class2/coul/long 10 10

neighbor 0.3 bin
neigh_modify delay 0 every 1

# ***** Define a few parameters *******

variable T equal 600
variable DT equal 0.4 # 0.4fs

# ******* Potential Function **********

velocity all create \${T} 12345 mom yes rot yes dist gaussian

#能量最小化
timestep \${DT}
min_style cg
minimize 1.0e-10 1.0e-10 10000 100000

# equilibration and thermalization

fix 1 all nvt temp \$T \$T 5
thermo 4000
thermo_style custom step temp
run 1000000
unfix 1
write_restart restart.nvt.first

fix 1 all npt temp \$T 450 5 iso 1 1 50
#dump 2 all atom 4000 dump.atom
thermo 4000
thermo_style custom step temp vol density
run 3000000
unfix 1
#undump 2

fix 1 all npt temp 450 450 0.005 iso 1 1 0.05
thermo 4000
thermo_style custom step temp vol density
run 1000000
unfix 1

fix 1 all npt temp 450 300 0.005 iso 1 1 0.05
thermo 4000
thermo_style custom step temp vol density
run 3000000
unfix 1
write_restart restart.npt.1

fix 1 all npt temp 300 300 0.005 iso 1 1 0.05
thermo 4000
thermo_style custom step temp vol density
run 1000000
unfix 1
write_restart restart.npt.2

What is the point of heating this system with variable cell?
You are just needlessly expanding it without the desire to run a simulation at that state.
Since it is much easier to expand a system than to condense it again, and particularly so with polymers (note the keyword “jamming”) and glasses, you are just needlessly complicating your job.
Also, your choice of time constants for thermostat and barostat are for the most part far too small.

Sorry, let me revise my thermostat and barostat. So do you suggest that I relax only in the context of the NVT ensemble?

It depends on what your ultimate goal is for the kind of equilibrated system that you want to get.
Please keep in mind that during equilibration it is not necessary to do a physically meaningful “production” simulation for as long as it gets you efficiently into or close enough to the equilibrium state.

My guess was, that you want to melt/relax your possibly too “regular” initial geometry into a more irregular one by running at higher temperature and the cool it down to the target temperature.
Because of the difficulties of shrinking a previously expanded system, particularly with entangled polymers and glasses, it is preferable to adjust the system to the desired target density and then try to manipulate only the kinetic energy. That is usually also more effectively done with a dissipative thermostat (like fix langevin, fix temp/csld, fix gld, fix gle etc.) rather than a nose-hoover thermostat.
only once the system is sufficiently relaxed and termalized, I would switch to the nose-hoover for final equilibration.

With respect to achieving the desired density and relaxing to the desired pressure state, you need to monitor the pressure and potential energy during equilibration and I would switch to a barostat only after the system has sufficiently equilibrated and the pressure and energy stabilized. You may not get a perfect match because force fields are always approximations.

Due to the difficulties of relaxing jammed/entangled systems, you may need to look into some monte carlo based methods to accelerate equilibration.

Thank you so much for the guide, I’ll explore it again