[gmx-users] Concrete pull code explanation needed

[Justin Lemkul]

On 7/13/17 8:11 PM, Du, Yu wrote:
> 
> 
> 
>> -----Original Messages-----
>> From: "Justin Lemkul" <jalemkul at vt.edu>
>> Sent Time: Thursday, July 13, 2017
>> To: gmx-users at gromacs.org
>> Cc:
>> Subject: Re: [gmx-users] Concrete pull code explanation needed
>>
>>
>>
>> On 7/11/17 8:23 PM, Du, Yu wrote:
>>>> On 7/10/17 11:19 PM, Du, Yu wrote:
>>>>> Dear Justin and gmx users,
>>>>>
>>>>>
>>>>> I have gone through mdp-option and Justin A. Lemkul's COM pulling tutorial serveral times.
>>>>>
>>>>>
>>>>> The following is Justin's pull code.
>>>>>
>>>>>
>>>>> ; Pull code
>>>>> pull                    = yes
>>>>> pull_ngroups            = 2
>>>>> pull_ncoords            = 1
>>>>> pull_group1_name        = Chain_B
>>>>> pull_group2_name        = Chain_A
>>>>> pull_coord1_type        = umbrella      ; harmonic biasing force
>>>>> pull_coord1_geometry    = distance      ; simple distance increase
>>>>> pull_coord1_groups= 1 2
>>>>> pull_coord1_dim         = N N Y
>>>>> pull_coord1_rate        = 0.01          ; 0.01 nm per ps = 10 nm per ns
>>>>> pull_coord1_k           = 1000          ; kJ mol^-1 nm^-2
>>>>> pull_coord1_start       = yes           ; define initial COM distance > 0
>>>>>
>>>>>
>>>>> My understanding lists as follows,
>>>>>
>>>>>
>>>>> Justin defines two pull groups, with `pull-ngroups = 2`, each of them has a name in the index.ndx generated by `gmx make_ndx -f npt.gro`and their names are defined by `pull_group1_name = Chain_B` and `pull_group2_name = Chain_A`.
>>>>>
>>>>>
>>>>> My question is about the definition of pulling coordination and the orientation of pulling force.
>>>>>
>>>>>
>>>>> 1) I learnt from [gmx-users] Change to umbrella sampling pull code，
>>>>> "You need: pull-coord1-groups = 1 2 otherwise the reaction coordinate is undefined, or otherwise defaults to the entire system, I can't remember which. -Justin"
>>>>>
>>>>>
>>>>> I know it's a plot :) in the input.pdb that proteins are placed exquisitely along the z-axis which is the same as the pulling coordinate but it makes pull code confused and here I need a concrete explanation.
>>>>> 1.The pull coordinate is the line that connects COM of group1 and group2 with `pull_coord1_groups= 1 2`.
>>>>> OR 2. The pull coordinate is the z axis with `pull_coord1_dim = N N Y`.
>>>>> Which is correct?
>>>>>
>>>>
>>>> The z-component of the vector connecting the COMs of the two groups.
>>>>
>>>>>
>>>>> 2)Then turn to the orientation of pulling forces.
>>>>> My understanding is that force1 acts on pull_group1, force2 acts on pull_group2 and the orientation of force 1 and 2 is opposite, both forces have a rate of 10nm per ns.
>>>>> Is my understanding and the below schematic draw right?
>>>>>
>>>>
>>>> There is one force.  It acts on the spring connecting the two groups.
>>>
>>>
>>> How does the spring connect the two groups? Does the spring link to the COM of the whole two groups?
>>
>> Yes.
>>
>>> How does Gromacs define the orientation of the force of pulling? Or by default is the pulling force just positive along the z axis with `pull_coord1_geometry = distance` and `pull_coord1_dim = N N Y`?
>>>
>>
>> This is not the default, but it is precisely what is specified by those .mdp
>> settings.
>>
> 
> 
> Could you please show which line specifies the orientation of the pulling force? Is it the positive `pull_coord1_rate`? So `pull_coord1_rate` can be either positive or negative?
> 

pull_coord1_dim or pull_coord1_vec, depending on the geometry used.

Pull rates can be positive (to increase distance between the groups over time), 
negative (to decrease the distance), or zero (to maintain a distance).  The rate 
has nothing to do with the direction.

>>>>
>>>>>
>>>>> Z-axis-0-------------5----------->-------positive-orientation------->-------------25------>
>>>>>       <--------Force1-------pull_group1~pull_group2------Force2---------->
>>>>>
>>>>>
>>>>> The last question is about the umbrella sampling.
>>>>> I learnt from [gmx-users] Re: doubt about your Umbrella Sampling tutorial that it's ok to remove the pores of Chain B during the US. But in longer simulation time and in the periodical box, will the COM of Chain B and A be affacted by the boundary? and then affact the calculation of US potential?
>>>>>
>>>>
>>>> The tutorial system won't work if you turn off the restraint.  Eventually the
>>>> system will rotate and the groups will cross periodic boundaries, which will
>>>> cause the chosen pull geometry to fail.  So the restraints serve a dual purpose:
>>>> (1) to mimic the stability of larger amyloid assemblies and (2) to reduce the
>>>> system size, as several hundred thousand atoms was not feasible for me at the time.
>>>
>>>
>>> So your point is that during US, we can't remove the Chain B's restrain. Right?
>>
>> *For this specific case* yes.
>>
>>> During US, is there any means to study the flexiblity of both groups? (i.e. there is no restrain on both groups except the umbrella potential between them and at the same time both groups will not cross periodic boundaries and will not affact the umbrella potential geometry, which is a necessary part of my study)
>>
>> Normally one should not apply any position restraints, except in niche cases
>> like mine.
> 
> 
> During US, Why did you use Chain B's position restrain? Why is your case specific?
> 

Please read the paper I linked from the tutorial and the many, many times I've 
explained this on the mailing list.

> 
> Again, the same question is that if I don't use Chain B's position restrain, the groups can cross the box boundaries, the US may be affacted by the periodical boundary after a long time of simulation. Right?
> 

Yes, in which case you have to use a sufficiently large box to prevent this problem.

-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

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