1% 2% $Id$ 3% 4% $Id$ 5\label{sec:cosmo} 6 7COSMO is the continuum solvation `COnductor-like Screening MOdel' 8of A. Klamt and G. Sch\"{u}\"{u}rmann to describe dielectric screening 9effects in solvents. 10 11\begin{enumerate} 12\item A. Klamt and G. Sch\"{u}\"{u}rmann, J.~Chem.~Soc.~Perkin Trans.~2, 1993 13799 (1993). 14\end{enumerate} 15 16The NWChem COSMO module implements algorithm for calculation of the 17energy for the following methods: 18\begin{enumerate} 19\item Restricted Hartree-Fock (RHF), 20\item Restricted open-shell Hartree-Fock (ROHF), 21\item Restricted Kohn-Sham DFT (DFT), 22\item Unrestricted Kohn-Sham DFT (ODFT), 23\end{enumerate} 24by determining the solvent reaction field self-consistently 25with the solute charge distribution from the respective methods. 26Note that COSMO for unrestricted Hartree-Fock (UHF) method 27can also be performed by invoking the DFT module with appropriate 28keywords. 29 30Correlation energy of solvent molecules may also be evaluated at 31\begin{enumerate} 32\item MP2, 33\item CCSD, 34\item CCSD+T(CCSD), 35\item CCSD(T), 36\end{enumerate} 37levels of theory. It is cautioned, 38however, that these correlated COSMO calculations determine 39the solvent reaction field using the HF charge distribution of 40the solute rather than the charge distribution of the correlation 41theory and are not entirely self consistent in that respect. 42In other words, these calculations assume that the correlation 43effect and solvation effect are largely additive, and the combination 44effect thereof is neglected. 45COSMO for MCSCF has not been implemented yet. 46 47In the current implementation the code 48calculates the gas-phase energy of the system followed by the 49solution-phase energy, and returns the electrostatic contribution 50to the solvation free energy. 51At the present gradients are calculated by finite 52difference of the energy. Known problems include that the code does not 53work with spherical basis functions. 54The code does not calculate the 55non-electrostatic contributions to the free energy, except for 56the cavitation/dispersion contribution to the solvation free energy, 57which is computed and printed. 58It should be noted that one must in general take into account 59the standard state correction besides the electrostatic 60and cavitation/dispersion contribution to the solvation free energy, 61when a comparison to experimental data is made. 62 63Invoking the COSMO solvation model is done by specifying the input 64COSMO input block with the input options as: 65\begin{verbatim} 66cosmo 67 [off] 68 [dielec <real dielec default 78.4>] 69 [radius <real atom1> 70 <real atom2> 71 . . . 72 <real atomN>] 73 [rsolv <real rsolv default 0.00>] 74 [iscren <integer iscren default 0>] 75 [minbem <integer minbem default 2>] 76 [maxbem <integer maxbem default 3>] 77 [ificos <integer ificos default 0>] 78 [lineq <integer lineq default 1>] 79end 80\end{verbatim} 81followed by the task directive specifying the wavefunction and 82type of calculation, e.g., \verb+task scf energy+, \verb+task mp2 energy+, 83\verb+task dft optimize+, etc. 84 85\verb+off+ can be used to turn off COSMO in a compound (multiple task) 86run. By default, once the COSMO solvation model has been defined it will 87be used in subsequent calculations. Add the keyword \verb+off+ if COSMO 88is not needed in subsequent calculations. 89 90\verb+Dielec+ is the value of the dielectric constant of the medium, 91with a default value of 78.4 (the dielectric constant for water). 92 93\verb+Radius+ is an array that specifies the radius of the spheres 94associated with each atom and that make up the molecule-shaped cavity. 95Default values are Van der Waals radii. Values are in units of angstroms. 96The codes uses the following Van der Waals radii by default: 97\begin{verbatim} 98 data vdwr(103) / 99 1 0.80,0.49,0.00,0.00,0.00,1.65,1.55,1.50,1.50,0.00, 100 2 2.30,1.70,2.05,2.10,1.85,1.80,1.80,0.00,2.80,2.75, 101 3 0.00,0.00,1.20,0.00,0.00,0.00,2.70,0.00,0.00,0.00, 102 4 0.00,0.00,0.00,1.90,1.90,0.00,0.00,0.00,0.00,1.55, 103 5 0.00,1.64,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00, 104 6 0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00, 105 7 0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00, 106 8 0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00, 107 9 0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00, 108 1 0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,1.65, 109 2 0.00,0.00,0.00/ 110\end{verbatim} 111with 0.0 values replaced by 1.80. Other radii can be used as well. 112See for examples: 113 114\begin{enumerate} 115\item E. V. Stefanovich and T. N. Truong, Chem.~Phys.~Lett. 244, 65 (1995). 116\item V. Barone, M. Cossi, and J. Tomasi, J.~Chem.~Phys. 107, 3210 (1997). 117\end{enumerate} 118 119\verb+Rsolv+ is a parameter used to define the solvent accessible 120surface. See the original reference of Klamt and Schuurmann for a 121description. The default value is 0.00 (in angstroms). 122 123\verb+Iscren+ is a flag to define the dielectric charge scaling option. 124``{\tt iscren 1}'' implies the original scaling from Klamt and Sch\"{u}\"{u}rmann, 125mainly ``$(\epsilon-1)/(\epsilon+1/2)$'', where $\epsilon$ is the dielectric constant. 126``{\tt iscren 0}'' implies the modified scaling suggested by Stefanovich and 127Truong, mainly ``$(\epsilon-1)/\epsilon$''. Default is to use the modified scaling. 128For high dielectric the difference between the scaling is not 129significant. 130 131The next three parameters define the tesselation of the unit sphere. 132The approach follows the original proposal by Klamt and Sch\"{u}\"{u}rmann. 133A very fine tesselation is generated from \verb+maxbem+ refining 134passes starting from either an octahedron or an icosahedron. The 135boundary elements created with the fine tesselation are condensed 136down to a coarser tesselation based on \verb+minbem+. The induced 137point charges from the polarization of the medium are assigned to 138the centers of the coarser tesselation. Default values are 139``{\tt minbem 2}'' and ``{\tt maxbem 3}''. The flag \verb+ificos+ serves to 140select the original tesselation, ``{\tt ificos 0}'' for an octahedron 141(default) and ``{\tt ificos 1}'' for an icoshedron. Starting from an icosahedron 142yields a somewhat finer tesselation that converges somewhat faster. 143Solvation energies are not really sensitive to this choice for 144sufficiently fine tesselations. 145 146The \verb+lineq+ parameter serves to select the numerical algorithm to solve 147the linear equations yielding the effective charges that represent 148the polarization of the medium. ``{\tt lineq 0}'' selects an iterative method 149(default), ``{\tt lineq 1}'' selects a dense matrix linear equation solver. 150For large molecules where the number of effective charges is large, 151the codes selects the iterative method. 152 153The following example is for a water molecule in `water', using 154the HF/6-31G** level of theory: 155 156\begin{verbatim} 157start 158echo 159 title "h2o" 160geometry 161o .0000000000 .0000000000 -.0486020332 162h .7545655371 .0000000000 .5243010666 163h -.7545655371 .0000000000 .5243010666 164end 165basis segment cartesian 166 o library 6-31g** 167 h library 6-31g** 168end 169cosmo 170 dielec 78.0 171 radius 1.40 172 1.16 173 1.16 174 rsolv 0.50 175 lineq 0 176end 177task scf energy 178\end{verbatim} 179