SPCE water density at the boundary

GROMACS version:
GROMACS modification: Yes/No
Greetings of the day!
I am running MD simulations of only the SPCE water box with a length of 20 x 20 x 20 angstroms at NVT ensembles. I calculated the density of water along the Z-axis. At the boundary (0 and 20 angstroms) the density (at first and last bins) is less than 1 g/cm^3. Please see the attached input file and density file, let me know where is wrong. I am expecting density at the boundaries also 1 g/cm^3. Any help is appreciated.
Thank you very much.
Best regards,
Nagaraju

Here is my input file:
########################################
title = Bulk Water (268 molecules) NVT equilibration
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 5000000 ;1 * 5000 = 5000 ps (5 ns)
dt = 0.001 ; 1 fs

; Output control
nstxout = 50 ; save coordinates every 0.05 ps (50 fs)
nstvout = 50 ; save velocities every 0.05 ps (50 fs)
nstenergy = 50 ; save energies every 0.05 ps (50 fs)
nstlog = 50 ; update log file every 0.05 ps (50 fs)
nstxout-compressed = 50 ; save compressed coordinates every 10.0 ps
compressed-x-grps = System ; save the whole system

; Bond parameters
continuation = yes ; restarting dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy

; Nonbonded settings
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet
rcoulomb = 0.9 ; short-range electrostatic cutoff (in nm)
rvdw = 0.9 ; short-range van der Waals cutoff (in nm)

; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme

; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT

; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = SOL ; two coupling groups - more accurate
tau_t = 0.1 ; time constant, in ps
ref_t = 300 ; reference temperature, one for each group, in K

; Pressure coupling is off
pcoupl = no ; no pressure coupling in NVT

; Periodic boundary conditions
pbc = xyz ; 3-D PBC

; Velocity generation
gen_vel = no ; Velocity generation is off