I’m currently trying to equilibrate a system of approximately 13 000 atoms using an NPT ensemble with the following fix: fix NPT all npt temp 300.0 300.0 1000.0 aniso 1.0 1.0 1000.0.
The system is a slit pore model composed of a solid silicate structure with water confined in between.
After 60 ns of simulation, the total energy still seems to be decreasing slightly. I calculated the mean and standard deviation of the energy over the last 5 ns, and used a moving average to smooth the data for the graphs.
My question is: are these fluctuations in pressure and energy expected due to the system size, or could it be a sign that the system is still not fully relaxed?
I would really appreciate any insights or suggestions.
This kind of question is very difficult to answer with confidence unless one has significant experience in that specific area of research. Thus, your chances to get meaningful help are best when talking to people that know and care about your research like your adviser(s), tutor(s), colleague(s), or collaborator(s).
Please also note that this is not really a question about LAMMPS, but about the science of your system and thus it is borderline off-topic.
I agree with what @akohlmey said—it’s difficult for us to say anything meaningful about a system we don’t know.
That being said, looking at your potential energy, it’s clear that your system exhibits a characteristic exponential relaxation, with a characteristic time Tc ~ 20 nanoseconds. Regardless of the system, when such relaxations occur, you typically need to wait at least 5 \times Tc for full equilibration. After that, you should see the potential energy reaching a plateau.
So, I would say no—your system is not fully relaxed.
You cannot determine whether a quantity is constant or changing just by calculating its mean and standard deviation. Here are five alternative checks:
Perform a linear regression of the quantity against time. If the slope is significant, then you should reject the hypothesis that the quantity has become constant against time.
Perform a normality test (such as Kolmogorov-Smirnov or Shapiro-Wilk) on the quantity values. If they do not appear normally distributed then you do not have a value fluctuating normally around their mean value. (The quantity could have a non-Gaussian stationary distribution – but that’s not so likely.)
Alternatively, you can calculate the kurtosis of the quantity – a normal distribution will have a kurtosis of 3 (excess kurtosis of 0), while a quantity normally distributed about a non-constant mean will be “platykurtic”, with a lower kurtosis (negative excess kurtosis) denoting far weaker distribution tails than expected.
Find the mean and standard error of the quantity in the last 5 ns, then find the mean and standard error of the quantity in the 5 ns before that. If the two means overlap within three standard errors then they are probably the same. But this test is harder than it seems, since the standard deviation is “too wide” while the usual formula, standard deviation divided by the square root of sample size, is “too narrow” – the number of independent samples is often smaller than the actual number of samples.
As a more robust alternative, calculate the kurtosis of the two 5ns (or other size) windows separately, and then calculate the kurtosis of the combined 10ns. If the combined kurtosis is less than each of the split kurtos…es, then the sample mean is probably still moving.
This sentence tickles some memories from long, long ago.
Confined systems can be very slow to equilibrate from MD alone.
The problem is that your potential energy hypersurface can be quite rugged and thus you would need large energy fluctuations for molecules to “jump” around. Those are unlikely for plain MD. You may fare better by combining your MD simulation with occasional MC displacement or rotation moves.
Thank you all srtee, @simongravelle, and @akohlmey for your help and your time.
Thank you srtee for this methods — as @simongravelle mentioned, I let it run for 100 ns; the energy is still decreasing, but it seems to be converging now. @akohlmey, I’m currently testing this method, only doing it on the water in the pore. I’ll share the results once I have them.
