Treatment of Solvation Effects with COSMO

The Conductor-like Screening Model [190] (COSMO) is a continuum solvation model (CSM), where the solute molecule forms a cavity within the dielectric continuum of permittivity ε 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,

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.

A-Matrix Setup:
The 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} ≈ 3.8, 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 [191]. 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 [192] the dielectric energy for the disturbed state can be written as

Iterative MP2 COSMO: For ab initio MP2 calculations within the CSM framework three alternatives can be found in the literature [193]. The first approach, often referred to as PTE, performs a normal MP2 energy calculation on the solvated HF wave function. The response of the solvent, 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 [194,195]. In the so-called PTD approach the vacuum MP2 density is used to calculate the reaction field. The third approach, often called PTED, is iterative so that the reaction field reflects the density of the first-order wave function. In contrast to the PTE approach the reaction field, i.e. the screening charges, change during the iterations until self consistency is reached. Gradients are available on the formally consistent PTE level only [196].

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 [192,197]. 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.

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) [198] has been implemented in ridft and dscf.

COSMO-RS (COSMO for Real Solvents) [199,200] 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(ε) = 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 (σ_{t} = q_{t}∕a_{t}). Using the information
about the surface polarity σ and the interaction energy functional, one can obtain the so-called
σ-potential (μ_{S}(σ;T)). This function gives a measure for the affinity of the system S
to a surface of polarity σ. 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:

Solvation effects on excited states using COSMO in ricc2:
The COSMO approach has been recently implemented into the ricc2 module of TURBOMOLE. It
is now possible to equilibrate the solvent charges for any excited state or the ground state as
mentioned in the MP2 section above. Using the methods CCS/CIS or ADC(2) the implementation
is complete, for CC2 or higher methods, however it still has to be proven if there are terms
missing The recent implementation contains contributions to the off-diagonal elements of
the one-electron density. Furthermore the energy contributions for non-equilibrated
states can be calculated. Non-equilibrated means in this sense, that the slow part of the
solvent charges (described by f(ε)) are still equilibrated with a given initial state, while
the fast electronic part of the solvent charges (described by f(n^{2})) are in equilibrium
with the target state. To handle this one has to do a macro iteration, like in MP2.
This macro iteration can be managed with the script ’cc2cosmo’ which is the same as
’mp2cosmo’ but using ricc2 instead of rimp2 or mpgrad. To set up the basic settings
one can reuse the cosmoprep module, note that one has to specify the refractive index
’refind’ when observing excited states. To specify the state to which the solvent charges
should be equilibrated one inserts the keyword ’cosmorel state=(x)’, where the ground
state (x) is used normally but can be replaced by any requested excited state. Make
sure to request relaxed properties for any desired state, otherwise the COSMO macro
iteration will not work in ricc2. The off-diagonal contributions mentioned above can be
switched off by setting the keyword ’nofast’ in $cosmo. A typical input might look
like:

$ricc2

adc(2)

$excitations

irrep=a’ multiplicity= 1 nexc= 1 npre= 1 nstart= 1

irrep=a" multiplicity= 1 nexc= 1 npre= 1 nstart= 1

exprop relaxed states=all

$response

fop relaxed

$cosmo

epsilon= 50.000

rsolv= 1.30

refind= 3.0000

cosmorel state=(a" 1)

# nofast

$cosmo_correlated

$cosmo_atoms

...

adc(2)

$excitations

irrep=a’ multiplicity= 1 nexc= 1 npre= 1 nstart= 1

irrep=a" multiplicity= 1 nexc= 1 npre= 1 nstart= 1

exprop relaxed states=all

$response

fop relaxed

$cosmo

epsilon= 50.000

rsolv= 1.30

refind= 3.0000

cosmorel state=(a" 1)

# nofast

$cosmo_correlated

$cosmo_atoms

...

This would deliver an excited state calculation for the lowest singlet A’ and A" excitations using
the ADC(2) method. The solvent charges are equilibrated to state 1^{1}A" and the non-equilibrium
energy contributions for the MP2 ground state and the 1^{1}A’ excited state are calculated
furthermore. All contributions to the one-electron density are included since the proper keyword is
commented out. Note: when doing solvent relaxations with the CCS/CIS model, no request of
relaxed ground–state properties are needed, since the relaxed ground state is identical to the HF
ground state. A short summary of the COSMO input is given at the beginning of the ricc2
output as well as a summary of the energy contributions is given at the end of the ricc2
output.