No S1 state for benzene when EOM-CCSD is used without symmetry

Hello Q-Chem Community,

I need to study the S1 state of benzene when the symmetry is broken. This is needed for PES calculation.
First I made some test when the symmetry is on, it uses D2h group, so S1 is described by B2u representation, [0,0,0,0,0,0,1,0], the excitation value is 4.8634 eV. When I turn off symmetry, then the lowest excitation energy is 6.3832 eV while I expected 4.8634 eV.

I tried to do the same test in Molpro and it worked. What am I doing wrong?
The inputs are below.
With symmetry

0 1
H    4.739433138    0.000000000    0.000000000
H  -4.739433138    0.000000000    0.000000000
H   2.369716569    4.104470040    0.000000000
H  -2.369716569    4.104470040    0.000000000
H  -2.369716569   -4.104470040    0.000000000
H   2.369716569   -4.104470040    0.000000000
C   2.704198094    0.000000000    0.000000000
C  -2.704198094    0.000000000    0.000000000
C   1.352099047    2.341903580    0.000000000
C  -1.352099047    2.341903580    0.000000000
C  -1.352099047   -2.341903580    0.000000000
C   1.352099047   -2.341903580    0.000000000

METHOD eom-ccsd
BASIS aug-cc-pvdz
EE_STATES [0,0,0,0,0,0,2,0]

Without symmetry, the molecule block is the same

METHOD eom-ccsd
BASIS aug-cc-pvdz

P.S. I am using the 5.1 version

This could happen due to the numeric nature of Davidson solver. I recommend to try the following, in this order:

  1. Use more guess states (EOM_GUESS_XXX keyword), e.g., 4 should do the trick.
  2. Use EOM_PRECONV_SINGLES (to find CIS-like solutions to be used as guesses in the EOM calculations (and use 4 guess vectors).
  3. Increase number of states (most expensive option, but also most reliable).
  4. Sometimes using tighter EOM_CONVERGENCE would help as well.

Thank you for your reply.
Did you mean EOM_NGUESS_SINGLES in the first point ?
There is no keyword EOM_GUESS_XXX.
I have tried all possible options including their combinations, but only 3 has worked so far.
Increasing number of states to 20 did the job. Now I am trying to figure out the minimum number of states
needed to obtain the desired state.

Another question, is there any way to get more digits in the output for
“Total energy = -231.33759266 a.u.”? The force constants are very sensitive to the number of s.f.
Is it possible to do it in the input? I am sorry if I missed it in the manual.
If the only option is recompiling, then it’s not for me, we have only binaries.