.. _solvation:
Implicit Solvation Models
=========================
One aproach to include solvent effects on your calculation is through the so-called "implicit solvent models", and here we will discuss two main options available in ORCA.
Overview
--------
In these models, the solute is placed in a cavity of roughly molecular shape. The solvent is described by a continuum that interacts with the charges on the cavity surface, which are in turn determined by the solute and the problem is solved iteratively (see :ref:`bib:[Cammi2005]` for a good review on the subject).
The general picture of the method can be depicted as:
.. image:: CPCM/surfaces.png
:align: center
:width: 450
where the "probe" is an ideal solvent molecule and the molecular cavities are defined from:
* "Van der Waals surface" (vdW), built from the Van der Waals radii (the default in ORCA);
* "Solvent-Accessible Surface" (SAS), defined by the center of the probe;
* or the "Solvent-Excluded Surface" (SES), obtained if you follow the inward part of the solvent probe sphere rolling over vdW surface.
The solvation energy is then decomposed in two main terms, electrostatic (:math:`\Delta G_{ENP}`) and cavity-dispersion (:math:`\Delta G_{CDS}`):
.. math::
\Delta G_{solv}^o = \Delta G_{ENP} + \Delta G_{CDS}
and the methods differ on how to compute these terms.
Conductor-like Polarizable Continuum Model (CPCM)
-------------------------------------------------
In the CPCM method :ref:`bib:[Cossi1998]`, the bulk solvent is treated as a conductor-like polarizable continuum and the main parameters to define the method are the **refractive index** and the **dielectric constant** of the medium.
The electrostatic contribution (:math:`\Delta G_{ENP}`) that arises from the interaction of the medium and the molecular surface charges is included in the SCF calculation - so that you even get "solvated" orbitals - and the cavity term (:math:`\Delta G_{CDS}`) can be obtained from more complicated schemes.
ORCA has a list of predefined solvents that can be called by using::
!CPCM(solvent)
where solvent can be any from:
.. table::
:align: center
=============== ============== =============
Solvent Name Dielec. Const. Refrac. Index
=============== ============== =============
Water 80.4 1.33
Acetone 20.7 1.359
Acetonitrile 36.6 1.344
Ammonia 22.4 1.33
Benzene 2.28 1.501
CCl4 2.24 1.466
CH2Cl2 9.08 1.424
Chloroform 4.9 1.45
Cyclohexane 2.02 1.425
DMF 38.3 1.430
DMSO 47.2 1.479
Ethanol 24.3 1.361
Hexane 1.89 1.375
Methanol 32.63 1.329
Octanol 10.3 1.421
Pyridine 12.5 1.510
THF 7.25 1.407
Toluene 2.4 1.497
=============== ============== =============
If you want to include any other solvent, one can always input these parameters manually using::
%CPCM EPSILON 80.4
REFRAC 1.33
END
More options and details can be found on the `ORCA manual `_. One example for the single-point calculation of an Aspirin molecule in water is::
!BP86 DEF2-SVP CPCM(WATER)
* XYZFILE 0 1 aspirin.xyz
.. note::
Both the analytic gradients and Hessians are available for the CPCM, so that can be used together with the OPT and FREQ keywords as well!
Then, before the SCF, a header with all the CPCM-related information is printed::
--------------------
CPCM SOLVATION MODEL
--------------------
CPCM parameters:
Epsilon ... 80.4000
Refrac ... 1.3300
Rsolv ... 1.3000
Surface type ... GAUSSIAN VDW
Epsilon function type ... CPCM
Radii:
Radius for C used is 3.8550 Bohr (= 2.0400 Ang.)
Radius for O used is 3.4469 Bohr (= 1.8240 Ang.)
Radius for H used is 2.4944 Bohr (= 1.3200 Ang.)
Calculating surface ... done! ( 0.0s)
GEPOL surface points ... 1160
GEPOL Volume ... 1403.2346
GEPOL Surface-area ... 757.5566
Calculating surface distance matrix ... done! ( 0.0s)
Performing Cholesky decomposition & store ... done! ( 0.1s)
Overall time for CPCM initialization ... 0.2s
after the SCF calculation ends, there will be a the summary of the CPCM-related information::
Components:
Nuclear Repulsion : 765.01084602 Eh 20817.00344 eV
Electronic Energy : -1413.26452682 Eh -38456.88288 eV
One Electron Energy: -2400.58284083 Eh -65323.18007 eV
Two Electron Energy: 987.31831400 Eh 26866.29718 eV
CPCM Dielectric : -0.01874113 Eh -0.50997 eV
[...]
CPCM Solvation Model Properties:
Surface-charge : -0.01502552
Charge-correction : -0.00002872 Eh -0.00078 eV
Free-energy (cav+disp) : This term is not implemented in the current solvation scheme
As one can see, **dielectric contribution** to the energy, the solute **net charge**, the **charge** and **cavity corrections** are printed separately (the latter in this case is not calculated). The `Charge-correction` term is not included in the SCF energy, and is just printed for information purposes.
The *FINAL SINGLE POINT ENERGY* now already includes all computed solvation terms. In terms of the equations shown above, :math:`\Delta G_{ENP}` corresponds to `CPCM Dielectric` and :math:`\Delta G_{CDS}` to the `Free-energy (cav+disp)`.
.. important::
There is an important difference with respect to previous ORCA versions. Now the cavity term is **not** calculated or added for regular CPCM by default.
.. note::
For more information on cavity terms and CPCM, please consult the `ORCA manual `_.
Combining DFT with post-HF methods
..................................
Since the HF density is of lower quality than those obtained from DFT, so is the CPCM correction obtained from HF calculations.
What one usually can do is to take the :math:`\Delta G_{ENP}` (*CPCM Dielectric*) and :math:`\Delta G_{CDS}` (*Free-energy (cav+disp)*) from some calculation like DFT and add that to the total energy obtained with that method, e.g. if you want to combine CPCM with DLPNO-CSSD.
Gaussian point charges
......................
The usual point charge scheme might lead to instabilities in the energy, e.g. if two points end up too close. In ORCA, we use a smeared Gaussian charge to obtain a smoother potential energy surface :ref:`bib:[Karplus1999]` :ref:`bib:[Neese2020]` instead. From ORCA 5, the default is::
%CPCM SURFACETYPE VDW_GAUSSIAN END
But a SES surface can also be chosen with::
%CPCM SURFACETYPE GEPOL_SES_GAUSSIAN END
For more detailed information, please consult the `ORCA manual `_
Universal Solvation Model (SMD)
-------------------------------
The SMD method :ref:`bib:[Truhlar2009]` can be thought as an improvement over the CPCM, since it uses the full solute electron density to compute the cavity-dispersion contribution instead of the area only.
This method requires more parameters, which makes it less flexible for unknown solvents, but ORCA has currently 179 solvents available (`ORCA manual `_)! In order to use it, add::
%CPCM SMD TRUE
SMDSOLVENT "SOLVENT"
END
to the input. The same example from before would be::
!BP86 DEF2-SVP
%CPCM SMD TRUE
SMDSOLVENT "WATER"
END
* XYZFILE 0 1 aspirin.xyz
The initial output is similar to CPCM, except that one now has the SMD descriptors::
SMD-CDS solvent descriptors:
Soln ... 1.3328
Soln25 ... 1.3323
Sola ... 0.0000
Solb ... 0.0000
Solg ... 0.0000
Solc ... 0.0000
Solh ... 0.0000
and the output before the the *TOTAL SCF ENERGY* header has::
CPCM Dielectric : -0.02952804 Eh -0.80350 eV
SMD CDS (Gcds) : 0.01088483 Eh 0.29619 eV
One can use this :math:`\Delta G_{CDS}` (*Gcds*) term, together with the :math:`\Delta G_{ENP}` (*CPCM Dielectric*) to correct energies from further calculations.
Example - Octanol/Water Partition Coefficient
---------------------------------------------
We will use SMD and try to predict the octanol/water partition coefficient (:math:`P_{o/w}`) for the drug Diazepam:
.. math::
P_{o/w} = \frac {[solute]_{octanol}} {[solute]_{water}}
That can be obtained from the relationship between the equilibrium constant and the Gibb's free energy difference between the solute in octanol and water :math:`\Delta G^o_{o/w} = \Delta G^o_{o} - \Delta G^o_{w}` as :ref:`bib:[Truhlar2009]`:
.. math::
logP_{o/w} = \frac {-\Delta G^o_{o/w}} {2.303RT}
First, we optimize both geometries in each solvent, for instance using::
!B3LYP DEF2-SVP OPT NUMFREQ D4
%CPCM SMD TRUE
SMDSOLVENT "1-OCTANOL"
END
* XYZFILE 0 1 diazepam.xyz
or::
!B3LYP DEF2-SVP OPT NUMFREQ D4
%CPCM SMD TRUE
SMDSOLVENT "WATER"
END
* XYZFILE 0 1 diazepam.xyz
to get the Gibbs free energy of each solvated compound and its geometry. Here is an image of the molecule with a VdW cavity surface in blue:
.. image:: CPCM/cavity.png
:align: center
:width: 400
After making the calculations, one obtains an energy difference of :math:`3.41 kcalmol^{-1}`, that corresponds to a :math:`logP_{o/w}=2.50`, which is not far from the experimental value of :math:`logP_{o/w}^{exp}=2.99`, considering the simplicity of the method.
.. note::
Since there are no analytic frequencies using the SMD model yet, we are using !NUMFREQ. It might take long if you do not parallelization here.
Changing from Gas to Liquid phase
---------------------------------
The equation given at the beginning for the solvation free energy was not complete, it actually has a third term:
.. math::
\Delta G_{solv}^o = \Delta G_{ENP} + \Delta G_{CDS} + \Delta G^o_{conc}
that arises when changing from gas phase at :math:`1 atm` and :math:`298 K` (which is how the :math:`G^o` is calculated when using FREQ by default) to a solution phase at :math:`1 molL^{-1}`, and that is :ref:`bib:[Cramer2004]`:
.. math::
\Delta G^o_{conc}= RTln(24.5) = 1.89 kcalmol^{-1}.
.. warning::
Do not forget to add this term to your :math:`G^o` when predicting solution thermodynamics. The absence of this correction can lead to large errors!
CPCM and COSMO
--------------
ORCA does **not** have any COSMO implementation, but one can use the COSMO epsilon function by using the CPCMC(solvent) keyword, with an extra "C"::
!B3LYP DEF2-SVP CPCMC(WATER)
Starting structures
-------------------
Aspirin
.......
::
C -2.64076 2.23326 0.00014
C -3.28499 1.00146 -0.00001
C -2.53276 -0.17323 -0.00006
C -1.24586 2.29692 0.00023
C -0.48687 1.12113 0.00013
C -1.12591 -0.13504 0.00002
C -0.44272 -1.46662 -0.00003
O -1.06358 -2.51853 -0.00017
O 0.90426 -1.44815 0.00010
O 0.94078 1.15087 0.00026
C 1.78023 2.27870 -0.00045
O 1.43411 3.44948 -0.00130
C 3.21351 1.83099 0.00005
H 3.28673 0.73990 -0.00017
H 3.70751 2.20906 0.89847
H 3.70821 2.20936 -0.89786
H 1.28154 -0.55226 0.00027
H -3.22416 3.15173 0.00019
H -4.37149 0.95196 -0.00007
H -0.79315 3.28276 0.00040
H -3.05624 -1.12937 -0.00015
Diazepam
........
::
C -2.76533 -0.90855 1.85461
C -4.10278 -0.72362 2.22052
C -5.02528 -0.30117 1.27344
C -4.62706 -0.08566 -0.04023
C -3.28690 -0.29486 -0.41749
C -2.32208 -0.67243 0.53842
Cl -6.66626 -0.05897 1.72353
N -0.93200 -0.82351 0.23569
C -0.24076 -0.01864 -0.67861
C -1.02949 1.08452 -1.37013
C -2.96608 -0.16369 -1.86911
N -1.95301 0.48903 -2.32972
H -0.32054 1.71528 -1.91810
H -1.52297 1.73752 -0.64155
O 0.97434 -0.13069 -0.86774
C -0.12886 -1.74693 1.02668
H 0.21402 -1.23663 1.93228
H 0.75039 -2.08063 0.46656
H -0.70989 -2.63829 1.28413
C -3.87195 -0.81596 -2.86776
C -4.45465 -2.06403 -2.61294
C -4.11329 -0.17301 -4.08765
C -5.29672 -2.64689 -3.56241
C -5.55029 -1.99264 -4.76807
C -4.95622 -0.75907 -5.03306
H -3.64497 0.78461 -4.30431
H -4.24920 -2.60250 -1.69208
H -5.74819 -3.61652 -3.36692
H -6.20254 -2.45023 -5.50807
H -5.14291 -0.25572 -5.97837
H -5.36297 0.23683 -0.77401
H -4.40857 -0.89624 3.24980
H -2.07050 -1.19831 2.63954