[gmx-users] PMF from pull code, unexpected results
> Michael: your mdp options indicated using US and my stance is that, > using >US, you are incorrect to say that "And you will get the SAME > result if you >do this calculation in one or in three dimensions > (pulldim NNY or pulldim >YYY)". this has nothing to do with the mdp file ... i was talking about a purely hypothetical case, if you use pull = constraints, rather than umbrella sampling you don't NEED to run the calculation, since, as i said you can figure out the result in your head ... with constraints the "average" force will always be the exact LJ force at that distance, and if you integrate that you will necessarily end up with the original LJ potential, no matter whether you do that in one or in three dimensions! if you only have two spherically symmetric particles, with an interaction energy that depends only on the distance, and no solvent then the "free energy" and the "energy" are identical (apart from the entropic contribution discussed in the Neumann paper) ... that means: if you COULD do this calculation with constraints (in practice you probably can't because there'd be too few DOFs and the temperature would be ill defined) you SHOULD get identical results, independent of pulldim ... unless the term in eqn 6.1 is somehow included explicitly, but i can't see that anywhere in the code ... > Nevertheless, I'm a little perturbed by > your call out for help from somebody else "in the know", so I'll leave > > you >at this point. Good luck. thanks, doesn't look as if there was much interest, I guess I'll have to go back to my statistical mechanics textbooks... mic -- gmx-users mailing listgmx-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
[gmx-users] PMF from pull code, unexpected results
Michael: your mdp options indicated using US and my stance is that, using US, you are incorrect to say that "And you will get the SAME result if you do this calculation in one or in three dimensions (pulldim NNY or pulldim YYY)". Nevertheless, I'm a little perturbed by your call out for help from somebody else "in the know", so I'll leave you at this point. Good luck. Chris. -- original message -- Re: PMF from pull code, unexpected results Michael Brunsteiner mbx0009 at yahoo.com Wed Feb 23 23:16:07 CET 2011 * Previous message: [gmx-users] widom insertion * Next message: [gmx-users] trr files transferability * Messages sorted by: [ date ] [ thread ] [ subject ] [ author ] chris.neale at utoronto.ca wrote: Looking at http://www.brunsteiner.net/tbut-pmf.gif you seem to be getting exactly what I would expect. The difference is the entropy term. Note that the spherical shell increases volume as r increases for pulldim YYY but this effect is absent for pulldim NNY. This is why you can correct as per an RDF. The "entropy term" (Eqn 6.1 in the manual) has already been the source of some confusion in this list ... If you make a "Gedankenexperiment" and consider the PMF between two atoms in vacuum that interact ONLY through a simple radial, for example a LJ, potential, and then calculate the PMF using the pull-code ... if you don't do umbrella sampling + WHAM but instead use constraints (pull = constraint) you can do the calculation in your head, to find that the resulting PMF is, of course, nothing else but the original LJ potential. And you will get the SAME result if you do this calculation in one or in three dimensions (pulldim NNY or pulldim YYY) so where does the entropy and the correction term in eqn 6 come in here?? In Section 6.1 of the manual it says: "But when calculating a PMF between two solutes in a solvent, for the purpose of simulating without solvent, the entropic contribution should be removed." I find this a bit confusing, first ... "simulating without solvent" ... why would that be my (or anybodies) purpose? and second: according to eqn 6.1 this "entropic contribution" is: PMF(r) = - (nc-1)kBT log(r) for this to make sense r would have to be dimensionless wouldn't it? I could divide r by the distance at which i arbitrarily chose my PMF to vanish, call it r_c, which would have the advantage that then the correction is zero at this distance ... but then this is just wild guess of mine ... anyway, that would mean that, if I call the PMF coming out of mdrun/g_wham PMF_wham, then the "true" PMF is: PMF(z) = PMF_wham (z) + (nc-1)kBT log(z/r_c) (the plus comes from the double minus, as in: "removing" something thats negative) is that correct? ... and if so, why does it, seemingly, not apply in the above Gedankenexperiment? Yours truly (and maybe some other people) would really appreciate if somebody in the know would clarify that! thanks in advance! Michael -- gmx-users mailing listgmx-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
[gmx-users] PMF from pull code, unexpected results
Looking at http://www.brunsteiner.net/tbut-pmf.gif you seem to be getting exactly what I would expect. The difference is the entropy term. Note that the spherical shell increases volume as r increases for pulldim YYY but this effect is absent for pulldim NNY. This is why you can correct as per an RDF. Please note that I didn;t check everything thoroughly or look at the other images so there could still be something weird going on, but I disagree that http://www.brunsteiner.net/tbut-pmf.gif shows anything strange. Chris. -- original message -- 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 = 10; 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
[gmx-users] PMF from pull code, unexpected results
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 = 10; 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 listgmx-users@gromacs.org http://li
