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.mi...@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 > > > > -- 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.
