Choosing Time Step and the Fictitious Mass

Having the cutoff and the box dimensions optimized we turn to the dynamics. First test is to run microcanonical MD simulation for 5-10 ps with default values of fictitious electron mass (400 a.u.) and time step (5 a.u.) and see how the energy components (fictitious electronic, ionic kinetic¹, Kohn-Sham, total) behave. Already at this stage the correction (i.e. decrease) of the time step may be required. You can assess that by looking at the root mean square deviation of the total energy from simulations with different time step (see Figure 8).

What can be expected from NVE simulation is that the fictitious kinetic energy of the electrons may seem to be conserved very tightly. This however may be only illusory, since not all degrees of freedom may be excited in the microcanonical regime. In order to verify this you should turn on the thermostat for the ions. See ``Choosing the Nosé-Hoover chain thermostat parameters'' section for discussion on how to do this. In addition, running VIBRATIONAL ANALYSIS may be helpful as well as calculating the power spectrum of the electronic degrees of freedom using fourier.x program (see section 10.3). Nevertheless, before going any further you have to make sure, that during NVE simulation all the energy components are conserved. Any irregular behavior here is most probably the sign of too large time step. You can also ``stabilize'' the CP-dynamics by increasing the fictitious mass, but that would also increase the ``drag'' and consequently the effective mass of the ions.

Figure 9: Fictitious kinetic energy of the electrons ( $E_{\mathrm {KINC}}$ ) with different pseudopotentials, electron mass ($\mu$ ) and time step ($\Delta {t}$ ) from the CPMD simulation with thermostat on the ions.

During MD run with Nosé-Hoover chain thermostat with coupling parameter properly chosen all degrees of freedom become excited after some time. If the electronic and ionic spectra yet overlap you can see that $E_{\mathrm {KINC}}$ drifts like on the Figure 9. Decreasing electron mass helps improving adiabaticity, however this might not always lead to the strict conservation. If this was the case you can try thermostating the electrons with another Nosé-Hoover chain, choosing the coupling frequency much higher than that for the ions, e.g. 10'000 cm$^{-1}$ or more. The protocol here would be to run initial simulation with the thermostat on the ions only and then restart from the wave function, coordinates, velocities and the ionic thermostat parameters with:

RESTART WAVEFUNCTION COORDINATES VELOCITIES NOSEP ACCUMULATORS

and continue the run with the target electronic fictitious energy equal to the mean value at the end of the previous run. Doing so may result in negative drift in the total, conserved energy $E_{\mathrm {HAM}}$ though.

What you can also learn from Figure 9 is that switching from Troullier-Martins to Vanderbilt pseudopotentials and decreasing the plane wave cutoff from 70 to 35 Ry improves stability of the electron dynamics in this specific case (cf. green and dark blue plots, that share the same electron mass and time step). This can be attributed to reduced drag, since less wavefunction coefficients are coupled to the ionic degrees of freedom.