[gmx-users] Free energy perturbation for larger molecule - big Coulomb solvation energies

[Krzysztof Makuch]

Hello!

I am trying to calculate solvation energy, using Free Energy Perturbation
for lutein in water and bilayer. Simulation runs smoothly, but Coulomb part
of calculated energy is extremly high for hydrophobic molecule
(C40H54(OH)2): above 40kcal/mol // 180 kJ/mol both in water and bilayer. I
am working with topology, with OPLS-AA charges.

To check possibility if the charges are somehow wrong I've parametrized
molecule with automatic server for charmm - cgenff. Parametrization of the
molecule in charmm clearly isn't perfect, but since what I want is only
rough estimation of possibility of wrong OPLS-AA charges it should be
sufficient. However, I am still getting even higher value of 300kJ/mol.
I've also checked methanol under same conditions in mdp file and obtained
results were reasonable.

Now, I suspect I've messed something in mdp FEP settings,could you please
look at it and give me a hint, what may be wrong with these?

Best,
Krzysiek

; Run control
integrator               = sd       ; Langevin dynamics
tinit                    = 0
dt                       = 0.002
nsteps                   = 1750000   ; 3.5 ns
nstcomm                  = 100
; Output control
nstxout                  = 500
nstvout                  = 500
nstfout                  = 0
nstlog                   = 500
nstenergy                = 500
nstxout-compressed       = 0
; Neighborsearching and short-range nonbonded interactions
cutoff-scheme            = verlet
nstlist                  = 20
ns_type                  = grid
pbc                      = xyz
rlist                    = 1.2
; Electrostatics
coulombtype              = PME
rcoulomb                 = 1.2
; van der Waals
vdwtype                  = cutoff
vdw-modifier             = potential-switch
rvdw-switch              = 1.0
rvdw                     = 1.2
; Apply long range dispersion corrections for Energy and Pressure
DispCorr                  = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing           = 0.1
; EWALD/PME/PPPM parameters
pme_order                = 5
ewald_rtol               = 1e-06
epsilon_surface          = 0
; Temperature coupling
; tcoupl is implicitly handled by the sd integrator
tc_grps                  = System
tau_t                    = 0.6
ref_t                    = 310
; Pressure coupling is on for NPT
Pcoupl                   = Parrinello-Rahman
pcoupltype         = isotropic
tau_p                    = 1.0
compressibility          = 4.5e-05
ref_p                    = 1.0
; Free energy control stuff
free_energy              = yes
init_lambda_state        = 0
delta_lambda             = 0
calc_lambda_neighbors    = 1        ; only immediate neighboring windows
; 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
vdw_lambdas              = 0    0    0    0    0    0    0    0    0
0    0    0    0    0    0    0    0    0    0    0    0
coul_lambdas             = 0    0.05    0.1    0.15    0.2    0.25
0.3    0.35    0.4    0.45    0.5    0.55    0.6    0.65    0.7    0.75
0.8    0.85    0.9    0.95    1
; We are not transforming any bonded or restrained interactions
bonded_lambdas           = 0    0    0    0    0    0    0    0    0
0    0    0    0    0    0    0    0    0    0    0    0
restraint_lambdas        = 0    0    0    0    0    0    0    0    0
0    0    0    0    0    0    0    0    0    0    0    0
; Masses are not changing (particle identities are the same at lambda = 0
and lambda = 1)
mass_lambdas             = 0    0    0    0    0    0    0    0    0
0    0    0    0    0    0    0    0    0    0    0    0
; Not doing simulated temperting here
temperature_lambdas      = 0    0    0    0    0    0    0    0    0
0    0    0    0    0    0    0    0    0    0    0    0
; 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
couple-moltype           = lutein  ; name of moleculetype to decouple
couple-lambda0           = vdw-q      ; all
couple-lambda1           = vdw     ; vdw
couple-intramol          = yes
nstdhdl                  = 100
; No velocities during EM
gen_vel                  = no
; options for bonds
constraints              = h-bonds  ; we only have C-H bonds here
; Type of constraint algorithm
constraint-algorithm     = lincs
; Do not constrain the starting configuration
continuation             = no
; Highest order in the expansion of the constraint coupling matrix
lincs-order              = 4

-- 
Jagiellonian University
Department of Computational Biophysics and Bioinformatics
tel.1: (12) 664 61 49
tel.2: (48) 664 086 049