# Umbrella sampling PMFs with additional restraints **URL:** \ **Category:** User discussions\ **Tags:** mdp-parameters, simulation-setup\ **Created:** [October 20, 2025, 8:29am UTC](https://gromacs.bioexcel.eu/t/umbrella-sampling-pmfs-with-additional-restraints/12815 "2025-10-20T08:29:55Z")\ **Posts on this page:** 1\ **Page:** 1
**Author:** ![hzhekova](https://avatars.discourse-cdn.com/v4/letter/h/58f4c7/32.png) [@hzhekova](https://gromacs.bioexcel.eu/u/hzhekova)\ **Post date:** [October 20, 2025, 8:29am UTC](https://gromacs.bioexcel.eu/t/umbrella-sampling-pmfs-with-additional-restraints/12815/1 "2025-10-20T08:29:55Z")
GROMACS version: 2022.6 for GPU GROMACS modification: No Hello! Apologies in advance for what will be a long post, but I have some questions that I couldn’t find clear answers to by myself and I was advised to post them here. I’m trying to run umbrella sampling (US) simulations in which I move together groups of counterions (1 anion, with either 2 or 1 Na) along a reaction coordinate using the pulling module. My system is a membrane protein (embedded in POPC bilayer, with 0.15 NaCl simulation solution) with a wide vestibule in which the ions are well hydrated and have a lot of space to move around. I have assessed the anion pathway in this wide vestibule with previous SILCS and MD simulations and it’s pretty much a straight line upwards that however is slightly tilted and involves changes along all three axes (x, y, and z) as opposed to the usual change only along z during vertical pulling in systems like the potassium ion channels. This pathway is my reaction coordinate. I generated 71 windows along this reaction coordinate by moving the group of ions along it at small increments for a total of ~3 nm upward displacement with respect to the ion binding site at the bottom of the vestibule (the protein starting coordinates are identical in all 71 windows, as there is very little protein backbone reorganization during the ion permeation in the vestibule, as I’ve seen from previous MD simulations). For the umbrella sampling simulations, I’ve followed the tutorial from the Lemkul lab ([Umbrella Sampling](http://www.mdtutorials.com/gmx/umbrella/05_pull.html)) with some modifications to the US mdp file. This is the mdp file I’m using. title = Umbrella pulling simulation define = -DPOSRES\_B ; Run parameters integrator = md dt = 0.002 tinit = 0 nsteps = 5000000 ; 10 ns nstcomm = 10 ; Output parameters nstxout-compressed = 5000 ; every 10 ps nstenergy = 5000 ; Bond parameters constraint\_algorithm = lincs constraints = h-bonds continuation = yes ; Single-range cutoff scheme cutoff-scheme = Verlet nstlist = 20 ns\_type = grid rlist = 1.2 vdwtype = Cut-off vdw-modifier = Force-switch rvdw\_switch = 1.0 rvdw = 1.2 coulombtype = PME rcoulomb = 1.2 ; Berendsen temperature coupling is on in two groups Tcoupl = Nose-Hoover tc\_grps = Protein Non-Protein tau\_t = 1.0 1.0 ref\_t = 310.15 310.15 ; Pressure coupling is on pcoupl = C-rescale pcoupltype = semiisotropic tau\_p = 5.0 compressibility = 4.5e-5 4.5e-5 ref\_p = 1.0 1.0 refcoord-scaling = all ; Pull code pull = yes pull\_ncoords = 1 ; only one reaction coordinate pull\_ngroups = 2 ; two groups defining one reaction coordinate pull\_group1\_name = Chain\_A ; 2Na-anion or 1Na-anion group pull\_group2\_name = Chain\_B ; a group of protein atoms defining the ion binding site at the bottom of the vestibule pull\_coord1\_type = umbrella ; harmonic potential pull\_coord1\_geometry = distance ; simple distance increase pull\_coord1\_dim = Y Y Y ; changes along the x, y, and z axes pull\_coord1\_groups = 1 2 pull\_coord1\_start = yes ; define initial COM distance \> 0 pull\_coord1\_rate = 0.0 ; restrain in place pull\_coord1\_k = 1000 ; kJ mol^-1 nm^-2 The calculations run, but I noticed that if I don’t restrain the Na ions from the ion group that I’m pulling, they eventually dissociate from the anion and exit the vestibule. Also, occasionally, Na and Cl ions from the external NaCl solution enter the vestibule and interact with the anion/Na from the ion group which remain in the vestibule. Without any additional restrains (other than the US pulling settings) I get very negative PMF energies and weird coordination space sampling issues, as you can see on the Figure below (on the left, the black line). In the Figure, on the right you can see zoomed in the same graph as on the left, but within -10 - 0 kcal/mol : ![PMF](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/0/016e0a67f52e35ce948b340a4d584c9b66b45d84.png) To avoid this, I introduced distance restrains for the Na ions in the anion-Na group (see the red and green lines in the Figure above, right) I’m pulling by adding this text to my topology file (the example is for 1Na-anion): ; restraints for distances [intermolecular\_interactions] [bonds] ; ai aj type b0 kb 7924 7928 6 0.33 500 This seems to be doing what it is supposed to, as I don’t see dissociation of the Na ions anymore and the Na-anion group remains together in all windows. I also tried to prevent Cl from coming in the vestibule by applying an inverse flat-bottom potential (see the blue line in the Figure above, right) by adding this to the Cl itp file: [position\_restraints] 1 2 5 -2.7 100 This seems to work ok as well, as I don’t see Cl ions coming in the vestibule in my trajectories and the plotted Cl density map does not show Cl density in the vestibule. I’m aware that some parts of the reaction space haven’t been sampled enough (see the histograms below) and I’ll try to fix this in the future: ![Histo](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/2/28a4d5251d731771f9d55b12e570a73573241c37.jpeg) My questions are as follows: 1. I imagine that the additional restrains I add contribute to the PMF and have to be removed somehow. Am I correct? If yes, any ideas how to subtract the energetic contributions from the PMF would be very appreciated. In some groups I have 2Na, in others - only 1. So, it’s not like a systematic error that would cancel if I look at energy differences (e.g. how much deeper the first pronounced energy minimum for the 2Na-anion case is compared to the 1Na-anion case). So, how do I subtract from the PMF the contribution from 1 or 2 Na-anion restrains? Also, any ideas how to remove the contribution from the Cl flat bottom potential? There it’s probably a systematic type of error that will cancel in the different simulations if I apply the flat bottom potential in the same way and I look at energy differences instead of absolute energies. 2. Can we over-sample in the US calculations? I understand that under-sampling is problematic, but with my 71 windows I have some areas of the reaction coordinate where multiple histograms overlap (see the histogram graphs above). Should I remove some of the overlapping windows? 3. I’ve updated Nose-Hoover and C-rescale, since I’ve used them before and they don’t seem to be leading to anything weird in standard MD simulations. But is this combination ok for PMF with additional position restrains? Thank you for your patience with this long post and for any suggestions! Let me know if you need any further information.