Hi
I am running umbrella sampling for a protein-organic polymer complex. After the SMD step (pull-geometry =direction), I observed that in the end frames, the ligand is going out of the PBC (My guess is due to periodicity/jump). When I start sampling, the nvt_umbrella runs fine. However, md_umbrella run stops by giving this error:
**Fatal error: Distance between pull groups 1 and 2 (7.770726 nm) is larger than 0.49 times the box size (7.770657). You might want to consider using “pull-geometry = direction-periodic” instead.
And if I change the pull-geometry to direction-periodic, it shows a fatal error but still runs:
Fatal error:
Can not have dynamic box while using pull geometry ‘direction-periodic’ (dim y)
I think that a nojump trjconv could solve this issue. However, I am not sure if this is the right step to proceed with; please let me know if this is the correct way to tackle this. If not, then what is the proper approach to tackle such an issue?
The window spacing is set to 0.1, and the pull code with pull-geometry=direction is as follows:
The gmx trjconv -nojump is a post processing unwrapping of the trajectory and as such won’t change anything here.
The problem is that once you pull too much (> half the box size) away from the reference point then the PBCs are ill defined (think about it: which is the new point closer to you? The old one or the new one given by the PBC images after moving more than half the box size? And how does this affect the pulling?). As suggested by the error, you can change the geometry to direction-periodic to avoid this problem. However, as per the manual and as the error you report says, you can’t have a dynamics box ( = barostat) in the pulling direction. To get rid of this problem you can turn off the barostat in the direction you are pulling (which in your case is along the y axis, as you set pull-coord1-vec = 0 1 0).
Thank you for the clarification and suggestion. It worked!! I switched the barostat off during the MD part of sampling. However, I have some more doubts:
Is it okay if I change the pull-geometry to direction-periodic during the sampling part and use the direction for the SMD step (where I generate the windows)?
From what I have gathered through the papers, NPT ensemble is generally used during sampling. Is turning the barostat off going to affect the outcome of the Umbrella sampling?
A technical point before getting through your questions. To avoid the error you were getting you can i) switch off the barostat, that is, run an NVT sim, as you did, or ii) you can switch it off partially, meaning that you set the compressibility to zero in the direction you are pulling but keep the barostat on. This means that you are still running an NPT simulation but the box can deform only in the directions where the compressibility is not zero, which can be achieved only with a semi-isotropic or an anisotropic pressure coupling. Back in the day I used to do this to generate sampling windows for lipid bilayers, where the pressure coupling is naturally semi-isotropic and I set to zero the compressibility along the z axis (as that was the direction in which I was pulling). I guess your case is a box of water with a protein-ligand system and a isotropic pressure coupling? Then, whether you should completely switch off the barostat is not a easy answer. However, if you have an isotropic coupling, this would mean that you have to switch to a semi-isotropic/anisotropic coupling and then generate the windows. In a situation like yours where the geometry is isotropic by construction this could induce box deformations as all your directions should be coupled together because you have the same phase (water) everywhere. So, I guess the NVT is the safest option with your box? Nevertheless, I would seriously try to understand if those distant sampling windows are really needed or not.
Regarding your questions:
No, that’s actually the opposite. The SMD has a distance that is changing in time, and that’s where I would expect the problem to appear as the reference position becomes ill defined after > 0.5 box length. So I would say that you need the direction-periodic for SMD, and then you can probably revert back to direction during the sampling, as the position is defined once and doesn’t change (too much, if restrained properly). So maybe you can use NVT to generate the windows and revert back to NPT + direction or distance during the sampling? The change of ensemble might not be a problem if you equilibrate the windows in NPT, hoping that the box rescaling won’t change much the windows positioning. Hard to predict without trying, though.
In general, I think yes. NPT and NVT are different ensembles and you are sampling different equilibrium distributions, and therefore different free energies. Whether this will give rise to different free energy profiles is something that cannot be answered in general, but in principle NVT and NPT are technically different ensembles.
Yeah, the conditions defined for NPT are according to isotropic pressure coupling. As per the explanation, I will try running the simulation using semi-isotropic coupling.
The ligand seems to be tightly bound with protein (it starts pulling away after _150 ps), so initially I kept nsteps = 1 ns. However, at 1 ns, the ligand gets pulled towards the extreme edge of the PBC (which I guess is problematic).
So, in the current run, together with changing geometry, I also changed the nsteps to 500 ps so that do not overestimate the results. Also, this time, I specifically indexed the binding site and treated it as a pull-coord-group1.
Will look into all the suggestions you made. Thanks a ton for the help!!