# 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 [Cammi2005] for a good review on the subject).

The general picture of the method can be depicted as:

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 (\(\Delta G_{ENP}\)) and cavity-dispersion (\(\Delta G_{CDS}\)):

and the methods differ on how to compute these terms.

## Conductor-like Polarizable Continuum Model (CPCM)#

In the CPCM method [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 (\(\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 (\(\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:

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, \(\Delta G_{ENP}\) corresponds to CPCM Dielectric and \(\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 \(\Delta G_{ENP}\) (*CPCM Dielectric*) and \(\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 [Karplus1999] [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 [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 \(\Delta G_{CDS}\) (*Gcds*) term, together with the \(\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 (\(P_{o/w}\)) for the drug Diazepam:

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 \(\Delta G^o_{o/w} = \Delta G^o_{o} - \Delta G^o_{w}\) as [Truhlar2009]:

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:

After making the calculations, one obtains an energy difference of \(3.41 kcalmol^{-1}\), that corresponds to a \(logP_{o/w}=2.50\), which is not far from the experimental value of \(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:

that arises when changing from gas phase at \(1 atm\) and \(298 K\) (which is how the \(G^o\) is calculated when using FREQ by default) to a solution phase at \(1 molL^{-1}\), and that is [Cramer2004]:

Warning

Do not forget to add this term to your \(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
```