Pitzer-Debye-Hückel long-range electrostatic correction
*******************************************************

Systems with charged species often demonstrate behavior that is not
easily captured with standard thermodynamic models meant for neutral
molecules.  A major reason for this that the chemical potential
requires an additional contribution due to the electrostatic
potential.  The electrostatic potential decays with r^{-1}, where r is
the distance between two charges.  This inverse dependence on r makes
the electrostatic potential a long-range term, especially as compared
to other intermolecular energetic contributions, which are typically
only relevant at short distances.  One successful model for adding
necessary electrostatic corrections for systems with charged species
is the Pitzer-Debye-Hückel (PDH) [1] model.  In the COSMO-RS package,
we use a modification of the PDH term that allows for mixed solvents.


Mixing rules and required property inputs for the PDH term
==========================================================

In the original publication, the PDH model is used only for pure
solvents.  To generalize to mixed solvent systems, we have included
mixing rules to estimate required parameters for the PDH model.  These
are discussed below in addition to required parameters.


Molecular weight
----------------

No additional input is needed for this property as this is directly
determined by the atomic composition of a molecule.  The average
molecular weight for the entire system, M, is given as a function of
all the individual molecular weights m_i :

   M = \frac{\sum\limits_i n_i m_i}{\sum\limits_i n_i}

where n_i is the number of moles of species i.


Density
-------

The density of each compound is also required for the calculation of
the PDH term.  For ionic species, one may simply use the density of
the combined salt.  The average density for the entire system, D, is
given as a function of all the individual densities d_i :

   D = \frac{\sum\limits_i n_i}{\sum\limits_i \frac{n_i}{d_i}}


Dielectric constant
-------------------

Additionally, the dielectric constant of each compound is required for
the PDH term.  These terms do not need to be experimentally accurate
numbers, but should be approximate.  The dielectric constant of a
combined salt can be used for the individual ionic species.  The
average dielectric constant for the entire system, \mathcal{E}, is
given as a function of all the individual dielectric constants
\epsilon_i :

   \mathcal{E} = \frac{\sum\limits_i n_i \frac{m_i
   \epsilon_i}{d_i}}{\sum\limits_i n_i \frac{m_i}{d_i}}


Derivation of the PDH term for general mixtures
===============================================

In [1] , the Debye-Hückel model is adapted to account for higher
concentrations of ions (the Debye-Hückel model works well for low
concentrations).  The authors choose from a variety of possible
models, ultimately deciding on one that is deemed most accurate
empirically.  The expression for the PDH excess energy (long range
electrostatic correction) is given below.

   \frac{G^{PDH}}{RT} = -4 \left({\frac{1000}{M}}\right)^{1/2} \frac{N
   A_x I_x}{\rho} \ln( 1+\rho I_x^{1/2} )

where I_x is the ionic strength of the system, which is a function of
the charge z_i of each species:

   I_x = \frac{1}{2} \frac{\sum\limits_i n_i z_i^2}{\sum\limits_i n_i}

and A_x is the Debye-Hückel parameter, defined as:

   A_x = \frac{1}{3} \frac{(2\pi N_A D e^6 )^{1/2}}{(4 \pi \epsilon_0
   \mathcal{E} kT)^{3/2}}

where e is the charge of the electron, N_A is Avogadro's number, k is
the Boltzmann constant, and \epsilon_0 is the permittivity of free
space.

Finally, \rho is the closest approach parameter.  The specifics of
this parameter are still under development, and for now the following
simple formula from [1] is used.

   \rho = 2150 \left( \frac{D}{ \mathcal{E} T } \right)^{1/2}

Finally, we can calculate the PDH contribution to the activity
coefficient as follows:

   \ln(\gamma_i^{PDH}) =& \frac{\delta}{\delta n_i} \left(
   \frac{G^{PDH}}{RT} \right)_{T,P} \\ =& -4
   \left({\frac{1000}{M}}\right)^{1/2} \frac{A_x I_x}{\rho} \Bigg[
   \ln( 1+\rho I_x^{1/2}) \bigg( \frac{M-m_i}{2M} + \frac{z_i^2}{2I}
   \\ & +  \frac{1}{2}\left(1-\frac{D}{d_i}\right) +\frac{3N}{2}
   \frac{m_i}{d_i \sum\limits_j n_j \frac{m_j}{d_j}}
   \left(1-\frac{\epsilon_i}{\mathcal{E}}\right) -\frac{N}{\rho}
   \frac{\delta\rho}{\delta n_i} \bigg) \\ & + \frac{1}{1+\rho
   I_x^{1/2}} \left( \frac{\rho I_x^{1/2}}{2N} \left(
   \frac{z_i^2}{2I}-1 \right) + I_x^{1/2} \frac{\delta\rho}{\delta
   n_i} \right) \Bigg]

Note that in the above, the analytical form of the term
\frac{\delta\rho}{\delta n_i} is not used to simplify the expression.


Tutorial on using the PDH correction
====================================

There is a tutorial on using the PDH correction that demonstrates the
use of this correction with the GUI.

-[ References ]-

[1] Pitzer, Kenneth S., and John M. Simonson. *Thermodynamics of
    multicomponent, miscible, ionic systems: theory and equations.*
    The Journal of Physical Chemistry 90 (1986): 3005-3009.
