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:
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.
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.
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.
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!
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?
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.
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)?
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?
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.
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?)
