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

[chris.neale at utoronto.ca]

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