Unexpected Energy Ordering in TDDFT Nonequilibrium PCM

Questions
#1

Hi all,

I have a question about running TDDFT calculations for charged species with nonequilibrium PCM.

I am trying to read the solvent reaction field from an initial charged product state and apply it to the excited state using TDDFT for S1. However, the CDFT product state is a broken-symmetry open-shell singlet, with some spin contamination.

My questions are:

  1. Is this level of spin contamination acceptable for generating the product-state PCM reaction field?
  2. Could the broken-symmetry CDFT reference, or the way I am using rf_ptss_save/read, lead to incorrect energy ordering?
  3. In my calculation, the TDDFT S1 energy is coming out lower than the CDFT ground-state energy, which seems unphysical. Could this be related to the saved/read reaction field setup?

Any insights would be greatly appreciated.

$molecule
0 1
H         -3.40839       -2.04574       -1.88829
C         -3.89279       -1.12149       -1.58909
C         -3.10104       -0.00218       -1.20761
C         -3.82661        1.17366       -0.86618
H         -3.29083        2.05879       -0.53757
C         -5.21195        1.22146       -0.90647
H         -5.71613        2.14434       -0.62944
C         -5.96406        0.10440       -1.28259
H         -7.04827        0.14450       -1.31051
C         -5.27809       -1.06602       -1.62121
H         -5.83430       -1.94974       -1.92507
C         -1.65556       -0.05596       -1.16384
C         -0.84180        1.11650       -1.05000
C          0.52800        1.06514       -0.98466
C          1.25836       -0.16475       -1.02666
C          0.44651       -1.33598       -1.15686
C         -0.92368       -1.28519       -1.22076
C          2.70211       -0.21906       -0.93417
C          3.47372        0.89003       -0.48915
C          4.85659        0.83513       -0.39934
C          5.55770       -0.32564       -0.73927
C          4.82470       -1.43306       -1.17632
C          3.44203       -1.38498       -1.27497
H         -1.31065        2.09706       -1.05180
H          1.06787        2.00744       -0.93838
H          0.91670       -2.31604       -1.16250
H         -1.46127       -2.22826       -1.27427
H          2.97243        1.80128       -0.17804
H          5.39777        1.70981       -0.04623
H          6.63971       -0.36632       -0.66496
H          5.34214       -2.34808       -1.45466
H          2.92131       -2.26162       -1.64781
N          0.27583        0.13371        4.43300
C          1.52609        0.86536        4.39911
C          1.74843        1.66921        5.68680
C          0.28732       -1.31395        4.49453
C          0.86220       -1.81878        5.82503
C         -0.98315        0.85068        4.45617
C         -1.69789        0.69710        5.80570
H          1.46873        1.54604        3.54179
H          2.33335        0.15149        4.23464
H          0.96832        2.41964        5.83049
H          2.70800        2.18164        5.59646
H          1.77977        1.00989        6.55644
H          0.91398       -1.65926        3.66467
H         -0.73230       -1.66973        4.34451
H          1.89910       -1.50332        5.95835
H          0.83286       -2.90979        5.80611
H          0.26562       -1.46303        6.66722
H         -1.60252        0.42536        3.65871
H         -0.78248        1.89882        4.23283
H         -1.94709       -0.34626        6.01001
H         -2.62562        1.27016        5.75618
H         -1.08333        1.08771        6.61901
$end

$rem
   BASIS  =  6-31G**
   GUI  =  2
   JOB_TYPE  =  SP
   METHOD  =  wB97X-D
   SCF_CONVERGENCE  =  8
   SOLVENT_METHOD  =  PCM
   CDFT = TRUE
   CDFT_PRINT = TRUE
   CDFT_BECKE_POP = TRUE
   MEM_TOTAL = 248000
   MEM_STATIC = 4000
$end

$pcm
   THEORY  CPCM
   heavypoints 590
   method swig
   radii bondi
   solver inversion
   rf_ptss_save  true
$end

$solvent
   DIELECTRIC  37.5
   OPTICALDIELECTRIC  1.8068
$end

$cdft
1
1 1 32
1
1 1 32 s
$end

@@@

$molecule
 read
$end

$rem
   JOB_TYPE = SP
   METHOD = wB97X-D
   BASIS = 6-31G**
   CIS_N_ROOTS  =  10
   CIS_SINGLETS = TRUE
   CIS_TRIPLETS = FALSE
   RPA = TRUE
   CIS_RELAXED_DENSITY = TRUE
   SCF_MAX_CYCLES  = 200
   SOLVENT_METHOD = PCM
   SCF_CONVERGENCE = 8
$end

$pcm
   Theory  CPCM
   heavypoints 590
   method swig
   radii bondi
   solver inversion
   rf_ptss_read     true  
$end

$solvent
   Dielectric        37.5
   OpticalDielectric 1.8068
$end
#2

My first thought was “shouldn’t you get a negative excitation energy with a CDFT reference state?”, because the reference is not a true minimum in LCAO space due to the constraints, so there should exist some orbital rotation that lowers the energy. However, I don’t get that, and I get an excitation energy of just less than 1 eV. I do see <S^2> = 1.0 for the ground state, which is typical of an open-shell singlet and you’re not wrong to worry about it. If that’s a concern, you could always try a Delta-SCF or Noodleman-style spin purification.