[gmx-users] Charging free energies are unaffected by couple-intramol

[Justin Lemkul]

On 9/8/15 10:24 AM, Maksim Mišin wrote:
> Dear Gromacs users,
>
>
> I was trying to compute charging free energies of simple solutes such as
> ions and dipoles using BAR.
>
> In the manual it says that setting mdp option couple-intramol to no makes
> decoupled state of the molecule correspond to the proper vacuum state
> without periodicity effects. However, two runs with this option on and off
> gave practically identical results:
>
> couple-intramol = no
> LJ particle:  -59.667  +-  0.037 kcal/mol
> Dipole:        -35.172  +-  0.032 kcal/mol
>
> couple-intramol = yes
> LJ particle:  -59.729  +-  0.037 kcal/mol
> Dipole:        -35.177  +-  0.032 kcal/mol
>
> If I understand everything correctly (and this has been mentioned recently
> on mail list:
> https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users/2015-July/099237.html)
> these energies should differ by the charging free energy in vacuum, which
> for these molecules is:
>
> LJ particle: 18.12  kcal/mol
> Dipole:        0.074 kcal/mol
>
> I have also taken look at tpr files produced by grompp using gmxdump.
> Except for different initial velocities tpr files are identical with either
> option, which is somewhat strange.
>
> Is this a bug or intended behaviour?
>

For a dipole, intended behavior, or at least, it depends on the value of nrexcl. 
  In a typical setup, atoms connected directly by a bond experience no nonbonded 
interactions.  This is the way all additive force fields work.  So here, the use 
of couple-intramol is irrelevant.

-Justin

>
> Parameters:
> LJ particle: sigma = 3.8 A, eps = 0.125 kcal/mol
> Dipole: same LJ parameters for both atoms with bond length = 2 A
> Water model is SPC/E; cubic box with a length of 2.5 nm was used for both
> solvent and vacuum simulations. Other free energy estimation methods such
> as TI and MBAR gave very similar results.
>
> Here is typical mdp file I used for these calculation:
>
> ; MD Run control
> integrator               = sd
> tinit                    = 0
> dt                       = 0.002
> nsteps                   = 2500000
> nstcomm                  = 100
>
>   ; Output control
> nstxout                  = 0
> nstvout                  = 0
> nstfout                  = 0
> nstlog                   = 500
> nstenergy                = 500
> nstxout-compressed       = 0
>
> ; Neighboursearch and short-range nondbonded interactions
> cutoff-scheme            = verlet
> nstlist                  = 40
> ns_type                  = grid
> pbc                      = xyz
> rlist                    = 1.2
>
> ; Electrostatics
>   coulombtype              = pme
> rcoulomb                 = 1.2
> epsilon-rf               = 73.5
> ; van der Waals
> vdwtype                  = cutoff
> vdw-modifier             = potential-switch
> rvdw-switch              = 1.0
> rvdw                     = 1.2
> DispCorr                 = EnerPres
> fourierspacing           = 0.15
> ; EWALD/PME/PPPM parameters
> pme_order                = 6
> ewald_rtol               = 1e-06
> epsilon_surface          = 0
>
> ; Temperature coupling
> tc_grps                  = system
> tau-t                    = 1.0
> ref-t                    = 298.15
>
> ; Presure coupling
> Pcoupl                   = berendsen
> tau_p                    = 1.0
> compressibility          = 4.5e-05
> ref_p                    = 1.0
>
>   ; Free energy control
> free_energy              = yes
> init_lambda_state        = 0
> delta_lambda             = 0
> calc_lambda_neighbors    = -1
> coul_lambdas      = 0.0 0.25 0.5 0.75 1.0
>
> ; Options for the decoupling
> sc-alpha                 = 0.5
> sc-coul                  = no
> sc-power                 = 1.0
> sc-sigma                 = 0.3
> couple-moltype           = MOL
> couple-lambda0           = vdw-q
> couple-lambda1           = vdw
> couple-intramol          = no  ; (or yes)
> nstdhdl                  = 100
>
> ; Velocity generation
> gen-vel                 = yes
> gen-temp                = 298.15
>
> ; Bond constraint
> constraints              = all-bonds
> constraint-algorithm     = lincs
> continuation             = no
> lincs-order              = 12
>
>
> Thank you,
> Maksim Mišin
>

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

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

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