Very high pressure in NVT Polycrystal system simulation

LAMMPS
,
1

Hello,

I am running a LAMMPS MD simulation for a grain boundary polycrystalline system (~62k atoms) with ML potential.

LAMMPS version is the lammps deepmd module

Workflow:

  • Gentle shakeout using nve/limit + langevin

  • Energy minimization (CG)

  • NVT heating (300 → 1000 K) to hold at 1000K to study gas atoms diffusion and segregation at grain boundaries.

However, I consistently observe very high pressures (~50,000–60,000 bar) after minimization and during NVT runs, while the simulation box remains fixed as set by NVT.

I also tried using the following workflow: Step 0: Energy minimization, Step 1: NPT equilibration (20 ps), Step 2: NVT production for segregation/diffusion, but in this case, when I reach step 1:NPT the grain boundaries start to dissolve, and the system starts to become amorphous rather than retaining the polycrystalline nature.

From my understanding, this may be due to internal stress since the box cannot relax under NVT conditions. So, I would like to request for guidance regarding the workflow that should be adopted in this case. Attached is the lammps script used.

Thank you for your help.
clear
units metal
dimension 3
boundary p p p
atom_style atomic

log log.gb_he_300_to_1000K_nvt

variable Tstart equal 300.0
variable Tprod equal 1000.0

read_data vac_HeGB_140atoms.lmp

mass 1 238.02891
mass 2 91.224
mass 3 4.002602

plugin load libdeepmd_lmp.so
pair_style deepmd home/frozen_model.pb
pair_coeff * *

neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes

timestep 0.001

group UZr type 1 2
group gas type 3

thermo 100
thermo_style custom step temp pe ke etotal press pxx pyy pzz vol lx ly lz
thermo_modify flush yes

restart 50000 restart.gb_he_nvt.*.bin

============================================================

STEP 1: Very gentle shake-out

============================================================

reset_timestep 0

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

thermo 10
thermo_style custom step temp pe ke etotal press pxx pyy pzz vol lx ly lz
thermo_modify flush yes

fix relax1 all nve/limit 0.005
fix damp all langevin {Tstart} {Tstart} 0.1 12345 zero yes

run 2000

unfix damp
unfix relax1

write_data after_gentle_shakeout.lmp

============================================================

STEP 2: Atomic minimization only

============================================================

reset_timestep 0

thermo 100
thermo_style custom step pe etotal press pxx pyy pzz vol lx ly lz
thermo_modify flush yes

min_style cg
minimize 1.0e-8 1.0e-8 5000 50000

write_data after_atomic_minimization_no_boxrelax.lmp
write_restart after_atomic_minimization_no_boxrelax.restart

============================================================

STEP 3: NVT heating from 300 K to 1000 K

============================================================

reset_timestep 0
timestep 0.001

velocity all create ${Tstart} 24680 mom yes rot yes dist gaussian

dump HEAT all custom 100 dump.nvt_heat_300_to_1000K.lammpstrj id type xu yu zu vx vy vz
dump_modify HEAT sort id

thermo 10
thermo_style custom step temp pe ke etotal press pxx pyy pzz vol lx ly lz
thermo_modify flush yes

fix heat all nvt temp {Tstart} {Tprod} 0.1

run 80000 # 80 ps

unfix heat
undump HEAT

write_data after_nvt_heating_1000K.lmp
write_restart after_nvt_heating_1000K.restart

============================================================

STEP 4: NVT equilibration hold at 1000 K

============================================================

reset_timestep 0
timestep 0.001

dump HOLD all custom 100 dump.nvt_hold_1000K.lammpstrj id type xu yu zu vx vy vz
dump_modify HOLD sort id

thermo 10
thermo_style custom step temp pe ke etotal press pxx pyy pzz vol lx ly lz
thermo_modify flush yes

fix hold all nvt temp {Tprod} {Tprod} 0.1

run 80000 # 80 ps

unfix hold
undump HOLD

write_data after_1000K_NVT_equilibrated.lmp
write_restart after_1000K_NVT_equilibrated.restart

============================================================

STEP 5: NVT production for He GB diffusion/segregation

============================================================

reset_timestep 0
timestep 0.001

dump PROD all custom 100 dump.he_gb_nvt_1000K.lammpstrj id type xu yu zu vx vy vz
dump_modify PROD sort id

compute msdGas gas msd com yes
compute msdUZr UZr msd com yes

fix MSD gas ave/time 100 10 1000 c_msdGas[1] c_msdGas[2] c_msdGas[3] c_msdGas[4] file He_MSD_1000K.txt

fix prod all nvt temp {Tprod} 0.1

thermo 10
thermo_style custom step temp pe ke etotal press pxx pyy pzz vol lx ly lz c_msdGas[4] c_msdUZr[4]
thermo_modify flush yes

run 1000000 # 1 ns

unfix MSD
unfix prod
undump PROD

write_data final_He_GB_1000K_NVT.lmp
write_restart final_He_GB_1000K_NVT.restart

2

I cannot comment on the DeepMD part. DeepMD is not part of LAMMPS and developed and maintained separately. Any questions related to the potential have to be directed at the DeepMD developers. They are off-topic in this forum, since this is about LAMMPS.

If the pressure in your simulation is high, this means that the system “wants” to expand. There are two possible reasons: 1) your density is too high for the chosen conditions, or 2) your potential does not represent the pressure for the chosen state of your system well. Generally, it is considered a mistake to assume that a geometry for a simulated system will match any experimental data exactly. Also, keep in mind that the density will change with temperature and that condensed systems are not very compressible, so small changes in density can result in very large changes in pressure.

It is up to you to determine which of the cases applies (it can also be both) and how to address them. There is no simple “do this, not that” kind of advice that can be given, especially since we don’t know which choices you have made and you probably need to consult with an expert (how about your adviser or supervisor?) on the type of system you are simulating what a reasonable choice would be.