[gmx-users] To get PMF using pull geometry of cylinder

[Justin A. Lemkul]

Lin, Dejun wrote:
> Hi all,
> 
> I want to do an umbrella sampling to calculate the potential of mean force of
> a micelle fusing to a lipid membrane. In the starting configuration, the
> membrane normal is set in z direction and the membrane is under the micelle
> (z(COM-membrane) < z(COM-micelle)). At the beginning, I used pull_geometry =
> distance but the sampling took to long to converge. And then I considered
> using cylinder instead and set up the pull code as: pull		= umbrella 
> pull_geometry	= cylinder pull_dim	= N N Y pull_start	= no pull_ngroups	= 1 
> pull_group0	= membrane pull_group1	= micelle pull_vec1	= 0.0 0.0 -1.0 pull_r1
> = 2.3 pull_r0         = 2.5 pull_init1	= 1.8 pull_rate1	= 0.0 pull_k1		= 1000
>  pull_nstxout	= 1000 pull_nstfout	= 1000
> 
> And the output files (px.xvg and pf.xvg) look like: (px.xvg) 0.0000		23.3397
> 4.90267 100.0000	        0.742426	        1.86002 200.0000		0.737842		1.86232
>  300.0000		0.762641		1.83985 400.0000		21.6814		1.86777 500.0000		0.771017
> 1.82482 600.0000		0.789036		1.83254 700.0000		0.793745		1.83503 800.0000
> 0.804732		1.84217 900.0000		0.80137		1.84166 1000.0000	0.817979		1.83546 
> 1100.0000	0.806272		1.83535 1200.0000	0.792355		1.80372 1300.0000	0.786839
> 1.8438 1400.0000	0.798686		1.82625 1500.0000	0.802859		1.84138 1600.0000
> 0.803844		1.84379
> 
> (pf.xvg) 100.0000		-480.158 200.0000		-498.547 300.0000		-318.839 400.0000
> -542.14 500.0000		-198.526 600.0000		-260.291 700.0000		-280.22 800.0000
> -337.385 900.0000		-333.317 1000.0000	-283.644 1100.0000	-282.761 1200.0000
> -29.7684 1300.0000	-350.409 1400.0000	-209.991 1500.0000	-331.029 1600.0000
> -350.309
> 
> And I visualized the trajectory in VMD and found that an "explosion" of the
> micelle (with all its residues scattered in the box). I then changed the
> pull_vec1	= 0.0 0.0 -1.0 to pull_vec1 = 0.0 0.0 1.0 and the "explosion" was
> gone and everything seemed to work. But I'm not quite sure if I'm setting up
> the pull-code right because if z(COM-membrane) < z(COM-micelle), then
> pull_vec1	= 0.0 0.0 -1.0 should be the right choice. Basically, I think my

The negative is incorrect.  If the micelle is "above" the membrane (i.e., it has 
a greater z-coordinate), then the reference vector should be positive with 
respect to the reference position.  It would appear from the pullx.xvg file that 
you've got some weird PBC crossing going on, since the reference positions 
oscillate dramatically at the beginning of the simulation.

> pull-code is restraining the COM of the micelle to the COM of a "cylinder" of
> the membrane at a distance of 1.8nm, with the vector pointing from
> COM(micelle) to  COM(cylinder) being (0,0,-1). So can anyone explain what
> caused the explosion or am I setting up the pull-code right? And what are the
> 2nd and 3rd column in the px.xvg output in this case?
> 

For such a simple restraint, I don't think there is any need for cylinder 
geometry.  You're not using different weighting.  Any of the other pull_geometry 
settings should work, with "distance" probably being the most intuitive.  I 
suspect it would be very easy to use the "position" setting as well.

The output in the pullx.xvg file should be labeled in the header.  The generic 
output is (x,y,z) of the reference group, followed by delta(x,y,z) for the 
restrained group, i.e. how far away the restrained group is along each 
coordinate axis.  If you're only restraining/pulling along the z-axis, then the 
only terms written are z and dz, since, in principle, x, y, dx, and dy are all 
fixed, with the last two terms being zero.

-Justin

-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================