Choosing an appropriate SMASS value for NVT AIMD simulations

[ahterpet]

Hi everyone!

I'm running NVT AIMD simulations on pristine tantalum and tungsten cells at several temperatures ranging from 300 K to 1800 K. I've observed that decreasing the Nosé–Hoover thermostat effective mass (SMASS), corresponding to shorter thermal fluctuation periods (e.g., ~20 and ~10 time steps), leads to smaller temperature fluctuations.

From my understanding, a smaller SMASS increases the coupling strength between the thermostat and the system, which should suppress temperature fluctuations. If that’s the case, why not simply reduce SMASS as much as possible to minimize these fluctuations? How should one determine an appropriate SMASS value for a given simulation? Are there other factors that might contribute to large temperature fluctuations in AIMD?

I've attached plots of temperature vs. simulation step for both Ta and W at 1800 K, showing the reduction in temperature fluctuation with decreasing SMASS. I’ve also included the corresponding project directories for the default SMASS (period of ~40 time steps).

Any insight or recommendations would be greatly appreciated!

[marie-therese.huebsch]

Hi, thank you for posting this very nice visualization.

Regarding the choice of SMASS, you can find the following advice on the Wiki

For SMASS≥0, a canonical ensemble is simulated using the algorithm of Nosé. The Nosé mass controls the frequency of the temperature oscillations during the simulation. For SMASS=0, a Nosé-mass corresponding to period of 40 time steps will be chosen. The Nosé-mass should be set such that the induced temperature fluctuation show approximately the same frequencies as the typical 'phonon'-frequencies for the specific system. For liquids something like 'phonon'-frequencies might be obtained from the spectrum of the velocity auto-correlation function. If the ionic frequencies differ by an order of magnitude from the frequencies of the induced temperature fluctuations, Nosé thermostat and ionic movement might decouple leading to a non-canonical ensemble. The frequency of the approximate temperature fluctuations induced by the Nosé-thermostat is written to the OUTCAR file.

In other words, you want the thermostat to be able to populate all phonon modes. This can be achieved by tuning SMASS to obtain temperature oscillations at similar frequencies as the characteristic phonon frequencies. You can obtain that for pristine tantalum and tungsten following Computing the phonon dispersion and DOS. Personally, I would also recommend using NHC. The standard NV thermostat is recovered as a special case of NHC with NHC_NCHAINS=0, but the time scale of the temperature oscillations can be set (much more transparently) using NHC_PERIOD instead of SMASS. If the thermostat fails to populate the phonon modes, you may loose ergodicity. That means you won't sample the correct ensemble.

Does this answer your question?

[ahterpet]

Thank you! I really appreciate the detailed and insightful response.

I had one follow-up question: is there a simpler or more heuristic way to estimate an appropriate SMASS value without explicitly computing the phonon modes or DOS? For example, can one make a reasonable approximation based on system size, atomic mass, temperature, or the Debye frequency?

Additionally, are there any checks or diagnostics I could perform on the AIMD trajectory to verify that the simulation is indeed sampling the correct canonical ensemble (e.g., by examining certain distributions or correlation functions)? In other words, how do I know if my simulation results are reasonable/physical?

[marie-therese.huebsch]

I had one follow-up question: is there a simpler or more heuristic way to estimate an appropriate SMASS value without explicitly computing the phonon modes or DOS?

Yes, the Debey frequency should be an indicator as well. Or for common materials, you can look up if there are phonon spectra in literature.

For example, can one make a reasonable approximation based on system size [...]?

Mind that in smaller systems you would need larger thermal fluctuations to sample the correct ensemble. Imagine a perfectly thermalized system very large system. If you zoom in and look at only a small subset of atoms, they will have a higher or lower instantaniouse temperature. The smaller the subset the higher the diviation.

Additionally, are there any checks or diagnostics I could perform on the AIMD trajectory to verify that the simulation is indeed sampling the correct canonical ensemble (e.g., by examining certain distributions or correlation functions)? In other words, how do I know if my simulation results are reasonable/physical?

Depends on the quantity of interest. There is really a platora of checks one can think of. Using the ICONST file you can monitor certain parameters and check if the distribution behaves as expected. Likewise you can check the distribution of the velocities (VELOCITY tag). Finally, the velocity autocorrelation function can be used to compute the phonon spectra. Those should be close to the one obtained using the force constant matrix.

Please note that this goes slightly beyond the scope here. This forum is for advice and technical support on using VASP, but for scientific advise we recommend consulting your peers, literature or your supervisor. Feel free to post questions regarding, e.g. scientific best practices, in from users for users. We also always take a look at questions in from users for users and answer if we feel it is within the area of our expertise, but other users may have much more experience in a certain area of application. So it is great if you can exchange your knowledge there.

[ahterpet]

Thanks again for the detailed response! Apologies if my follow-up went a bit beyond the scope of this forum. I’ve been stuck on this issue for a while, but your suggestions have given me much clearer direction moving forward.

[marie-therese.huebsch]

Glad I could help :) Feel free to ask before you get stuck for too long. We can always move it to the other forum if appropriate.