Facing issue with spin-contamination

Hi everyone,

I am optimizing benzene using the SF-TDDFT method and I am encountering a spin-contamination issue for excited state 3. Is there any way to reduce or remove this spin contamination?
Using 5.4 version Q-Chem,

Input:

$molecule
0 3
  H         1.9585035166    1.5031400334   -0.0000157462
  C         1.1001149685    0.8451022363   -0.0000193933
  C         1.2819562868   -0.5301429871   -0.0000314124
  H         2.2809984962   -0.9445603477   -0.0000378001
  C         0.1818011255   -1.3751354500   -0.0000356229
  H         0.3241201502   -2.4474529701   -0.0000443085
  C        -1.1001007300   -0.8451473631   -0.0000279667
  H        -1.9584748809   -1.5032017595   -0.0000276636
  C        -1.2819349136    0.5300852147   -0.0000195873
  H        -2.2809813311    0.9444960887   -0.0000201224
  C        -0.1817723386    1.3750857557   -0.0000140019
  H        -0.3240958718    2.4474017619   -0.0000021974
$end

$rem
   BASIS  =  6-31G*
   JOB_TYPE  =  Optimization
   METHOD  =  BHHLYP
   SYMMETRY_IGNORE = TRUE
   MAX_CIS_CYCLES = 500
   MAX_SCF_CYCLES = 500
   THRESH = 14
   SPIN_FLIP = true
   CIS_N_ROOTS = 8
   CIS_STATE_DERIVATIVE = 1
   sts_mom = true
   UNRESTRICTED = FALSE
   SCF_CONVERGENCE  =  8
$end

and the output is:

        SF-DFT Excitation Energies              

(The first “excited” state might be the ground state)

Excited state 1: excitation energy (eV) = -5.0417
Total energy for state 1: -232.09385324 au
<S**2> : 0.0145
S( 2) → S( 1) amplitude = 0.9976 alpha

Excited state 2: excitation energy (eV) = 0.4614
Total energy for state 2: -231.89161772 au
<S**2> : 2.0243
D( 20) → S( 1) amplitude = 0.6681
S( 2) → V( 1) amplitude = 0.7403 alpha

Excited state 3: excitation energy (eV) = 0.6288
Total energy for state 3: -231.88546432 au
<S**2> : 1.0104
S( 1) → S( 1) amplitude = 0.7757 alpha
S( 2) → S( 2) amplitude = 0.6288 alpha

Excited state 4: excitation energy (eV) = 1.0717
Total energy for state 4: -231.86918786 au
<S**2> : 1.0142
S( 1) → S( 1) amplitude = -0.6225 alpha
S( 2) → S( 2) amplitude = 0.7716 alpha

Excited state 5: excitation energy (eV) = 1.6278
Total energy for state 5: -231.84875283 au
<S**2> : 0.0574
D( 20) → S( 1) amplitude = 0.7320
S( 2) → V( 1) amplitude = -0.6643 alpha

Excited state 6: excitation energy (eV) = 3.3823
Total energy for state 6: -231.78427658 au
<S**2> : 1.0177
S( 2) → V( 2) amplitude = 0.9919 alpha

I also tried to give keywords for ROHF (unrestricted = false) but I am getting error: TDDFT based on ROKS is not available

The SF-TDDFT method is know to have this problem. See, e.g.,

ROHF-based TDDFT is not available. You can use spin-adapted spin-flip TDDFT,
https://doi.org/10.1063/1.4937571
which is rigorously free of spin contamination.