[gmx-users] question regarding to PMF calculation between 2 plates

[Zhe Wu]

Hi Chris,

Thank you so much for the idea! It is a good idea of using umbrella
sampling technique. But I am not sure whether it is just a matter of
averaging 3 PMFs as the final result.

To be physical, we need the 2 plates to be parallel to each other at
any time, right? With 3 pull groups (or more groups), indeed we
defined 2 plates. But in simulation, independent variation of these 3
distances will make angles between the plates. The angles could be
large when the umbrella is relatively small, as a balance between
distribution overlap and # of windows. For sure, their average
behavior should be parallel, but I don't know whether simple averaged
PMFs is physical. Of course, trying 4 or 5 pull groups as tests will
be helpful to examine this idea.

Meanwhile, we still need to deal with the rigidity of the plates: to
make sure they are flat during the simulation. So as I can see, it is
inevitable to use some sort restraint in the system anyway. That is a
reason why I used position restraint (PR) as mentioned in my last
email. As long as PR works in the NPT, I can just simply use the
potential obtained, followed by free energy perturbation scheme, get
the PMF. The problem is, the PR doesn't work in NPT in my test cases.

Again, thank you very much for helping me with your idea! I really
appreciate it. I am now testing even smaller time steps to see whether
PR will work (as suggested in discussion from Pawan, Florian and Mark
in a few emails after my post). Will keep you updated for any new
progress. Sure, further ideas will be really appreciated.

Thanks,
Zhe


On Sun, Nov 21, 2010 at 9:46 AM,  <chris.neale at utoronto.ca> wrote:
> Zhe, here's what comes to my mind. It certainly requires a whole bunch of
> testing, but if you really must calculate this then you could give it a try.
>
> In order to keep 2 plates planar, it seems to me that the distance between
> the planes must be controlled at 3 or more points (3 points defining a
> plane). Ideally, you can use 3 pull groups, each with a different reference
> group arranged in an equilateral triangle in one plane that each control the
> distance to a group in a similar triangle in the other plane. I would then
> do wham separately on each of the 3 pairs to get 3 evaluations of the PMF
> that I would average together. You could then redo this calculation for 4
> points in a square, 5 points, etc. to quantify any systematic drift in the
> resulting answer as you increase the number of restraints and find the
> asymptote.
>
> I'm not sure if the newest pull code can handle multiple reference groups
> though. If not, then you might need to triangulate distances and do the
> pulling in XYZ with many more than 3 groups.
>
> Finally, if you want to get fancy and collate all of the data together from
> many pull groups prior to putting it into wham, then you'll need to
> rigorously ensure that you get the math right. I would personally use only
> one restrained distance set per wham calculation.
>
> Since I don't do anything like this myself, I can't add much more. Good
> luck.
>
> Chris.
>
> -- original message --
>
> Dear GROMACS users,
>
> I have a question for calculating the potential of mean force between
> 2 plates (i.g., graphene) in water. It is not technically as simple as
> it seem. I found this question has been asked many times (I sum them
> up in the bottom of this note, so maybe it will be helpful for someone
> with similar question in future) but I cannot see a clear answer
> completely solving the problem. Information in the manual is quite
> limited on this topic as well.
>
> There are two ways I can see to do the calculation: 1. umbrella
> sampling; 2. free energy perturbation (JACS, 2005, 127, 3556). Both of
> them, however, need to deal with 2 questions first in order to get any
> meaningful result: 1. keep the plates to be really plates (flat
> surface); 2. inhibit plate rotation. There are 3 ways that come into
> my mind to solve these 2 problems, but none of them works as I tried.
> (Similar idea has been suggested in the previous posts regarding to
> the problem.)
>
> I am using a simple NPT system in the free energy perturbation scheme,
> because it is physically more meaningful and simpler to test than
> umbrella potential. But problem comes from the 3 method I tried:
> 1) Merely constructing a really complicated itp file with bond, angle,
> dihedral, imp dihedral cannot easily maintain a planar plate. Rotation
> of the plates is not easy to deal with, which I still don't have any
> idea yet. (Someone suggested other methods rather than this one in the
> previous post.)
> 2) Freezing all plates will make the plates do what they should do,
> but there is problem with the pressure coupling according to the
> manual. GMX 4.0.x wouldn't work as specified in the manual. GMX 3.3.3
> works as it seems but the system is very easy to blow up.
> 3) Use position RESTRAINT on all atoms on the plates in XYZ direction
> (basically fixing them in the initial position) as suggested in
> previous posts. As an example:
> -----------------------------------------------------------------------------------------------
> [ position_restraints ]
> ; ai funct fc
>  1  1    95000  95000   95000
>  2  1    95000  95000   95000
>  3  1    95000  95000   95000
>  4  1    95000  95000   95000
> .... every atom in the plates
> ------------------------------------------------------------------------------------------------
> However, or small force constant, the plate won't be maintained. For
> large force constant (the above example), the system will easily blow
> up even with time step of 1fs. For sure, once we have too large
> vibration in the system, an even smaller time step is required. Then
> the efficiency is going to be really really low. (To provide further
> details: I minimize the system with both steep and cg first and run
> the 400 ps NVT(Nose) run before the NPT (Parrinello, isotropic or
> semiisotropic) run. The test system is fairly small: 570 waters
> (SPC/E) with 2 plates of graphene (60 atoms each, amber parameters).
> So there should not be problem with the equilibration. I tuned the
> force constant in the restraint to see whether the plate change or
> not.)
>
> I will really appreciate it if someone can suggest some way out or
> where my problem is in the simulation. Is it the right way of
> restraining the system? Is there better way of doing this? It will be
> great if discussions coming along. Thank you in advanced for your
> help! (There has been tons of paper doing this with some other
> package, I am sure GROMACS can do it. I am just not the person who can
> figure this out.)
>
> Thanks,
> Zhe
>
>
> Links for the previous representative post related:
> http://lists.gromacs.org/pipermail/gmx-users/2010-February/048615.html
> http://lists.gromacs.org/pipermail/gmx-users/2010-January/048326.html
> http://lists.gromacs.org/pipermail/gmx-users/2009-December/047696.html
>
>
> --
> gmx-users mailing list    gmx-users at 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 thewww interface
> or send it to gmx-users-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>