Hi,
I am trying to perform a single-point (SP) calculation for a ground-state out-of-equilibrium condition using B3LYP-D3/6-311G** with non-equilibrium PCM implicit solvation, with the state specified via Marcus charges. However, I do not see any indication in the output file that the non-equilibrium solvation condition has been applied. I am sharing my input file below and would appreciate any guidance on whether this setup is correct, or what the recommended approach is for this type of calculation in Q-Chem.
$molecule
0 1
– OPP
-1 2
H -2.3947597 -2.6409220 -0.2454001
C -3.0518069 -1.7876196 -0.1277040
C -2.5125863 -0.4667684 -0.1623793
C -3.4659280 0.5866937 -0.0228681
H -3.1378048 1.6191429 -0.0553140
C -4.8162386 0.3352052 0.1475658
H -5.4984939 1.1751432 0.2496615
C -5.3112673 -0.9697263 0.1886599
H -6.3705717 -1.1609066 0.3223051
C -4.4032394 -2.0228000 0.0428483
H -4.7605915 -3.0491406 0.0605752
C -1.1053100 -0.2148207 -0.2990558
C -0.5567698 1.1104637 -0.3888651
C 0.7860462 1.3486905 -0.4907584
C 1.7545655 0.2955413 -0.5189614
C 1.2112801 -1.0264213 -0.4933254
C -0.1296088 -1.2684151 -0.3778886
C 3.1795211 0.5520356 -0.6221959
C 3.7294912 1.8318580 -0.3632027
C 5.0893338 2.0753459 -0.4635083
C 5.9749575 1.0595061 -0.8232920
C 5.4595600 -0.2102991 -1.0825247
C 4.1004707 -0.4610785 -0.9865479
H -1.2215805 1.9668608 -0.3977569
H 1.1162003 2.3785291 -0.5786257
H 1.8812596 -1.8793355 -0.5183401
H -0.4526656 -2.3026609 -0.3423568
H 3.0794702 2.6419450 -0.0522110
H 5.4662808 3.0716106 -0.2494934
H 7.0398354 1.2526379 -0.8999823
H 6.1274663 -1.0166352 -1.3726637
H 3.7401883 -1.4559559 -1.2228686
– TEA
1 2
N -0.5493130 -0.1222658 2.7468223
C 0.8886422 -0.1638791 2.7153696
C 1.5679838 0.5397708 3.8926644
C -1.2936050 -1.3517403 2.9464583
C -1.0600668 -1.9865012 4.3205058
C -1.2451067 1.1436980 2.8729565
C -1.7697236 1.4268815 4.2829139
H 1.2162543 0.3025379 1.7719636
H 1.1884373 -1.2109927 2.6697941
H 1.3265395 1.6039954 3.9253812
H 2.6490548 0.4443409 3.7727555
H 1.2827434 0.0846894 4.8435882
H -0.9864791 -2.0489919 2.1634966
H -2.3500338 -1.1294744 2.7960813
H -0.0076690 -2.2375020 4.4696957
H -1.6390550 -2.9106524 4.3777440
H -1.3806515 -1.3285182 5.1297210
H -2.0887181 1.1132129 2.1758088
H -0.5662601 1.9305654 2.5446317
H -2.5062600 0.6828014 4.5932736
H -2.2605376 2.4025735 4.2764386
H -0.9624457 1.4507261 5.0172042
$end
$rem
BASIS = 6-31G(d,p)
GUI = 2
JOB_TYPE = SP
METHOD = CAMB3LYP
SCF_CONVERGENCE = 8
SOLVENT_METHOD = PCM
SCF_GUESS = FRAGMO
SCF_PRINT_FRGM = TRUE
MEM_TOTAL = 248000
MEM_STATIC = 4000
$end
$pcm
heavypoints 590
STATESPECIFIC MARCUS
$end
$solvent
DIELECTRIC 37.5
OPTICALDIELECTRIC 1.8068
$end
The nonequilibrium PCM formalism has to do with excited states, not with ground-state optimization. See Section 5 of this review:
https://doi.org/10.1002/wcms.1519
Thank you so much for the paper. It was really useful. However, I have a question about calculating a single-point energy under non-equilibrium solvation conditions for an electron-transfer system. In my case, the reactant state (R) is OPP + TEA*, and the product state (P) is OPP*⁻ + TEA⁺. For the non-equilibrium PCM calculation at point 2 (i.e., the vertical energy of the product state evaluated at the reactant geometry), my understanding is that one must explicitly specify the initial state so that the solvent polarization corresponds to the reactant configuration. In this context, I am using FRAGMO to enforce the product-state charge distribution.
For non-equilibrium PCM, the second job is typically expected to use SCF_GUESS = READ in order to carry over the solvent polarization from the initial (reactant) calculation. However, when using FRAGMO, the second step instead requires SCF_GUESS = FRAGMO. Because of this, I am unsure whether the non-equilibrium calculation is actually reading and preserving the solvent polarization from the first job, or whether the FRAGMO guess effectively overrides the solvent response associated with the initial state. I would appreciate guidance on the correct and internally consistent way to combine FRAGMO with non-equilibrium PCM for this type of vertical electron-transfer calculation. I am sharing my input file:
$molecule
0 1
H -2.3947597 -2.6409220 -0.2454001
C -3.0518069 -1.7876196 -0.1277040
C -2.5125863 -0.4667684 -0.1623793
C -3.4659280 0.5866937 -0.0228681
H -3.1378048 1.6191429 -0.0553140
C -4.8162386 0.3352052 0.1475658
H -5.4984939 1.1751432 0.2496615
C -5.3112673 -0.9697263 0.1886599
H -6.3705717 -1.1609066 0.3223051
C -4.4032394 -2.0228000 0.0428483
H -4.7605915 -3.0491406 0.0605752
C -1.1053100 -0.2148207 -0.2990558
C -0.5567698 1.1104637 -0.3888651
C 0.7860462 1.3486905 -0.4907584
C 1.7545655 0.2955413 -0.5189614
C 1.2112801 -1.0264213 -0.4933254
C -0.1296088 -1.2684151 -0.3778886
C 3.1795211 0.5520356 -0.6221959
C 3.7294912 1.8318580 -0.3632027
C 5.0893338 2.0753459 -0.4635083
C 5.9749575 1.0595061 -0.8232920
C 5.4595600 -0.2102991 -1.0825247
C 4.1004707 -0.4610785 -0.9865479
H -1.2215805 1.9668608 -0.3977569
H 1.1162003 2.3785291 -0.5786257
H 1.8812596 -1.8793355 -0.5183401
H -0.4526656 -2.3026609 -0.3423568
H 3.0794702 2.6419450 -0.0522110
H 5.4662808 3.0716106 -0.2494934
H 7.0398354 1.2526379 -0.8999823
H 6.1274663 -1.0166352 -1.3726637
H 3.7401883 -1.4559559 -1.2228686
N -0.5493130 -0.1222658 2.7468223
C 0.8886422 -0.1638791 2.7153696
C 1.5679838 0.5397708 3.8926644
C -1.2936050 -1.3517403 2.9464583
C -1.0600668 -1.9865012 4.3205058
C -1.2451067 1.1436980 2.8729565
C -1.7697236 1.4268815 4.2829139
H 1.2162543 0.3025379 1.7719636
H 1.1884373 -1.2109927 2.6697941
H 1.3265395 1.6039954 3.9253812
H 2.6490548 0.4443409 3.7727555
H 1.2827434 0.0846894 4.8435882
H -0.9864791 -2.0489919 2.1634966
H -2.3500338 -1.1294744 2.7960813
H -0.0076690 -2.2375020 4.4696957
H -1.6390550 -2.9106524 4.3777440
H -1.3806515 -1.3285182 5.1297210
H -2.0887181 1.1132129 2.1758088
H -0.5662601 1.9305654 2.5446317
H -2.5062600 0.6828014 4.5932736
H -2.2605376 2.4025735 4.2764386
H -0.9624457 1.4507261 5.0172042
$end
$rem
JOB_TYPE = SP
METHOD = CAMB3LYP
BASIS = 6-31G(d,p)
CIS_N_ROOTS = 3
CIS_SINGLETS = TRUE
CIS_TRIPLETS = FALSE
RPA = TRUE
CIS_STATE_DERIV = 1
SOLVENT_METHOD = PCM
SCF_CONVERGENCE = 8
$end
$pcm
HEAVYPOINTS 590
NONEQUILIBRIUM
PrintLevel 1
$end
$solvent
Dielectric 37.5
OpticalDielectric 1.8068
$end
@@@
$molecule
0 1
– OPP
-1 2
H -2.3947597 -2.6409220 -0.2454001
C -3.0518069 -1.7876196 -0.1277040
C -2.5125863 -0.4667684 -0.1623793
C -3.4659280 0.5866937 -0.0228681
H -3.1378048 1.6191429 -0.0553140
C -4.8162386 0.3352052 0.1475658
H -5.4984939 1.1751432 0.2496615
C -5.3112673 -0.9697263 0.1886599
H -6.3705717 -1.1609066 0.3223051
C -4.4032394 -2.0228000 0.0428483
H -4.7605915 -3.0491406 0.0605752
C -1.1053100 -0.2148207 -0.2990558
C -0.5567698 1.1104637 -0.3888651
C 0.7860462 1.3486905 -0.4907584
C 1.7545655 0.2955413 -0.5189614
C 1.2112801 -1.0264213 -0.4933254
C -0.1296088 -1.2684151 -0.3778886
C 3.1795211 0.5520356 -0.6221959
C 3.7294912 1.8318580 -0.3632027
C 5.0893338 2.0753459 -0.4635083
C 5.9749575 1.0595061 -0.8232920
C 5.4595600 -0.2102991 -1.0825247
C 4.1004707 -0.4610785 -0.9865479
H -1.2215805 1.9668608 -0.3977569
H 1.1162003 2.3785291 -0.5786257
H 1.8812596 -1.8793355 -0.5183401
H -0.4526656 -2.3026609 -0.3423568
H 3.0794702 2.6419450 -0.0522110
H 5.4662808 3.0716106 -0.2494934
H 7.0398354 1.2526379 -0.8999823
H 6.1274663 -1.0166352 -1.3726637
H 3.7401883 -1.4559559 -1.2228686
– TEA
1 2
N -0.5493130 -0.1222658 2.7468223
C 0.8886422 -0.1638791 2.7153696
C 1.5679838 0.5397708 3.8926644
C -1.2936050 -1.3517403 2.9464583
C -1.0600668 -1.9865012 4.3205058
C -1.2451067 1.1436980 2.8729565
C -1.7697236 1.4268815 4.2829139
H 1.2162543 0.3025379 1.7719636
H 1.1884373 -1.2109927 2.6697941
H 1.3265395 1.6039954 3.9253812
H 2.6490548 0.4443409 3.7727555
H 1.2827434 0.0846894 4.8435882
H -0.9864791 -2.0489919 2.1634966
H -2.3500338 -1.1294744 2.7960813
H -0.0076690 -2.2375020 4.4696957
H -1.6390550 -2.9106524 4.3777440
H -1.3806515 -1.3285182 5.1297210
H -2.0887181 1.1132129 2.1758088
H -0.5662601 1.9305654 2.5446317
H -2.5062600 0.6828014 4.5932736
H -2.2605376 2.4025735 4.2764386
H -0.9624457 1.4507261 5.0172042
$end
$rem
BASIS = 6-31G(d,p)
GUI = 2
JOB_TYPE = SP
METHOD = CAMB3LYP
SCF_CONVERGENCE = 8
SOLVENT_METHOD = PCM
SCF_GUESS = FRAGMO
SCF_PRINT_FRGM = TRUE
MEM_TOTAL = 248000
MEM_STATIC = 4000
$end
$pcm
HEAVYPOINTS 590
STATESPECIFIC MARCUS
PrintLevel 1
$end
$solvent
Dielectric 37.5
OpticalDielectric 1.8068
$end
SCF_GUESS = READ does not actually read the solvent reaction field in Q-Chem’s implementation, just the MOs (density matrix) from the previous job. Nonequilibrium PCM means using two different dielectric constants, but applied to the converged SCF density (as described in the review mentioned above).
For the record, you should note that SCF_GUESS = FRAGMO does not “enforce” any kind of charge constraints. It’s possible that, given a particular form for the fragment densities, this guess might allow you to find an SCF solution (relaxed orbitals) with similar character, but that’s not guaranteed and you should always check. SCF_ALGORITHM = DIIS (which is the default) may not be the best choice for such cases.