From NWChem

COSMO Solvation Model

COSMO is the continuum solvation `COnductor-like Screening MOdel' of A. Klamt and G. Schüürmann to describe dielectric screening effects in solvents^[1].

The NWChem COSMO module implements algorithm for calculation of the energy for the following methods:

1. Restricted Hartree-Fock (RHF),
2. Restricted open-shell Hartree-Fock (ROHF),
3. Restricted Kohn-Sham DFT (DFT),
4. Unrestricted Kohn-Sham DFT (ODFT),

by determining the solvent reaction field self-consistently with the solute charge distribution from the respective methods. Note that COSMO for unrestricted Hartree-Fock (UHF) method can also be performed by invoking the DFT module with appropriate keywords.

Correlation energy of solvent molecules may also be evaluated at

1. MP2,
2. CCSD,
3. CCSD+T(CCSD),
4. CCSD(T),

levels of theory. It is cautioned, however, that these correlated COSMO calculations determine the solvent reaction field using the HF charge distribution of the solute rather than the charge distribution of the correlation theory and are not entirely self consistent in that respect. In other words, these calculations assume that the correlation effect and solvation effect are largely additive, and the combination effect thereof is neglected. COSMO for MCSCF has not been implemented yet.

In the current implementation the code calculates the gas-phase energy of the system followed by the solution-phase energy, and returns the electrostatic contribution to the solvation free energy. At the present gradients are calculated by finite difference of the energy. Known problems include that the code does not work with spherical basis functions. The code does not calculate the non-electrostatic contributions to the free energy, except for the cavitation/dispersion contribution to the solvation free energy, which is computed and printed. It should be noted that one must in general take into account the standard state correction besides the electrostatic and cavitation/dispersion contribution to the solvation free energy, when a comparison to experimental data is made.

Invoking the COSMO solvation model is done by specifying the input COSMO input block with the input options as:

cosmo
  [off]
  [dielec  <real dielec default 78.4>]
  [radius  <real atom1>
           <real atom2>
      . . .
           <real atomN>]
  [rsolv   <real rsolv default 0.5>]
  [iscren  <integer iscren default 0>]
  [ifscrn  <integer ifscrn default 1>]
  [minbem  <integer minbem default 2>]
  [maxbem  <integer maxbem default 3>]
  [ificos  <integer ificos default 0>]
  [lineq   <integer lineq default 1>]
  [do_gasphase <logical do_gasphase default True>]
end

followed by the task directive specifying the wavefunction and type of calculation, e.g., "task scf energy", "task mp2 energy", "task dft optimize", etc.

"off' can be used to turn off COSMO in a compound (multiple task) run. By default, once the COSMO solvation model has been defined it will be used in subsequent calculations. Add the keyword "off" if COSMO is not needed in subsequent calculations.

"Dielec" is the value of the dielectric constant of the medium, with a default value of 78.4 (the dielectric constant for water).

"Radius" is an array that specifies the radius of the spheres associated with each atom and that make up the molecule-shaped cavity. Default values are Van der Waals radii. Values are in units of angstroms. The codes uses the following Van der Waals radii by default:

Default radii provided by Andreas Klamt (Cosmologic)

vdw radii: 1.17 (+/- 0.02) * Bondi radius^[2]

optimal vdw radii for H, C, N, O, F, S, Cl, Br, I^[3]

for heavy elements: 1.17*1.9

     data (vander(i),i=1,102)
    1 / 1.300,1.638,1.404,1.053,2.0475,2.00,
    2   1.830,1.720,1.720,1.8018,1.755,1.638,
    3   1.404,2.457,2.106,2.160,2.05,2.223,
    4   2.223,2.223,2.223,2.223,2.223,2.223,
    5   2.223,2.223,2.223,2.223,2.223,2.223,
    6   2.223,2.223,2.223,2.223,2.160,2.223,
    7   2.223,2.223,2.223,2.223,2.223,2.223,
    8   2.223,2.223,2.223,2.223,2.223,2.223,
    9   2.223,2.223,2.223,2.223,2.320,2.223,
    1   2.223,2.223,2.223,2.223,2.223,2.223,
    2   2.223,2.223,2.223,2.223,2.223,2.223,
    3   2.223,2.223,2.223,2.223,2.223,2.223,
    4   2.223,2.223,2.223,2.223,2.223,2.223,
    5   2.223,2.223,2.223,2.223,2.223,2.223,
    6   2.223,2.223,2.223,2.223,2.223,2.223,
    7   2.223,2.223,2.223,2.223,2.223,2.223,
    7   2.223,2.223,2.223,2.223,2.223,2.223/

For examples see Stefanovich et al.^[4] and Barone et al.^[5]

"Rsolv" is a parameter used to define the solvent accessible surface. See the original reference of Klamt and Schuurmann for a description. The default value is 0.5 (in angstroms).

"Iscren" is a flag to define the dielectric charge scaling option. "iscren 1" implies the original scaling from Klamt and Schüürmann, mainly "(ε − 1) / (ε + 1 / 2)", where ε is the dielectric constant. "iscren 0" implies the modified scaling suggested by Stefanovich and Truong^[6], mainly "(ε − 1) / ε". Default is to use the modified scaling. For high dielectric the difference between the scaling is not significant.

"ifscrn" controls the correction of the COSMO charges. The charges on the cavity surface have to sum to the charge of the solvated system as a consequence of Gauss's law. In practice this condition is only satisfied approximately, and corrections have been suggested^[7] (see Eq.23). However, introducing these corrections leads to inconsistencies in the energy expression that cause inaccuracies in the gradient evaluation. For this reason the default "ifscrn 1" does not apply any corrections to the COSMO charges. "ifscrn 2" does correct the charges by adding fractional charges so that Gauss's law is exactly satisfied. In practice "ifscrn 1" works better in geometry optimizations as correct gradients are used, but the total energy is affected. On the other hand "ifscrn 2" gives better total energies, but it can cause problems in geometry optimizations as the gradients are not strictly correct.

The next three parameters define the tesselation of the unit sphere. The approach follows the original proposal by Klamt and Schüürmann. A very fine tesselation is generated from "maxbem" refining passes starting from either an octahedron or an icosahedron. The boundary elements created with the fine tesselation are condensed down to a coarser tesselation based on "minbem". The induced point charges from the polarization of the medium are assigned to the centers of the coarser tesselation. Default values are "minbem 2" and "maxbem 3". The flag +ificos+ serves to select the original tesselation, "ificos 0" for an octahedron (default) and "ificos 1" for an icoshedron. Starting from an icosahedron yields a somewhat finer tesselation that converges somewhat faster. Solvation energies are not really sensitive to this choice for sufficiently fine tesselations.

The "lineq" parameter serves to select the numerical algorithm to solve the linear equations yielding the effective charges that represent the polarization of the medium. "lineq 0" selects an iterative method (default), "lineq 1" selects a dense matrix linear equation solver. For large molecules where the number of effective charges is large, the codes selects the iterative method.

"do_gasphase" is a flag to control whether the calculation of the solvation energy is preceded by a gas phase calculation. The default is to always perform a gas phase calculation first and then calculate the solvation starting from the converged gas phase electron density. However, in geometry optimizations this approach can double the cost. In such a case setting "do_gasphase false" suppresses the gas phase calculations and only the solvated system calculations are performed. This option needs to be used with care as in some cases starting the COSMO solvation from an unconverged electron density can generate unphysical charges that lock the calculation into strange electron distributions.

The following example is for a water molecule in `water', using the HF/6-31G** level of theory:

start
echo
 title "h2o"
geometry
 o                  .0000000000         .0000000000        -.0486020332
 h                  .7545655371         .0000000000         .5243010666
 h                 -.7545655371         .0000000000         .5243010666
end
basis segment cartesian
 o library 6-31g**
 h library 6-31g**
end
cosmo
 dielec 78.0
 radius 1.40
        1.16
        1.16
 rsolv  0.50
 lineq  0
end
task scf energy

Instead of listing COSMO radii parameters in the input, the former can now be loaded using an external file through the following directive (placed outside the cosmo block)

set cosmo:map cosmo.par

The format for such file (named as cosmo.par in the above case) consists of the atom name (as found in geometry block) followed by the radii. The file HAS TO BE PLACED IN THE PERMANENT DIRECTORY. In the case of the water example shown above it can take the following form

O 1.40
H 1.16

The input file in this case is

start
echo
 title "h2o"
geometry
 o                  .0000000000         .0000000000        -.0486020332
 h                  .7545655371         .0000000000         .5243010666
 h                 -.7545655371         .0000000000         .5243010666
end
basis segment cartesian
 o library 6-31g**
 h library 6-31g**
end
cosmo
 dielec 78.0
 rsolv  0.50
 lineq  0
end

set cosmo:map cosmo.par

task scf energy

References

1. ↑ Klamt, A; Schuurmann, G (1993). "COSMO: A new approach to dielectric screening in solvents with explicit expressions for the screening energy and its gradient". Journal of the Chemical Society, Perkin Transactions 2: 799-805. doi:10.1039/P29930000799. 
2. ↑ A. Bondi (1964). "van der Waals volums and radii". Journal of Physical Chemistry 68: 441-451. doi:10.1021/j100785a001. 
3. ↑ A. Klamt, V. Jonas (1998). "Refinement and parametrization of COSMO-RS". Journal of physical chemistry A 102: 5074-5085. doi:10.1021/jp980017s. 
4. ↑ E. V. Stefanovich, T. N. Truong (1995). "Optimized atomic radii for quantum dielectric continuum solvation models". Chemical Physics Letters 244: 65-74. doi:10.1016/0009-2614(95)00898-E. 
5. ↑ V. Barone, M. Cossi (1997). "A new definition of cavities for the computation of solvation free energies by the polarizable continuum model". Journal of Chemical Physics 107: 3210-3221. doi:10.1063/1.474671. 
6. ↑ E. V. Stefanovich, T. N. Truong (1995). "Optimized atomic radii for quantum dielectric continuum solvation models". Chemical Physics Letters 244: 65-74. doi:10.1016/0009-2614(95)00898-E. 
7. ↑ MA Aguilar, FJ Olivares del Valle (1993). "Nonequilibrium solvation: An ab initio quantummechanical method in the continuum cavity model approximation". Journal of Chemical Physics 98: 7375-7384. doi:10.1063/1.464728.