25.2.14 Keywords for COSMO

The Conductor-like Screening Model (Cosmo) is a continuum solvation model, where the solute molecule forms a cavity within the dielectric continuum of permittivity epsilon that represents the solvent. A brief description of the method is given in chapter 21.2. The model is currently implemented for SCF energy and gradient calculations (dscf/ridft and grad/rdgrad), MP2 energy calculations (rimp2 and mpgrad) and MP2 gradients (rimp2), vibrational frequencies using aoforce and response calculations with escf. The ricc2 implementation is described in section 21.2.5.

For simple HF or DFT single point calculations or optimizations with standard settings, we recommend to add the $cosmo keyword to the control file and to skip the rest of this section.

Please note: due to improvements in the ${\bf A}$ matrix and cavity setup the Cosmo energies and gradients may differ from older versions (5.7 and older). The use_old_amat option can be used to calculate energies (not gradients) using the old cavity algorithm of Turbomole 5.7.
The basic Cosmo settings are defined in the $cosmo and the $cosmo_atoms block.
Example with default values:

$cosmo
  epsilon=infinity
  refind= 1.3
  nppa= 1082
  nspa=   92
  disex= 10.0000
  rsolv= 1.30
  routf= 0.85
  cavity closed
  ampran= 0.1D-04
  phsran= 0.0
# the following options are not used by default
  solvent=not set
  gauss
  nleb= 3
  klamt
  isorad
  allocate_nps= 140 
  use_old_amat 
  use_contcav
  no_oc
  ccf
  force_cosmowrite
  name= out.cosmo

epsilon=: defines a finite permittivity used for scaling of the screening charges. If the option ion is added to the same input line, the scaling factor for ions $f(\varepsilon) = \frac{\varepsilon-1}{\varepsilon + x}$ with $x=0$ will be used. Alternatively the x value can be set by adding ion=x with x as a real value.
refind=: refractive index used for the calculation of vertical excitations and num. frequencies (the default 1.3 will be used if not set explicitly)
solvent=: requests to use the values for epsilon and refind for a given solvent named from the cosmosolvent list. E.g., for water, the values epsilon=80.10 and refind=1.3334 are set. Note that additional settings for epsilon and refind (see above) after solvent will overwrite these values.
allocate_nps=: skips the Cosmo segment statistics run and allocates memory for the given number of segments.
no_oc: skips the outlying charge correction.
use_old_amat: uses ${\bf A}$ matrix setup of TURBOMOLE 5.7
use_contcav: in case of disjunct cavities only the largest contiguous cavity will be used and the smaller one(s) neglected. This makes sense if an unwanted inner cavity has been constructed e.g. in the case of fullerenes. Default is to use all cavities.
point_charges: Allows the COSMO response to point charges. By default it is disables, so that point charges interacts only with the molecule. Note that when point charges are used together with COSMO, the COSMO surface is still constructed considering only the atoms.
gauss: switches to Gaussian Charge Model COSMO-GCM which uses Gaussian charges and Lebedev grids for the cavity construction. Recommended for geometry optimizations and vibrational frequency calculations as well as response properties. See also section 21.2.
nleb=: selects the grid size for the Lebedev cavity construction, valid values are 1 (small) to 9 (very large). Default is nleb= 3 if gauss option is used.
klamt: switches back to the pre-V8.0 default for the cavity construction and the usage of point charges.
isorad: switches on the radii based isosurface cavity (also named FINE cavity), this option is equivalent to adding $cosmo_isorad keyword, see section 21.2.
ccf: writes out compressed cosmo file format (ccf) instead of the default plain text cosmo format. Setting a name of the resulting cosmo/ccf file using cosmo_out file= including the .cosmo or .ccf appendix overrules this option.
force_cosmowrite: enforces generation of a cosmo file even in case a calculation does not converge.
name=: sets the name of the resulting cosmo or ccf file. This is an alternative to setting the $cosmo_out keyword. In case both this option and the keyword is set for the name, the cosmo_out settings are used.

All other parameters affect the generation of the surface and the construction of the $\mathbf{A}$ matrix:

nppa=: number of basis grid points per atom
(allowed values: $i=10 \times 3^k \times 4^l+2=12, 32, 42, 92...$)
nspa=: number of segments per atom
(allowed values: $i=10 \times 3^k \times 4^l+2=12, 32, 42, 92...$)
disex=: distance threshold for A matrix elements (Ångstrom)
rsolv=: distance to outer solvent sphere for cavity construction (Ångstrom)
routf=: factor for outer cavity construction in the outlying charge correction
cavity closed: pave intersection seams with segments
cavity open: leave untidy seams between atoms
ampran=: amplitude of the cavity de-symmetrization
phsran=: phase of the cavity de-symmetrization

If the $cosmo keyword is given without further specifications the default parameter are used (recommended). For the generation of the cavity, Cosmo also requires the definition of atomic radii. User defined values can be provided in Ångstrom units in the data group $cosmo_atoms, e.g. for a water molecule:

$cosmo_atoms
# radii in Angstrom units
o  1                                                   \
   radius=  1.7200
h  2-3                                                 \
   radius=  1.3000

If this section is missing in the control file, the default values defined in the radii.cosmo file (located in $TURBODIR/parameter) are used. A user defined value supersedes this defaults.

$cosmo and $cosmo_atoms can be set interactively with the Cosmo input program cosmoprep after the usual generation of the TURBOMOLE input.

The Cosmo energies and total charges are listed in the result section. E.g.:

  SCREENING CHARGE:
    cosmo      :  -0.003925
    correction :   0.003644
    total      :  -0.000282
  ENERGIES [a.u.]:
    Total energy            =      -76.0296831863
    Total energy + OC corr. =      -76.0297567835
    Dielectric energy       =       -0.0118029468
    Diel. energy + OC corr. =       -0.0118765440
    The following value is included for downward compatibility
    Total energy corrected  =      -76.0297199849

The dielectric energy of the system is already included in the total energy. OC corr denotes the outlying charge correction. The last energy entry gives the total outlying charge corrected energy in the old definition used in TURBOMOLE 5.7 and older versions.

The Cosmo result file, which contains the segment information, energies, and settings, can be set using: $cosmo_out file=

Alternatively set the name in the $cosmo section.

Isodensity Cavity:

This option can be used in HF/DFT single point calculations only. The $cosmo_isodens section defines the settings for the density based cavity setup (see also chapter 21.2). If the $cosmo_isodens keyword is given without suboptions, a scaled iosodensity cavity with default settings will be created. Possible options are:

$cosmo_isodens: activates the density based cavity setup. The default values of nspa and nsph are changed to 162 and 92, respectively. This values are superseded by the user defined nspa value of the $cosmo section. By default the scaled density method is used. The atom type dependent density values are read from the radii.cosmo file (located in $TURBODIR/parameter).
dx=: spacing of the marching tetrahedron grid in Å (default: 0.3Å).
all_dens=: use one isodensity value for all atom types (value in a.u.)

The outlying charge correction will be performed with a radii based outer cavity. Therfore, and for the smoothing of the density changes in the intersection areas of the scaled density method, radii are needed as for the standard Cosmo cavity. Please note: The isodensity cavity will be constructed only once at the beginning of the SCF calculation. The density constructed from the initial mos will be used (file mos or alpha/beta in case of unrestricted calculations). Because the mos of an initial guess do not provide a good density for the cavity construction, it is necessary to provide mos of a converged SCF calculation (e.g. a Cosmo calculation with standard cavity). We recommend the following three steps: perform a standard Cosmo calculation, add the isodensity options afterwards, and start the calculation a second time.

Radii based Isosurface Cavity:

The $cosmo_isorad section defines the radii defined isosurface cavity construction. Alternatively add the isorad option directly to the $cosmo keyword. The option uses the algorithm of the isodensity cavity construction but the objective function used depends on the cosmo radii instead of the electron density. The default values of nspa and nsph are changed to 162 and 92, respectively. This values are superseded by the user defined nspa value of the $cosmo section. The resulting surface exhibits smoother intersection seams and the segment areas are less diverse than the ones of the standard radii bases cavity construction.

$cosmo_isorad
dx=: spacing of the marching tetrahedron grid in Å (default: 0.3Å).

COSMO in MP2 Calculations:

The iterative Cosmo PTED scheme (see chapter 21.2) can be used with the mp2cosmo and cc2cosmo scripts. Options are explained in the help message (mp2cosmo -h or cc2cosmo -h). Both MP2 modules ricc2 and mpgrad can be utilized. The control file can be prepared by a normal Cosmo SCF input followed by a ricc2 or mpgrad input.

COSMO-MP2 geometry optimizations are currently only possible for PTE-COSMO-MP2 with the ricc2 module. For this you need to set in addition the (normal) COSMO input and (vacuum-like) MP2 gradient inputs the PTE option in the data group $reaction_field:

$ricc2
  geoopt model=mp2
$reaction_field
  PTE

COSMO in Numerical Frequency Calculations:

NumForce can handle two types of Cosmo frequency calculations. The first uses the normal relaxed Cosmo energy and gradient. It can be performed with a standard dscf or ridft Cosmo input without further settings. This is the right method to calculate a Hessian for optimizations. The second type, which uses the approach described in chapter 21.2, is implemented for ridft only. The input is the same as in the first case but NumForce has to be called with the -cosmo option. If no solvent refractive index refind=real is given in the $cosmo section of the control file the program uses the default (1.3).

COSMO in vertical excitations and polarizabilities:

Cosmo is implemented in escf and will be switched on automatically by the Cosmo keywords of the underlying SCF calculation. The refractive index, used for the fast term screening of the vertical excitations, needs to be defined in the cosmo section of control file (refind=real).

COSMO in CC2 and ADC(2) calculations:

For the calculation of ground-state energy at COSMO-CC2, vertical excitation energy at COSMO-CC2 and COSMO-ADC(2) and excited-state analytic gradient at COSMO-ADC(2), the post-SCF reaction-field scheme will be switched on by the $reaction_field keyword as:

 $reaction_field
   post-SCF
   ccs-like

For details on the $reaction_field data group see Sec. 25.2.13

DCCOSMO-RS:

The DCOSMO-RS model (see chapter 21.2) has been implemented for restricted and unrestricted DFT and HF energy calculations and gradients (programs: dscf/ridft and grad/rdgrad). In addition to the Cosmo settings defined at the beginning of this section, the $dcosmo_rs keyword has to be set.

$dcosmo_rs file=: activates the DCOSMO-RS method. The file defined in this option contains the DCOSMO-RS $\sigma$-potential and related data (examples can be found in the defaut potentials in the $TURBODIR/parameter directory).

If the potential file cannot be found in the local directory of the calculation, it will be searched in the $TURBODIR/parameter directory. The following $\sigma$-potential files for pure solvents at $25\,{}^{\circ}\mathrm{C}$ are implemented in the current Turbomole distribution (see parameter subdirectory):

h2o_25.pot

ethanol_25.pot

methanol_25.pot

thf_25.pot

propanone_25.pot

chcl3_25.pot

ccl4_25.pot

acetonitrile_25.pot

nitromethane_25.pot

dimethylsulfoxide_25.pot

diethylether_25.pot

hexane_25.pot

cyclohexane_25.pot

benzene_25.pot

toluene_25.pot

aniline_25.pot

The DCOSMO-RS energies and total charges are listed in the Cosmo section of the output:

  SCREENING CHARGE:
    cosmo      :  -0.012321
    correction :   0.011808
    total      :  -0.000513
 (correction on the COSMO level)
  ENERGIES [a.u.]:
    Total energy            =      -76.4841708454
    Outlying charge corr. (COSMO)    =    -0.0006542315
    Outlying charge corr. (DCOSMO-RS)=    -0.0011042856
    Combinatorial contribution of the solute =    -0.0017627889
    (at inf. dil. in the mixture/pure solvent. Not included in the total energy above)

The outlying charge correction cannot be defined straight forward like in the normal Cosmo model. Therefore, the output shows two corrections that can be added to the Total energy. The first one is the correction on the Cosmo level (COSMO) and the second is the difference of the DCOSMO-RS dielectric energy calculated from the corrected and the uncorrected Cosmo charges, respectively (DCOSMO-RS). The charges are corrected on the Cosmo level only. The Total energy includes the $E_{diel,RS}$ defined in section 21.2.2. Additionally the combinatorial contribution at infinute dilution of the COSMO-RS model is given in the output. The use of this energy makes sense if the molecule under consideration is different than the used solvent or not component of the solvent mixture, respectively. To be consistent one should only compare energies containing the same contributions, i.e. same outlying charge correction and with or without combinatorial contribution. Please note: the COSMO-RS contribution of the DCOSMO-RS energy depends on the reference state and the COSMO-RS parameterization (used in the calculation of the chosen COSMO-RS potential). Therefore, the DCOSMO-RS energies should not be used in a comparision with the gas phase energy, i.e. the calculation of solvation energies.