[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