Asymmetric distribution of ions on two sides of charged surface

GROMACS version: Gromacs 2023.1
GROMACS modification: No
Hello everyone,
I am simulating the diffused double layer of a negatively charged surface(neutralized by K+) in a water-NaCl solution. Even though a long production simulation was performed, the four peaks (8.5A, 10A, 27.5A, 29A) of Na+ density (black curves) are still not equivalent. The third peak at 27.5A was not obvious either. What could be the cause that the peaks at two equivalent places on two sides of the surface are not at the same height?

Simulation system:
The surface is 1515 nm2 large and 1 nm thick. The box is 1515*10nm3. I put the surface 1nm from the bottom of the box.

MDP file:
; Bond parameters
continuation = yes ; Restarting after NPT
constraint_algorithm = lincs ; holonomic constraints
constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy

; Neighborsearching
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = 1.4 ; short-range electrostatic cutoff (in nm)
rvdw = 1.4 ; short-range van der Waals cutoff (in nm)

; 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
tc-grps = System ; Treat the whole system as one
bd_fric = 0 ; Brownian dynamics friction coefficient. When bd-fric is 0, the friction coefficient for each particle is calculated as mass/ tau-t.
tau_t = 1 ; time constant, in ps
ref_t = 300 ; reference temperature, one for each group, in K

image

I would say that the peaks look quite good to me, even if they are definitely not identical. There are two possible explanations I can think of:

  1. Simulation time. You say that a “long simulation was performed”, but you do not say how long. Are we talking microseconds?
  2. The surface is not completely symmetric and/or fluctuates slightly during the simulation.

I would also guess that the number of ions might make a difference. It seems like there are quite a few K+ ions, but not that many Na+ and Cl- ions in the system. More ions would almost certainly help making the distribution more symmetric.