Dear all, I would like to calculate the potential of mean force between two molecules in aqueous solution using the pull code. For a start i performed a number of calulations with a couple of very simple model systems with settings loosely based on the example in the gmx tutorial section (see mdp file included at the bottom). My box is 8x8x8 nm**3, pmf is calculated from some r_min up to approx 3.8 nm. Performing simulations with two t-buntanol molecules in implicit solvent, anlyzing the output (x and f*.xvg files) with g_wham, i get rather unexpected results:
The PMFs look very differnt depending on whether I use "pulldim Y Y Y" or "pulldim "N N Y": http://www.brunsteiner.net/tbut-pmf.gif Since we look at two neutral molecules with a permanent dipole moment one would expect the PMF to be negative up to a distance at which the two molecules come into contact, but with "pulldim Y Y Y" the PMF is positive, i.e., repulsive, throughout. With "pulldim N N Y" I constrain the two molecules in two dimensions (by freezing the central atom in the x and y dimensions only, BTW any ideas how else i could achieve that?) One could argue that a combination of frozen atoms and pull code might be problematic, but i freeze and pull in different dimensions, so this should be OK, and, more importantly: with "pulldim N N Y" i get results that look much more reasonable than with "pulldim "Y Y Y" and NO additional constraints. With "pulldim Y Y Y" the RDF (option nolog in g_wahm) looks, in fact, like a NON-normalized rdf: http://www.brunsteiner.net/tbut-rdf.gif and, interestingly if i normalize it myself (by dividing through z**2 before taking the logarithm) the resulting PMF looks much more reasonable. Although what I see suggests that with "pulldim Y Y Y" the RDF just doesn't normalized properly, this issue seems to be more involved: Looking at the average forces and positions (the content of the f*xvg and x*xvg files) http://www.brunsteiner.net/tbut-f.gif http://www.brunsteiner.net/tbut-z.gif suggests that there's already something wrong in the mdrun output meaning that the problem is further upstream and not connected to anything done by g_wham. I repeated the calculations with an even simpler system (2 single atoms that only interact via a LJ potential) to get qualitatively the same results: http://www.brunsteiner.net/2at-pmf.gif http://www.brunsteiner.net/2at-rdf.gif http://www.brunsteiner.net/2at-z.gif http://www.brunsteiner.net/2at-f.gif A few more remarks: 1) Convergence does not seem to be an issue here as i extended some of the calculations to include 35+ windows, 2 ns each, and the PMF remains the same nearly quantitatively. 2) The length of my reaction coordinates is always shorter than half the box length. 3) I've compared calculations cut-off vs PME, to get very similar results. 4) If I use "pulldim Y Y Y" with no additional constraints but with "comm_mode = Angular" i get results somewhere inbetween the two cases above. my specs: Gromacs version: 4.5.3 Linux 2.6.35-23-generic Ubuntu x86_64 gcc version 4.4.5 I am not sure whether i overlooked something in my input, or whether there's a problem with code. I'd be grateful for any ideas/suggestions! cheers, Michael ============= mdp file: integrator = sd ; this is better than md/NHT for systems with very few DOFs tau_t = 2.0 2.0 ref_t = 290 290 ; a bit lower than 298 since sd with default parameters ; has problems controlling the temperature. tc_grps = p1 p2 ; also tried System here, no difference ; dt = 0.002 tinit = 0 nsteps = 100000 ; played with that, results seem to have converged nstxout = 1000 nstvout = 0 nstfout = 1000 nstenergy = 100 ; constraint_algorithm = lincs constraints = hbonds continuation = no ; comm_mode = Linear nstcomm = 1 pbc = xyz ; also tried "pbc = no" same result nstlist = 10 ns_type = grid rlist = 4.0 rcoulomb = 4.0 rvdw = 4.0 ; coulombtype = Shift vdwtype = Shift epsilon_r = 80 ; pull = umbrella pull_geometry = distance pull_dim = N N Y ; or "Y Y Y" pull_start = yes pull_ngroups = 1 pull_group0 = p1 pull_group1 = p2 pull_rate1 = 0.0 pull_k1 = 500 ; i let this number vary depending on the pull distance; ; also played with the numbers, same result pull_nstxout = 100 pull_nstfout = 100 pull_pbcatom0 = 1 ; my system geometry is so that the molecules never leave the box pull_pbcatom1 = 16 ; or cross a cell boundary. ; freezegrps = c1 c2 ; with pulldim Y Y Y these lines are freezedim = Y Y N Y Y N ; removed or commented out ; energygrps = p1 p2 -- gmx-users mailing list gmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
