21.2 Treatment of Solvation Effects with Cosmo

The Conductor-like Screening Model [280] (Cosmo) is a continuum solvation model (CSM), where the solute molecule forms a cavity within the dielectric continuum of permittivity \(\varepsilon\) that represents the solvent. The charge distribution of the solute polarizes the dielectric medium. The response of the medium is described by the generation of screening charges on the cavity surface.

CSMs usually require the solution of the rather complicated boundary conditions for a dielectric in order to obtain the screening charges. Cosmo instead uses the much simpler boundary condition of vanishing electrostatic potential for a conductor,

\[\begin{align*} {\bf \Phi}^{tot} = 0. \end{align*}\]

This represents an electrostatically ideal solvent with \(\varepsilon=\infty\). The vector of the total electrostatic potential on the cavity surface segments is determined by the solute potential \({\bf \Phi}^{sol}\), which consist of the electronic and the nuclear part, and the vector of the screening charges \(\bf q\), \[\begin{align*} {\bf \Phi}^{tot} = {\bf \Phi}^{sol}+ {\bf Aq} = 0. \end{align*}\] \(\bf A\) is the Coulomb matrix of the screening charge interactions. For a conductor, the boundary condition \({\bf \Phi}^{tot}=0\) defines the screening charges as

\[\begin{align*} {\bf q } & = -{\bf A}^{-1}{\bf \Phi}^{sol}. \end{align*}\]

To take into account the finite permittivity of real solvents, the screening charges are scaled by a factor.

\[\begin{align*} f(\varepsilon) &= \frac{\varepsilon-1}{\varepsilon+\frac{1}{2}}\\ {\bf q}^{\star} &= f(\varepsilon){\bf q } \end{align*}\]

The deviation between the Cosmo approximation and the exact solution is rather small. For strong dielectrics like water it is less than 1%, while for non-polar solvents with \(\varepsilon\approx 2\) it may reach 10% of the total screening effects. However, for weak dielectrics, screening effects are small and the absolute error therefore typically amounts to less than one kcal/mol.

As shown in [354] ions can be described more accurately by using the scaling factor \(f(\varepsilon) = \frac{\varepsilon-1}{\varepsilon + 0}\). This also leads to results (e.g. by comparing solvation free energies) more or less identical to IEFPCM or SS(V)PE. This is possible since TURBOMOLE version 7.1 by adding the keyword ions to the epsilon option in the $cosmo section (see 25.2.14).

The dielectric energy, i.e. the free electrostatic energy gained by the solvation process, is half of the solute-solvent interaction energy.

\[\begin{align*} E_{diel} = \frac{1}{2}f(\varepsilon){\bf q}^{\dagger}{\bf \Phi}^{sol} \end{align*}\]

The total free energy of the solvated molecule is the sum of the energy of the isolated system calculated with the solvated wave function and the dielectric energy

\[\begin{align*} E=E(\Psi^{solv})+E_{diel}. \end{align*}\]

A Cosmo energy calculation starts with the construction of the cavity surface grid. Within the SCF procedure, the screening charges are calculated in every cycle and the potential generated by these charges is included into the Hamiltonian. This ensures a variational optimization of both the molecular orbitals and the screening charges and allows for the evaluation of analytic gradients.

How to switch on COSMO:

To enable COSMO with the default dielectric constant \(\varepsilon=\infty\) it is sufficient to add the keyword $cosmo to the control file.

For a specific solvent either add the \(\varepsilon\) value (and refractive index if needed for optical properties), or use the solvent=solvent-name as option from the list of available solvents. There are about 100 solvents available, see the file parameter/cosmosolvents in your TURBOMOLE  installation.

For more options and to use non-default settings please see section 25.2.14 for a list of keywords and options.

Gaussian Charge Model COSMO-GCM:

In order to generate smooth and reliable potential energy surfaces, the Gaussian charge model (GCM) for COSMO is highly recommended. Especially geometry optimizations and vibrational frequency calculations using analytical 2nd derivatives with aoforce  show less numercial errors and fewer imaginary frequencies. COSMO-GCM constructs the Coulomb matrix \(\bf A\) using Gaussian charges and Lebedev grids. This leads to consistent analytical gradients as well as geometry Hessians. COSMO-GCM can be enabled using an addition to the COSMO keyword:

$cosmo
  gauss
  nleb= 3

Lebedev grids between 1 (small) and 9 (very large) can be used. Default settings use nleb= 3, if this information is omitted in the control file, this grid is chosen. COSMO-GCM can also be enabled using the cosmoprep module.

Radii based Cavity Construction:

In order to ensure a sufficiently accurate and efficient segmentation of the molecular shaped cavity the COSMO implementation uses a double grid approach and segments of hexagonal, pentagonal, and triangular shape. The cavity construction starts with a union of spheres of radii \(R_{i}+RSOLV\) for all atoms \(i\). In order to avoid problems with symmetric species, the cavity construction uses de-symmetrized coordinates. The coordinates are slightly distorted with a co-sinus function of amplitude AMPRAN and a phase shift PHSRAN. Initially a basis grid with NPPA segments per atom is projected onto atomic spheres of radii \(R_{i}+RSOLV\). In order to avoid the generation of points in the problematic intersections, all remaining points, which are not in the interior of another sphere, are projected downwards onto the radius \(R_{i}\). In the next step a segment grid of NSPH segments per H atom and NSPA segments for the other atoms is projected onto the surface defined by \(R_{i}\). The basis grid points are associated to the nearest segment grid centers and the segment coordinates are re-defined as the center of area of their associated basis grid points, while the segment area is the sum of the basis grid areas. Segments without basis grid points are discarded. In order to ensure nearest neighbor association for the new centers, this procedure is repeated once. At the end of the cavity construction the intersection seams of the spheres are paved with individual segments, which do not hold associated basis grid points.

Density based Cavity Construction:

Instead of using atom specific radii the cavity can be defined by the electron density. In such an isodensity cavity construction one can use the same density value for all atoms types or the so-called scaled isodensity values. In the later approach different densities are used for the different atom types. The algorithm implemented in Turbomole uses a marching tetrahedron algorithm for the density based cavity construction. In order to assure a smooth density change in the intersection seams of atoms with different isodensity specification, this areas are smoothened by a radii based procedure.

Radii based Isosurface Cavity:

A cavity construction algorithm based on the triangulation of an iso-surface is available as an alternative to the radii or density based construction. It overcomes deficiencies which have become apparent for the original COSMO standard cavity, especially in concave regions of the molecular shaped cavity. The new construction, called FINE Cavity, is described in details in [355].

To enable the new radii based isosurface cavity the keyword $cosmo_isorad has to be added to the control file.

\({\bf A}\)-Matrix Setup:

The \({\bf A}\) matrix elements are calculated as the sum of the contributions of the associated basis grid points of the segments \(k\) and \(l\) if their distance is below a certain threshold, the centers of the segments are used otherwise. For all segments that do not have associated basis grid points, i.e. intersection seam segments, the segment centers are used. The diagonal elements \(A_{kk}\) that represent the self-energy of the segment are calculated via the basis grid points contributions, or by using the segment area \(A_{kk} \approx 3.8 \sqrt{a_k}\), if no associated basis grid points exist.

Outlying charge correction:

The part of the electron density reaching outside the cavity causes an inconsistency that can be compensated by the "outlying charge correction". This correction will be performed at the end of a converged SCF or an iterative MP2 calculation and uses an outer surface for the estimation of the energy and charge correction [356]. The outer surface is constructed by an outward projection of the spherical part of the surface onto the radius \(R_{i}+ROUTF*RSOLV\). It is recommended to use the corrected values.

Numerical Frequency Calculation:

The calculation of harmonic frequencies raises the problem of non-equilibrium solvation in the Cosmo framework, because the molecular vibrations are on a time scale that do not allow a re-orientation of the solvent molecules. Therefore, the total response of the continuum is split into a fast contribution, described by the electronic polarization, and a slow term related to the orientational relaxation. As can be shown [357] the dielectric energy for the disturbed state can be written as \[\begin{align*} E^d_{diel}=\frac{1}{2}f(\varepsilon){\bf q}({\bf P}^0){\bf \Phi}({\bf P}^0) +\frac{1}{2}f(n^2){\bf q}({\bf P}^\Delta){\bf \Phi}({\bf P}^\Delta) + f(\varepsilon){\bf q}({\bf P}^0){\bf \Phi}({\bf P}^\Delta), \end{align*}\] where \({\bf P}^\Delta\) denotes the density difference between the distorted state and the initial state with density \({\bf P}^0\). The interaction is composed of three contributions: the initial state dielectric energy, the interaction of the potential difference with the initial state charges, and the electronic screening energy that results from the density difference. The energy expression can be used to derive the correspondent gradients, which can be applied in a numerical frequency calculation. Because the Cosmo cavity changes for every distorted geometry the initial state potential has to be mapped onto the new cavity in every step. The mapped potential of a segment of the new cavity is calculated from the distance-weighted potentials of all segments of the old cavity that fulfill a certain distance criterion. The mapped initial state screening charges are re-calculated from the new potential.

21.2.1 Vertical excitations and Polarizabilities for TDDFT, TDA and RPA:

The escf program accounts for the Cosmo contribution to the excitation energies and polarizabilities. The Cosmo settings have to be defined for the underlying Cosmo dscf or ridft calculation. In case of the excitation energies the solvent response will be divided into the so-called slow and fast term [357, 358]. The screening function of the fast term depends on the refractive index of the solvent which can be defined in the input. If only the Cosmo influence on the ground state should be taken into account we recommend to perform a normal Cosmo calculation (dscf or ridft) and to switch off Cosmo (i.e. deactivate $cosmo) before the escf calculation.

21.2.2 The Direct COSMO-RS method (DCOSMO-RS):

In order to go beyond the pure electrostatic model a self consistent implementation of the COSMO-RS model the so-called "Direct COSMO-RS" (DCOSMO-RS) [359] has been implemented in ridft and dscf.

COSMO-RS (COSMO for Real Solvents) [360, 361] is a predictive method for the calculation of thermodynamic properties of fluids that uses a statistical thermodynamics approach based on the results of COSMO SCF calculations for molecules embedded in an electric conductor, i.e. using \(f(\varepsilon)=1\). The liquid can be imagined as a dense packing of molecules in the perfect conductor (the reference state). For the statistical thermodynamic procedure this system is broken down to an ensemble of pair wise interacting surface segments. The interactions can be expressed in terms of surface descriptors. e.g. the screening charge per segment area (\(\sigma_t = q_t / a_t\)). Using the information about the surface polarity \(\sigma\) and the interaction energy functional, one can obtain the so-called \(\sigma\)-potential (\(\mu_S(\sigma;T)\)). This function gives a measure for the affinity of the system \(S\) to a surface of polarity \(\sigma\). The system \(S\) might be a mixture or a pure solvent at a given temperature \(T\). Because the parabolic part of the potential can be described well by the Cosmo model, we subtract this portion form the COSMO-RS potential:

\[\begin{align*} \tilde{\mu}_S(\sigma;T) = \mu_S(\sigma;T)-(1-f(\varepsilon))c_0\sigma^2. \end{align*}\]

The parameter \(c_0\) can be obtained from the curvature of a COSMO-RS \(\sigma\)-potential of a nonpolar substance, e.g. hexane.

Thus, the remaining part of the chemical potential of a compound \(i\) with mole fraction \(x_i\) in the mixture \(S\)i can be expressed as:

\[\begin{align*} \mu^i \cong \sum^m_{t=1} f_{pol} a_t \tilde{\mu}_S(\sigma;T)+\mu^i_{C,S}+kT\ln(x_i). \end{align*}\]

where the combinatorial term \(\mu^i_{C,S}\) accounts for effects due to the size and shape differences of the molecules in the mixture and \(a_t\) denotes the area of segment \(t\). The \(kT\ln(x_i)\) can be skipped for infinite dilution. The factor \(f_{pol}\) has been introduced to account for the missing solute-solvent back polarization. The default value is one in the current implementation. The free energy gained by the solvation process in the DCOSMO-RS framework is the sum of the dielectric energy of the Cosmo model and the chemical potential described above:

\[\begin{align*} E_{diel,RS}= \frac{1}{2} f(\varepsilon) {\bf q}^{\dagger}{\bf \Phi}^{sol}+ \mu^i = E_{diel}+\mu^i \end{align*}\]

From the above expression the solvent operator \(\hat{V}^{RS}\) can be derived by functional derivative with respect to the electron density:

\[\begin{align*} \hat{V}^{RS} = - \sum^m_{t=1} \frac{f(\varepsilon) q_t+q^{{\Delta}RS}_t}{\lvert {\bf r}_t -{\bf r} \rvert}= \hat{V}^{cos}-\sum^m_{t=1} \frac{q^{{\Delta}RS}_t}{\lvert {\bf r}_t -{\bf r} \rvert}. \end{align*}\]

Thus, the solvation influence of the COSMO-RS model can be viewed as a correction of the Cosmo screening charges \(q_t\). The additional charges denoted as \(q^{{\Delta}RS}_t\) can be obtained from \({\bf q}^{{\Delta}RS} = -{\bf A}^{-1}{\bf \Phi}^{{\Delta}RS}\), where the potential \({\bf \Phi}^{{\Delta}RS}\) arises from the chemical potential of the solute in the solvent:

\[\begin{align*} \phi^{{\Delta}RS}_t = a_t \left( \frac{\delta \tilde{\mu}_S}{\delta q}\right)_{q=q_t}. \end{align*}\]

In order to get a simple and differentiable representation of the COSMO-RS \(\sigma\)-potential \(\mu_S(\sigma;T)\), we use equally spaced cubic splines.

An approximate gradient of the method has been implemented. DCOSMO-RS can be used in SCF energy and gradient calculations (geometry optimizations) with dscf, ridft, grad, and rdgrad. Please regard the restriction of the DCOSMO-RS energy explained in the keyword section 25.2.14. Because the DCOSMO-RS contribution can be considered as a slow term contribution in vertical excitations it does not have to be taken into account in response calculations. For the calculation of vertical excitation energies it is recommended to use the mos of a DCOSMO-RS calculation in a Cosmo response calculation (see above).

21.2.3 COSMO-MP2 with mpgrad

For MP2 calculations within the CSM framework three alternatives can be found in the literature [362]. The first approach, often referred to as PTE, performs a normal MP2 energy calculation on the solvated HF wave function. The electrostatic potential due to the response of the solvent to the solute, also called reaction field, is still on the HF level. It is the only of the three approaches that is formally consistent in the sense of second-order perturbation theory [363, 364]. In the so-called PTD approach the vacuum MP2 density is used to calculate the reaction field. The PTD approach is not implemented in TURBOMOLE. The third approach, often called PTED, is iterative and equilibrated the reaction field with the (orbital-relaxed) MP2 density so that the reaction field reflects the correlated density. In contrast to the PTE approach, the reaction field, i.e. the screening charges, change in the PTED calculation during the iterations until self consistency is reached.

The implementation in mpgrad is limited to PTE and PTED. Gradients are only available on the formally consistent PTE level [365] with radii based cavity construction, i.e., not for PTED and not in combination with an isodensity ($cosmo_isodens) or radii based isosurface ($cosmo_isorad) cavity.

COSMO-MP2 calculations with mpgrad do not require any special input beyond what is needed for COSMO-HF. It will be switched on automatically if the ($cosmo) data group is found in the control file. By default a PTE-COSMO-MP2 energy and gradient calculation will be done when mpgrad is invoked after a COSMO-HF calculation. Geometry optimizations for PTE-COSMO-MP2 can be done with jobex.

Self-consistent PTED-COSMO-MP2 calculations can be carried with the script mp2cosmo. Use mp2cosmo -h for available options.

21.2.4 COSMO in combination with the ricc2 program

The following table shows the avaliability of COSMO-CCS (or COSMO-CIS), COSMO-MP2, COSMO-ADC(2), and COSMO-CC2 with the PTED and post-SCF coupling schemes for different functionalities in the ricc2 program.

The combination of COSMO-CC2 and COSMO-ADC(2) with SCS and SOS are for some functionalities also available in the ricc2 program. Furthermore, the calculations of all of the above listed features are possible for open-shell systems and also with SMP (OpenMP) and MPI parallelization.

COSMO-MP2 and COSMO-CC2 for ground-state calculations

Post-SCF coupling scheme

Ground state energy, orbital-unrelaxed and relaxed densities and properties, and gradients are available for COSMO-CC2 within the post-SCF[283] reaction-field scheme with a CCS-like approximation for the density that is used to determine a correlation contribution to the reaction field[283]. In the post-SCF scheme, the ground-state reaction field is first determined self-consistently at the COSMO-HF level. In the subsequent correlation treatment, the correlation effect to the reaction field is calculated from the correlation contribution to the unrelaxed density and included in the equations for the ground-state wavefunction parameters. In order to set the input file for the calculation of ground-state gradient, energy and relaxed properties the following data groups must be included in the control file.

$reaction_field
  post-SCF
  ccs-like
$cosmo
  epsilon=   50.000
  rsolv= 1.30
$cosmo_atoms
...

For MP2, the post-SCF coupling scheme with the singles-only (CCS-like) coupling density is equivalent to the PTE scheme. Ground-state energies, orbital-relaxed densities and one-electron properties, and gradients are available for COSMO-MP2 with the PTE (or post-SCF) coupling scheme.

PTED coupling scheme

Ground-state PTED-COSMO-MP2 and PTED-COSMO-CC2 calculations can be done with the script cc2cosmo. Before starting the script, the control file has to contain the settings needed for the $cosmo and $ricc2 data groups. Additional input in the $reaction_field data group can be specified, but is not needed if the default values set by cc2cosmo are sufficient:

$reaction_field
   PTED
   econv=1.d-$econv        # by default $econv=6
   qmaxconv=1.d-$qmconv    # by default $qmconv=5
   qrmsconv=1.d-$qconv     # by default $qconv=6

Note, that by default the reaction field is equilibrated with the (orbital-relaxed and correlated) ground-state density. Orbital-relaxed (and for CC2 also unrelaxed) one-electron densities are one-electron properties are obtained as by-product of the PTED calculation. Also ground-state gradients are available through the usual input for the ricc2 program. With the setting

$ricc2
   geoopt model=mp2

ground-state geometry optimizations for PTED-COSMO-MP2 and PTED-COSMO-CC2 can be done with

   jobex -level cc2cosmo

Note, that depending on the convergence threshold for the geometry, the convergence thresholds for the PTED self-consistent reaction field cycle (options econv, qmconv, and qconv for script cc2cosmo) need to be tightened.

21.2.5 Solvation effects on excited states using COSMO in ricc2:

The COSMO approach has been recently implemented into the ricc2 module of TURBOMOLE for excitation energies with CC2 and ADC(2). The ADC(2) method has been implmented in combination with the iterative PTED and the more economical post-SCF reaction field schemes. CC2 is currently only available in combination with the post-SCF reaction field scheme.

Iterative COSMO-ADC(2) within the PTED(LR) scheme

With the PTED reaction field (RF) scheme it is possible to equilibrate the solvent charges for the ground state at MP2 or any excited state. The implementation of the PTED scheme for electronically excited states in the ricc2 program goes back to Ref. [281], where an extended PTED scheme which includes also linear response (LR) and non-equilibrium corrections for the non-equilibrated states has been presented for CCS (or CIS) and ADC(2) excitation energies. Non-equilibrated means in this sense, that the slow part of the solvent charges (described by \(f(\varepsilon)\)) are equilibrated with a given initial state, while the fast electronic part of the solvent charges (described by \(f(n^2)\)) are in equilibrium with a target state different from the equilibrated initial state. This approach has later been extended to CC2 excitation energies and excited-state gradients.[366].

Calculations for this extended PTED scheme require similar as PTED-MP2 self-consistency iterations over a reference Hartree-Fock and an excitation energy calculation with ricc2. These macro iterations can be managed with the script cc2cosmo which is similar to mp2cosmo but uses ricc2 instead of mpgrad and can handle a gradient calculation after convergence of the SCRF iterations.

The input for calculations with the PTED(LR) scheme has to contain:

• The data groups $ricc2 and $excitations with input for electronic structure method (CC2, ADC(2), or CCS) and for gradients/geometry optimizations in addition the option geoopt in $ricc2.

• The data group $cosmo with the input for the COSMO model including an input for the refractive index (either via solvent or the refind option). The latter can be skipped if the fast LR contributions are switched off (vide infra).

• The data group $reaction_field with the option scrf state for defining the state for which (the slow part of) the reaction field should be equilibrated. (The short-hand notation (s1) is in this case not recognized.) If not specified, the reaction field will be equilibrated for the ground state. The linear response or “fast” contributions mentioned above can be switched off by setting the option ’nofast’ in $reaction_field.

A typical input might look like:

$ricc2
  adc(2)
$excitations
  irrep=a'  multiplicity=1  nexc=2  npre=4  nstart=8
  irrep=a"  multiplicity=1  nexc=2  npre=4  nstart=8
$cosmo
  solvent=acetonitrile
$reaction_field
  scrf state=(a" 1)
  # nofast

This would request an excitation energy calculation for the lowest two singlet A’ and A" excitations using the ADC(2) method. The solvent charges are equilibrated for state \(1^1A"\) and the non-equilibrium energy contributions for the MP2 ground state and the other three excited states are calculated in addition to the energy of the \(1^1A"\) reference state. One-electron properties for all states will be obtained as side result of the calculation.

When planning PTED calculations one should be aware of that:

• Due to the SCRF cycles and the need to compute in each iteration the (orbital-relaxed correlated) ground- or excited-state density, these calculations take a factor of 3-20 longer than a respective gas phase calculation.

• The state-specific contributions can within the variational iterative optimization lead to a change in the energetic ordering of the states including a flip between the ground and an excited state. The SCRF iterations might then not converge or converge to unphysical solutions.

• The ricc2 program updates with each invocation the charges on the cavity surface which are stored on a file and used for the next iteration. As a consequence, ricc2 can not be restarted without first rerunning the Hartree-Fock calculation.

COSMO-ADC(2) within the post-SCF scheme:

COSMO-ADC(2) with the post-SCF reaction–field scheme enables the calculations of vertical excitations energies,transition moments (TM),excited state first-order properties, and analytic gradients.[282] Unlike the PTED scheme, the COMSO-ADC(2) method with post-SCF does not need macro iterations over the Hartree-Fock and the correlated calcualtion to obtain vertical excitation energies, meaning that after solving the ground–state reaction field self-consistently at the COSMO-HF level the program computes linear response reaction-field contributions to excitation energies without coupling excited-state or correlated contributions back to the ground–state HF.

COSMO-CC2 within the post-SCF scheme:

The implementation of COSMO-CC2 within the post-SCF reaction field enables the calculations of vertical excitation energies, transition moments and Faraday \(\cal{B}\) terms for the simulation of UV-Vis and magnetic circular dichroism (MCD) absorption spectrum.[262]

A typical input for COSMO-CC(2) and COSMO-ADC(2) any excited-state calculations with the post-SCF scheme

include following data groups in addition to the typical COSMO data groups:

$reaction_field
  post-SCF
  ccs-like
$cosmo
  solvent=acetone

The input for other data groups like $ricc2, $excitations, etc. is the same as without COSMO.