This is an old revision of the document!


The COSMO model

The Conductor-like Screening Model (COSMO) (A. Klamt and G. Schüürmann, J. Chem. Soc. Perkin Trans. II 799-805 (1993)) is currently available for HF (RHF,UHF) and DFT (RKS,UKS) energy calculations and the corresponding gradients.

The COSMO model is invoked by the COSMO card:

COSMO[,option$_1$=value$_1$, option$_2$=value$_2, \ldots$]

where option can be

  • NPPA size of the underlying basis grid. The value must satisfy: $value = 10 \times 3^k \times 4^l +2$ (default = 1082; type integer).
  • NSPA number of segments for non hydrogen atoms. The value must satisfy: $values = 10 \times 3^k \times 4^l +2$ (default = 92; type integer).
  • CAVITY the intersection seams of the molecular surface are closed (1) or open (0) (default = 1; type integer).
  • EPSILON dielectric permittivity (default = -1.d0, which means $\epsilon=\infty$; type real)
  • DISEX distance criteria for the $\bf A$-matrix setup. Short range interactions (segment centre distances $<$ DISEX $\times$ mean atomic diameter) are calculated using the underlying basis grid. Long range interactions are calculated via the segment centres (default = 10.d0; type float).
  • ROUTF factor used for outer cavity construction. The radii of the outer cavity are defined as: $r^{out}_i=r_i+ROUTF \times RSOLV$ (default = 0.85d0; type float)
  • PHSRAN phase offset of coordinate randomization (default = 0.d0; type float)
  • AMPRAN amplitude factor of coordinate randomization (default = 1.0d-5; type float)
  • RSOLV additional radius for cavity construction (default = -1d0, the optimized H radius is used; type float).
  • MAXNPS maximal number of surface segments (default = -1, will be estimated; type integer).

It is recommended to change the default values for problematic cases only.

By default the program uses optimized radii if existent and 1.17$\times$vdW radius else. The optimized radii [Å] are: H=1.30, C=2.00, N=1.83, O=1.72, F=1.72, S=2.16, Cl=2.05, Br=2.16, I=2.32. Own proposals can be given directly subsequent to the cosmo card:
RAD,symbol, radius

where the radius has to be given in Å.

Example:
cosmo rad,O,1.72 rad,H,1.3 Output file:
The COSMO output file will be written after every converged SCF calculation. The segment charges and potentials are corrected by the outlying charge correction. For the total charges and energies corrected and uncorrected values are given. The normal output file contains uncorrected values only. It is recommended to use the corrected values from the output file.

Optimizations:
It is recommended to use optimizer that operates with gradients exclusively. Line search techniques that use energies tends to fail, because of the energy discontinuities which may occur due to a reorganization of the segments after a geometry step. For the same reasons numerical gradients are not recommended.

COSMO is a continuum solvation model, in which the solvent is represented as a dielectric continuum of permittivity $\epsilon$. The solute molecule is placed in a cavity inside the continuum. The response of the continuum due to the charge distribution of the solute is described by the generation of a screening charge distribution on the cavity surface. This charge distribution can be calculated by solving the boundary equation of vanishing electrostatic potential on the surface of a conductor. After a discretization of the cavity surface into sufficiently small segments, the vector of the screening charges on the surface segments is

$${\bf q^{*}} = -{\bf A^{-1}} \Phi$$

where $\Phi$ is the vector of the potential due to the solute charge distribution on the segments, and $\bf A$ is the interaction matrix of the screening charges on the segments. This solution is exact for an electric conductor. For finite dielectrics the true dielectric screening charges can be approximated very well by scaling the charge density of a conductor with $f(\epsilon)$.

$${\bf q} = f(\epsilon) {\bf q^{*}} ; \quad f(\epsilon) = (\epsilon - 1)/(\epsilon + 0.5)$$

In every SCF step the screening charges $\bf q$ have to be generated from the potential $\Phi$, and then added to the Hamiltonian as external point charges. The total energy of the system is $$E_{tot} = E_{0} + E_{diel}; \quad E_{diel} = \frac{1}{2} \Phi {\bf q}$$ where $E_{0}$ is the bare self-energy of the system and $E_{diel}$ the dielectric energy.
Cavity construction:
First a surface of mutually excluding spheres of radius $R_i+rsolv$ is constructed, where the $R_i$ are the radii of the atoms, defined as element specific radii and $rsolv$ is some radius representing a typical maximum curvature of a solvent molecular surface. $rsolv$ should not be misinterpreted as a mean solvent radius, nor modified for different solvents. Every atomic sphere is represented by an underlying basis grid of nppa points per full atom. Basis grid points which intersect a sphere of a different atom are neglected. In a second step the remainder of the basis grid points are projected to the surface defined by the radii $R_i$. As a third step of the cavity construction the remaining basis grid points are gathered to segments, which are the areas of constant screening charges in the numerical solution. Finally, the intersection seams between the atoms are filled with additional segments.
Now the A-matrix can be set up. The matrix elements will be calculated from the basis grid points of the segments for close and medium segment distances (governed by the disex value), or using the segment centres for large segment distances.
Outlying charge correction:
The non vanishing electron density outside the cavity causes an error that can be corrected by the outlying charge correction. This correction uses the potential on the so called outer surface (defined by the radii $R_i+\mbox{rsolv} \times \mbox{routf}$) to estimate a correction term for the screening charges and the energies (A. Klamt and V. Jonas, J. Chem. Phys. 105, 9972-9981(1996)). The correction will be performed once at the end of a converged SCF calculation. All corrected values can be found in the COSMO output file.

CUBE,filename,iflag,$n_1$,$n_2$,$n_3$

calls a module which dumps the values of various properties on a spatial parallelopipedal grid to an external file. The purpose is to allow plotting of orbitals, densities and other quantities by external programs. The format of the file is intended to be the same as that produced by other programs.

  • filename is the unix path name of the file to be written, and its specification is mandatory.
  • iflag If iflag is negative (default), a formatted file will be written, otherwise unformatted fortran i/o will be used.
  • $n_1$,$n_2$,$n_3$ specify the number of grid points in each of three dimensions. If not specified, sensible defaults are chosen.

By default, the last density computed is evaluated on the grid, and written to filename. This behaviour can be modified by one or more of the following subcommands.

STEP,[stepx],[stepy], [stepz]
stepx,stepy, stepz specify the point spacing in each of three axis directions. By default, the value of stepx,stepy, stepz is determined by the number of grid points, the Bragg radii of the atoms, and some related parameters.

DENSITY,[density-source]
GRADIENT,[density-source]
LAPLACIAN,[density-source]

Compute the density and, optionally, its gradient and laplacian. $<$density-source$>$ may be a record number containing the required density, and may contain further qualification, such as set or state number, in the usual way. By default, the last computed density is taken.

Example:

{multi;wf,spin=0,symmetry=1;state,2;dm;orbital,2140.2}
{cube,multi_1;density,2140.2,state=1.1}   ! cube for state 1
{cube,multi_2;density,2140.2,state=2.1}   ! cube for state 2
{mrci;wf,spin=0,symmetry=1;state,2;dm,3000.2}   !save densities in record 3000.2
{cube,mrci_1;density,3000.2,state=1.1}    ! cube for state 1
{cube,mrci_2;density,3000.2,state=2.1}    ! cube for state 2

ORBITAL,[orbital-list],[RECORD=orbital-source]

$<$orbital-list$>$ is a list of one or more orbital numbers of the form number.symmetry or keywords chosen from HOMO, LUMO, OCC (all occupied orbitals), ALL. If nothing is specified, the default is HOMO. $<$orbital-source$>$ may be a record number containing the required density, and may contain further qualification, such as set number, in the usual way. By default, the last computed orbitals are taken.

Note that the CUBE file format precludes simultaneous orbital and density dumps, but that this may be achieved in the GOPENMOL format (see GOPENMOL --- calculate grids for visualization in gOpenMol).

AXIS,x,y,z

x,y,z specify the unnormalised direction cosines of one of the three axes defining the grid. Up to three AXIS commands can be given, but none is required. Axes need not be orthogonal. By default, the first axis is the cartesian $x$, the second is orthogonal to the first and to the cartesian $z$, and the third is orthogonal to the first two.

Based on the direction of the coordinate axes, a parallelopiped (in the usual case of orthogonal axes, a cuboid) is constructed to contain the molecule completely. The atoms are assumed to be spherical, with an extent proportional to their Bragg radii, and the constant of proportionality can be changed from the default value using

BRAGG,scale

After the parallelopiped has been constructed, the grid is laid out with equal spacing to cover it using the number of points specified on the CUBE command.

ORIGIN,x,y,z

x,y,z specify the centroid of the grid. It is usually not necessary to use this option, since the default should suffice for most purposes.

TITLE,title

Set a user defined title in the cube file.

DESCRIPTION,description

Set a user defined description in the cube file.

The formatted cube file contains the following records

  • (A) job title.
  • (A) brief description of the file contents.
  • (I5,3F12.6) number of atoms, coordinates of grid origin (bohr).
  • (I5,3F12.6) number of grid points $n_1$, step vector for first grid dimension.
  • (I5,3F12.6) number of grid points $n_2$, step vector for second grid dimension.
  • (I5,3F12.6) number of grid points $n_3$, step vector for third grid dimension.
  • (I5,4F12.6) atomic number, charge and coordinates; one such record for each atom.
  • (6E13.5) $n_1 \times n_2$ records of length $n_3$ containing the values of the density or orbital at each grid point. With default choice of the axes, the inner loop will be over $n_3$ points along the cartesian $z$ axis, and the outer loop over $n_1$ points along the cartesian $x$ axis.

In the case of a number of orbitals $m$, the record length is $m \times n_3$, with the data for a single grid point grouped together. In the case of the density gradient, there is first a record of length $n_3$ containing the density, then one of length $3n_3$ containing the gradient, with the three cartesian components contiguous. For the laplacian, there is a further record of length $n_3$.

GOPENMOL,filename,iflag,$n_1$,$n_2$,$n_3$

The syntax and sub-options are exactly the same as for CUBE, except that the files produced are in a format that can be used directly in the gOpenMol visualization program. The following should be noted.

  • Only the base name (up to the last ’.’) in filename is used, and is appended by different suffices to create several different files:
    • .crd A CHARMm CRD-format file containing the coordinates is always produced, and may be used in the invocation of gOpenMol: rungOpenMol -ifilename.crd
    • _density.plt If DENSITY is given, then the file filename_density.plt is produced and contains the density grid in gOpenMol internal format.
    • _orbital_number.symmetry.plt If ORBITAL is given, then for each orbital number.symmetry specified, the file filename_orbital_number.symmetry.plt is produced and contains the orbital grid in gOpenMol internal format.
  • The default is not to produce any orbitals or densities, and so only the atomic coordinates are dumped.
  • The default is to use unformatted binary files, and this should not normally be changed.
  • The ORIGIN and AXIS commands should not be used.
  • If INTERACT is given in the input, when all the grids have been calculated, an attempt is made to start gOpenMol by executing the Unix command rungOpenMol. If rungOpenMol is not in $PATH, then nothing happens. Otherwise, gOpenMol should start and display the molecule. Any .plt files produced can be added to the display by following the Plot;Contour menu item. The name of the Unix command may be changed from the default rungOpenMol by specifying it as the first argument to the INTERACT directive. By default, gOpenMol is not started, and this is equivalent to giving the command BATCH.