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) 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 :
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:
My questions are as follows:
- 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.
- 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?
- 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.