Using implicit solvent model in binding/hydration free energy simulations

User discussions
1

GROMACS version: 5.1.3
GROMACS modification: No

Hi,

I’m trying to use implicit solvent model (GB) in running FEP simulations for absolute binding free energy. As a test, I started with running FEP simulations to compute absolute hydration free energies of some small molecules (e.g. Ethanol). I input the GBSA parameters for GAFF atom types into the methanol topology. Then in the mdp file for production run I specified no PBC condition and no pressure coupling (NVT). For implicit solvent model parameters I specified as following:

; SOLVENT MODEL
implicit_solvent = GBSA
gb_algorithm     = OBC
nstgbradii       = 1
gb_epsilon_solvent = 80.0
sa_algorithm    = Ace-approximation
gb_obc_alpha    =   1    ; for OBC2
gb_obc_beta     =   0.78   ;for OBC2
gb_obc_gamma    =   4.85   ; for OBC2

For free energy control part:

; Free energy control stuff
free_energy              = yes
init_lambda_state        = 0
delta_lambda             = 0
calc_lambda_neighbors    = -1       ; all lambda points written out (for MBAR)
couple-moltype           = Ethanol  ; name of moleculetype to decouple
couple-lambda0           = vdw-q    ; All interactions on at lambda=0
couple-lambda1           = none     ; turn off everything
couple-intramol          = no
; Vectors of lambda specified here
; Each combination is an index that is retrieved from init_lambda_state for each simulation
; init_lambda_state        0    1    2    3    4    5    6    7    8    9    10   11   12   13   14   15   16   17   18   19   20
fep_lambdas              = 0.00 0.20 0.40 0.60 0.80 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
vdw_lambdas              = 0.00 0.00 0.00 0.00 0.00 0.00 0.05 0.10 0.20 0.30 0.40 0.50 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
; Options for the decoupling
sc-alpha                 = 0.5
sc-coul                  = no       ; linear interpolation of Coulomb (none in this case)
sc-power                 = 1
sc-sigma                 = 0.3
nstdhdl                  = 1000

Other parameters are not shown here but can be posted if it is helpful.
Then I ran my simulation without any warnings/errors. But in the outcome dhdl file, all of the columns are 0.0 for any lambda ensemble. It seems I’m simulating my system in the vacuum or the GB parameters are not dependent on the lambda values. Could you please give me some suggestions on how to fix it? Thank you!

2

The combination of implicit solvent and the free energy code has probably never been tested and was never intended for use. Implicit solvent features have been removed from GROMACS due to limited use and poor compatibility.