[gmx-users] solvation free energy, electrostatic transformation

[Justin Lemkul]

On 7/27/15 9:55 AM, Daniele Veclani wrote:
> I'm trying to calculate the  solvation free energy of a molecule (M).
>
> I have done a VdW.  transformation.
>
> I have done also a  electrostatic transformation. But I get an incorrect
> value (large overestimation of DG solvation).
>

So presumably you have an experimental value of DGsolv?  How different are the 
values?

> I use gromacs 5.0.4 with gromos 54a7 force fields and my .mdp file (for
> lambda 00) is:
> integrator               = sd
> tinit                    = 0
> dt                       = 0.002
> nsteps                   = 2500000
> comm_mode                = angular
> nstcomm                  = 100
> nstxout                  = 500
> nstvout                  = 500
> nstfout                  = 0
> nstlog                   = 500
> nstenergy                = 500
> nstxout-compressed       = 0
> cutoff-scheme            = group
> nstlist                  = 0
> ns_type                  = simple
> pbc                      = no
> rlist                    = 0
> coulombtype              = cutoff
> rcoulomb                 = 0
> epsilon_r                = 1
> vdwtype                  = cutoff
> rvdw                     = 0
> DispCorr                  = no
> fourierspacing           = 0.12
> pme_order                = 6
> ewald_rtol               = 1e-06
> epsilon_surface          = 0
> tcoupl                   = berendsen
> tc_grps                  = system
> tau_t                    = 0.1
> ref_t                    = 300
> Pcoupl                   = no
> tau_p                    = 1.0
> compressibility          = 4.5e-05
> ref_p                    = 1.0
> free_energy              = yes
> init_lambda_state        = 0
> delta_lambda             = 0
> calc_lambda_neighbors    = 1
> vdw_lambdas              = 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
> 1.00 1.00
> coul_lambdas             = 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80
> 0.90 1.00
> bonded_lambdas           = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> 0.00 0.00
> restraint_lambdas        = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> 0.00 0.00
> mass_lambdas             = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> 0.00 0.00
> temperature_lambdas      = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> 0.00 0.00
> sc-alpha                 = 0.5
> sc-coul                  = no
> sc-power                 = 1.0
> sc-sigma                 = 0.3
> couple-moltype           = QUI
> couple-lambda0           = vdw-q
>
> couple-lambda1           = vdw
> couple-intramol          = no
> nstdhdl                  = 10
> gen_vel                  = no
> constraints              = all-bonds  ; we only have C-H bonds here
> constraint-algorithm     = lincs
> continuation             = no
> lincs-order              = 12
>
>
> What's wrong?
>

If you're getting a result that does not agree with a known quantity, without 
any errors in the run or the analysis, it means your topology for that molecule 
is bad and needs to be better parametrized.

-Justin

-- 
==================================================

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

==================================================