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

Justin A. Lemkul jalemkul at vt.edu
Thu May 26 19:09:54 CEST 2011

```
Lin, Dejun wrote:
> Ah...Sorry if I confused you. Actually what I meant is if I set the pull code, e.g., as:
> pull_geometry = position
> pull_dim = N N Y
> pull_rate1 = 0.1
> pull_init1 = 3.0
> pull_vec1 = 0 0 1
> (where pull rate is NOT zero)
> Then the position = pull_init + time*pull_rate*pull_vec1 and since  pull_vec1 is z-component-positive (0, 0, 1), the position between the reference and the pull group increases monotonically and that's not what I want. I want the distance to decrease so the pull_vec1 should be (0,0,-1). Am I right?
>

If you want the two species to approach rather than separate, do not change the
pull_vec, change the pull_rate so that it is less than zero.  A negative
pull_rate will decrease the separation over time.

Note that pull_dim is irrelevant for pull_geometry=position; only pull_vec1 is
required.

-Justin

> Thanks,
> 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: Thursday, May 26, 2011 11:26 AM
> To: Discussion list for GROMACS users
> Subject: Re: [gmx-users] To get PMF using pull geometry of cylinder
>
> 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
>
> ========================================
> --
> gmx-users mailing list    gmx-users at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-request at gromacs.org.
>

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

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

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

```