I’m trying to run grand canonical monte carlo with particle exchanges only in order to test the response of my protein to different levels of water chemical potential. I have ran 1ns NVT and 10 ns NPT without a problem but when I try running muVT I get an “out of range atoms - cannot compute PPPM error” as well as some shake warnings - “Shake determinant < 0”. I can’t work out why the gcmc is causing these errors. I have tried altering the neighbor list update to delay 0 every 1 check yes but it has not helped. I am using colvars to restrain the rotation of the molecule but the error still occurs even with colvars off. Any advice would be useful.
I have attached the input file for the muVT simulation as well as the output and log file for that run. I have also attached the in.minimisation file where I set up the system if that would be of interest and the colvars file.
A simulation box with a protein in it is a very inhomogeneous environment, on the length scale of one water molecule, and thus the “chemical potential” of a water molecule will be a very jagged function of position. You could easily “explode” your protein inserting a water molecule just about anywhere inside it.
My initial guess is that this is physically and mathematically valid behaviour, and you would have to think carefully about why you are applying GCMC in this situation, instead of controlling the water’s “chemical potential” (again, not too well-defined in such heterogeneous all system) with a more conventional control like a barostat.
The GCMC implementation in LAMMPS is very primitive and thus has no configurational bias that will determine if a possible insertion location is suitable. This can lead, especially for extreme values of the chemical potential, to insertion of molecules where they have very high potential energy due to close contacts. This can lead to the shake issue and the PPPM out of range (which is caused by unusually fast moving atoms). There is very little that can be done for such cases. It is possible to reduce the impact of close contacts through using soft-core potentials (suffix /soft) with a suitable setting so that it does not impact the overall structure and dynamics too much. This will specifically avoid cases where any insertion leads to divergences (and thus forces that cannot be represented as floating point numbers).
It is probably best to look around for other MD/MC codes that would have a more advanced GCMC implementation.
Blockquote
instead of controlling the water’s “chemical potential” (again, not too well-defined in such heterogeneous all system) with a more conventional control like a barostat.
If I use a barostat to control the water content do you know of a good way to connect this to an experimentally measurable quantity such as the humidity?
no – but either your protein is fully immersed in water (in which case the isothermal compressibility of water is very small and thus a barostat would have minimal effect unless pressure values were absurdly high or low) or it isn’t (in which case you can just simulate the protein “in vacuo” with the simulation box having some amount of water corresponding to the appropriate volumetric concentration of water vapour – but that situation is completely outside the calibration data for MD force fields ).