Hello, I am a beginner in the field of materials simulation and have just started learning LAMMPS.Recently, I’ve been encountering an energy explosion problem while conducting a simulation.This is a simulation that includes water molecules and upper and lower copper atomic walls. The simulation requires non-periodic boundaries.The scenario involves randomly filling a large number of TIP4P water molecules under a non-periodic boundary condition.After filling, neither using the CG algorithm for energy minimization nor using NVT for stabilization worked; the Thermo output showed complete pressure and potential energy maximization.The s-boundary doesn’t cause an energy explosion, however it changes the initial geometry I set. I tried changing it back to the f-boundary using change_box, but the energy explosion occurred. how can i solve this problem？Below is my inscript
units metal
boundary p p s
atom_style full
read_restart cavitation_restart.2000000.bin
pair_style hybrid lj/cut/tip4p/cut 1 2 1 1 0.1546 12.0 eam/alloy
bond_style harmonic
angle_style harmonic
pair_coeff 1 1 lj/cut/tip4p/cut 0.008035 3.1589
pair_coeff 2 2 lj/cut/tip4p/cut 0.0000 1.0
pair_coeff 1 3 lj/cut/tip4p/cut 0.064658 2.72808
pair_coeff 2 3 lj/cut/tip4p/cut 0.0 0.0
pair_coeff * * eam/alloy Cu.eam.alloy NULL NULL Cu
bond_coeff 1 19.5 0.9572
angle_coeff 1 2.38 104.52
change_box all boundary p p f
change_box all z final 0.0 94.5 remap units box
fix 1 all nvt temp 300 300 0.1
timestep 0.01
thermo_style custom step temp press vol
thermo 200
run 1000
change_box all z final 0.0 200 units box
group Cu type 3 # Cu atoms (type 3)
group water type 1 2 # water atoms (types 1 & 2)
region bot block INF INF INF INF 0 3 units box
region top block INF INF INF INF 92 95 units box
group gbot_region region bot
group gtop_region region top
group gbot intersect Cu gbot_region
group gtop intersect Cu gtop_region
group mobile subtract all gbot gtop
fix freeze_bot gbot setforce 0.0 0.0 0.0
velocity gbot set 0.0 0.0 0.0
fix fnvt gtop nvt temp 300.0 300.0 0.1
variable P0 equal 1200.0
variable conv equal 6.241509e-7
variable A equal lx*ly
variable Ntop equal count(gtop)
print “Ntop (number of atoms in gtop) = ${Ntop}”
variable F0 equal v_P0v_convv_A/v_Ntop
variable omega equal 2.0PI5.0e-3
variable Fz equal v_F0sin(v_omegatime)
apply force to top wall
fix fUS gtop addforce 0.0 0.0 v_Fz
timestep 0.001 # production timestep
compute zwater water chunk/atom bin/1d z lower 0.5 units box
fix zden water ave/chunk 100 10 1000 zwater density/mass file rho_water_z.profile
compute zall all chunk/atom bin/1d z lower 0.5 units box
fix zdenall all ave/chunk 100 10 1000 zall density/mass file rho_all_z.profile
---------------- output / trajectory ----------------
thermo_style custom step time temp press pe ke etotal
thermo 1000
dump snap all custom 1000 snap.lammpstrj id type mol q x y z ix iy iz
dump_modify snap first yes
run 100000
write final data
write_data 3system_from_restart.data
And this is what I meant by using change to change the s-boundary to the f-boundary, because I would encounter the problem of hydrogen atom loss when using the f-boundary.
