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

[Lin, Dejun]

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 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?

Thanks,
Dejun