Why do the results of NPT dynamics change so much with the heating rate?

GROMACS 2022.4:
GROMACS modification: No
Hello everyone,

I’m working on the heating kinetics of crystalline materials under NPT conditions, aiming to reach a target temperature of 550K. I’m using the same structural model and topology file for the experiment. The variables such as step size, the heat bath, the pressure bath, and corresponding time constants are kept constant. The only variable I’m changing is the rate of temperature increase. I’ve set the heating to occur at five different rates: 2ps, 5ps, 10ps, 20ps, and 40ps, to observe the time it takes for a 5° increase.

However, the outcomes have shown two distinct behaviors: a phase change occurring at the 10ps and 20ps rates, and no phase change at the 2ps, 5ps, and 40ps rates. I’ve prepared the mdp files and the results files for these five scenarios, displaying the relationships between temperature/density and time, as illustrated in the attached figure.

Could anyone shed some light on why this might be happening? Any advice would be greatly appreciated.
at40ps.mdp (1.2 KB)

at2ps.mdp (1.1 KB)

at5ps.mdp (1.1 KB)

at10ps.mdp (1.2 KB)

at20ps.mdp (1.2 KB)

Are you sure the difference is actually from the differences in heating rate, or might there be random differences? If you repeat the simulations (try first with repeating the ones with heating rates of 5K/5ps and 5K/40ps) five times, do you still not see this phase change behaviour? If you repeat the 5K/20ps simulation five times, do you always see the same phase change? Also, of the plots you posted, the 5K/10ps and 5K/20ps plots look suspiciously similar. Are they really from two different simulations with different heating rates?

Dear MagnusL,
I deeply regret the significant error I made during the data importation process. To clarify, I mistakenly imported the results for a 10ps simulation as though they were for 20ps, which made the two sets of data appear deceptively similar. I have since corrected this mistake by importing the accurate data and updating the erroneous figure. Consequently, the phase change is now observed exclusively at the 10ps rate, with no phase changes at the 2ps, 5ps, 20ps, and 40ps rates. This outcome is still perplexing. Following your advice, I intend to run additional simulations to better understand the results.

Regarding your query about “random differences,” could you provide some insights on how to reduce this randomness as much as possible?

Thank you very much.

Regarding “random differences”, provided that you start from different TPR files (with different random seeds, for generating velocities or for v-rescale (ld-seed)) the trajectories and simulations will be different. Therefore you should. in general, run multiple replicas to ensure that your observations are relevant and reproducible.

In your case, observing phase transitions will almost certainly be probability based, with higher probability at higher temperatures. The longer you simulate at a higher temperature the higher the probability that you will observe a phase transition. So, with enough replicas at each temperature, I would expect more observations of phase transitions at 5K/40ps. At least, I would be very careful stating that the heating rate itself has a large influence, based on one simulation (5K/10ps) with a moderate heating rate. Are you sure your phase transition is not just due to some system instability? It happens very quickly, so it might not be related to the rising temperature at all. I presume you have run a sufficiently long NPT equilibration at the starting temperature, right?