Umbrella sampling of Lipid Protein system using Martini FF

GROMACS version: 2020.4
GROMACS modification: No

I am calculating potential of mean force (PMF) for a protein-lipid bilayer system using Coarse Grained Martini force field and Gromacs 2020.4. In steered MD simulation I applied force to get a reaction coordinates in z-direction, where protein moved closer to the lipid bilayer and reached to the bilayer center (using pull code of Gromacs).
I generated a series of initial configurations, each corresponding to a location where the protein was harmonically restrained at decreasing center-of­-mass (COM) distance from the lipid molecule (from 3.9 nm to 0.5 nm) using an umbrella biasing potential following Dr. Justin A. Lemkul tutorial.
For Umbrella sampling (US), I generated a total of 37 different starting configurations along the z-reaction coordinate. Each configurations were independently equilibrated followed by the US for 10ns using following pull code:

; Pull code
pull = yes
pull_ncoords = 1 ; only one reaction coordinate
pull_ngroups = 2 ; two groups defining one reaction coordinate
pull_group1_name = PROTEIN
pull_group2_name = MEMBRANE
pull_coord1_type = umbrella ; harmonic potential
pull_coord1_geometry = distance ; simple distance increase
pull_coord1_dim = N N Y ; pull along z
pull_coord1_groups = 1 2 ; groups 1 and 2 define the reaction coordinate
pull_coord1_start = yes ; define initial COM distance > 0
pull_coord1_rate = 0.0 ; 0.001 nm per ps = 1 nm per ns
pull_coord1_k = 1000 ; kJ mol^-1 nm^-2
pull_coord1_init = 0.2

pull_group1_pbcatom = 20
pull-pbc-ref-prev-step-com = yes

This is generating 37 trajectories, each of 10 ns. To use in WHAM I am using the last 5ns from each trajectory. Everything seems to be running fine, however, I am getting histogram, which seems not to be correct (see attachment please), as each window is not overlapping.
I have three questions:
Q.1 Why my COM distance is increased form 0.5 nm to over 1nm during production run despite using pulling force?
Q.2 Why windows in my histogram is not overlapping? If I use small force constant, the COM starts increasing.
Q.3 How can I resolve it?

Many thanks for taking time and reading the message. Highly appreciated. Looking forward for the solution.


The lack of overlapping between each window means that you did not reach convergence. A option is to expend the run and check for convergence.

Concerning your second issue, according to mdp file you define the distance as COM of starting conformation at t=0 + 0.2 nm)

pull_coord1_start = yes => add the COM distance of the starting conformation to pull-coord1-init
pull_coord1_init = 0.2

I am not sure that it is what you want. Using this setting, to have to sure what is the distance at t=0. One possible alternative is to use pull_coord1_start = no and define pull_coord1_init = selected distance.

Kind regards