[gmx-users] ligand moving out during umbrella sampling

[Justin Lemkul]

On 5/15/17 2:45 AM, abhisek Mondal wrote:
> On Thu, May 11, 2017 at 9:03 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>
>>
>>
>> On 5/11/17 9:21 AM, abhisek Mondal wrote:
>>
>>> On Thu, May 11, 2017 at 11:02 AM, abhisek Mondal <abhisek.mndl at gmail.com>
>>> wrote:
>>>
>>> Hi,
>>>>
>>>> Thank you for the explanation. It really cleared some concepts. But I'm
>>>> still having my ligand moving in this step. I have modified the code as:
>>>> ; Pull code
>>>> pull                    = umbrella
>>>> pull_ngroups            = 1
>>>> pull_group0             = Protein_chain_A
>>>> pull_group1             = ACO
>>>> pull_geometry           = direction ; simple distance increase
>>>> pull_dim           = Y Y Y         ; not to allow ligand move along other
>>>> dir
>>>> pull_rate1         = 0.0
>>>> pull_k1           = 1000           ; kJ mol^-1 nm^-2
>>>> pull_start       = yes           ; define initial COM distance > 0
>>>> pull_vec1               = 0 0 -1
>>>>
>>>>
>> Note that with "direction" geometry, only pull_vec1 is acting.  pull_dim
>> is ignored.
>>
>> The ligand was previously moving along x,y direction when I was using
>>>> pull_dim  = N N Y. So I changed it to Y in all direction and provided 0
>>>> as
>>>> vector  and pull_rate1=0.0, so that it does not move much. But at the end
>>>> of a 10ns run, I see that the ligand is still moving as it was earlier.
>>>>
>>>> It shows me:
>>>  Pull group  natoms  pbc atom  distance at start     reference at t=0
>>>        0      1132    936665
>>>        1        59      1618  -1.555                -1.555
>>> Is it ok withe negative value ? Anyway this setup is not working.
>>>
>>>
>> Again you're trying to just apply the restraint to one dimension and it
>> looks to be fairly arbitrary.  I already suggested using the vector
>> connecting the ligand COM with the binding site residues' COM and using
>> that as pull_vec1.  Draw it out.  It makes a lot more sense than trying to
>> restrain only along one axis, which as I have said before, makes no sense
>> in this case.
>>
>> Thank you for such detailed suggestion.
> I followed on as per your suggestion. Calculated COM of protein and Ligand.
> Calculated protein-lig vector (using COM) to be used for pulling (as
> pull_vec1).
> Pulling also achieved successfully.
> But after pulling, when I performed the brief npt_umbrella run with
> pull_rate1=0, I found the ligand is moving little bit. Could not understand
> what I have mistaken this time.
> So I moved on for umbrella sampling mdrun step, with pull_rate1=0 and
> pull_vec1=as determined from COM calculations. Despite I found that the
> ligand is moving vigorously and got pulled away probably.
> I'm putting npt_umbrella.mdp,md_umbrella.mdp,conf140.gro(used for
> npt_umbrella), npt140.gro,md_umbrella run video in the following link:
> https://drive.google.com/drive/folders/0B6O-L5Y7BiGJQ1FIc2tIRFE2dE0
>
>
> Am I still missing some drastic steps ? Please suggest me if I do. I'm
> totally lost here in this regard.
>

Again, don't restrain the protein.  I've said this multiple times.

The pull setup looks reasonable.  Maybe you just need a stronger force constant, 
or you should not use the COM of the whole protein, instead the COM of a few 
important residues (use an index group with gmx traj -ox -com to get its 
coordinates).  If the specified vector is off, so too will be the resulting 
biasing potential.

-Justin

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

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

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