Use implicit solvent in Langevin simulation

GROMACS version: 2016-2018
GROMACS modification: Yes/No

Dear GMX users,

I know the implicit solvent function has been removed from the newer version, I recently need to run some Langevin simulations with implicit solvent.
I have Na+ and Cl- seperated at 2 nm in a vacuum box. For comparsion, I changed the “gb_epsilon_solvent” and extracted the force on Na atom as a function of the Na-Cl distance, but the relation seems to be unaffected by the choice of epsilon.
This doesn’t seem right, so I am asking if there is something I missed in my mdp file or how can I look at the different implicit solvent screening effect.

For reference, here are the mdp parameters I used:

#################################### INPUT ####################################
ld_seed = -1 ; Use random seed
tinit = 0 ; Initial time is 0
integrator = sd ; Langevin thermostat
dt = 0.002 ; Timestep (ps)
nsteps = 25000 ; Simulation duration (timesteps)
nstcomm = 250 ; Center of mass motion removal interval
comm_mode = linear ; Center of mass motion removal mode

ref_t = 300 ; System temperature (K)
gen_vel = yes ; Initialize velocities from Maxwellian distribution
tau_t = 2.0 ; Thermostat time constant (ps)
tc_grps = system ; Apply thermostat to complete system

############################## IMPLICIT SOLVENT ###############################
implicit_solvent = GBSA ; Generalized Born implicit solvent
gb_algorithm = OBC
rgbradii = 0 ; Cutoff for Born radii calculation (A)
gb_epsilon_solvent = 10 ; dielectric constant of implicit solvent
nstgbradii = 0
gb_dielectric_offset = 0.009
gb_saltconc = 0 ; Salt Concentration
gb_obc_alpha = 1
gb_obc_beta = 0.8
gb_obc_gamma = 4.85
sa_algorithm = Ace-approximation
sa_surface_tension = 2.25936

cutoff_scheme = group ; Method of managing neighbor lists
pbc = no ; Periodic boundary conditions disabled
coulombtype = cut-off ; Calculate coulomb interactions using cutoff
rcoulomb = 0 ; Coulomb cutoff of infinity
vdw_type = cut-off ; Calculate van der Waals interactions using cutoff
rvdw = 0 ; Van der Waals cutoff of infinity
rlist = 0 ; Neighbor list cutoff
nstlist = 0 ; Do not update neighbor list

nstlog = 50 ; Log output interval (timesteps)
nstenergy = 50 ; Energy output interval (timesteps)
nstcalcenergy = 50 ; Energy calculation interval (timesteps)
nstxout = 50 ; Trajectory output interval (timesteps)
nstvout = 50 ; Velocity outout interval (timesteps)
nstfout = 50 ; Force output interval (timesteps)

Thank you for any suggestion!

Shi