[gmx-users] How to define pull-init for two pull groups

[Justin Lemkul]

On 1/13/14, 7:40 PM, Rini Gupta wrote:
> Dear gmx-users,
>
> I want to calculate Potential of Mean Force (PMF) of a lipid in a bilayer
> by umbrella sampling method using coarse grain molecular dynamics
> simulations. My simulation system consists of
> 76 (PEPC and GDPE lipids mixture), 2495 water molecules and 16 other
> molecules which are acting as a surface attached to this bilayer system.
>
>
> In order to generate starting structures for the pulling simulations, after
> minimization, I have equilibrated the system for 25 ns under NVT ensemble
> followed by 25 ns NPT run with position restraints applied to  both lipids
> and surface molecules. Following this, I have ran production simulations
> under NPT ensemble for another 25 ns. During production phase, position
> restraints were removed from all lipids except surface molecules.
>    I want to calculate PMF for pulling out two lipids simultaneously one in
> top and one in bottom layer of this system.
>   Following Justin tutorial, I understand that I need to define two pull
> groups and one reference group for my problem. For this, I have chosen all
> carbon's of lipid bilayer as reference group  and Phosphate atom of two
> lipids which are placed in top and bottom layers respectively,
> as pull group1=phosup and pull group2 =phoslow
> My problem is how to define pull-init for two pull groups:
> I have to directly place the initial coordinates of two selected lipids as
> pull-init or first I need to assume that  bilayer center is shifted to
> (0,0,0 ) and then put initial cordinates as pull-init.
> Please let me know if anything is wrong here.
> Following second approach, I have used the following grompp.mdp file for
> pulling simulation:
>
> ; Pull code
> pull            = umbrella
> pull_geometry   = position  ;
> pull_dim        = N N Y
> pull_start      = yes       ; define initial COM distance > 0
> pull_ngroups    = 2
> pull_group0     = cbilayer  ; center of carbon's of bilayer
> pull_group1     = phosup   ; phosphate atom of 16th PEPC in top
> pull_pbcatom1   = 0
> pull_vec1       = 0.0 0.0 -1.0
> pull_init1      = 0.058 -1.88226 -1.91169
> pull_rate1      = 0.01      ; 0.01 nm per ps = 10 nm per ns
> pull_k1         = 1000      ; kJ mol^-1 nm^-2
> pull_group2     = phoslow   ; phosphate atom of 56th PEPC in bottom
> pull_pbcatom2   = 0
> pull_vec2       = 0.0 0.0 1.0
> pull_init2      = 1.778 1.20774 2.20831
> pull_rate2      = 0.01      ; 0.01 nm per ps = 10 nm per ns
> pull_k2         = 1000      ; kJ mol^-1 nm^-2
>

One problem is that you are using both "pull_start = yes" and specifying 
non-zero values for pull_init[12].  What you are doing is specifying the 
reference distance as the initial COM distance + pull_init[12].  If you are 
using "pull_start = yes," set pull_init[12] to a zero vector.  If you manually 
specify the actual COM distance in pull_init[12], then set "pull_start = no." 
The trajectories indicate that you are probably getting something close to what 
you want to achieve, but likely the forces and resulting output distances are 
going to be somewhat screwy.

>
> After this, following tutorial , I extracted the separate trajectory frames
> using trjconv and then in order to measure the COM distances between two
> pull groups and a reference group
> I have created two groups.txt files and as a result I got two
> summary_distances.dat files corresponding to two pulling lipids with
> respect to a reference group.
> However, looking at the contents of both files I do not observe the any
> definite decreasing or increasing COM distances between pull group and
> reference.
>   But, on observing the snapshots,
> I found that pull groups are moving in right direction.
> If I am pulling a lipid across bilayer,
> these COM distance between pull group and reference group should first
> decrease and then increase right?
> Please tell me what I am doing wrong here.
>
>   How to take some of these configurations as starting window structures for
> next step of pulling simulations?
>

That depends on what you define as a reaction coordinate.  If you are interested 
in sampling states during which both lipids are traversing the bilayer, then you 
need configurations that are appropriate for both lipids along the reaction 
coordinate.  There's nothing really special about it, except that you likely 
have two target distances to simultaneously satisfy.

-Justin

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

Justin A. Lemkul, Ph.D.
Postdoctoral Fellow

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

jalemkul at outerbanks.umaryland.edu | (410) 706-7441

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