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

Zhe Wu zephyrwuz at gmail.com
Fri Nov 19 23:06:30 CET 2010

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

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.)


Links for the previous representative post related:

More information about the gromacs.org_gmx-users mailing list