Awh mdp options

GROMACS version: 2021.5
GROMACS modification: No


I am trying the seemingly simple task of pulling a molecule through a bilayer using the pull-code + awh, but I am confused about some of the options, and how they influence each other.

The initial setup I have is the molecule outside of the bilayer, and “above” the bilayer with reference to the z-axis. My aim is to pull it to the middle of the bilayer, a distance of approx. 5 nm. I am using the cylinder pull geometry, here is my mdp-file:

pull-group1-name = bilayer
pull-group2-name = MOL

pull-coord1-groups = 1 2
pull-coord1-geometry = cylinder
pull-coord1-vec = 0.0 0.0 -1.0
pull-cylinder-r = 1.0
pull-coord1-type = external-potential
pull-coord1-potential-provider = AWH

awh = yes ; AWH on.
awh-nstout = 500000
awh-nbias = 1
awh1-ndim = 1
awh1-dim1-coord-index = 1
awh1-dim1-start = 7
awh1-dim1-end = 2
awh1-dim1-force-constant = 100000
awh1-dim1-diffusion = 5e-3

With this setup I was intending to enable a simulation where the molecule starts from its initial position (at z=7, pull for 5 nm (until z=2), so along the vector 0, 0, -1.) However, grompp tells me that the distance at start (I presume between the molecule and the bilayer) is -1.627 nm, and the awh gives a fatal error saying that: For the non-periodic pull coordinates awh1-dim1-start (7.000000) cannot be larger than awh1-dim1-end (2.000000).

Please advise how I need to modify the mdp-file to acheive what I want.

Kind regards

Anyone, please?

Why do you used 0,0,-1 as pull vector and not 0,0,1 ? If the initial value of the pull coordinate is -1.627, it seems to me that your molecule is above the bilayer in the starting configuration and not below.

But if you change to 0,0,1 as vector, you starting pull value is still outside of the interval. So either you need to run a short steered MD run to move your molecule to 2 or beyond, or you need to change the interval in the AWH parameters.

Thanks Berk, much appreciated!

The molecule is indeed above the bilayer in the starting configuration, and my reasoning for using 0,0,-1 was to pull it in the negative direction along z.
But I am clearly not understanding 100% how this works. When I run grompp, the distance that is reported changes sign (from -1.627 to 1.627) when the pull vector instead is 0,0,1. This does not make sense to me, why would the initial distance be affected by this vector?

Sedondly, how should the awh interval be specified for pulling from above to below the bilayer?


If you use -1 as z-component of the vector, then a distance of 1.627 above the bilayer changes to a location of -1.627.

The pull coordinate value is given as the z-position with as reference, or zero-, point the middle of the bilayer. If you want to pull it through, you can set the interval to e.g. -1.8 to 1.8.

Thanks again, I think I understand now. Will do some runs to see if the outcome is as intended.

So, still something I am not clear about, now about the values of awh-dim1-start and -end. Assuming a pull-vec of 0,0,-1, and (therefore) an initial position of -1.627, I put
awh-dim1-start = -1.8 and awh-dim1-end to 3.2 (I want to pull about 5 nm in total, approximately to the middle of the bilayer). Grompp works fine, but after about 100000 steps I get this from gromacs:

Simulation instability detected:
Coordinate 1 of an AWH bias has a value -1.872609 which is more than 10 sigma
out of the AWH range of [-1.800000, 3.200000].

How should the start and end awh interval be specificied in this case?
I really appreciate your help!

Hello again!

Unfortunately still struggling with this, any input would be appreciated. How should the values for awh-dim1-start and awh-dim1-end be picked in relation to the intial position of my molecule above the lipid bilayer, with the intention of using awh to sample passage across the bilayer?

Your start and end look reasonable to me. I don’t see how the molecule can end up so much beyond the start, unless you system is unstable.

Are you still using a force constant of 100000? What is the mass of your molecule?

Thanks again Berk for your reply! The mass of my molecule is ~1000 Da, and the force constant is 10000.

Over the weekend I was working more with this, testing also different pulling geometries. Before it was always cylinder. I have now also tried distance and direction, and I have additional questions as a consequence of that.

  1. In cases like this (pulling from one side of a bilayer to the other), what should be group1 and group2 in the pull-code? I was always putting the bilayer as group 1 and the permeant as group2, but does it matter (other than the sign of the distance being reversed)?

  2. With pull-geometry=direction, (and the bilayer as group2) the distance as reported by grompp is -3.069 nm, but with pull-geometry = cylinder is is the ~1.6 nm value from above. This makes sense since the COM of the bilayer is different for these two cases, but given the different values, how should one think when specifiying the awh boundaries to make sure the entire passage is sampled? I suppose I am trying to understand still if the distances and positions that grompp prints and uses are directly comparable or not to the awh start and end positions?

  3. With pull-geometry = direction, and the awh start and end put to -3.1 and 3.1, I obtained the following PMF:

When I visualise the trajectory in VMD, it seems as if my molecule is pulled to about the center of the bilayer, possibly a little bit further. How does this PMF match with that trajectory, since the PMF is flat from about -1 and onwards. I was expecting a more “classical” PMF that was changing as a function of the reaction coordinate inside the bilayer?

Apologies for the ton of questions, but I would very much like to get going with AWH for permeation studies.


With geometry=cylinder, the cylinder group should be group 1.

The PMF you show is the result of a simulation where the reaction coordinate never reached the flat part. You would need to run longer.

But cylinder is the right choice. I don’t understand what happens with the distance ending up outside the interval. You did not set any more pull and AWH parameters than what you show in your first post, right?

Since I have no clue what is going wrong, could you change the pull direction to 0,0,1 and exchange the start and end values and their sign? I don’t think this will make a difference, but maybe there is a bug somewhere.

This is the full pull and awh parts of my mdp-file.

pull-group1-name = bilayer
pull-group2-name = MOL

pull-coord1-groups = 1 2
pull-coord1-geometry = cylinder
pull-coord1-vec = 0.0 0.0 1.0
pull-group1-pbcatom = 20100
pull-pbc-ref-prev-step-com = yes
pull-coord1-k = 1000
pull-coord1-start = yes
pull-cylinder-r = 1.5
pull-coord1-type = external-potential
pull-coord1-potential-provider = AWH

awh = yes
awh-nstout = 50000
awh-nbias = 1
awh1-ndim = 1
awh1-dim1-coord-index = 1
awh1-dim1-start = 3.1
awh1-dim1-end = -3.1
awh1-dim1-force-constant = 10000
awh1-dim1-diffusion = 5e-5
awh1-dim1-coord-provider = pull

Here is the output that I get with these settings:

Pull group natoms pbc atom distance at start reference at t=0
1 40200 20100
2 131 40266 1.661 nm 1.661 nm

Program: gmx grompp, version 2021.5
Source file: src/gromacs/applied_forces/awh/read_params.cpp (line 994)

Fatal error:
For the non-periodic pull coordinates awh1-dim1-start (3.100000) cannot be
larger than awh1-dim1-end (-3.100000)

For more information and tips for troubleshooting, please check the GROMACS
website at Common Errors — GROMACS webpage documentation

I also attach a screen-shot of the initial position of my molecule and the bilayer (lipid headgroups as spheres):

You need to swap start and end, or swap their signs, which is equivalent in this case.

You don’t need the pull-pbc-ref-prev-step-com = yes setting. I don’t think this causes issues, but it might be good to try without.

Thanks, I will try and report back!
Before starting with this, is there some setting that I can change/add to the mdp-file that allows me to quickly assess if the simulation is exploring the entire awh-range as intended? (I am thinking similarly to specifying a high pull-rate when not combining with awh?)


If you run gmx awh -more you get the coordinate distribution.