I intend to do moving piston based nemd shock simulation for an explosive material. For this I am keeping a void portion at one edge. While doing so in NVE ensemble pressure value of the relaxed structure increases. When I do the simulation without void the pressure stays around 0, whereas with around 30-40A void size in X direction, the pressure increases upto 4000-5000 bar. Why is this happening? and How can I get rid of this?
The input script is as follows:
units real
atom_style full
boundary p p p
bond_style harmonic
angle_style harmonic
dihedral_style harmonic
improper_style cvff
read_data trial.dat
pair_style lj/cut/coul/long 9.5
include paircoefficient.txt
kspace_style pppm 1.0e-6
pair_modify mix arithmetic
special_bonds amber
special_bonds lj 0.0 0.0 0.5 coul 0.0 0.0 0.5
neighbor 2 bin
neigh_modify every 10 delay 0 check no
#============= minimize ================================================
EQUILIBRATION
timestep 0.1
#velocity all create 300 12345 mom yes rot no
#min_style cg
#min_modify line quadratic
#minimize 1e-40 1e-40 1000000 1000000
#write_data modeltkx50_14.data
#write_restart 50.restart
variable tstp equal 0.1
timestep ${tstp}
variable tdamp equal “v_tstp * 100”
velocity all create 300 12345 mom yes rot no
fix 1 all nve
fix 2 all langevin 300 300 ${tdamp} 90429997
compute ch_id all chunk/atom bin/1d x lower 0.05 units reduced
fix ave_ch_d all ave/chunk 1 10 100 ch_id density/mass file density.out
fix ave_ch_v all ave/chunk 1 10 100 ch_id vx file velocity.out
fix ave_ch_t all ave/chunk 1 10 100 ch_id temp file temperature.out
compute temperature all temp/chunk ch_id temp com yes
fix ave_ch_t2 all ave/time 1 10 100 c_temperature[1] file temperature2.out mode vector
set number of bins
variable nbins index 20
variable fraction equal 1.0/v_nbins
define bins as chunks
compute cchunk all chunk/atom bin/1d x lower ${fraction} units reduced
compute stress all stress/atom temperature
apply conversion to pressure early since we have no variable style for processing chunks
variable press atom -(c_stress[1]+c_stress[2]+c_stress[3])/(3.0vol{fraction})
variable shear atom -(c_stress[4])/(vol*{fraction})
compute binpress all reduce/chunk cchunk sum v_press
compute binshear all reduce/chunk cchunk sum v_shear
fix avg all ave/time 1 10 100 c_binpress mode vector file ave_stress.dat
fix avg_shear all ave/time 1 10 100 c_binshear mode vector file ave_shear.dat
thermo 1000
thermo_style custom step lx ly lz press pxx pyy pzz pe temp pxy pyz pxz density
dump 1.50 all custom 1000 dump1.mixed.lammpstrj id mol type q xu yu zu ix iy iz
run 80000
write_data modeltkx50_20.data
write_restart 50.restart
unfix 2