MECP using penalty-constrained algorithm

Aarzoo · August 9, 2024, 10:22am

Hi all,
I calculated MECP using a penalty-constrained optimization algorithm between S0 and S1 states. I used the SF-TDDFT method for this.
Here is my input:
$molecule
0 3
coordinates
$end

$rem
BASIS = cc-pVDZ
METHOD = wB97X-D
JOB_TYPE = optimization
MECP_OPT = true
mecp_methods = penalty_function
MAX_CIS_CYCLES = 400
SCF_CONVERGENCE = 8
MAX_SCF_CYCLES = 400
GEOM_OPT_MAX_CYCLES = 300
THRESH = 14
sts_mom = true
UNRESTRICTED = true
SPIN_FLIP = true
mecp_state1 = [0,1]
mecp_state2 = [0,2]
cis_s2_thresh = 120
CIS_N_ROOTS = 4
SYMMETRY_IGNORE = true
SYMMETRY = false
SOLVENT_METHOD = PCM
$end

$solvent
DIELECTRIC 32.613000
OPTICALDIELECTRIC 1.765709
$end

and In output I got:

at the last optimization cycle:

Excited state 1: excitation energy (eV) = 0.4876
Total energy for state 1: -899.00611584 au
<S**2> : 0.1121
S( 2) → S( 1) amplitude = 0.9704 alpha
S( 2) → V( 2) amplitude = 0.1743 alpha

Excited state 2: excitation energy (eV) = 0.5069
Total energy for state 2: -899.00540911 au
<S**2> : 1.0957
S( 1) → S( 1) amplitude = 0.9881 alpha

Excited state 3: excitation energy (eV) = 0.8563
Total energy for state 3: -898.99256568 au
<S**2> : 0.9926
D( 67) → S( 2) amplitude = 0.1521
S( 2) → S( 2) amplitude = 0.9652 alpha

Excited state 4: excitation energy (eV) = 2.7267
Total energy for state 4: -898.92383098 au
<S**2> : 1.0458
D( 65) → S( 1) amplitude = -0.2151
D( 66) → S( 1) amplitude = -0.8870
D( 67) → S( 1) amplitude = -0.3176
D( 69) → S( 1) amplitude = 0.1715

and

– Optimization of the CI objective function (Cycle # 2) –
CI_Sigma (Lagrange multiplier) = 14.000
CI_Alpha (stability parameter) = 0.020 a.u.
Lower state energy (state 1) = -899.0061158375 a.u.
Upper state energy (state 2) = -899.0054091106 a.u.
Energy gap = 0.000707 a.u. (0.019 eV)

Value of F (objective function) = -899.00542478 a.u.
Value of sigma*G (penalty) = +0.00033769 a.u.

 Convergence criteria                        value     thresh   cnvgd?
 Component of F^x parallel      to G^x:  -0.00000606 0.00030000  YES
 Component of F^x perpendicular to G^x:  +0.00017731 0.00030000  YES
 Change in the value of F:               -0.00000013 0.00000100  YES
 Energy gap:                             +0.00070673 0.00100000  YES

I want to ask here that In this it is considering E.S-1 and E.S-2 as S0 and S1 states respectively but E.S-2 S2 value is 1.09. I have given CIS_S2_THRESH = 120 which is why it is considered it singlet state but normally can we say that the calculated results are correct? This is the energy gap between S0 and S1 not between S0 and T1?

kaushik · August 9, 2024, 8:26pm

It seems like your SF target states are quite spin contaminated, which is not uncommon. You may refer to the last paragraph in Section 9.8.3 of the Manual (9.8.3 Minimum-Energy Crossing Pointsâ£ 9.8 Nonadiabatic Couplings and Optimization of Minimum-Energy Crossing Points â£ Chapter 9 Exploring Potential Energy Surfaces: Searches for Critical Points and Molecular Dynamics â£ Q-Chem 6.2 Userâs Manual) and consider the SA-SF approach for verification.

Aarzoo · August 10, 2024, 4:18am

Thanks @kaushik for your suggestion. I tried the SA-SF approach but I think this will be on ROHF reference state.
And I am doing Unrestricted = true.
I am getting this error “SA-SF-CIS/SA-SF-TDDFT only supports restricted open-shell reference state”

jherbert · August 10, 2024, 3:40pm

That’s correct, SA-SF requires ROHF reference state, but you can at least try to match the spin-contaminated states with proper spin eigenstates and see if you can determine what’s going on.

Aarzoo · August 12, 2024, 6:30am

Hi @jherbert and @kaushik
I tried the SA-SF method and I gave these inputs:

$molecule
0 3
Already optimized coordinates (MECP using penalty-constrained optimization by SF-TDDFT method)
$end

$rem
BASIS = cc-pVDZ
METHOD = wB97X-D
JOB_TYPE = SP
SCF_CONVERGENCE = 8
UNRESTRICTED = false
SASF_CIS = 1
CIS_N_ROOTS = 5
SYMMETRY_IGNORE = true
SYMMETRY = false
SOLVENT_METHOD = PCM
$end

$solvent
DIELECTRIC 32.613000
OPTICALDIELECTRIC 1.765709
$end
and as output, I got:

     SA-SF-RPA Excitation Energies

Excited state 1: excitation energy (eV) = 2.8808
Total energy for state 1: -898.914099921657
Singlet
D( 69) → S( 1) amplitude = 0.3593
S( 2) → S( 1) amplitude = -0.9138
Excited state 2: excitation energy (eV) = 3.9907
Total energy for state 2: -898.873312581605
Singlet
D( 69) → S( 2) amplitude = -0.2508
S( 1) → S( 1) amplitude = -0.6730
S( 2) → S( 2) amplitude = 0.6730
Excited state 3: excitation energy (eV) = 5.2111
Total energy for state 3: -898.828465709689
Singlet
D( 63) → S( 1) amplitude = -0.2879
D( 65) → S( 1) amplitude = 0.4122
D( 66) → S( 1) amplitude = 0.4556
D( 68) → S( 1) amplitude = -0.5429
D( 69) → S( 2) amplitude = 0.1501
S( 1) → V( 2) amplitude = -0.3915
Excited state 4: excitation energy (eV) = 5.8726
Total energy for state 4: -898.804152906251
Singlet
D( 65) → S( 1) amplitude = -0.5919
D( 66) → S( 1) amplitude = 0.7301
D( 68) → S( 1) amplitude = 0.2699
Excited state 5: excitation energy (eV) = 6.0833
Total energy for state 5: -898.796412367205
Singlet
D( 63) → S( 1) amplitude = 0.2849
D( 65) → S( 1) amplitude = -0.2726
D( 68) → S( 1) amplitude = -0.5680
D( 69) → S( 1) amplitude = -0.1602
D( 69) → S( 2) amplitude = -0.5420
S( 1) → V( 2) amplitude = 0.1815
S( 2) → V( 3) amplitude = 0.2421

Here, it shows that E.S-2 and 3 are singlet states.

With one more input, it tried:
$molecule
0 3
Already optimized coordinates (MECP using penalty-constrained optimization by SF-TDDFT method)
$end

$rem
BASIS = cc-pVDZ
METHOD = wB97X-D
JOB_TYPE = SP
MECP_OPT = true
mecp_methods = penalty_function
MAX_CIS_CYCLES = 400
SCF_CONVERGENCE = 8
MAX_SCF_CYCLES = 400
THRESH = 14
UNRESTRICTED = false
SASF_CIS = 1
mecp_state1 = [0,1]
mecp_state2 = [0,2]
CIS_N_ROOTS = 4
SYMMETRY_IGNORE = true
SYMMETRY = false
SOLVENT_METHOD = PCM
$end

$solvent
DIELECTRIC 32.613000
OPTICALDIELECTRIC 1.765709
$end

and as output, I got:

     SA-SF-RPA Excitation Energies

Excited state 1: excitation energy (eV) = 2.8808
Total energy for state 1: -898.914094717520
Singlet
D( 69) → S( 1) amplitude = -0.3593
S( 2) → S( 1) amplitude = 0.9138
Excited state 2: excitation energy (eV) = 3.9907
Total energy for state 2: -898.873307332235
Singlet
D( 69) → S( 2) amplitude = -0.2508
S( 1) → S( 1) amplitude = -0.6730
S( 2) → S( 2) amplitude = 0.6730
Excited state 3: excitation energy (eV) = 5.2111
Total energy for state 3: -898.828460126063
Singlet
D( 63) → S( 1) amplitude = 0.2879
D( 65) → S( 1) amplitude = -0.4122
D( 66) → S( 1) amplitude = -0.4556
D( 68) → S( 1) amplitude = 0.5429
D( 69) → S( 2) amplitude = -0.1501
S( 1) → V( 2) amplitude = 0.3915
Excited state 4: excitation energy (eV) = 5.8726
Total energy for state 4: -898.804147613736
Singlet
D( 65) → S( 1) amplitude = -0.5919
D( 66) → S( 1) amplitude = 0.7302
D( 68) → S( 1) amplitude = 0.2699

I want to confirm that I have given the inputs to check the status of the E.S. states. Is that correct?