S1 Geom optimization based on TDDFT and delta SCF approach

Hello, I am comparing S₁ geometry optimizations for a conjugated π-system (OPP) using:

  • TDDFT (CAM-B3LYP/6-31G(d,p), CPCM solvent)

  • ΔSCF with SGM (same functional, basis, and solvent)
    Both calculations correspond to the lowest singlet excited state, which TDDFT indicates is predominantly a HOMO→LUMO π→π* excitation. Here is what I observe:

  • The TDDFT S₁ optimization shows noticeable structural relaxation relative to S₀ (bond-length alternation changes and larger overall distortion).

  • The ΔSCF S₁ optimization results in much smaller geometry changes, even though:

    • The occupation was enforced correctly (HOMO→LUMO),
    • ⟨S²⟩ ≈ 0 (no spin contamination),
    • SCF converged cleanly using SGM.

The vertical excitation energies from TDDFT and ΔSCF are reasonably close (3.4 eV and 3.3 eV), but the optimized geometries differ significantly in terms of relaxation magnitude. My questions: Is it expected that ΔSCF may underestimate structural relaxation compared to TDDFT for local π→π* excitations? And what would be the correct approach to handle this situation?
I am sharing both the input file:
TDDFT:
$molecule
0 1
H 1.3292676 -1.6836893 -3.0809952
C 0.7356205 -0.9500791 -3.6169987
C 0.0000000 0.0000000 -2.8999730
C -0.7356205 0.9500791 -3.6169987
H -1.3292676 1.6836893 -3.0809952
C -0.7357354 0.9503581 -5.0068087
H -1.3188870 1.6920649 -5.5435948
C 0.0000000 0.0000000 -5.7078619
H 0.0000000 0.0000000 -6.7931470
C 0.7357354 -0.9503581 -5.0068087
H 1.3188870 -1.6920649 -5.5435948
C 0.0000000 0.0000000 -1.4156563
C -0.0005202 1.1978303 -0.6937517
C 0.0005202 1.1978303 0.6937517
C 0.0000000 0.0000000 1.4156563
C -0.0005202 -1.1978303 0.6937517
C 0.0005202 -1.1978303 -0.6937517
C 0.0000000 0.0000000 2.8999730
C 0.7356205 0.9500791 3.6169987
C 0.7357354 0.9503581 5.0068087
C 0.0000000 0.0000000 5.7078619
C -0.7357354 -0.9503581 5.0068087
C -0.7356205 -0.9500791 3.6169987
H 0.0109345 2.1444854 -1.2243742
H -0.0109345 2.1444854 1.2243742
H 0.0109345 -2.1444854 1.2243742
H -0.0109345 -2.1444854 -1.2243742
H 1.3292676 1.6836893 3.0809952
H 1.3188870 1.6920649 5.5435948
H 0.0000000 0.0000000 6.7931470
H -1.3188870 -1.6920649 5.5435948
H -1.3292676 -1.6836893 3.0809952
$end

$rem
BASIS = 6-31G(d,p)
GUI = 2
JOB_TYPE = Optimization
GEOM_OPT_MAX_CYCLES = 200
METHOD = CAMB3LYP
SCF_CONVERGENCE = 8
SCF_MAX_CYCLES = 100
CIS_N_ROOTS = 3
CIS_SINGLETS = TRUE
CIS_TRIPLETS = FALSE
RPA = TRUE
CIS_STATE_DERIV = 1
SOLVENT_METHOD = PCM
SYMMETRY = FALSE
SYM_IGNORE = TRUE
MEM_TOTAL = 248000
MEM_STATIC = 4000
$end

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

$solvent
DIELECTRIC 37.5
OPTICALDIELECTRIC 1.8068
$end
delSCF:
$molecule
0 1
H 1.3292676 -1.6836893 -3.0809952
C 0.7356205 -0.9500791 -3.6169987
C 0.0000000 0.0000000 -2.8999730
C -0.7356205 0.9500791 -3.6169987
H -1.3292676 1.6836893 -3.0809952
C -0.7357354 0.9503581 -5.0068087
H -1.3188870 1.6920649 -5.5435948
C 0.0000000 0.0000000 -5.7078619
H 0.0000000 0.0000000 -6.7931470
C 0.7357354 -0.9503581 -5.0068087
H 1.3188870 -1.6920649 -5.5435948
C 0.0000000 0.0000000 -1.4156563
C -0.0005202 1.1978303 -0.6937517
C 0.0005202 1.1978303 0.6937517
C 0.0000000 0.0000000 1.4156563
C -0.0005202 -1.1978303 0.6937517
C 0.0005202 -1.1978303 -0.6937517
C 0.0000000 0.0000000 2.8999730
C 0.7356205 0.9500791 3.6169987
C 0.7357354 0.9503581 5.0068087
C 0.0000000 0.0000000 5.7078619
C -0.7357354 -0.9503581 5.0068087
C -0.7356205 -0.9500791 3.6169987
H 0.0109345 2.1444854 -1.2243742
H -0.0109345 2.1444854 1.2243742
H 0.0109345 -2.1444854 1.2243742
H -0.0109345 -2.1444854 -1.2243742
H 1.3292676 1.6836893 3.0809952
H 1.3188870 1.6920649 5.5435948
H 0.0000000 0.0000000 6.7931470
H -1.3188870 -1.6920649 5.5435948
H -1.3292676 -1.6836893 3.0809952
$end

$rem
JOB_TYPE SP
METHOD CAMB3LYP
BASIS 6-31G(d,p)
SCF_ALGORITHM DIIS
SCF_CONVERGENCE 8
MAX_SCF_CYCLES 200
SYM_IGNORE TRUE
SOLVENT_METHOD PCM
MEM_TOTAL 248000
MEM_STATIC 4000
$end

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

$solvent
DIELECTRIC 37.5
OPTICALDIELECTRIC 1.8068
$end

@@@

$molecule
read
$end

$rem
JOB_TYPE OPTIMIZATION
METHOD CAMB3LYP
BASIS 6-31G(d,p)
UNRESTRICTED TRUE
SCF_ALGORITHM SGM_LS
SCF_CONVERGENCE 6
SCF_GUESS READ
MAX_SCF_CYCLES 200
SYM_IGNORE TRUE
SOLVENT_METHOD PCM
MEM_TOTAL 248000
MEM_STATIC 4000
$end

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

$solvent
DIELECTRIC 37.5
OPTICALDIELECTRIC 1.8068
$end

$occupied
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 62
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61
$end

If you’re running an unrestricted (spin-symmetry broken) DeltaSCF calculation the <S^2> should be roughly 1. Are you sure you’re landing in the state you want? 7.8.2 Approximate Spin Purification‣ 7.8 Restricted Open-Shell and Δ SCF Methods ‣ Chapter 7 Open-Shell and Excited-State Methods ‣ Q-Chem 6.4 User’s Manual

Hi, I actually tried using OPSING = TRUE for approximate spin purification. But still I ended up getting <S^2> = 0.000000136. I am sharing my input file:
$comment
OPP Ground State (S0) in MeCN PCM
$end

$molecule
0 1
H 1.3292676 -1.6836893 -3.0809952
C 0.7356205 -0.9500791 -3.6169987
C 0.0000000 0.0000000 -2.8999730
C -0.7356205 0.9500791 -3.6169987
H -1.3292676 1.6836893 -3.0809952
C -0.7357354 0.9503581 -5.0068087
H -1.3188870 1.6920649 -5.5435948
C 0.0000000 0.0000000 -5.7078619
H 0.0000000 0.0000000 -6.7931470
C 0.7357354 -0.9503581 -5.0068087
H 1.3188870 -1.6920649 -5.5435948
C 0.0000000 0.0000000 -1.4156563
C -0.0005202 1.1978303 -0.6937517
C 0.0005202 1.1978303 0.6937517
C 0.0000000 0.0000000 1.4156563
C -0.0005202 -1.1978303 0.6937517
C 0.0005202 -1.1978303 -0.6937517
C 0.0000000 0.0000000 2.8999730
C 0.7356205 0.9500791 3.6169987
C 0.7357354 0.9503581 5.0068087
C 0.0000000 0.0000000 5.7078619
C -0.7357354 -0.9503581 5.0068087
C -0.7356205 -0.9500791 3.6169987
H 0.0109345 2.1444854 -1.2243742
H -0.0109345 2.1444854 1.2243742
H 0.0109345 -2.1444854 1.2243742
H -0.0109345 -2.1444854 -1.2243742
H 1.3292676 1.6836893 3.0809952
H 1.3188870 1.6920649 5.5435948
H 0.0000000 0.0000000 6.7931470
H -1.3188870 -1.6920649 5.5435948
H -1.3292676 -1.6836893 3.0809952
$end

$rem
JOB_TYPE SP
METHOD CAMB3LYP
BASIS 6-31G(d,p)
SCF_ALGORITHM DIIS
SCF_CONVERGENCE 8
MAX_SCF_CYCLES 200
SYM_IGNORE TRUE
UNRESTRICTED TRUE
SOLVENT_METHOD PCM
MEM_TOTAL 248000
MEM_STATIC 4000
$end

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

$solvent
DIELECTRIC 37.5
OPTICALDIELECTRIC 1.8068
$end

@@@

$molecule
read
$end

$comment
OPP S1 Delta-SCF (HOMO->LUMO) in MeCN PCM
$end

$rem
JOB_TYPE SP
METHOD CAMB3LYP
BASIS 6-31G(d,p)
GEN_SCFMAN true
UNRESTRICTED TRUE
OPSING TRUE
SCF_ALGORITHM SGM_LS
SCF_GUESS READ
SCF_CONVERGENCE 6
MAX_SCF_CYCLES 200
SYM_IGNORE TRUE
SOLVENT_METHOD PCM
MEM_TOTAL 248000
MEM_STATIC 4000
$end

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

$solvent
DIELECTRIC 37.5
OPTICALDIELECTRIC 1.8068
$end

$occupied
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 62
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61
$end
Can anyone help me to solve this issue?

I don’t think the open-shell singlet (OPSING or Yamaguchi spin purification) is working in this job, for reasons that aren’t immediately clear. Please run the sample jobs OpSing1.in and OpSing2.in, you will notice that in each case Q-Chem runs a singlet and a triplet calculation at each geometry, which I don’t see when I run your input. It may be that you’ve included some other input options that are incompatible with OPSING=TRUE so I’d encourage you to strip down your input file until you isolate the problem, then let us know what it is.

Thank you so much for your suggestions — I will try them. In the meantime, I used the MOM approach to prevent variational collapse. With MOM, I obtained ⟨S²⟩ ≈ 1 and the excitation energy is comparatively higher. I was wondering whether the MOM approach also performs any kind of spin purification?
$comment
OPP Ground State (S0) in MeCN PCM
$end

$molecule
0 1
H 1.3292676 -1.6836893 -3.0809952
C 0.7356205 -0.9500791 -3.6169987
C 0.0000000 0.0000000 -2.8999730
C -0.7356205 0.9500791 -3.6169987
H -1.3292676 1.6836893 -3.0809952
C -0.7357354 0.9503581 -5.0068087
H -1.3188870 1.6920649 -5.5435948
C 0.0000000 0.0000000 -5.7078619
H 0.0000000 0.0000000 -6.7931470
C 0.7357354 -0.9503581 -5.0068087
H 1.3188870 -1.6920649 -5.5435948
C 0.0000000 0.0000000 -1.4156563
C -0.0005202 1.1978303 -0.6937517
C 0.0005202 1.1978303 0.6937517
C 0.0000000 0.0000000 1.4156563
C -0.0005202 -1.1978303 0.6937517
C 0.0005202 -1.1978303 -0.6937517
C 0.0000000 0.0000000 2.8999730
C 0.7356205 0.9500791 3.6169987
C 0.7357354 0.9503581 5.0068087
C 0.0000000 0.0000000 5.7078619
C -0.7357354 -0.9503581 5.0068087
C -0.7356205 -0.9500791 3.6169987
H 0.0109345 2.1444854 -1.2243742
H -0.0109345 2.1444854 1.2243742
H 0.0109345 -2.1444854 1.2243742
H -0.0109345 -2.1444854 -1.2243742
H 1.3292676 1.6836893 3.0809952
H 1.3188870 1.6920649 5.5435948
H 0.0000000 0.0000000 6.7931470
H -1.3188870 -1.6920649 5.5435948
H -1.3292676 -1.6836893 3.0809952
$end

$rem
JOB_TYPE SP
METHOD CAMB3LYP
BASIS 6-31G(d,p)
SCF_ALGORITHM DIIS
SCF_CONVERGENCE 8
MAX_SCF_CYCLES 200
SYM_IGNORE TRUE
UNRESTRICTED TRUE
SOLVENT_METHOD PCM
MEM_TOTAL 248000
MEM_STATIC 4000
$end

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

$solvent
DIELECTRIC 37.5
OPTICALDIELECTRIC 1.8068
$end

@@@

$molecule
read
$end

$comment
OPP S1 Delta-SCF (HOMO->LUMO) in MeCN PCM
$end

$rem
JOB_TYPE SP
METHOD CAMB3LYP
BASIS 6-31G(d,p)
GEN_SCFMAN true
UNRESTRICTED TRUE
MOM_START TRUE
MOM_METHOD MOM
SCF_ALGORITHM SGM_LS
SCF_GUESS READ
SCF_CONVERGENCE 6
MAX_SCF_CYCLES 200
SYM_IGNORE TRUE
SOLVENT_METHOD PCM
MEM_TOTAL 248000
MEM_STATIC 4000
$end

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

$solvent
DIELECTRIC 37.5
OPTICALDIELECTRIC 1.8068
$end

$occupied
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 62
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61
$end

MOM is simply a convergence algorithm, it does not perform spin purification although you could do it yourself by also running the triplet calculation and using Yamaguchi’s formula, Eq. (8) in this paper:
https://doi.org/10.1021/acs.jctc.0c00502
(This is the same formula used in the OPSING procedure.)

With <S^2> = 1, you are likely getting the open-shell singlet state. Within a single-determinant formalism, it is an approximately equal mixture of singlet and triplet.

Hi, I was able to perform ΔSCF excited-state calculations using the OPSING approach. In the gas-phase calculation, the job completed successfully and it did complete both singlet and triplet state energy. However, when I include PCM (CPCM) in the calculation, the triplet energy is not being computed or reported. The calculation completes, but I do not see a stable triplet solution in the output. I am adding both the input files:
OPSING in Gas Phase:
$comment
OPP Ground State (S0) in GAs
$end

$molecule
0 1
H 1.3292676 -1.6836893 -3.0809952
C 0.7356205 -0.9500791 -3.6169987
C 0.0000000 0.0000000 -2.8999730
C -0.7356205 0.9500791 -3.6169987
H -1.3292676 1.6836893 -3.0809952
C -0.7357354 0.9503581 -5.0068087
H -1.3188870 1.6920649 -5.5435948
C 0.0000000 0.0000000 -5.7078619
H 0.0000000 0.0000000 -6.7931470
C 0.7357354 -0.9503581 -5.0068087
H 1.3188870 -1.6920649 -5.5435948
C 0.0000000 0.0000000 -1.4156563
C -0.0005202 1.1978303 -0.6937517
C 0.0005202 1.1978303 0.6937517
C 0.0000000 0.0000000 1.4156563
C -0.0005202 -1.1978303 0.6937517
C 0.0005202 -1.1978303 -0.6937517
C 0.0000000 0.0000000 2.8999730
C 0.7356205 0.9500791 3.6169987
C 0.7357354 0.9503581 5.0068087
C 0.0000000 0.0000000 5.7078619
C -0.7357354 -0.9503581 5.0068087
C -0.7356205 -0.9500791 3.6169987
H 0.0109345 2.1444854 -1.2243742
H -0.0109345 2.1444854 1.2243742
H 0.0109345 -2.1444854 1.2243742
H -0.0109345 -2.1444854 -1.2243742
H 1.3292676 1.6836893 3.0809952
H 1.3188870 1.6920649 5.5435948
H 0.0000000 0.0000000 6.7931470
H -1.3188870 -1.6920649 5.5435948
H -1.3292676 -1.6836893 3.0809952
$end

$rem
JOB_TYPE SP
METHOD CAMB3LYP
BASIS 6-31G(d,p)
SCF_ALGORITHM DIIS
SCF_CONVERGENCE 8
MAX_SCF_CYCLES 200
SYM_IGNORE TRUE
UNRESTRICTED TRUE
MEM_TOTAL 248000
MEM_STATIC 4000
$end

@@@

$molecule
read
$end

$comment
OPP S1 Delta-SCF (HOMO->LUMO) in GAS
$end

$rem
JOB_TYPE SP
METHOD CAMB3LYP
BASIS 6-31G(d,p)
UNRESTRICTED TRUE
OPSING TRUE
SCF_ALGORITHM GDM
SCF_GUESS READ
SCF_CONVERGENCE 6
MAX_SCF_CYCLES 200
SYM_IGNORE TRUE
MEM_TOTAL 248000
MEM_STATIC 4000
$end

$occupied
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 62
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61
$end
OPSING with CPCM:
$comment
OPP Ground State (S0) in MeCN PCM
$end

$molecule
0 1
H 1.3292676 -1.6836893 -3.0809952
C 0.7356205 -0.9500791 -3.6169987
C 0.0000000 0.0000000 -2.8999730
C -0.7356205 0.9500791 -3.6169987
H -1.3292676 1.6836893 -3.0809952
C -0.7357354 0.9503581 -5.0068087
H -1.3188870 1.6920649 -5.5435948
C 0.0000000 0.0000000 -5.7078619
H 0.0000000 0.0000000 -6.7931470
C 0.7357354 -0.9503581 -5.0068087
H 1.3188870 -1.6920649 -5.5435948
C 0.0000000 0.0000000 -1.4156563
C -0.0005202 1.1978303 -0.6937517
C 0.0005202 1.1978303 0.6937517
C 0.0000000 0.0000000 1.4156563
C -0.0005202 -1.1978303 0.6937517
C 0.0005202 -1.1978303 -0.6937517
C 0.0000000 0.0000000 2.8999730
C 0.7356205 0.9500791 3.6169987
C 0.7357354 0.9503581 5.0068087
C 0.0000000 0.0000000 5.7078619
C -0.7357354 -0.9503581 5.0068087
C -0.7356205 -0.9500791 3.6169987
H 0.0109345 2.1444854 -1.2243742
H -0.0109345 2.1444854 1.2243742
H 0.0109345 -2.1444854 1.2243742
H -0.0109345 -2.1444854 -1.2243742
H 1.3292676 1.6836893 3.0809952
H 1.3188870 1.6920649 5.5435948
H 0.0000000 0.0000000 6.7931470
H -1.3188870 -1.6920649 5.5435948
H -1.3292676 -1.6836893 3.0809952
$end

$rem
JOB_TYPE SP
METHOD CAMB3LYP
BASIS 6-31G(d,p)
SCF_ALGORITHM DIIS
SCF_CONVERGENCE 8
MAX_SCF_CYCLES 200
SYM_IGNORE TRUE
UNRESTRICTED TRUE
SOLVENT_METHOD PCM
MEM_TOTAL 248000
MEM_STATIC 4000
$end

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

$solvent
DIELECTRIC 37.5
OPTICALDIELECTRIC 1.8068
$end

@@@

$molecule
read
$end

$comment
OPP S1 Delta-SCF (HOMO->LUMO) in MeCN PCM
$end

$rem
JOB_TYPE SP
METHOD CAMB3LYP
BASIS 6-31G(d,p)
UNRESTRICTED TRUE
OPSING TRUE
SCF_ALGORITHM GDM
SCF_GUESS READ
SCF_CONVERGENCE 6
MAX_SCF_CYCLES 200
SOLVENT_METHOD PCM
SYM_IGNORE TRUE
MEM_TOTAL 248000
MEM_STATIC 4000
$end

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

$solvent
DIELECTRIC 37.5
OPTICALDIELECTRIC 1.8068
$end

$occupied
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 62
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61
$end

What version of Q-Chem are you running? When I modify the testing job OpSing1.in that ships with the code, to set SOLVENT_METHOD=PCM, it seems to work.

Update: when I run your input with Q-Chem v. 6.4, the singlet part completes normally but the triplet calculation (which does run) converges to NaN with a large DIIS error. Not sure what the issue is.

Update #2: I made two changes to make this work successfully. First, remove the HeavyPoints keyword from $pcm section because that dense of a grid is overkill. Second, use the default SCF algorithm rather than GDM. I suspect this is GDM leading the system into some bad solution.

Thank you so much. It worked.

(post deleted by author)

Hello, I performed a ΔSCF calculation with OPSING = TRUE (CAM-B3LYP/6-31G(d,p), CPCM solvent) and it ran successfully with both singlet and triplet calculation.
From the output I obtained:

  • Ground-state (S₀) SCF energy: –693.97978272 Hartree
  • Open-shell singlet energy (OPSING): –694.1180913962160 Hartree

This implies that the S₁ energy is lower than the S₀ energy, which physically this does not seem reasonable.
I am sharing input file. Can you please help to resolve this issue.
$comment
OPP Ground State (S0) in MeCN PCM
$end

$molecule
0 1
H 1.3292676 -1.6836893 -3.0809952
C 0.7356205 -0.9500791 -3.6169987
C 0.0000000 0.0000000 -2.8999730
C -0.7356205 0.9500791 -3.6169987
H -1.3292676 1.6836893 -3.0809952
C -0.7357354 0.9503581 -5.0068087
H -1.3188870 1.6920649 -5.5435948
C 0.0000000 0.0000000 -5.7078619
H 0.0000000 0.0000000 -6.7931470
C 0.7357354 -0.9503581 -5.0068087
H 1.3188870 -1.6920649 -5.5435948
C 0.0000000 0.0000000 -1.4156563
C -0.0005202 1.1978303 -0.6937517
C 0.0005202 1.1978303 0.6937517
C 0.0000000 0.0000000 1.4156563
C -0.0005202 -1.1978303 0.6937517
C 0.0005202 -1.1978303 -0.6937517
C 0.0000000 0.0000000 2.8999730
C 0.7356205 0.9500791 3.6169987
C 0.7357354 0.9503581 5.0068087
C 0.0000000 0.0000000 5.7078619
C -0.7357354 -0.9503581 5.0068087
C -0.7356205 -0.9500791 3.6169987
H 0.0109345 2.1444854 -1.2243742
H -0.0109345 2.1444854 1.2243742
H 0.0109345 -2.1444854 1.2243742
H -0.0109345 -2.1444854 -1.2243742
H 1.3292676 1.6836893 3.0809952
H 1.3188870 1.6920649 5.5435948
H 0.0000000 0.0000000 6.7931470
H -1.3188870 -1.6920649 5.5435948
H -1.3292676 -1.6836893 3.0809952
$end

$rem
JOB_TYPE SP
METHOD CAMB3LYP
BASIS 6-31G(d,p)
SCF_ALGORITHM DIIS
SCF_CONVERGENCE 8
MAX_SCF_CYCLES 200
SYM_IGNORE TRUE
UNRESTRICTED TRUE
SOLVENT_METHOD PCM
MEM_TOTAL 248000
MEM_STATIC 4000
$end

$pcm
THEORY CPCM
method swig
radii bondi
solver inversion
$end

$solvent
DIELECTRIC 37.5
OPTICALDIELECTRIC 1.8068
$end

@@@

$molecule
read
$end

$comment
OPP S1 Delta-SCF (HOMO->LUMO) in MeCN PCM
$end

$rem
JOB_TYPE SP
METHOD CAMB3LYP
BASIS 6-31G(d,p)
UNRESTRICTED TRUE
OPSING TRUE
SCF_ALGORITHM DIIS
SCF_GUESS READ
SCF_CONVERGENCE 6
MAX_SCF_CYCLES 200
SOLVENT_METHOD PCM
SYM_IGNORE TRUE
MEM_TOTAL 248000
MEM_STATIC 4000
$end

$pcm
THEORY CPCM
method swig
radii bondi
solver inversion
$end

$solvent
DIELECTRIC 37.5
OPTICALDIELECTRIC 1.8068
$end

$occupied
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 62
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61
$end

In general the spin-purified singlet lives in a different universe so I’m not sure it’s appropriate to compare it to the broken-symmetry result.

That said, when I run your input I don’t get a broken-symmetry result. The state with E = -693.9798 Ha also has <S^2> = 0.

I think that you need to consult with your advisor about this. This forum is for help running Q-Chem, we don’t have the capacity to answer questions about the science.

I tried your calculation with the new driver.

$molecule
  0  1
  H  1.3292676 -1.6836893 -3.0809952
  C  0.7356205 -0.9500791 -3.6169987
  C  0.0000000  0.0000000 -2.8999730
  C -0.7356205  0.9500791 -3.6169987
  H -1.3292676  1.6836893 -3.0809952
  C -0.7357354  0.9503581 -5.0068087
  H -1.3188870  1.6920649 -5.5435948
  C  0.0000000  0.0000000 -5.7078619
  H  0.0000000  0.0000000 -6.7931470
  C  0.7357354 -0.9503581 -5.0068087
  H  1.3188870 -1.6920649 -5.5435948
  C  0.0000000  0.0000000 -1.4156563
  C -0.0005202  1.1978303 -0.6937517
  C  0.0005202  1.1978303  0.6937517
  C  0.0000000  0.0000000  1.4156563
  C -0.0005202 -1.1978303  0.6937517
  C  0.0005202 -1.1978303 -0.6937517
  C  0.0000000  0.0000000  2.8999730
  C  0.7356205  0.9500791  3.6169987
  C  0.7357354  0.9503581  5.0068087
  C  0.0000000  0.0000000  5.7078619
  C -0.7357354 -0.9503581  5.0068087
  C -0.7356205 -0.9500791  3.6169987
  H  0.0109345  2.1444854 -1.2243742
  H -0.0109345  2.1444854  1.2243742
  H  0.0109345 -2.1444854  1.2243742
  H -0.0109345 -2.1444854 -1.2243742
  H  1.3292676  1.6836893  3.0809952
  H  1.3188870  1.6920649  5.5435948
  H  0.0000000  0.0000000  6.7931470
  H -1.3188870 -1.6920649  5.5435948
  H -1.3292676 -1.6836893  3.0809952
$end

$rem
  JOB_TYPE        SP
  METHOD          CAMB3LYP
  BASIS           6-31G(d,p)
  SCF_ALGORITHM   DIIS
  SCF_CONVERGENCE 8
  MAX_SCF_CYCLES  200
  SYM_IGNORE      TRUE
  MEM_TOTAL       128000
  MEM_STATIC      4000
  DELTA_SCF       TRUE
$end

$delta_scf
  singlet unrestricted
  triplet restricted
  singlet_SOMO_print both
$end

$plots
  grid_spacing 0.35
  grid_range   4.0
$end

Gives me

***************************************************************
*                                                             *
*              RESULTS OF DeltaSCF CALCULATION                *
*                                                             *
***************************************************************

 Energies
   Ground state                            =   -693.970595 Ha
   Unrestricted spin symmetry-broken state =   -693.807293 Ha
   Restricted open-shell triplet state     =   -693.828433 Ha

 <S^2> values
   Unrestricted spin symmetry-broken state =      1.007453 Ha
   Restricted open-shell triplet state     =      2.000000 Ha

 Excitation (singlet / triplet) energies
   Unrestricted spin symmetry-broken state =      0.163302 Ha /       4.443685 eV
   Restricted open-shell triplet state     =      0.142162 Ha /       3.868433 eV
   Spin-projected singlet state            =      0.184760 Ha /       5.027577 eV

 Transition dipole moment between ground state and unrestricted spin-symmetry broken state:
        < Phi_0 | mu_x | Phi_i^a >         =   0.000000 a.u. /   0.000000 Debye
        < Phi_0 | mu_y | Phi_i^a >         =  -0.000000 a.u. /  -0.000000 Debye
        < Phi_0 | mu_z | Phi_i^a >         =   2.531280 a.u. /   6.433871 Debye

And per the orbital overlaps this is the HOMO → LUMO transition. Not too sure what that lower state your getting is. I’m not sure if the SCF calculations with the new driver are aware of PCM but I know for sure they don’t work with geometry optimizations (yet). You can also try ROKS with

$delta_scf
  singlet restricted
  singlet_SOMO_print both
$end

Ok. in your input you’re getting

--- Corrected Open Shell Singlet Energy ---


    Singlet energy =     -693.9798018798941
    Triplet energy =     -693.8415123635707
    Open-shell singlet energy =     -694.1180913962176

 SCF   energy =  -694.11809140
 Total energy =  -694.11809140

because your unrestricted mixed-state calculation is collapsing to the ground state, so the energy computed with the OPSING procedure is nonsense. You need to prevent it from collapsing with MOM, SGM, STEP, or something like that. Add

MOM_START 1

to your $rem variables.