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

Justin A. Lemkul jalemkul at vt.edu
Thu May 26 17:26:59 CEST 2011

```
Lin, Dejun wrote:
> Thank you Justin. But I'm still confused about the pull_vec1. Is the vector
> pointing from the reference or from the pull group for all the geometry
> settings? For instance, if I use pull_geometry = position, pull_dim = N N Y,
> pull_rate1 = 0.1, pull_init1 = 3.0, pull_vec1 = 0 0 1, the position =
> pull_init + time*pull_rate*pull_vec1 according to the manual and it's
> actually pulling the pull group away from the reference because the vector
> has a positive z component.
>

The reference group serves as the reference :)  Thus, the vector connecting your
reference and pull group is (0,0,1) since the pull group has a larger
z-coordinate.  You're right about the position equation, but you've got a
pull_rate1 of zero, so the added term drops out and position = pull_init1.

> BTW, I did consider using position geometry but can I get the PMF using
> pull_rate NOT equal to zero? The reason I turned to cylinder is that with
> "distance" geometry the membrane concaved as the micelle was restrained to a
> very short distance (any distance below the critical distance where the
> fusion was supposed to occur) without fusing and this would persist for some
> hundred nanoseconds.
>

To obtain the PMF, the pull rate should be zero, using several (many)
independent simulations along the desired reaction coordinate. See the tutorial:

http://www.gromacs.org/Documentation/Tutorials#Pull_Code_and_Umbrella_Sampling

-Justin

> Thanks again for your help! Dejun
>
> ________________________________________ From: gmx-users-bounces at gromacs.org
> [gmx-users-bounces at gromacs.org] On Behalf Of Justin A. Lemkul
> [jalemkul at vt.edu] Sent: Wednesday, May 25, 2011 7:44 PM To: Discussion list
> for GROMACS users Subject: Re: [gmx-users] To get PMF using pull geometry of
> cylinder
>
> 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
>
> ======================================== -- gmx-users mailing list
> gmx-users at gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at
> don't post (un)subscribe requests to the list. Use the www interface or send
> it to gmx-users-request at gromacs.org. Can't post? Read
> http://www.gromacs.org/Support/Mailing_Lists
>

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

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

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

```