Solvent effects: COSMO

You can study chemistry in solution, as contrasted to the gas phase, with the implementation in ADF [66] of the Conductor like Screening Model (COSMO) of solvation [67-69]. The energy derivatives can also be calculated, so geometry optimization, harmonic frequencies, et cetera are available within this model.

The COSMO model is a dielectric model in which the solute molecule is embedded in a molecule-shaped cavity surrounded by a dielectric medium with given dielectric constant e. Energy-related terms are computed for a conductor first, then scaled by the function

f(ε) = (ε-1)/(ε+x)   (2.1.1)

The empirical scaling factor x is specified in the input data block for the solvation key. The block key SOLVATION turns the solvation calculation on. In most cases default values are available for the involved parameters.

The COSMO routine has never been tested with non-gas phase fragments. It is designed to correct the "total" energy (wrt fragments gas phase) for solvation. The COSMO energy of the fragments is never taken into account, not even for the atoms, and no information is passed in regarding the COSMO energy when the fragments are defined.

If a calculations was done on a fragment, then the wavefunction obtained would be optimal for the fragment in solution, but not optimal for gas phase. The energies with the gas phase Hamiltonian would be higher, and the apparent solvation contribution to bonding would also be higher. The net point is that normally the COSMO procedure reports the energy of Esolv(AB), but to get the solvation energy, you need to subtract the E(AB, solv) from E(AB,gas) because the wavefunction changes (unless you are doing it post-SCF.) For this you need to have the same reference fragments in each case A(g) and B(g).

SOLVATION
 {Surf Esurf {NOKEEP}}
 {Div {Ndiv=3} {Min=0.5} {OFAC=0.8}}
 {NOASS}
 {Solv {Name=solvent}{Eps=78.4}{Del=1.4}{Rad=1.4}{Emp=0.0}{Cav0=1.321}{Cav1=0.0067639}}
 {RADII
 name1=value1
 name2=value2
  ...
 subend }
 {Charged {Method=meth} {Conv=1e-6} {Omega=1.0} (Iter=300} {Corr} }
 {C-Mat How {SCF} tol=1e-10 }
 {DISC {SC=0.01} {LEG=4} {TOL=0.1} }
 {SCF {When} {How}
 {LPRT }
End

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. Follows a description of the subkeys

Surf

Esurf must be Wsurf, Asurf, Esurf or Klamt. Four different cavity types are available.
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 (the default) 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.
These 3 surfaces are constructed with the GEPOL93 algorithm [70]
The SES surface is the default in ADF.
The fourth surface option is Klamt as described in [67]. It excludes the cusp regions also.

The optional parameter NOKEEP controls surface creation during calculation of frequencies by numerical differentiation. By default, the surface is constructed only once at the central geometry and is used for the rest of the calculation. If the NOKEEP is specified then ADF will construct a new surface at each displaced geometry. The NOKEEP option was the default in ADF2005 and earlier versions but it was found to cause problems. Since ADF2006 one needs to specify SURF NOKEEP to get the same behavior.
The actual construction of the surface involves a few technical parameters controlled with the subkey

DIV

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 Ndiv=3
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 Min=0.5
Ofac is a maximum allowed overlap of new created spheres, in the construction procedure. Default Ofac=0.8

NOASS

By default all new spheres that are created in the surface-construction are assigned to atoms, for the purpose of gradient computations (geometry optimization). Specifying the noass subkey turns this off. It has no argument.

Solv

Solvent details. Eps specifies the dielectric constant (the default relates to water). In ADF 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.

Emp addresses the empirical scaling factor x in the formula above.
Other options specify a linear parameterization of non-electrostatic terms as a function of surface area:

E_non-elst = f(ε) ∗ (CAV0 + CAV1∗area)   (2.1.2)

In order to construct the surface you have to specify the atomic ('Van der Waals') radii. There are three ways of doing this. In the first method you append 'R=value' to the atomic coordinates record, in the ATOMS key block. This would look like, for instance

C 1 2 3 CC CCO CCOH f=C.dz R=2.0

It assigns a radius of 2.0 to the Carbon atom.
In the second method you apply the same format, but specify a symbol (identifier) rather than a value

C 1 2 3 CC CCO CCOH f=C.dz R=C-sp3

The identifiers must be defined in the (optional) RADII subkey block in the Solvation data block (see next).
In the third method, you don't modify the Atoms block at all. In this case, the RADII subkey must be used and the 'identifiers' in it must be exactly the atom type names in the Atoms block.

Radii

This subkey is block type. Its data block (if the subkey is used) must terminate with a record subend. In the Radii data block you give a list of identifiers and values

SOLVATION
  ...
 Radii
 name1=value1
 name2=value2
  ...
 Subend
  ...
End

The values are the radii of the atomic spheres, in the same units of length as used in the Atoms block (angstrom or bohr). The names specify to which atoms these values apply. As discussed for the Solv subkey this depends on the Atoms block. If in the specification of atomic coordinates you have used the 'R=' construct to assign radii, with identifiers rather than values for the R-value, these identifiers must be defined in the Radii sub block. If no 'R=' construct was applied in the Atoms block, you must use the atom type names as they occurred in the Atoms data block. Referring to the example given in the Solv subkey discussion, you might have

...
 Radii
  C-sp3=2.0
  ...
 Subend
...

A simple atom type reference might look like

...
 Radii
  C=2.0
  ...
 Subend
...

When no radius specified a default value is used. The default value for an atom is the corresponding Van der Waals radius from the MM3 method by Allinger divided by 1.2.

This concludes the discussion of the Radii subkey.

CHARGED

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 Aq=-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 specifies the equation-solving algorithm. Meth=INVER requests direct inversion. Meth=GAUS calls for the Gauss-Seidel iterative method. Meth=Jacobi activates another standard iterative procedure. The latter two methods require a positive-definite matrix (which may fail to be the case in an actual calculation) and can be used with a relaxation technique, controlled by the relaxation parameter OMEGA (1.0=no relaxation).
Meth=CONJ (default) 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 methods.

Some of the molecular electronic charge distribution may be located outside the cavity. This affects the assumptions underlying the COSMO equations. Specifying the CORR option to the CHARGED subkey constrains the computed solvent surface charges to add up to the negative of the molecular charge.

C-MATRIX

- How: For the potential f we need the Coulomb interaction between the charges q and the molecular electronic density (and nuclei). Three methods are available, specified by the first option to the C-Matrix subkey.
a) EXACT: compute the straightforward Coulomb potential due to the charge q in each point of the molecular numerical integration grid and integrate against the electronic charge density. This is, in principle, exact but may have inaccuracies when the numerical integration points are very close to the positions of a charge q. To remedy this, the point charges q can be 'smeared out' and represented by a disc, see the next subkey DISC.
b) FIT: same as EXACT, but the q-potentials are now integrated not against the exact electronic charge density, but against the (much cheaper-to-compute) fitted density. The same DISC considerations apply.
c) POT: evaluate the molecular potential at the position of the charge q and multiply against the q-strength. Since the molecular Coulomb potential is computed from the fit density, any difference in results between the FIT and the POT approach should be attributed to the DISC issue.
POT is the default, because it is faster, and is only inadequate if the fit density is very inaccurate, which would be a problem anyway.
- SCF: If you specify this option, the computation of the Coulomb interaction matrix (between electrons and surface charges) is carried out during the SCF procedure, but this turns out to hamper the SCF convergence behavior. Therefore: not recommended. IF you use it, the program will switch to one of the other 3 methods, as given by the 'How' option, as soon as the SCF convergence error drops below
TOL: (applies only to the SCF option, which is not recommended).

DISC

Applies only when the C-matrix method is EXACT or FIT. Note, however, that the default for the C-matrix method is POT, in which case the DISC subkey has no meaning. The DISC key lets the program replace the point charges q by a solid uniformly charged spherical surface disc whenever the numerical integration accuracy requires so, i.e. for those charges that are close to numerical integration points.
Options:
SC defines a shrinking factor, by which the actual disc radius used is reduced from its 'normal' value: an inscribed disc in the triangular surface partitions that define the distribution of surface charges, see the subkey DIV.
LEG gives the polynomial expansion order of the disc potentials. The Legendre expansion converges rapidly and the default should be adequate.

TOL is a tolerance parameter to control the accuracy of the disc potential evaluations.

SCF

In COSMO calculations you can include the surface charges in the Fock operator self-consistently, i.e. by recomputing the charges q at every SCF cycle and include them in the equations, or in a perturbational manner, i.e. post-SCF. This is controlled with the first option. The When option must be either VAR or PERT, for variational and perturbational, respectively. Default is VAR.
The second (HOW) option applies only to the WHEN=VAR case and may affect the speed of SCF convergence. The COSMO calculation implies a considerable increase in CPU time! Values for HOW:
- ALL: This includes it in all SCF cycles (except for the first SCF cycle, which is gas-phase)
- LAST: This lets the program first converge the SCF completely without any solvent effects. Thereafter, the COSMO is turned on, hopefully converging in fewer cycles now, to compensate for the 'double' SCF effort.
- TOL=0.1 (or another value) is an in-between approach: converge the gas-phase SCF until the SCF error is below TOL, then turn on COSMO.

LPRT

This is a debug switch and triggers a lot more output related to the cavity construction etc.