Hello everyone. Following the tutorial from Prof Lemkul, I performed a SMD of a protein-ligand system followed by Umbrella Sampling. The simulations ran correctly, and I followed with the WHAM analysis.
As you can see, there is a spike just before the 3nm, corresponding with the lack of a window sampling the 2.9nm mark of the simulation. So, I checked the distances file I used to select the windows, and ran the window number 80, where the distance from the COM of the ligand and the protein was at 2.933nm.
The issue I have is that after running that specific window, a new histogram appears, but not at the 2.9nm mark but rather just before it. I went to check the pullx.xvg file from window 80 and noticed that the position is nowhere near the 2.933 nm that the summary_distances.dat file says the ligand should be at that window:
Am I doing something wrong? Did I misunderstand something? I used the same script given in the tutorial to generate the summary_distances.dat file, only replacing the names of the groups for the ones relevant in my simulation.
When your are using SMD to generate starting configurations you are pulling one group from the other with a “spring”. If the spring is not stiff enough (depending on the free energy surface) the actual distance will not be what you expect. The actual distance is what you see in the pullx.xvg file - from the SMD simulation. If that is nowhere near the reference distance (the distance you are trying to keep using a spring-like potential) then your pull force constant is probably too low or you are pulling too fast.
The tutorials are meant to work as they are. You cannot necessarily use the same settings with other systems or force fields.
Thank you for your help! I did adjust the spring constant according to literature for my system, but left the pull rate the same as the tutorial. Since at first I just wanted the rupture force, I tested several slower rates and saw minimal changes, thus I used the faster rate. I will run the SMD again with a slower rate and try again.
Once again, thank you for your time and help!
Pato
If you want to be really sure that the pulled group is where you expect it to be (according to the pulling reference position) you can set pull-coord1-type = constraint during the SMD (not during the umbrella simulations). Then you don’t need to set a force constant. But this can definitely cause problems if you pull too quickly.
