[gmx-users] Obtaining PMF for change in domain position

Abhi Acharya abhi117acharya at gmail.com
Sun Jan 4 18:48:56 CET 2015

Hello Gmx Users,

As suggested previously in this thread, I tried to use enforced rotations
on my system. First I calculated the amount of domain rotation and
translation required to obtain conformation 2 from conformation 1 using
'measure rotation' tool in chimera. I received an output as below:

  Matrix rotation and translation
    -0.96887254  -0.01658882  -0.24700365  59.55727674
     0.17648032  -0.74599360  -0.64214347  85.34692780
    -0.17361074  -0.66574646   0.72570032  11.00502990
  Axis  -0.11353457  -0.35303290   0.92869676     (Axis vector for rotation
  Axis point  28.32087779  46.35125870   0.00000000     (point on the axis)
  Rotation angle (degrees) 174.03353949
  Shift along axis -26.67174721  (in Angstroms)

For generating the windows for umbrella sampling, I used the pull code for
generating the translation and enforced rotations for the required rotation
simultaneously. For pull code I used the axis vector as the pull vector and
a pulling rate of -0.0002667 nm/ps for a 10 ns simulation. For enforced
rotations, I used the same Axis vector as the rotation vector and a
rotation rate of 0.0174 deg/ps.  I have provided the relevant parameter
sections from the .mdp file below.

Although I can see the expected movement of the domain of interest, but in
the final conformation it fails to reach the intended final position; falls
short by ~21 Angstroms. Visual analysis revealed that although the domain
rotated correctly, the translation of the domain did not work effectively,
so I suspect the pull code didn't work as I expected. What could be the
possible reasons for such a behavior ? Maybe I am doing some thing wrong.

The link to image containing the superposition of the initial, final
obtained and final expected conformations is provided in the link for

The .mdp portion is as follows:

; Pull type: no, umbrella, constraint or constant-force
pull                     = umbrella
; Pull geometry: distance, direction, cylinder or position
pull_geometry            = direction
; Select components for the pull vector. default: Y Y Y
pull_dim                 = Y Y Y
; Cylinder radius for dynamic reaction force groups (nm)
pull_r1                  = 1
; Switch from r1 to r0 in case of dynamic reaction force
pull_r0                  = 1.5
pull_constr_tol          = 1e-06
pull_start               = no
pull_nstxout             = 10
pull_nstfout             = 1
; Number of pull groups
pull_ngroups             = 1
; Group name, weight (default all 1), vector, init, rate (nm/ps),
pull_group0              = KH-domain
pull_weights0            =
pull_pbcatom0            = 0
pull_group1              = GD-domain
pull_weights1            =
pull_pbcatom1            = 0
pull_vec1                = -0.11353457 -0.35303290 0.92869676
pull_init1               = 0.0
pull_rate1               = -0.0002667
pull_k1                  = 1000
pull_kB1                 = 1000

; Enforced rotation: No or Yes
rotation                 = Yes
; Output frequency for angle, torque and rotation potential energy for the
whole group
rot_nstrout              = 100
; Output frequency for per-slab data (angles, torques and slab centers)
rot_nstsout              = 1000
; Number of rotation groups
rot_ngroups              = 1
; Rotation group name
rot_group0               = GD-C-alpha
; Rotation potential. Can be iso, iso-pf, pm, pm-pf, rm, rm-pf, rm2,
rm2-pf, flex, flex-t, flex2, flex2-t
rot_type0                = flex
; Use mass-weighting of the rotation group positions
rot_massw0               = no
; Rotation vector, will get normalized
rot_vec0                 = -0.11353457 -0.35303290 0.92869676
; Pivot point for the potentials iso, pm, rm, and rm2 (nm)
rot_pivot0               = 2.832087779  4.635125870   0.00000000
; Rotation rate (degree/ps) and force constant (kJ/(mol*nm^2))
rot_rate0                = -0.0174
rot_k0                   = 400
; Slab distance for flexible axis rotation (nm)
rot_slab_dist0           = 1.5
; Minimum value of Gaussian function for the force to be evaluated (for
flex* potentials)
rot_min_gauss0           = 0.001
; Value of additive constant epsilon' (nm^2) for rm2* and flex2* potentials
rot_eps0                 = 0
; Fitting method to determine angle of rotation group (rmsd, norm, or
rot_fit_method0          = norm
; For fit type 'potential', nr. of angles around the reference for which
the pot. is evaluated
rot_potfit_nsteps0       = 21
; For fit type 'potential', distance in degrees between two consecutive
rot_potfit_step0         = 0.25

Thanks in advance.

Abhishek Acharya
Shasara Research Foundation

On Tue, Dec 30, 2014 at 5:15 PM, Abhi Acharya <abhi117acharya at gmail.com>

> Abhishek Acharya
> Shasara Research Foundation
> On Tue, Dec 30, 2014 at 10:42 AM, Justin Lemkul <jalemkul at vt.edu> wrote:
>> On 12/29/14 11:34 AM, Abhi Acharya wrote:
>>> On Mon, Dec 29, 2014 at 8:59 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>>>> On 12/29/14 6:57 AM, Abhi Acharya wrote:
>>>>  Hello GROMACS Users,
>>>>> This is a problem I am facing for the first time. Kindly guide be to
>>>>> the
>>>>> best options.
>>>>> I have a protein which has two large domains connected by a flexible
>>>>> linker
>>>>> peptide (~10 aa). The two domains seem to interact with each other and
>>>>> have
>>>>> been crystallized in three different conformations. I want to calculate
>>>>> the
>>>>> change in binding energy of the two domains wrt change in their
>>>>> relative
>>>>> position, i.e. keeping position of one domain constant what is change
>>>>> in
>>>>> binding energy as the other domain moves from conformation 1, through
>>>>> conformation 2 to finally, conformation 3. What is the best way to do
>>>>> so?
>>>>>  You need to describe what these three conformations are.  Do they
>>>> involve
>>>> rotations of the domains with respect to one another, or is it a simple
>>>> linear distance that varies?
>>> Yes, there are significant rotations+translations along different axes
>>> involved which makes it challenging. Is there a way to model such
>>> movements
>>> using the pull code?
>> Perhaps the enforced rotation options can be useful here.
> I had a brief look at these options; seem to be appropriate for my
> problem. I will try them first
>>>> At first, I considered using umbrella sampling to address the problem.
>>>> Here
>>>>> too, I was not sure how to write a pull code for such a complex
>>>>> reaction
>>>>> coordinate. Also, even if I can generate the conformations, what would
>>>>> be
>>>>> the ideal number of windows that would be needed to correctly generate
>>>>> a
>>>>> PMF ? My feeling is it would be a large number, but again it just a
>>>>> guess.
>>>>> Is there is a better way than this?
>>>>>  Hard to know a priori.  You may need some trial and error here, but
>>>> with
>>>> umbrella sampling it's trivial to just go back and add windows, since
>>>> the
>>>> simulations are still all independent of one another.
>>> My concern is that since the domain movement from initial to final
>>> conformation is of about 10 nm, the number of windows required maybe too
>>> large. I read that the windows need to be spaced close enough so that
>>> each
>>> one samples some part of the next window. My system is also very large
>>> (200K atoms), so it seems to be computationally very expensive.
>> The number of windows needed is largely a function of the force constant,
>> e.g. how strong the biasing potential is.  But that's also a function of
>> how strong the interactions are that are intrinsic to the system.
>>  I was just now thinking to reduce the problem to  conducting umbrella
>>> sampling for each conformation, simply using COM distance of the two
>>> domains as the reaction coordinate to obtain the PMF in each case. From
>>> the
>>> PMF, the binding energy of the domains in each conformation can be
>>> obtain.
>>> This will allow me to circumvent the complications introduced by domain
>>> rotations and maybe reduce the number of simulations required a bit .
>>> Will
>>> this be a correct approach?
>> Possibly, but COM distance can be a degenerate measure in the case of
>> rotation.  The PLUMED plugin may be useful here rather than standard
>> Gromacs options, otherwise try to find an example in the literature and
>> follow a known procedure.
> I was suggested the same by someone else too, just thought to try options
> in gromacs before I try something unknown to me. I would try PLUMED in case
> the enforced rotations don't work.
> Thank you for the suggestions, always appreciated.
> Regards,
> Abhishek

More information about the gromacs.org_gmx-users mailing list