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

[Rini Gupta]

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


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?


Thanks and Regards,
Dr. Rini Gupta
Post Doctoral Fellow
UBC, Vancouver


PS. Full .mdp file is attached here.