Lipid bilayer membrane splits into two after nvt equilibration?

GROMACS version: 5.14
GROMACS modification: Yes/No

I am running lipid bilayer simulation with gram negative bacterial membrane. After running EM, my system looks exactly how I am expecting it to look, less ordered etc. However, after I run nvt equilibration for various lengths (100ps-2ns) I always get same result and that is that each of the layer splits from each other while water wall becomes smaller. I have tried modifying thermostat, running simulations for longer…but nothing helped. Anyone has any suggestions?

My nvt file looks like this:

title       = NVT equilibration
define      = -DPOSRES  ; position restrain the protein
; Run parameters
integrator  = md        ; leap-frog integrator
nsteps      = 15000     ; 2 * 50000 = 100 ps
dt          = 0.002     ; 2 fs
; Output control
nstxout     = 100       ; save coordinates every 0.2 ps
nstvout     = 100       ; save velocities every 0.2 ps
nstenergy   = 100       ; save energies every 0.2 ps
nstlog      = 100       ; update log file every 0.2 ps
; Bond parameters
continuation            = no            ; first dynamics run
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
ns_type     = grid      ; search neighboring grid cels
nstlist     = 5         ; 10 fs
rlist       = 1.2       ; short-range neighborlist cutoff (in nm)
rcoulomb    = 1.2       ; short-range electrostatic cutoff (in nm)
rvdw        = 1.2       ; 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
tcoupl      = V-rescale                         ; modified Berendsen thermostat
tc-grps     = MEMB  SOL_ION      ; two coupling groups - more accurate
tau_t       = 0.1   0.1                 ; time constant, in ps
ref_t       = 298   298                 ; 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
; Dispersion correction
DispCorr    = EnerPres  ; account for cut-off vdW scheme
; Velocity generation
gen_vel     = yes       ; assign velocities from Maxwell distribution
gen_temp    = 298       ; temperature for Maxwell distribution
gen_seed    = -1        ; generate a random seed
; COM motion removal
; These options remove motion of the protein/bilayer relative to the solvent/ions
nstcomm     = 1
comm-mode   = Linear

comm-grps = MEMB SOL_ION