TURBOMOLE manual
v8.0

10.9 Reaction field calculation (COSMO, PE)

The ricc2 program offers a number of reaction field coupling schemes and densities to include polarizable environments (COSMO or PE) in the calculations of energies and spectra. (Currently not yet available for linear and quadratic response properties.)

For responsive environments the embedding potential, also known as reaction field potential, is a functional of the QM systems’ electron density. For COSMO[280, 281, 282]: \[\begin{align} \hat{G}(\mathbf{D}) = - f\sum_{\kappa\lambda} \hat{V}_\kappa [{\cal A}^{-1} ]_{\kappa\lambda} \sum_{pq} V_{\lambda,pq} D_{pq} ~, \end{align}\](10.26) and for PE[283, 284, 285, 286]: \[\begin{align} \hat{G}(\mathbf{D}) = - \sum_{uv} \hat{\varepsilon}_u R_{uv} \sum_{pq} \varepsilon_{v,pq} D_{pq} ~. \end{align}\](10.27)

10.9.1 PTE/frozen solvent calculations

The “Perturbation on the Energy” (PTE) coupling scheme for correlated methods uses the embedding potential that has been determined in a preceeding self-consistent reaction field (SCRF) Hartree-Fock calculation, \(\hat{G}(\mathbf{D}^{\mathrm{HF}})\), and determines the parameters of the correlated wavefunction then in the usual way for the modified Hamiltonian \(\hat{H}^{vac}+\hat{G}(\mathbf{D}^{\mathrm{HF}})\).

In the ricc2 program, the PTE scheme implies for excitation energy and response calculations the frozen solvent approximation, i.e., also the response of \(\hat{G}\) to excitations in or perturbations of the QM system are neglected.

Although the PTE scheme can be applied with all methods to all properties that do not involve orbital-relaxed densities or orbital-relaxed derivatives, it is mainly interesting for ground-state energies, in particular for MP2 as PTE-MP2 is consistent with the trunction of the SCRF-FCI energy in second order. The change of the reaction field is, however, fully included for MP2 (orbital-relaxed) first-order properties and gradients which are available for the PTE scheme with COSMO and PE.

Required input data:

In addition to the input needed for the vacuum ricc2 calculation and the input needed to specify the environment (COSMO or PE) the only additional input for the PTE scheme is:

$reaction_field
   PTE

Availability:

Ground- and excited-state energies and orbital-unrelaxed properties with all methods available in ricc2. In addition orbital-relaxed first-order properties and gradients with COSMO and PE for MP2.

Limitations:

No orbital-relaxed properties for any other method than MP2 and Hartree-Fock (which is for closed-shell RHF and UHF identical with ground-state CCS).

10.9.2 The post-SCF coupling scheme

The post-SCF coupling includes at the correlated level correlation contributions to the ground-state reaction-field potential and determines the correlated wavefunction for \(\hat{H}^{vac}+\hat{G}(\mathbf{D}^{corr})\). Similar as for the PTE scheme, the calculations start from a converged SCRF Hartree-Fock calculation as reference. The scheme is non-iterative in the sense that correlation (and excitation) contributions to the reaction-field potential are not feeded back into the Hartree-Fock calculation.

In the ricc2 program the post-SCF scheme includes for excitation energies and spectra linear response (LR) contributions that are consistent with the respective ground-state Lagrangian with the additional approximation that the second derivative of the Lagrangian with respect to the Lagrange multipliers is neglected[283, 285]. Unlike the ground-state reaction field, the LR contributions are only evaluated with the fast (electronic) response of the environment. For COSMO this is done by using the screening factor \(f(n^2)\), where \(n\) is the refractive index, for PE the fast part is the (full) contribution from the environment’s polarizability, which models its electronic polarizability.

The coupling density, from which the correlation contributions to the reaction field is calculated, is for second-order methods by default truncated to avoid higher-order contributions that are not fully consistent with the second-order approximations and would break the computational efficiency. Available options for the coupling density are:

CCS-like:

the CC2 ground-state reaction field potential and the polarization energy are calculated from a singles-only approximation of the one-electron density \(D_{pq}^{\texttt{CCS-like}} = \langle \mathrm{HF}+\bar{T}_1| \exp(-T_1) a_p^\dagger a_q \exp(T_1)|\mathrm{HF}\rangle\): \[\begin{align} L^{\mathrm{CCS-like}}_{\mathrm{CC2}} = L^{vac}_{\mathrm{CC2}} + \langle \bar{T}_2 |[\hat{G}(\mathbf{D}^{\mathrm{HF}}),\hat{T}_2]|\mathrm{HF}\rangle + \tfrac{1}{2} \mathrm{Tr}\Big(\mathbf{D}^{\texttt{CCS-like}} \mathbf{G}(\mathbf{D}^{\texttt{CCS-like}})\Big) \end{align}\](10.28) The response contributions at the CC2 level are derived from this Lagrangian with the approximation that its second derivative with respect to \(\bar{T}_1\) is neglected (vide supra). This approximation was originaly proposed in the context of the polarizable embedding (PE) as PERI-CC2 in Ref.  and latter in Refs.  generalized to COSMO and ADC(2).

For MP2, the post-SCF scheme with the CCS-like density is equivalent to PTE-MP2. The Jacobian for CIS(D\(_\infty\)) is derived from that of CC2 by replacing the CC2 ground-state amplitudes by those from first-order perturbation theory. For ADC(2) the polarization (linear response) contributions to the vacuum-like secular matrix are derived by taking only the hermitian part of the CIS(D\(_\infty\)). The CIS(D) perturbative correction is in the post-SCF/CCS-like scheme related to the CIS(D\(_\infty\)) Jacobian in the same way as in the vacuum case. For CCS, the CCS-like density is not an approximation and the CCS excitation energies within the post-SCF coupling scheme are equivalent to CIS excitation energies with the linear response coupling scheme.

HF-T1sim:

the CC2 ground-state reaction field potential is calculated from the transition density between the Hartree-Fock and the CC wavefunction \(D_{pq}^{\texttt{HF-T1sim}} = \langle \mathrm{HF}| a_p^\dagger a_q \exp(T_1)|\mathrm{HF}\rangle\) and the CC2 Lagrangian defined as: \[\begin{align} L^{\mathrm{CCS-like}}_{\mathrm{CC2}} = L^{vac}_{\mathrm{CC2}} + \langle \bar{T}_2 |[\hat{G}(\mathbf{D}^{\mathrm{HF}}),\hat{T}_2]|\mathrm{HF}\rangle + \tfrac{1}{2} \mathrm{Tr}\Big(\mathbf{D}^{\texttt{CCS-like}} \mathbf{G}(\mathbf{D}^{\texttt{HF-T1sim}})\Big) \end{align}\](10.29) This approximation has been proposed in Refs. [287, 288] in the context of the polarizable embedding as sPERI-CC2. As part of the generalizations reported in Ref.  it can also be used with COSMO. For CC2, this approximation avoids the need to determine the ground-state Lagrange multipliers if only ground-state and excitation energies are requested, at the price of an increased asymmetry in the expression for the polarization energy.

CC2, ADC(2), CIS(D\(_\infty\)), and CIS(D) calculations with the post-SCF coupling scheme use per default the CCS-like approximation for the SCRF coupling contributions.

full:

includes for second-order methods all contributions needed to get excitation energies for singly excited states correct through second order, counting the correlation (fluctuation) contribution to the coulomb interaction with the environment as first-order as it is done for coulomb interaction between the electrons in the QM system.[281] It is—with the post-SCF coupling scheme—only available for the ADC(2) method as part of the work reported in Ref. [281]. For the MP2 ground state this approximation is equivalent to PTE-MP2 and, thus, also the ground-state amplitudes for ADC(2) are determined from first-order PTE-MP theory. For CCS, this approximation is equivalent to the CCS-like density since CCS does not include any second-order terms.

Perturbative state-specific corrections (ptSS)

In the post-SCF coupling scheme the reaction field is per construction equilibrated for the ground state, and only fully at the Hartree-Fock level while correlation effects are included only approximately, although the post-SCF scheme converges in the Full CI or Full CC limit for the ground-state to the equilibrium PTED Full CI limit, if no approximations are applied to the coupling density.

The perturbative state-specific corrections try to approximate the differences between the (non-iterative) post-SCF calculation and a state-specific equilibrium or non-equilibrium PTED calculation (vide infra) done at the same level and the same coupling density for a given reference state, which might be the ground or an excited state.

equilibrium case:

For the specified reference state, the equilibrium state-specific correction is calculated as the difference between the energy obtained in the post-SCF coupling scheme and the PTED free energy Lagrangian (for the same coupling density), evaluated with the wavefunction parameters obtaind from the post-SCF calculation. For CCS and MP2 or ADC(2), where no correlation effects to ground-state reaction field appear in the post-SCF scheme, the PTED Lagrangian is evaluated as: \[\begin{align} {\cal G}^{I,eq} & = {\cal G}^{eq}_{\mathrm{HF}} + E_{GS,corr} + \omega_I - \tfrac{1}{2} \mathrm{Tr}\Big(\mathbf{D}^{\mathrm{HF}} \mathbf{G}(\mathbf{D}^{\mathrm{HF}})\Big) \\ & - \mathrm{Tr}\Big(\big(\mathbf{D}^{I}-\mathbf{D}^{\mathrm{HF}}\big) \mathbf{G}(\mathbf{D}^{\mathrm{HF}})\Big) + \tfrac{1}{2} \mathrm{Tr}\Big(\mathbf{D}^{I} \mathbf{G}(\mathbf{D}^{I})\Big) \end{align}\](10.30) where \({\cal G}^{eq}_{\mathrm{HF}}\) is the SCRF-Hartree-Fock free energy, \(E_{GS,corr}\) and \(\omega_I\) are, respectively, the ground-state correlation and excitation energy for state \(I\) in the post-SCF coupling scheme. The sum of these three contributions are the uncorrected post-SCF total energies. \(\tfrac{1}{2} \mathrm{Tr}\Big(\mathbf{D}^{\mathrm{HF}} \mathbf{G}(\mathbf{D}^{\mathrm{HF}})\Big)\) the embedding energy (i.e. the interaction energy between QM system and environment plus the energy needed polarize the environment) at the HF level, and \(\tfrac{1}{2} \mathrm{Tr}\Big(\mathbf{D}^{I} \mathbf{G}(\mathbf{D}^{I})\Big)\) the same for the correlated and/or excited state. \(\mathrm{Tr}\Big(\big(\mathbf{D}^{I}-\mathbf{D}^{\mathrm{HF}}\big) \mathbf{G}(\mathbf{D}^{\mathrm{HF}})\Big)\) is the correlation contribution to the interaction energy between the correlated and/or excited wavefunction with the Hartree-Fock reaction field. These corrections are evaluated are evaluated[281] with the orbital-relaxed density \(\mathbf{D}^I\) for state \(I\).

non-equilibrium case:

For other states \(J\) than the reference state \(I\), the non-equilibrium free energy expression[289, 281] accounts for the fact the nuclear and the electronic relaxation and thus is polarization take place at different time scales. The assumption is that only the fast electronic polarization or reaction field \(\hat{G}^{fast}\) can follow perturbations at optical frequencies or electronic transitions while the slow nuclear polarization and reaction field \(\hat{G}^{slow} = \hat{G}-\hat{G}^{fast}\) remains equilibrated with a long living reference state \(I\). The free energy for a short living state \(J\) is under this assumption obtained as the expectation value for \(\hat{H}^{QM} + \hat{G}^{slow}(D^I) + \hat{G}^{fast}(D^J)\) with the wavefunction of state \(J\) plus the energies needed to induce the slow and fast components of the reaction field \(\hat{G}^{slow}(D^I)\) and \(\hat{G}^{fast}(D^J)\). For CCS and MP2 or ADC(2), where no correlation effects to the ground-state reaction field appear in the post-SCF scheme, the non-equilibrium free energy for a state \(J\) is evaluated as: \[\begin{align} {\cal G}^{J,neq} & = {\cal G}^{eq}_{\mathrm{HF}} + E_{GS,corr} + \omega_J - \tfrac{1}{2} \mathrm{Tr}\Big(\mathbf{D}^{\mathrm{HF}} \mathbf{G}(\mathbf{D}^{\mathrm{HF}})\Big) \\ & - \mathrm{Tr}\Big(\big(\mathbf{D}^{J}-\mathbf{D}^{\mathrm{HF}}\big) \mathbf{G}(\mathbf{D}^{\mathrm{HF}})\Big) + \mathrm{Tr}\Big(\mathbf{D}^J \mathbf{G}(\mathbf{D}^I)\Big) - \tfrac{1}{2} \mathrm{Tr}\Big(\mathbf{D}^{I} \mathbf{G}(\mathbf{D}^{I})\Big) \\ & + \tfrac{1}{2} \mathrm{Tr}\Big(\big(\mathbf{D}^J-\mathbf{D}^I\big)\big(\mathbf{G}^{fast}(\mathbf{D}^J-\mathbf{D}^I)\big)\Big) \end{align}\](10.31) As in the equilibrium case, the correction terms are evaluated with the orbital-relaxed densities.

For COSMO the fast part of the reaction field calculated with the screening factor \(f(n^2)\), where \(n\) is the refractive index, for PE the fast part of the reaction field is the (full) contribution from the environments polarizability.

Since the state-specific correction terms are evaluated with the orbital-relaxed densities, they are only printed out after the calculation of these densities and one-electron properties. The state-specific corrections effect the excitation energies. For CCS and ADC(2), where no correlation contribution to ground-state reaction field are included, the corrected energies for the transitions \(J \gets I\) become: \[\begin{align} \Delta {\cal G}_{J,I(eq)} & = {\cal G}^{J,neq} - {\cal G}^{I,eq} \\ & = \omega_J - \omega_I - \mathrm{Tr}\Big(\big(\mathbf{D}^{J}-\mathbf{D}^{I}\big) \big(\mathbf{G}(\mathbf{D}^{\mathrm{HF}}) - \mathbf{G}(\mathbf{D}^I)\big)\Big) \\ & + \tfrac{1}{2} \mathrm{Tr}\Big(\big(\mathbf{D}^J-\mathbf{D}^I\big)\big(\mathbf{G}^{fast}(\mathbf{D}^J-\mathbf{D}^I)\big)\Big) \end{align}\](10.32–10.33) The first correction term accounts for the fact the reaction field is equilibrated at the Hartree-Fock level instead for correlated and/or excited state \(I\), the second correction term for the fact that the fast (electronic) contribution to the reaction field should adapt instanteously to the transition to state \(J\).

Required input data:

For the post-SCF coupling scheme include in the control file:

$reaction_field
   post-SCF
   CCS-like

Availability:

The available functionalities depend on the method and approximation for the coupling density:

CCS-like

available for most energies and spectra

ground-state energies

CCS (equivalent to SCRF-HF), MP2 (equivalent to PTE-MP2), CC2

ground-state properties and gradients

CCS, MP2, CC2

excitation energies

CCS, CIS(D), CIS(D\(_\infty\)), ADC(2), CC2

excited-state properties and gradients

CCS, ADC(2)

one-photon TM

CCS, ADC(2), CC2

two-photon TM

CCS, CC2

induced TM

CCS, CC2

ptSS

MP2, CCS, ADC(2)

HF-T1sim

only developed for CC2 and 2c-CC2[288]

ground-state energies

CC2

excitation energies

CC2

one-photon TM

CC2

full

only developed for ADC(2)

excitation energies

ADC(2)

excited-state properties

ADC(2)

one-photon TM

ADC(2)

ptSS

MP2 (identical with CCS-like case), ADC(2)

Limitations:

  • Not yet available for linear and quadratic response properties

  • The “full” coupling density is not yet available for LT-SOS-ADC(2) and not yet for relativistic 2c-ADC(2) calculations.

  • ptSS corrections are only available for COSMO and not yet available for CC2 and not for relativistic 2-component calculations.

  • Ground- and excited-state properties and gradients are not available for relativistic 2-component calculations.

  • Gradients are only available for the uncorrected post-SCF energies.

10.9.3 The PTED coupling scheme

In the PTED coupling scheme, the reaction field potential is determined fully self-consistently with the density of a reference state which can be the ground or an excited state. This is achieved by macroiterations over the Hartree-Fock and the correlated/excitation calculations:

  1. Iteration 0:

    1. SCRF-Hartree-Fock calculation (dscf)

    2. correlated/response calculations with ricc2

    3. calculate RF potential for correlated ground- or excited-state density

  2. Iteration 1–n:

    1. Hartree-Fock calculation with fixed RF potential (dscf)

    2. correlated/response calculations with fixed RF potential (ricc2)

    3. calculate RF potential for correlated ground- or excited-state density

  3. energy and RF potential converged? If not, continue with step 10.9, if yes exit SCRF loop.

  4. Final cycle: If requested or needed, a final cycle is done to compute with the converged reaction field additional excitation energies and properties as e.g. gradients

These SCRF macroiterations can be done for COSMO with the script cc2cosmo and for PE with the script pecc2.

The PTED scheme implemented in ricc2 orginates to large parts from the work in Ref.  on PTED-COSMO-CCS and PTED-COSMO-ADC(2). It includes for excitation energies by default linear response contributions with the “full” coupling density. These can be switched off by setting the option nofast in $reaction_field.

The SCRF iteractions make the free energy Lagrangian for the specified reference state stationary. After convergence, the final equilibrium free energy for the refernce state is then obtained as \[\begin{align} {\cal G}^{I,eq} = {\cal G}^{neq}_{\mathrm{HF}} + E_{GS,corr} + \omega_I \end{align}\](10.34) without further corrections. \({\cal G}^{neq}_{\mathrm{HF}}\) is non-equilibrium Hartree-Fock free energy for the reaction field \(\mathbf{G}(\mathbf{D}^I)\). No further corrections are needed for the equilibrium case.

For those states for which the reaction field is not equilibrated non-equilibrium corrections[289, 281] are calculated and included in the final results for total and transition energies that take into account that the fast part of the reaction field should adapt instanteously to electronic transition from state \(I\) to state \(J\): \[\begin{align} {\cal G}^{J,neq} & = {\cal G}^{neq}_{\mathrm{HF}} + E_{GS,corr} + \omega_J + \tfrac{1}{2} \mathrm{Tr}\Big(\big(\mathbf{D}^J-\mathbf{D}^I\big)\big(\mathbf{G}^{fast}(\mathbf{D}^J-\mathbf{D}^I)\big)\Big) \end{align}\](10.35) Since the correction is evaluated from orbital-relaxed densities \(\mathbf{D}^J\) and \(\mathbf{D}^I\), the corrected total and transition energies are only evaluated and printed out after the calculation of the densities and one-electron properties.

The non-equilibrium correction effects the excitation energies: \[\begin{align} \Delta {\cal G}_{J,I(eq)} & = {\cal G}^{J,neq} - {\cal G}^{I,eq} = \omega_J - \omega_I + \tfrac{1}{2} \mathrm{Tr}\Big(\big(\mathbf{D}^J-\mathbf{D}^I\big)\big(\mathbf{G}^{fast}(\mathbf{D}^J-\mathbf{D}^I)\big)\Big) \end{align}\](10.36)

Required input data:

For the PTED coupling scheme include in the control file:

$reaction_field
   scrf state=(bu 1)

All further input options are set by the PTED scripts cc2cosmo and pecc2.

Availability:

The PTED coupling scheme is available for those methods (and states) for which (orbital-relaxed) one-electron densities are available:

ground state

MP2, CC2

excited states

CCS, ADC(2), CC2

Limitations:

  • Because of the macroiterations for the solution of SCRF equations, PTED calculations are not compatible with a loop over different methods within ricc2 (not counting MP2 when done as ground-state for ADC(2) or the automatic MP2 and CCS steps for the generation of start guesses).

  • Not available for relativistic 2-component calculations.

Gradients for the PTED coupling scheme

Gradients for the PTED coupling scheme are available for MP2 ground-state, ADC(2) excited-state and CC2 ground- and excited-state energies as described in Ref. . For geometry optimizations with jobex the latter as to be used the level option set to

cc2cosmo

for PTED-COSMO

pecc2

for PTED-PE

For geometry optimizations with the PTED coupling scheme it is important that the threshold for the convergence of the PTED cycles is tight enough to avoid numerical noise in the gradients. For options to control the convergence of the PTED cycle, call the PTED scripts with the option -h.

For COSMO it is strongly recommended to use Gaussian Charge Model for geometry optimizations and a sufficiently large grid to minimize numerical problems related to the incomplete rotational invariance of the cavity representation.