Corrected Linear Response

State-specific correction to solvent polarization response:
the corrected Linear Response approach
Ciro A. Guido and Stefano Caprasecca

Download: This report can be downloaded  here [DOI: 10.13140/RG.2.1.1903.7845]
Keywords
: Corrected Linear Response, cLR, PisaLR

These notes explain the correct procedure to obtain excitation energies (both absorption and emission) in solution by using the state-specific corrected linear response (cLR) approach of Caricato et al. [1], as developed in this group some time ago.

The method is already available in the G09 version of the Gaussian suite of programs [2].
The importance of such kind of simulations is due to the fact that for excitations involving a large density rearrangement, the Linear Response scheme is insufficient [3] because it does not account for the density-dependent relaxation of the solvent polarization. 

Please note that this is here explained for a continuum solvent described with the Polarizable Continuum Model.
It can however be used also with other polarizable models for the environment, such as the polarizable MM (MMPol) developed in this group. [4]

A summary of the theory is briefly detailed below. For further details, we refer to Refs. [1], [3], [4], [5].

A step-by-step guide on the usage in Gaussian follows.

Usage in Gaussian
Test
Schematic illustration of the steps necessary to compute absorption and fluorescence.
  1. Geometry optimization of the ground state
    Start by optimizing the ground state structure and calculating the frequencies, in the presence of the solvent:

    #p B3LYP/6-311G(d) SCRF=(IEFPCM,SOLVENT=ACETONITRILE) OPT FREQ
  2. Vertical excitation (absorption)
    The absorption calculation consists in a TD single-point excitation from the optimized ground state geometry in the presence of the solvent:

    #p CAM-B3LYP/6-311G(d) SCRF=(PISALR,SOLVENT=ACETONITRILE)
       TD=(NSTATES=5,ROOT=1,NONEQ) IOP(10/74=10)

    Here we have specified:
    (a) to use the cLR approach (PISALR);
    (b) the state of interest (e.g., ROOT=1);
    (c) that it is a non-equilibrium solvation (NONEQ). This indicates that we are considering a vertical process, where the dynamical response of the solvent immediately equilibrates with the excited state density, while the inertial part is not in equilibrium with it, and remains in equilibrium with the ground state;
    (d) the extra input (IOP) is needed when cLR is used with non-equilibrium solvation.
    Please note that, being a state-specific approach, the calculation should be repeated for each state you are interested in, changing the state of interest.
    In the output, apart from the standard LR results, the cLR energies will appear as follows:

    PCM Corrected Linear-Response - Non-equilibrium solvation:
    Corrected transition energy = 4.6316 eV
    Total energy after correction = -192.868770368 a.u.
  3.  Geometry optimization of the excited state
    As for step (1), a geometry optimization for the excited state of interest (here, state 1) is needed:

    #p CAM-B3LYP/6-311G(d) SCRF=(IEFPCM,SOLVENT=ACETONITRILE)
       TD(NSTATES=5,ROOT=1,EQSOLV) OPT FREQ=NUM

    Where:
    (a) specifying the OPT keyword when TD is active, and a certain state is selected, means that the geometry optimization refers to that particular excited state;
    (b) the keyword EQSOLV specifies that the solvent is now in equilibrium with the excited state
    (c) analytical frequencies are only available for the ground state. If they need to be computed for the excited state, numerical ones (much slower) must be used instead (FREQ=NUM)

  4. Determination of the reaction field of the excited state
    Once the excited state structure is optimized, the emission simulation requires a two-step calculation (here, steps 4 and 5). The two steps can be executed at once with one input file.
    In step 4, the inertial part of the reaction field due to the excited state is determined. This needs to be saved and then used, in the second calculation, for the final emission.
    (a) In TD, it is specified that the solvation is in equilibrium with the excited state (EQSOLV);
    (b) in the SCRF keywords, READ specifies that, at the end of the input section, extra keywords need to be read;
    (c) the extra keywords are written after the molecular structure. NONEQ=SAVE specifies that the inertial part of the solvation needs to be saved, to be used again in the following calculation
    (d) note that the checkpoint will contain all such information
  5. Vertical de-excitation (emission)
    In this step, the inertial reaction field is read (NONEQ=LOAD).Both the geometry and the guess for the electron density are read from the previous checkpoint (GUESS=READGEOM=CHECKPOINT)
    The following input file includes both steps 4 and 5.

    %chk=STEP4.chk
    
    #p CAM-B3LYP/6-311G(d) TD(NSTATES=5,ROOT=1,EQSOLV)
       SCRF=(PISALR,SOLVENT=ACETONITRILE,READ) IOp(10/74=20)
    
    TITLE SECTION: EMISSION FROM S1 IN ACN, STEP 4
    [charge] [multiplicity]
    [geometry]
    NONEQ = SAVE
    
    -- Link1 --
    %oldchk=STEP4.chk
    %chk=STEP5.chk
    
    #p CAM-B3LYP/6-311G(d) GUESS=READ GEOM=CHECKPOINT
       SCRF(SOLVENT=ACETONITRILE,READ)
    
    TITLE SECTION: EMISSION FROM S1 IN ACN, STEP 5
    [charge] [multiplicity]
    NONEQ=LOAD NONEQPARTITION=DYNAMICINERTIAL

    The de-excitation energy can be calculated in the following way:
    (a) take the total energy after correction:

    PCM State-Specific 1st Order Perturbation Theory - Equilibrium solvation.
    Corrected transition energy               =    2.5307          eV
    Total energy after correction             = -1032.54876403     a.u.

    (b) compare it with the ground state energy of the last calculation, i.e., step 5:

    SCF Done:  E(RCAM-B3LYP) =  -1032.64119112     A.U. after   11 cycles

    (c) the difference yields the corrected de-excitation energy in Hartree, which can be converted to eV:

    -1032.54876403 - (-1032.64119112) a.u. = 0.0924271 a.u.
    0.0924271 a.u. * 27.211399 eV/a.u. = 2.5151 eV
Theory overview
During the past decade, two principal formulations have been proposed and tested to couple the solute QM description with continuum models when describing the environment response following a vertical excitation/de-excitation. The two approaches, called Linear Response (LR) and State-specific (SS) models, respectively, can be illustrated by ideally reformulating the excitation in solution as a two-step process: first, the molecule, in its ground state, in equilibrium with the solvent, is excited to an excited state i in the presence of a solvent polarization frozen to the initial ground state (\omega^0_{0i}). In the second step, the dynamic component of the solvent polarization rearranges to equilibrate with the excite state charge density of the solute.
Based on the Rayleigh-Schrödinger perturbation theory, Cammi et al. [6] derived comparable expressions for both frameworks, namely:
\omega^{LR}_{0i} =\omega^0_{0i} + \left< i^{(0)} \right| \hat{V} \left| 0 \right> \left< 0 \right| \hat{Q}_{dyn} \left| i^{(0)} \right> =\omega^0_{0i} + R_{dyn}^{LR} (P^T_{0i})
\omega^{SS}_{0i} = \omega^0_{0i} +\frac{1}{2} \left[ \left< i^{(0)} \right| \hat{V} \left| i^{(0)} \right> - \left<0\right|\hat{V}\left|0\right>\right] \cdot
 \left[ \left< i^{(0)} \right| \hat{Q}_{dyn} \left| i^{(0)} \right>- \left<0\right| \hat{Q}_{dyn} \left|0\right> \right]
=\omega^0_{0i} +R_{dyn}^{SS} (P^\Delta_{0i})
where \hat{V} is the molecular electrostatic potential operator, \hat{Q}_{dyn} is the dynamical apparent charge operator. \left|i^{(0)}\right> is the i-th electronic state obtained in the presence of the fixed reaction field due to the ground state \left|0\right>. Note that this is just an explanatory picture, which does not correspond to the actual implementation.
In the LR formulation, the response R_{dyn} of the solvent dynamic polarization to the excitation is computed from the transition density (P^T_{0i}), while in the SS approach the same polarization is determined by the difference of the electron densities of the initial and final states, (P^\Delta_{0i}).
The LR formulation is generally well-suited for transitions characterized by a small change of the electron density, since the electrostatic interaction between the solvent and the excited state density is not very different from that of the ground state. By contrast, for excitations involving a large density rearrangement, the LR scheme is insufficient because the density-dependent relaxation of the solvent polarization is not included. In those cases, the difference between LR and SS schemes can therefore become significant. [3]
In the framework of TDDFT and CIS electronic structure methods, different state-specific formulations of TD with PCM have been suggested: the corrected linear response (cLR, also known as PisaLR [1]), and its self-consistent and variational generalization: the Vertical Excitation Method (VEM [4]). The cLR can be considered as the first cycle of the iterative process of VEM.
References
  1. M. Caricato, B. Mennucci, J. Tomasi, F. Ingrosso, R. Cammi, S. Corni & G. Scalmani
    J. Chem. Phys. 124, 124520 (2006) [DOI: 10.1063/1.2183309]
  2. Gaussian 09, M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr, J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, Ö. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski, D. J. Fox, Gaussian, Inc., Wallingford CT, 2009.
  3. C. A. Guido, D. Jacquemin, C. Adamo & B. Mennucci
    J. Chem. Theory Comput. 11, 5782 (2015) [DOI: 10.1021/acs.jctc.5b00679]
  4. A. V. Marenich, C. J. Cramer, D. G. Truhlar, C. A. Guido, B. Mennucci, G. Scalmani & M. J. Frisch
    Chem. Sci. 2, 2143 (2011) [DOI: 10.1039/C1SC00313E]
  5. C. Curutchet, A. Muñoz-Losa, S. Monti, J. Kongsted, G.D. Scholes & B. Mennucci
    J. Chem. Theory Comput. 5, 1838 (2009) [DOI: 10.1021/ct9001366]
  6. R. Cammi, S. Corni, B. Mennucci & J Tomasi
    J. Chem. Phys. 122, 104513 (2005) [DOI: 10.1063/1.1867373]