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 -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
