# Umbrella sampling of drug over whole membrane **URL:** \ **Category:** User discussions\ **Tags:** umbrella-sampling\ **Created:** [December 24, 2023, 4:47pm UTC](https://gromacs.bioexcel.eu/t/umbrella-sampling-of-drug-over-whole-membrane/7915 "2023-12-24T16:47:36Z")\ **Posts on this page:** 7\ **Page:** 1
**Author:** ![daniel\_tims](https://avatars.discourse-cdn.com/v4/letter/d/e9c0ed/32.png) [@daniel\_tims](https://gromacs.bioexcel.eu/u/daniel_tims)\ **Post date:** [December 24, 2023, 4:47pm UTC](https://gromacs.bioexcel.eu/t/umbrella-sampling-of-drug-over-whole-membrane/7915/1 "2023-12-24T16:47:36Z")
GROMACS version: 2022.2 GROMACS modification: Yes/No Here post your question Hi, I would like to get a PMF of drug translocation over the bilayer. To do so I have performed Umbrella sampling simulations of a drug over the entire bilayer. First I did steered MD to collect the initial coordinates of my drug in the membrane, and afterward performed MD for all US windows. I have followed the procedure described in several papers and set a reaction coordinate as COM of drug and membrane. When I construct my PMF for each leaflet with wham method it looks fine. However, when I try to reconstruct that over the whole bilayer (both leaflets), it is messed up. I suspect that the reaction coordinate is set as a distance between COM of my groups, and when I try to join both leaflets that distance overlaps. Can anyone suggest how I could create one profile considering both membrane sides? Is it actually a problem with defining reaction coordinate or maybe the assembly itself in WHAM? Any hints are appreciated. Here is my one-sided PMF. ![pmf](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/c/c6379a2967fdbb7db8c3b1a02ee753c2f797e7e2.png) My pull code: ;Pull code pull = yes pull\_ncoords = 1 pull\_ngroups = 2 pull\_group2\_name = drug pull\_group1\_name = POPC pull\_coord1\_type = umbrella pull\_coord1\_geometry = distance pull\_coord1\_dim = N N Y pull-pbc-ref-prev-step-com=yes pull\_coord1\_groups = 1 2 pull\_coord1\_start = yes pull\_coord1\_rate = 0.0 pull\_coord1\_k = 1500 ---
**Author:** ![milosz.wieczor](https://avatars.discourse-cdn.com/v4/letter/m/b19c9b/32.png) [@milosz.wieczor](https://gromacs.bioexcel.eu/u/milosz.wieczor)\ **Post date:** [December 28, 2023, 6:26pm UTC](https://gromacs.bioexcel.eu/t/umbrella-sampling-of-drug-over-whole-membrane/7915/2 "2023-12-28T18:26:38Z")
Is the membrane symmetric? If so, then it makes no physical difference on which side of the membrane your drug is, so you can just duplicate the raw PMF and add a minus sign to one copy’s X coordinate. If you want to e.g. assess convergence (which is always a good idea), you can run WHAM twice, once per every half of the dataset, and again reverse the sign in one of them. For future free energy calculations, you might want to use `direction` instead of `distance`, and supply a vector pointing in the Z direction (`0 0 1`). ---
**Author:** ![daniel\_tims](https://avatars.discourse-cdn.com/v4/letter/d/e9c0ed/32.png) [@daniel\_tims](https://gromacs.bioexcel.eu/u/daniel_tims)\ **Post date:** [December 29, 2023, 9:46am UTC](https://gromacs.bioexcel.eu/t/umbrella-sampling-of-drug-over-whole-membrane/7915/3 "2023-12-29T09:46:26Z")
Thanks a lot. I tried distance but the particle moved too fast. Now I adjusted the coord rate and it looks fine. The membrane is slightly asymmetric, so I will try to connect PMF of each leaflet. One more thing, is it fine to use my previous pull files (pullx and pullf from distance geometry) and modify only the .tpr files (with direction) for the wham analysis of both leaflets? Maybe I should additionally -rerun the trajectories with new tpr? Or rather these calculations cannot be reliable anymore. Thanks. ---
**Author:** ![milosz.wieczor](https://avatars.discourse-cdn.com/v4/letter/m/b19c9b/32.png) [@milosz.wieczor](https://gromacs.bioexcel.eu/u/milosz.wieczor)\ **Post date:** [December 29, 2023, 2:39pm UTC](https://gromacs.bioexcel.eu/t/umbrella-sampling-of-drug-over-whole-membrane/7915/4 "2023-12-29T14:39:22Z")
If the membrane is asymmetric, then you can still do this in one shot - as you say, rerun the trajectories with a .tpr that will give you signed values, and perform reweighting according to the original US weights. (In fact, reweighting is so general that you can perform it with any new collective variable/CV as long as it was sufficiently sampled in the original set of simulations.) If you ran your US using the distance CV \xi with an external potential in i -th window V\_i(\xi) = \frac12k\_i (\xi - \xi\_i^0)^2 and obtained an original PMF F(\xi), then every frame you got has an (unnormalized) statistical weight w = \exp(\beta(V(\xi) - F(\xi))). What you can then do is calculate this per-frame weight, bin the normalized weights along your new CV, and voila - you have the probability density along the new CV. (This procedure can be made binless with techniques like KDE, but that’s a technicality.) Then you do standard Boltzmann inversion, i.e. F(\xi') = -k\_BT \ln(p(\xi')), and you have your new corrected PMF. That’s the fully formal version, and it probably takes 10-15 lines of Python code to implement. Splitting the inputs based on the rerun will be an approximation, perhaps good enough if you don’t need a very precise description of the barrier region. ---
**Author:** ![daniel\_tims](https://avatars.discourse-cdn.com/v4/letter/d/e9c0ed/32.png) [@daniel\_tims](https://gromacs.bioexcel.eu/u/daniel_tims)\ **Post date:** [December 29, 2023, 11:00pm UTC](https://gromacs.bioexcel.eu/t/umbrella-sampling-of-drug-over-whole-membrane/7915/5 "2023-12-29T23:00:27Z")
Milosz, thanks a lot for your time and help. You’re right, I’ll write a script to comprehensively reweight all windows. ---
**Author:** ![milosz.wieczor](https://avatars.discourse-cdn.com/v4/letter/m/b19c9b/32.png) [@milosz.wieczor](https://gromacs.bioexcel.eu/u/milosz.wieczor)\ **Post date:** [December 30, 2023, 12:01pm UTC](https://gromacs.bioexcel.eu/t/umbrella-sampling-of-drug-over-whole-membrane/7915/6 "2023-12-30T12:01:43Z")
You’re welcome! There’s a part of it that always feels like reinventing the wheel but then you get to learn the most from doing it on your own. ---
**Author:** ![EOrfanakis](https://avatars.discourse-cdn.com/v4/letter/e/7cd45c/32.png) [@EOrfanakis](https://gromacs.bioexcel.eu/u/EOrfanakis)\ **Post date:** [March 26, 2025, 4:42pm UTC](https://gromacs.bioexcel.eu/t/umbrella-sampling-of-drug-over-whole-membrane/7915/7 "2025-03-26T16:42:06Z")
Why your initial energy has a value different from zero? PMF diagrams are about relative energy differences, did you do it on purpose?