Thin films atoms move across Z direction

Hi ,

I used below gro.mdp to run polystyrene (PS-top),sulphonated PS (PSS-bottom) thin films.
I could see once visualized PS goes to down as in the below visual.
Please let me know how can I prevent (atoms moving in Z direction) this by modifying the below gro.mdp file.

; Time step in ps
dt = 0.002

; Number of steps to run
nsteps = 600000000

; Remove center of mass motion every 100 steps
comm-mode = Linear
nstcomm = 100

; Write positions and velocities to the trr file every 5000 steps
; Output frequency for coords (x), velocities (v) and forces (f)
nstxout = 500000
nstvout = 500000

; Write to the log and energy file every 5000 steps
nstlog = 500000
nstenergy = 500000

; Calculate energies every 100 steps
nstcalcenergy = 100

; Write positions to three decimal places to xtc file every 5000 steps
nstxout-compressed = 500000
compressed-x-precision = 1000

; Use the Verlet cutoff scheme with grid neighbor searching; you should basically always use this when you can (the only time you wouldn’t would be if you were using tabulated potentials, buckingham interactions, or energy group exclusions; see manual for more detail)
cutoff-scheme = Verlet
ns_type = grid

; Periodic in all three dimensions
pbc = xyz

; Minimum cutoffs in nm for neighbor list, electrostatics, and LJ interactions. mdrun will tune rcoulomb for improved performance.
rlist = 1.2
rcoulomb = 1.2
rvdw = 1.2

; Use the PME algorithm (an alternative to PPPM) (this is very likely the electrostatics algorithm you want to use)
coulombtype = PME

; Smoothly shift vdw interactions to zero at the cutoff
vdw-type = Cut-off
vdw-modifier = Potential-shift-Verlet

; GENERATE VELOCITIES FOR STARTUP RUN =
gen-vel = yes
gen-temp = 500
gen-seed = 173529

; Use the v-rescale thermostat with a time constant of 0.1 ps and temperature of 600. V-rescale correctly samples the canonical distribution. Nose-hoover would be another good choice.
tcoupl = V-rescale
tc-grps = System
tau_t = 0.1
ref_t = 500

; Don’t use pressure coupling. For NPT this should be set. For equilibration, Berendsen is best since it is more stable. For production runs, Parrinello-Rahman is best since it samples the correct ensemble.
pcoupl = no
pcoupltype = Isotropic
tau-p = 1
refcoord-scaling = No

; Use the LINCS algorithm to turn all bonds with hydrogen atoms into constraints. This allows a larger timestep. (can also use shake if needed, but lincs is faster)
constraint_algorithm = lincs
constraints = H-bonds

; Do not generate initial velocities; rather this run is a continuation from a previous run. Alternately, set gen_vel to yes and give it a temperature to sample the maxwell distribution at, and set continuation to no.
;gen_vel = no
;continuation = yes

Position restraints may provide a potential solution, but it’s essential to thoroughly evaluate if they are suitable for your specific project and research goals. You can find more information about different types of restraints here.