COSMO: Conductor like Screening Model and the Solvation-key

You can study chemistry in solution, as contrasted to the gas phase, with the implementation in BAND of the Conductor like Screening Model (COSMO) of solvation.[36]

In the COSMO model all solvents are roughly the same, and approximated by an enveloping metal sheet. One explicit dependency on the solvent is that the solvation energy is scaled by

\[f(\epsilon) = \frac{\epsilon-1}{\epsilon+\chi}\]

and this depends on the dielectric constant of the solvent, and an empirical factor \(\chi\). The other is that the shape of the surface is influenced by the Rad parameter, see below.

SOLVATION (block type)

   {Surf [Delley|Wsurf|Asurf|Esurf|Klamt]}
   {Solv {Name=solvent} {Eps=78.4} {Rad=1.4} {Del=1.4} {Emp=0.0}}
   {Div {Ndiv=3} {NFdiv=1} {Min=0.5} {OFAC=0.8} {leb1=23} {leb2=29} {rleb=1.5}}
   {Charge {Method=meth} {Conv=1e-10} (Iter=1000} {NoCorr}}
   {Cvec [Exact|FitPot]}
   {SCF [Var|Pert]}

Presence of the SOLVATION key block triggers the solvent calculation and does not require additional data. With subkeys you can customize various aspects of the model, for instance to specify the type of solute. None of the subkeys is obligatory

See also the subkey StaticCosmoSurface in the GEOOPT block and the subkey StaticCosmoSurface in the FREQUENCIES block.

Below follows a description of the subkeys:

Surf [Delley|Wsurf|Asurf|Esurf|Klamt]

(Default: Delley) Five different cavity types are available. The Wsurf, Asurf, and Esurf surfaces are constructed with the GEPOL93 algorithm.[38]

Wsurf triggers the Van der Waals surface (VdW), which consists of the union of all atomic spheres.
Asurf gives the Solvent-Accessible-Surface (SAS). This is similar to VdW but consists of the path traced by the center of a spherical solvent molecule rolling about the VdW surface or, equivalently, a VdW surface created by atomic spheres to which the solvent radius has been added. These two surface types contain cusps at the intersection of spheres.
Esurf gives the Solvent-Excluding-Surface (SES), which consists of the path traced by the surface of a spherical solvent molecule rolling about the VdW surface. Primarily, this consists of the VdW surface but in the regions where the spheres would intersect, the concave part of the solvent sphere replaces the cusp. This SES surface is the default.
The fourth surface option is Klamt as described in [36]. It excludes the cusp regions also.
The fifth surface is the so called Delley surface.[39]

Solvent details.

Del is the value of Klamt’s delta_sol parameter, only relevant in case of Klamt surface.
Emp addresses the empirical scaling factor x for the energy scaling.
Name, Eps, Rad

Eps specifies the dielectric constant (the default relates to water). An infinite value for Eps is chosen if Eps is specified to be lower than 1.0. Rad specifies the radius of the (rigid sphere) solvent molecules, in Angstrom. Instead of specifying Eps and Rad one can specify a solvent name or formula after ‘name=’. The following table lists names and formulas that are recognized with the corresponding values for Eps and Rad. The names and formulas are case-insensitive.

Name Formula Eps Rad
AceticAcid CH3COOH 6.19 2.83
Acetone CH3COCH3 20.7 3.08
Acetonitrile CH3CN 37.5 2.76
Ammonia NH3 16.9 2.24
Aniline C6H5NH2 6.8 3.31
Benzene C6H6 2.3 3.28
BenzylAlcohol C6H5CH2OH 13.1 3.45
Bromoform CHBr3 4.3 3.26
Butanol C4H9OH 17.5 3.31
isoButanol (CH3)2CHCH2OH 17.9 3.33
tertButanol (CH3)3COH 12.4 3.35
CarbonDisulfide CS2 2.6 2.88
CarbonTetrachloride CCl4 2.2 3.37
Chloroform CHCl3 4.8 3.17
Cyclohexane C6H12 2 3.5
Cyclohexanone C6H10O 15 3.46
Dichlorobenzene C6H4Cl2 9.8 3.54
DiethylEther (CH3CH2)2O 4.34 3.46
Dioxane C4H8O2 2.2 3.24
DMFA (CH3)2NCHO 37 3.13
DMSO (CH3)2SO 46.7 3.04
Ethanol CH3CH2OH 24.55 2.85
EthylAcetate CH3COOCH2CH3 6.02 3.39
Dichloroethane ClCH2CH2Cl 10.66 3.15
EthyleneGlycol HOCH2CH2OH 37.7 2.81
Formamide HCONH2 109.5 2.51
FormicAcid HCOOH 58.5 2.47
Glycerol C3H8O3 42.5 3.07
HexamethylPhosphoramide C6H18N3OP 43.3 4.1
Hexane C6H14 1.88 3.74
Hydrazine N2H4 51.7 2.33
Methanol CH3OH 32.6 2.53
MethylEthylKetone CH3CH2COCH3 18.5 3.3
Dichloromethane CH2Cl2 8.9 2.94
Methylformamide HCONHCH3 182.4 2.86
Methypyrrolidinone C5H9NO 33 3.36
Nitrobenzene C6H5NO2 34.8 3.44
Nitrogen N2 1.45 2.36
Nitromethane CH3NO2 35.87 2.77
PhosphorylChloride POCl3 13.9 3.33
IsoPropanol (CH3)2CHOH 19.9 3.12
Pyridine C5H5N 12.4 3.18
Sulfolane C4H8SO2 43.3 3.35
Tetrahydrofuran C4H8O 7.58 3.18
Toluene C6H5CH3 2.38 3.48
Triethylamine (CH3CH2)3N 2.44 3.81
TrifluoroaceticAcid CF3COOH 42.1 3.12
Water H2O 78.39 1.93
Ndiv, NFdiv
Ndiv controls how fine the spheres, that in fact describe the surface, are partitioned in small surface triangles, each containing one point charge to represent the polarization of the cavity surface. Default division level for triangles Ndiv=3. Default final division level for triangles NFdiv=1 (NFdiv \(\leq\) Ndiv).
(Default: 0.5) Min specifies the size, in Angstrom, of the smallest sphere that may be constructed by the SES surface. For VdW and SAS surfaces it has no meaning.
(Default: 0.8) Ofac is a maximum allowed overlap of new created spheres, in the construction procedure.
leb1, leb2, rleb
For the Delley type of construction one needs to set the variables leb1 (Default: 23), leb2 (Default: 29), and rleb (Default: 1.5) to set the number of surface points. If the cavity radius of an atom is lower than rleb use leb1, otherwise use leb2. These values can be changed: using a higher value for leb1 and leb2 gives more surface points (maximal value leb1, leb2 is 29). A value of 23 means 194 surface points in case of a single atom, and 29 means 302 surface points in case of a single atom Typically one could use leb1 for the surface point of H, and leb2 for the surface points of other elements.

In the Radii data block you give a list of atom types and values:

  . . .
     H 0.7
     C 1.3
     . . .
  . . .

The values are the radii of the atomic spheres, in the default unit of length (Angstrom or Bohr). If not specified the default values are those by Allinger.[40]

This addresses the determination of the (point) charges that model the cavity surface polarization. In COSMO calculations you compute the surface point charges q by solving the equation \(A \cdot q = -f\), where f is the molecular potential at the location of the surface charges q and A is the self-interaction matrix of the charges. The number of charges can be substantial and the matrix A hence very large. A direct method, i.e. inversion of A, may be very cumbersome or even impossible due to memory limitations, in which case you have to resort to an iterative method. Meth (Default: INVER) specifies the equation-solving algorithm. The option INVER requests direct inversion. The option CONJ uses the preconditioned biconjugate gradient method. This is guaranteed to converge and does not require huge amounts of memory. CONV and ITER are the convergence criterion and the maximum number of iterations for the iterative method. In periodic calculations the unit cell needs to be neutral. Therefore, by default the COSMO charges are forced to sum to zero. When NOCORR is present, this constraint is not applied.

(Default: Exact) There are two ways to calculate the molecule-(COSMO)surface interaction. EXACT: Consider the interactions between the electron density and the potential due to the COSMO points. FITPOT: Use the COSMO point charges in the potential of the molecule (obtained by density fitting).

The EXACT option is exact on paper, but in reality numerical integration is used to get the energy. The standard grid is formally not suitable to handle the cusps at the COSMO points. Still this option works well because the singularity is softened by the smallness of the electronic density at the COSMO points. The EXACT option is much faster when doing a geometry optimization.

(Default: VAR) The option VAR is to apply the external potential of the COSMO surface during the SCF. To apply COSMO only post SCF you can use the PERT option.

Additional keys for periodic systems

PeriodicSolvation (block-type)
For the simulation of periodic structures additional information might be necessary. Therefore, the PeriodicSolvation block can be used.
   {RemovePointsWithNegativeZ [True|False]}
   {Nstar integer}
(Default: false) This option, expecting a boolean value (true or false), handles whether a COSMO surface is constructed on both sides of a surface. If one is only interested in the solvation effect on the upper side of a surface then this option should be set to true.
(Default: 4) This option, expecting an integer number (>2), handles the accuracy for the construction of the COMSO surface. The larger the given number the more accurate the construction.

General remarks: The accuracy of the result and the calculation time is influenced by the screening radius SCREENING%RMADEL (see Screening). If the calculation does take too long, defining a smaller radius does help. But: too small radii, especially smaller than the lattice constants, will give unphysical results.