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.
SOLVATION
{Surf Esurf}
{Div {Ndiv=3} {Min=0.5} {OFAC=0.8}}
{NOASS}
{Solv {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 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
Details
attributes of the solvent. Eps specifies the dielectric constant (the default
relates to water).
Rad specifies the radius of the (rigid sphere) solvent molecules, in angstrom.
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:
Enon-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
...
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.
