I am planning to use umbrella sampling to calculate the energy profile of a molecule’s permeation through an interface surfactant membrane. However, I’m uncertain about how to establish an appropriate pulling pathway to obtain a reliable energy profile.
This is my box (with water, surfactant-cyan, small molecule-green):
When I run the initial pulling simulation, I pulled the molecule COM upwards using
pull_coord1_type = umbrella
pull_coord1_geometry = direction-periodic
pull_coord1_groups = 0 1
pull_coord1_dim = N N Y
pull_coord1_vec = 0.0 0.0 1.0
However, I observed that the molecule’s movement in the x/y dimensions was not constrained, resulting in a trajectory where the molecule wound its way upward. I believe this is not the only pathway for the molecule to exit the solution.
If I use pull_coord1_geometry = cylinder or restrict the x/y dimensions of the COM, this would confine the molecule to exit the solution from a specific region of the interface. Given that the surfactants do not align uniformly, I suspect that a single PMF simulation may not provide an accurate energy profile for the permeation process.
In my experience with these systems, boxes with membranes tend to be treated as homogeneous along a plane (usually xy) and heterogeneous along an axis (usually z). In this way you can use z as a collective variable and sample some quantity with respect to it, like in your case the free energy would be with respect to the normal of the membrane (z, usually). Implicitly, this means that the position on xy along which you permeate is not relevant, that is, the free energy profile is independent of WHERE you go through in the bilayer, it just depends on the depth. As such, in my understanding it is to an extent mandatory to diffuse laterally as being independent from the position is clearly an approximation, and in principle you would want to sample decently the whole xy plane for a given z to have a good statistical description of that specific slab. Since lateral diffusion can be slow, you can actually combine several profiles with molecules at the same depth but different xy positions to sample better and faster.
The effect might be more pronounced especially for large inhomogeneous membranes, where you can have local distribution of lipids with different local permeabilities. Consider also that you are neglecting the inner freedom of the solute, as you are treating it as a point localized at a certain depth. The larger and the more flexible the solute, the worse the approximation becomes, as other solute-related degrees of freedom become relevant (e.g. orientation).
Ideally when running your umbrella sampling simulations, whilst restraining your molecules [z] movements, the molecule should diffuse freely in the x and y plane, you do not need to define the pathway through the membrane. Your umbrella sampling simulations should be long enough for your molecule to explore the [xy] plane at each sampled [z] level.
Depending on the density of your surfactant membrane, exploring the [xy] plane can be very slow and if your small molecule is large enough, your pulling simulation can cause disruption to your membrane that can take a long while to relax. You may want to look at an alternative way to place your molecule in the membrane at the different [z] levels, a series of short free energy perturbation simulations can be used to make your small molecule ‘appear’ in the membrane without the disruption that pulling may cause. If you find this to be the case, I’m happy to advise further.
That’s true, fully explore the xy plane will take a long time. I’m much interested in the alternative method you metioned that could make molecule appear in the membrane, could you tell me more about that? (like how to place the molecule at a specific region without cause any conflicts?)
Thanks for your helpful reply! Your comments have clarified my question.
I don’t quite understand this claim. In the simulation, the solute is not treated as a rigid object. Could you tell me why the inner freedom is neglected? Is this a result of using the COM as the collective variable?
Sorry I didn’t express properly what I mean. You are considering only the position of the center of mass of your solute. However, think about a complex solute, like a long flexible molecule that is polar (like a lipid, for example). Then, the orientation of the solute itself will be very relevant for the partition/permeation properties, and all these degrees of freedom should be properly explored if you want a well converged free energy profile. The issue is explored in literature, like here for example.
I don’t understand you argumentation for not using geometry=cylinder. I would say your arguments are exactly the reason to use geometry=cylinder. What matters is precisely where your molecule resides with respect to the local center of the membrane.
Thanks, Hess. I’m conducting this simulation to understand the global energy profile of a molecule passing through the surfactant membrane, rather than focusing on the energy profile of a specific pathway (is this achievable?). I’m uncertain whether using geometry=cylinder can effectively capture global information and provide a comprehensive understanding within just one simulation.
Yes, cylinder does that. A localized COM is required with a deformable membrane. Note that the cylinder is not fixed in space, it moves in the x-y plane along the with molecule.
If I remove the definition of pull-group1-pbcatom and pull-group2-pbcatom, there would be a error occur:
ERROR 1 [file pull.mdp]:
When the maximum distance from a pull group reference atom to other atoms
in the group is larger than 0.5 times half the box size a centrally
placed atom should be chosen as pbcatom. Pull group 2 is larger than that
and does not have a specific atom selected as reference atom.
I also discovered statements in the mdp document that mentioned,
pull-group1-pbcatom: This parameter is not used with pull-coord1-geometry cylinder.
So I’m confused and don’t know how to fix the error. Could you please let me know how to set mdp parameters correctly to run the pull code?
p.s. I didn’t find much detailed tutorial about pull code (particularly with geometry=cylinder) online, I’ll appreciate if you could share any helpful resources.
Hi hess, thanks for your previous guide on pull=cylinder, now I’m facing another problem:
If I’d like to calculate the free energy profile of such a process, e.g., a small molecule transiting across the membrane into and out of a liquid phase.
However, when the pull starts, the simulation collapsed and showed a lot of LINCS warning:
Step 405, time 0.81 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 1588.394653, max 132962.640625 (between atoms 18299 and 18300)
bonds that rotated more than 30 degrees:
In another case, if I put the small molecule close to the top membrane, the pull simulation could run normally.
I’ll appreciate if you can tell me why the first simulation does not work and how to set the parameters appropriately. Thank you.
Are the mdp options above all the pull options you provide? Does the simulation run normally when you start from the same state without pulling active?
The only suggestion I have it to include both membranes in pull group 1 and maybe also the liquid. That should improve your results. You may be lucky in that that also solves your crash.
