Hello all,
Something unusual is happening in my simulation. I’m modeling a supercell of the polymer Kapton, which I intend to bombard with oxygen atoms. The process begins by finding the equilibrium structure of Kapton, after which I extend the supercell along the z-axis. This creates a vacuum region above the Kapton layer, where I plan to introduce oxygen atoms with a defined velocity directed toward the polymer surface.
I’ve followed this same procedure with several different Kapton structures, always using the same code to generate the initial configurations. Before introducing the oxygen atoms, I apply a standard heating and quenching method to stabilize the system’s temperature — a process that has consistently worked in all previous cases.
However, with my latest supercell — the one I actually want to use — one polymer chain at the bottom of the supercell begins to decompose during the NVT ensemble phase of temperature stabilization. I suspected this might be related to changes in system density: by extending the z-dimension, Kapton might be experiencing unwanted compression in the x and y directions. To address this, I performed a few initial steps in the NPT ensemble before switching to NVT. Interestingly, no decomposition occurs during the NPT phase, but as soon as I switch to NVT, the same chain at the bottom decomposes again.
Here is my code
# Example of relaxation for the monomer
package kokkos
#-------------------------VARIABLES---------------------------#
variable dt equal 0.1 #fs
# lower z surface (look at surface using ovito,you will understand)
variable z_surface equal 142
# size supercell in z direction
variable z_max_supercell equal 700
# z at surface of Kapton
variable z_max equal 156
# Atoms with z between INF and z_frozen will be frozen
variable z_frozen equal 0
# Introduction AO zone :
variable zhi equal $(v_z_max + 80)
variable zlo equal $(v_z_max + 60)
# H&Q steps (increase step for better temperature stabilization or decrease T_step, check temperature constant in NVE ensemble)
variable T_heating_step equal 30.0 # delta T increase at each cycle during heating
variable T_quenching_step equal 20.0 # delta T decrease at each cycle during quenching
variable heating_step equal 3500 # nb steps for increasing T of T_heating_step
variable quenching_step equal 2500 # nb steps for decreasing T of T_quenching_step
variable init_T equal 100
variable inter_T equal 450
variable final_T equal 300
#-------------------------------------------------------------#
units real
boundary p p f
atom_style charge
read_data ./data.huge_supercell
# Define ReaxFF potential
pair_style reaxff/kk lmp_control
pair_coeff * * "/stck/dfernand/ATOX/2-LAMMPS/0-ReactiveForceFields/CHONSSiNaFZr_Rahnamoun_van-Duin.ff.txt" C H O N
# Neighbor list settings
neighbor 3 bin #radius cutoof interaction
neigh_modify every 1 delay 0 check yes
timestep ${dt}
# Charge equilibration
fix 2 all qeq/reaxff/kk 1 0.0 10.0 1.0e-6 reaxff # Charge equilibration
# Output settings
dump myDump all custom 300 dump.stabilization id type x y zu vx vy vz
thermo 20 # Print thermodynamic data every N steps
thermo_style custom step vol temp press density etotal
fix 3 all npt/kk temp ${init_T} ${init_T} $(100*dt) x 1.0 1.0 $(500*dt) y 1.0 1.0 $(500*dt)
run 3500
unfix 3
write_data data.after_npt
# Avoid loosing atoms in z direction as no peroidic boundary conditions in z
fix wall_zlo all wall/reflect zlo EDGE
fix wall_zhi all wall/reflect zhi EDGE
# Avoid that the polymer translatates in z direction
region bloc block INF INF INF INF INF ${z_frozen}
group frozen dynamic all region bloc every 3000
fix freeze frozen setforce 0 0 0
velocity frozen set 0 0 0
#------------------------------ T stabilization ----------------------------------#
#°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°° HEATING °°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°#
variable T equal $(15 + v_init_T)
thermo_style custom step vol temp press density etotal v_T
label heat_loop
fix 3 all nvt/kk temp ${T} $(v_T_heating_step + v_T) $(100*dt)
run ${heating_step}
fix 3 all nvt/kk temp $(v_T_heating_step + v_T) $(v_T_heating_step + v_T) $(100*dt)
run ${heating_step}
# Avoid that the polymer translatates in z direction
region bloc block INF INF INF INF INF ${z_frozen}
group frozen dynamic all region bloc every 3000
fix freeze frozen setforce 0 0 0
velocity frozen set 0 0 0
variable T equal $(v_T_heating_step + v_T)
if "$(v_T) < $(v_inter_T)" then "jump SELF heat_loop"
write_data data.after_heating
#°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°#