Why does the temperature fluctuate severely during the AIMD using an NVT ensemble?

Question

I am trying to raise the temperature of my model and keep the temperature at 300K, with an NVT ensemble and Nose-Hoover thermostat.
I am doing this job with VASP, here is the INCAR file.

DFT+U Calculation
LDAU   = .TRUE.        (Activate DFT+U)
LDATYPE=  2            (Dudarev; only U-J matters)
LDAUL  =  -1 -1 -1 -1 2 -1        (Orbitals for each species)
LDAUU  =  0 0 0 0  6  0       (U for each species)
LDAUJ  =  0 0 0 0 0 0        (J for each species)
LMAXMIX=  4            (Mixing cut-off; 4-d, 6-f)
 
Global Parameters
ISTART =  1            (Read existing wavefunction; if there)
ISPIN  =  2            (Non-Spin polarised DFT)
#MAGMOM = 4*0 4*0 48*0 24*1 1*1
LREAL  = Auto       (Projection operators: automatic)
ENCUT  =  520        (Cut-off energy for plane wave basis set, in eV)
PREC   =  Normal       (Precision level)
LWAVE  = .FALSE.        (Write WAVECAR or not)
LCHARG = .FALSE.        (Write CHGCAR or not)
ADDGRID= .TRUE.        (Increase grid; helps GGA convergence)
# LVTOT  = .TRUE.      (Write total electrostatic potential into LOCPOT or not)
# LVHAR  = .TRUE.      (Write ionic + Hartree electrostatic potential into LOCPOT or not)
# NELECT =             (No. of electrons: charged cells; be careful)
# LPLANE = .TRUE.      (Real space distribution; supercells)
NCORE   = 8           (Max is no. nodes; dont set for hybrids)
 
Electronic Relaxation
ISMEAR =  0            (Gaussian smearing; metals:1)
SIGMA  =  0.05         (Smearing value in eV; metals:0.2)
NELM   =  90           (Max electronic SCF steps)
NELMIN =  6            (Min electronic SCF steps)
EDIFF  =  1E-06        (SCF energy convergence; in eV)
PREC=Low
ISYM=0
MAXMIX = 40  
IWAVPR=11
############################# MD setting #####################################
# canonic (Nose) MD with XDATCAR updated every 10 steps
IBRION = 0
MDALGO = 2                     ä switch to select thermostat
SMASS =  1                     # Nose mass
ISIF = 2                       # this tag selects the ensemble in combination with the thermostat
NSW =1000
POTIM = 1
TEBEG = 300
TEEND = 300
NBLOCK = 4
##############################################################################

Here is the screenshot of my model:

My question is:
Why does the temperature fluctuate so much during AIMD?

Here is an result after I run 5ps AIMD:

Answer 1

The primary reason is that your simulation time is too short. To thermally equilibrate a system, the simulation time must be much longer than the period of the slowest phonon mode, so that even the slowest vibration has taken place multiple times. Taking a conservative estimate that the lowest phonon mode of your system is $100 \mathrm{~cm}^{-1}$, we see that the period of this mode is about $300 \mathrm{~fs}$ (note that this correlates very well with the oscillation period of the total energy!), so that you must perform the AIMD for at least a few picoseconds before the slowest vibration damps out.

Another, probably secondary, reason is that any small system will have notable thermal fluctuations even when already equilibrated. Your system is not particularly large, so this inherent thermal fluctuation may also contribute somewhat to the observed fluctuations.