I am simulating the energy profile of two identical flat hexagonal particles approaching each other in a parallel face-to-face way using the umbrella sampling method. Each particle has an edge length of 3 nm and a thickness of 1.7 nm.
To make them approach in a parallel way in the Z direction, I add harmonic position restraints on the heavy atoms of the upper particle(x y z 1000 1000 0). Similarly, to keep the below particle at the bottom, (x y z 1000 1000 1000) position restraints were added to its heavy atoms.
For the pushing process, I used a rate of 0.001 nm/ps and a force constant of 1000 kj/(mol nm²), and pushed for 4 nm. The pullf.xvg and pullx.xvg looks fine, but the energy profile figure I got from PMF looks very bumpy although the histogram figure looks okay except for the far left and the far right. Here shows two energy profiles with two different force constants settings. The less bumpy one uses 8000 kj/(mol nm²) (when they are far) and 40000 kj/(mol nm²) (when they are close), and the more bumpy one uses 1000 for the whole pushing process. I upload one MDP file for simulating every window.
I have two questions:
Will the position restraints added on heavy atoms affect the PMF energy profile?
Do I need to use more windows (I split 2.6 nm into about 100 windows) to smooth out the peaks that appear on the energy profile? Or it could be other reasons why the energy profiles are so bumpy?
If you need to use restraints when generating the input configurations for the umbrella sampling simulaitons, I don’t see any obvious problems. However, make sure that you do not use them during the umbrella sampling. You might need a longer equilibration stage (the -b option to gmx wham) in each umbrella sampling window if you were using restraints when generating the input configurations.
I would recommend more windows. There are regions where neighbouring histograms have very little overlap. I would suspect that those regions cause the noisy PMF.
A force constant of 40,000 kJ mol^-1 nm^-2 sounds quite high. I would think that you could manage with 8,000 or 10,000 kJ mol^-1 nm^-2 for the whole region (perhaps 5,000 kJ mol^-1 nm^-2 at > 2 nm). However, it is clear that 1,000 kJ mol^-1 nm^-2 is not enough, especially when the groups are close.
The reason why I still use restraints during the umbrella sampling for every window is that I want to make sure the two flat particles are parallel to each other and also to the xy plane, especially when they are close. If I remove restraints during the umbrella sampling, both of them will rotate along the centers of mass even though the COMs have a fixed distance. Do you know if there is a better way to keep the particle below parallel to the xy plane, and keep the upper particle parallel to another particle?
I will try to create more windows and use moderate force constants (10000 and 5000 kJ mol^-1 nm^-2). What do you mean by longer equilibration stage (the -b option to gmx wham)? Every window was simulated for 5 ns. Should I simulate longer, for example, 10 ns, and when I do the gmx wham analysis, I use -b 5ns to analyze the 5-10ns period instead of the default 0-10ns?
If what you want to measure is them exactly restrained like that, then I guess you should keep the restraints. You shouldn’t need the Z restraints in for the reference group either. The question is how to compare it to anything observable.
With longer equilibration time, I mean that if you change any parameters (including restraints) between the conformation generation and the umbrella simulations you might need to equilibrate the umbrella windows longer. In my opinion you should always use the -b option to discard the beginning of each simulation. I can’t say how much. It’s probably a bit overkill to discard as much as 5 ns out of a 10 ns simulation. The best way to check is usually to find a time block that gives a reasonably reliable estimate and then analyse using a sliding window to find the time when the simulations are reasonably stable, at least without a trend in the PMF values as you start analysing later. E.g.
If you find out that it seems stable after, e.g., 1500 ps, you can then use -b 1500 and analyse until the end of the simulation (5000 ps in your case currently, but might need to be more).
Thanks for the reply. I want a specific conformation between them because I plan to convert the shape of the energy profile to some parameters, and then transfer them to build a larger coarse-grained model for this flat particle.
I am now sliding windows as suggested by you. How do I know if I find a reasonably reliable estimate after using ‘gmx wham -it tpr-files.dat -if pullf-files.dat -b 0 -e 2500 -o profile-0-25.xvg -hist histo-0-25.xvg -unit kCal’?
I have another issue about using ‘gmx wham’. This command finished very fast when I did umbrella sampling for two small molecules (less than 100 atoms), but for these large flat particles (one particle has around 3000 atoms), it sometimes took hours. Most of the time, the terminal shows like below. Is it normal?
'Evaluating only 487 of 15200 expressions.
1) Maximum change 1.383798e+00
100) Maximum change 1.293938e+00
200) Maximum change 2.127842e-01
300) Maximum change 6.408089e-02
Unless you’ve got experimental data to compare to you won’t know if the estimate is reliable. If you wonder about the comparisons of sliding windows, then you can often plot all the resulting PMFs (with different start and end times) and see if they seem stable or if there is a trend (anywhere in the PMFs) when you increase the equilibration time. With trend I mean that, e.g., a peak or a trough or the whole PMF is changing systematically as you change the equilibration time.
The molecule size should not matter to gmx wham. But a high number of umbrellas might make it slower. I’m a bit surprised it’s that slow, though.
Hi Magnus, I did three tests in the last few days.
Sliding windows by every 5ns based on my previous umbrella sampling simulations based on your suggestion. Sorry that I forgot to set 4 nm as the zero point for a better comparison. It seems that 500-3000 or 1000-3500 has a more similar trend to the default one. What do you think of this sliding test?
Adding 39 windows between 2 and 4 nm. More windows did not smooth out the curve even though the added windows produced more overlapping in the histogram.
Changing restraints during umbrella sampling
(1) xyz restraint on heavy atoms of the lower particle and xy restraint on heavy atoms of the upper particle. The PMF curve is the same as the one in the previous post.
(2) Removing xyz restraint on heavy atoms of the lower particle and xy restraint on heavy atoms of the upper particle. Basically removing all restraints of the two particles. From the trajectory of every window, both particles tile and rotate along their COM.
(3) Only removing z restraint on heavy atoms of the lower particle. So both the lower particle and upper particle have the xy restraint.
It seems that removing restraints makes the energy profile curve smoother and the histogram distributes more evenly with fewer restraints. Combined with the finding from point 2, adding more windows did not improve the bumpy curve, but removing the restraint of z did.
Thanks for the detailed report. I’ll go through your findings step by step.
It is not surprising that the “Default” that contains all values are reasonably similar to those 1000-3500 and 1500-4000 as they are the in the middle of the series. However, you see the trend of the PMF rising with longer equilibration time and you don’t see them converging. So, I would conclude that discarding the first 2500 ps is not enough. I would recommend extending the simulations to at least 8 ns and continuing the series. What you might observe is that 2000-4500 and 2500-5000 mostly look shifted (whereas the others have different slopes) in the region > 1.6 nm. This indicates that it is the region < 1.6 nm that would require more sampling (and longer equilibration time).
I agree that the extra 39 windows did not seem to reduce the noise of the profile.
As you conclude, I would remove the Z restraints on the lower particle. But I would also carefully check the trajectories so that you sample the behaviour that you are interested in.
