Polarization & Polarization Models
Background:
- Stone's //Intermolecular Forces// for a fuller and more mathematical description of polarization.
Polarization is a classical phenomena. It is the energy lowering by the self-consistent polarization of the atoms/molecules by the permanent and induced moments. Consider just two interacting molecules: The permanent moments on one (if non-zero) will create a field at the other and polarize it, thereby altering the electornic charge density (we work in the Born-Oppenheimer approximation so the nuclei are not distorted in this process, though, in reality they are) and consequently the multipole moments of the partner. This process happens both ways, so we end up with a mutual self-consistent polarization. For details please see Stone's //Intermolecular Forces//.
The many-body nature of the polarization is manifest in the following example taken from Stone's book: Consider a polarizable atom inbetween two identical dipoles. If considered as the sum over pairs of interactions, the central atom will always be polarized as it will always experience a non-zero dipole field $F$ from each of the dipoles leading to a polarization energy $2 \times -\frac{1}{2}\alpha F^2$, where $\alpha$ is the polarizability of the central atom. However now consider the following cases. If both dipoles point in the same direction, the atom experiences a large non-zero field $2F$ and gets polarized with polarization energy $4 \times -\frac{1}{2} \alpha F^2$. But if the dipoles point in opposite directions the field at the atom is exactly zero so we don't get any polarization. To know the polarization state of the system we therefore need to know the interactions of all three units.
Points to note:
- From a pair-wise approach we get a qualitatively different polarization energy compared with the many-body approach.
- The many-body polarization can be (in this example) twice as large as the pair-wise summations would suggest. This would lead to stronger intermolecular binding and consequently to shorter intermolecular bonds.
- At short separations the molecular charge densities overlap.
- Consequently the electric fields calculated from molecular multipoles needs to be modified to include charge-density overlap effects.
Likewise, the small $R$ mathematical divergence of the multipole expansion needs to eliminated.
- These two effects are often taken care of by a //damping function//.
- We typically use the Tang-Toennies damping functions. These are incomplete Gamma functions:
$f_n(\beta R) = 1 - \exp(-\beta R) \sum_{k=0}^{n}\frac{(\beta R)^{k}}{k!}$
- These are known to work very well for the dispersion expansion.
But they require knowledge of the damping parameters (which may not be constants!) $\beta$.
- These parameters could be distance, angular and site-site dependent. Though, at present, most groups assume they are dependent on the pair of interacting molecules only.
Misquitta & Stone (2008) had argued that we should have $\beta = \sqrt{2 I_{\rm A}} + \sqrt{2 I_{\rm B}}$, where $I_{\rm A/B}$ are the vertical ionization potentials of the two molecules. The argument was based on the asymptotic form of the potential arising from the charge-density.
While Misquitta & Stone have shown that this choice works very well for the dispersion expansion, there is evidence that it is unsuitable for the polarization expansion.
- The reason for this is subtle and will be described below.
Consequently, we need to determine $\beta$ for the polarization expansion.
It does matter what value $\beta$ takes with the choice effecting the many-body energies. Once again, the reason for this is subtle.
Basically it comes down to the Charge-Transfer energy: the polarization expansion has typically been constructed to match the induction energy. But this is wrong as the induction energy consists of a charge-transfer and true polarization component. All that the polarization model should match is the polarization component of the induction. However, this separation is not made in standard versions of SAPT and is the subject of much research. The most promising separation method is a new one based on R-SAPT that is described in Charge-Transfer.
Quantities needed for the polarization:
- Distributed polarizabilities calculated using the WSM method.
- Distributed multipole moments calculated using the GDMA method.
- A suitable damping model.
The first two, the polarizabilities and moments, are readily calculated using the CamCASP program. This part of the calculation will be described elsewhere. The last, the damping model, is the main focus of this page.
Damping model
//Why damp the polarization model?// The polarization expansion is a multipole expansion in powers of $1/R$ and therefore diverges as $R \rightarrow 0$, i.e., as molecular separation decreases. While we will generally not see small values of $R$, if the system is subject to pressure or if it exhibits a good deal of strong hydrogen bonds or ionic bonds, intermolecular separations can get quite a lot smaller than would be expected from equilibrium geometries. This is primarily due to the many-body nature of the polarization. This is well described in Stone's book and won't be elaborated upon here.
Note
Insert water cluster example illustrating this phenomenon.
As mentioned above, the polarization expansion needs to be damped to
suppress the mathematical singularity in the $1/R$ expansion, and
- to account for charge-density overlap effects.
But the damping needs to be chosen so that the polarization expansion reproduces only the polarization part of the induction energy. This is very important. Let us denote by ${\rm Pol}^{(n)}$ the polarization energy at order $n$ as arising from the polarization expansion. Likewise, let ${\rm POL}^{(n)}$ denote the polarization part of the induction energy at order $n$, i.e., this is the energy from SAPT. And ${\rm CT}^{(n)}$ is the charge-transfer energy at order $n$. We have: $$ E^{(n)}_{\rm ind} = {\rm CT}^{(n)} + {\rm POL}^{(n)}.$$ That is, at each order, starting from $n=2$, the induction energy consistes of CT and POL parts. What we need to do with our damping is to determine it in such a way that: $$ {\rm Pol}^{(n)} = {\rm POL}^{(n)}.$$ If ${\rm Pol}^{(n)}$ included any part of the CT then we would have spurious CT components included in the many-body polarization energy. This could lead to non-sense results and overbinding.
Consequently, the first step towards defining the damping is the separation of ${\rm CT}^{(2)}$ and ${\rm POL}^{(2)}$. This is described in Charge-Transfer. At this stage we make this separation at second-order only.
Once the ${\rm CT}^{(2)}$ and ${\rm POL}^{(2)}$ have been separated, the damping can be obtained by a standard fitting procedure. Note that we will typically need a large number of dimer configurations to make a good fit to define the $\beta$ parameter. Also, if we wish to make a more complex model with damping parameters that depend on site pairs, then the data required needs to be commensurately larger. The steps are as follows:
- Define your polarization model: rank, multipoles. A good model will contain WSM polarizabilities to rank 2 or 3 and distributed multipoles to rank 4.
Using this model, evaluate ${\rm Pol}^{(2)}$ for at the dimer configurations. You will need to specify a damping parameter $\beta$ as **Orient** cannot yet optimize the damping parameters.
Optimize $\beta$ (by hand or by using an external optimizer) to get the best match between ${\rm Pol}^{(2)}$ and ${\rm POL}^{(2)}$.
A good choice for $\beta$ seems to be 75% of the $\beta = \sqrt{2 I_{\rm A}} + \sqrt{2 I_{\rm B}}$ value.
Typically $\beta$ needs to be close to $1.0$ for pairs of electronegative atoms, close to $2.0$ for pairs of hydrogen atoms, and close to $1.5$ for a mix of the two. These are only rules of thumb.
Once the damping parameters have been obtained, we can proceed to use the polarization model to calculate the full, infinite-order CT and POL.
Instructions for using **Orient** to calculate energies can be found here.
Inifinte-Order POl & CT
To get the inifinite-order polarization energy all we need do is calculate ${\rm Pol}^{(2-\infty)}$ usind the polarization model with the damping parameters obtained in the previous section. This is done by iterating the polarization model to convergence.
As described in the section on Charge-Transfer, the infinite order CT is defined as: $$ {\rm CT}^{(2-\infty)} = E^{(2-\infty)}_{\rm ind} - {\rm Pol}^{(2-\infty)} \\ \approx E^{(2)}_{\rm ind} + \delta^{\rm HF}_{\rm int} - {\rm Pol}^{(2-\infty)}$$ where $\delta^{\rm HF}_{\rm int}$ is the SAPT approximation to the higher-order interation energy which is known to be mostly higher-order induction.
Instructions for the infinite-order calculations with Orient are here.