Question about non-conserved energy in AIMD simulation

[zbwang]

Dear VASP Developers,

I have performed an AIMD simulation on the material Li10SiP2S12 using the PBE functional. The simulation was run in the NVT ensemble with a Nosé-Hoover thermostat.

According to my understanding (and the VASP wiki), the E column in the OSZICAR file for an NVT simulation represents the total conserved energy of the extended system. This value should include: (1) The ionic potential energy (part of E0); (2) The ionic kinetic energy (EK); (3) The thermostat's potential energy (SP); (4) The thermostat's kinetic energy (SK); Therefore, this total energy E should be conserved throughout the entire simulation.

However, in my calculation, the value for E keeps increasing linearly with simulation time, as shown in the attached plot. What I find puzzling is that while E is clearly drifting, the free energy term (F in the OSZICAR) appears to be stable and does not show a similar drift. To rule out a physical change in the system, I have already confirmed that the material did not undergo a phase transition and maintained its crystallinity throughout the run.

Could you please help me understand why the total energy E is not being conserved in my simulation? I have attached my INCAR file and the plots showing the trends for all relevant quantities from the OSZICAR.

Thank you for your help.

INCAR file:
ALGO = VeryFast
ENCUT = 400
EDIFF = 1e-04
IBRION = 0
ICHARG = 2
ISIF = 2
ISMEAR = 0
ISPIN = 1
ISYM = 0
KBLOCK = 100
LCHARG = False
LREAL = Auto
LWAVE = False
MAXMIX = 20
NBLOCK = 1
NELM = 100
NELMIN = 4
NSW = 26000
POTIM = 2.0
PREC = Low
SIGMA = 0.05
MDALGO = 2 ; SMASS = 0
TEBEG = 800
TEEND = 800
NPAR = 16

[marie-therese.huebsch]

Hi,

Short answer: your E drift is almost certainly numerical, not physical. In your INCAR a few settings systematically bias the forces/energies at every MD step even though the F looks steady. Tightening the electronic/plane-wave settings and using a smaller time step probably removes the linear ramp. Here is a list of things to try:

• In this first testing phase, you can set ISIF=0 since you don't use the stress tensor.

• Improve electronic convergence at each MD step. With ALGO=VeryFast and EDIFF=1e-4, the SCF is deliberately cheap. A small, systematic SCF residual biases the potential energy/forces in one direction, so the extended energy accumulates a linear drift.

• Tighten settings for the basis and FFT grids (PREC/ENCUT/LREAL).

You can test the above settings on a single structure (without running an MD). To this end, look at how the forces change w.r.t. improved settings. Try ALGO=Normal and EDIFF=1e-6, PREC=Normal, and increasing ENCUT.

• Time step (POTIM) may be too large at high T.

• Tune thermostat coupling / mass (SMASS).

The latter two can be tested in an MD run. Run for a couple of hundred steps and check if you still observe a linear drift.

Could you please try the above and let me know if it solves the issue?

[zbwang]

Thank you very much for your detailed suggestions!

I have now tried the tighter settings (ALGO=Normal, EDIFF=1e-6, PREC=Normal, and ENCUT=600). As you can see from the newly attached figures, the linear drift in the total energy (E (ETOTAL)) has been largely removed, and the simulation looks much more stable.

This leads me to a follow-up question. In my original (bad) simulation, the conserved total energy (E (ETOTAL)) showed a clear linear drift. However, the physical energy of the system (F (ETOTEN)) still appeared to be fluctuating correctly around an average value (similar to the orange line in the new, good plot). I am confused about this: Why was the physical energy F able to fluctuate normally around its average, as expected for an NVT simulation, even when the total conserved energy E was clearly non-conserved and drifting due to the numerical errors?

[marie-therese.huebsch]

Hi, I am glad you could improve the simulation with tighter settings.

However I have to say, there still seems to be a drift. Could you also plot T as a function of MD steps? After all in NVT the ensemble average of T should be conserved. Could be that the drift is related to the thermalization time.

The second point about the "free energy" F. This is really a bit confusing. This is not the free energy of the ion+electron system. It is actually the total energy of the electrons (summing over DFT eigenvalues and accounting for double counting) plus a term that arizes due to the fictitious electronic temperature introduced as smearing. This is required to perform the electronic minimization, but has nio connection to the temperature of the MD run.

Does this resolve the confusion?

[zbwang]

Thanks for the explanation! That clears up the confusion about the F term.

I’ve attached the plot of Temperature vs. Simulation Time below. It looks like the temperature is conserved around 800 K.

[marie-therese.huebsch]

Looks good :) I am glad I could help.