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

Maksim Mišin maksim.misin at strath.ac.uk
Tue Sep 8 16:28:09 CEST 2015


Sorry, forgot to specify charge.
For LJ particle it was +1 e, dipole had +1 one atom and -1 e on another.

On Tue, Sep 8, 2015 at 3:24 PM Maksim Mišin <maksim.misin at strath.ac.uk>
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?
>
>
> 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
>
>
>
>


More information about the gromacs.org_gmx-users mailing list