Polarization catastrophe
C9926: Problems and Issues of Molecular Modeling
Rahul Deb
CEITEC MUNI
2021
Thole, B. Th. Chemical Physics 59.3 (1981): 341-350.
Goal: To accurately calculate the effects of polarizability on electrostatic interaction (stabilization
of a charged or a dipolar molecule surrounded by a polarizable solvent)
Applequist et al. used isotropic atom polarizabilities, which is defined as the ratio of the
induced dipole moment of an atom to the electric field, that produces this dipole moment
Birge showed that the large anisotropy can be removed by inclusion of the effect of electron
repulsion, which essentially decreases the polarizability component in the direction of the bond
This model predicts mean polarizabilities very well but that the predicted anisotropy is too large
Birge used 14 parameters to describe the atomic polarizabilities in different chemical
environments and their tendencies to resist excessive polarizability along the bonds
However, the transferability of these parameters to other classes of molecules is bad
History
The atomic polarizability tensor remains isotropic
Thole’s approach
Molecule is considered as an arrangement of N atoms
each of which has a polarizability.
Induced dipole moment μp at atom p can be calculated
as a function of the applied electric field, E, at atom p
α is the atomic polarizability tensor of atom p and
Tp is the dipole field tensor
rpq the distance between atoms p and q, and
x, y and z are the cartesian components of the vector
connecting atoms p and q
This is in contrast with Birge’s model, which modifies atomic polarizability tensor keeping the field tensor unchanged
Point dipole -> smeared-out dipoles
Problem
Infinite molecular polarizabilities
For a diatomic molecule AB, the values of the polarizabilities parallel
and perpendicular to the bond axis can be derived as follows:
Head Tail
When r approaches (4αAαB)^1/6 , αII goes to infinity
Two induced dipoles
Bond
When the point dipole interaction atoms come closer than a certain limit, the polarizability becomes infinite,
which is unreal Polarization catastrophe
a dipole is thought to be built up from two infinitesimally shifted monopoles, which in this case will not be
point charges, but two charge distributions
The corresponding monopole-monopole interaction of some “well-behaved” charge distribution of
density ρ is
Solution -> the modification of the dipole field
tensor, T, so that it does not behave as the third
power of the distance, r, at small r
Thus, the unit point charge giving the point dipole equations can be considered “smeared out”
the field at u for a spherically symmetric charge
density is equal to that of the charge contained in the
sphere with radius u, concentrated at the center
And ρ is finite when u approaches zero
The polarizabilities of a wide range of chemically different molecules were calculated using for
each atom a polarizability independent of its chemical environment Model1:
point dipole interaction,
Model2: anisotropic atom point dipole interaction model,
Model3: modified dipole interaction.
When two atoms approach each other from infinity the interaction between their dipoles induced by an
external field results in two functions:
which govern the polarizability along the bond and perpendicular to it.
Along the bond is amplified and perpendicular is damped, keeping the mean polarizability approximately
conserved.
When the atoms penetrate, the interaction is effectively damped with respect to the point dipole case.
When two atoms touch, a new body is formed with an elongated shape giving anisotropic polarizability,
as in classical electrostatics.
When the atoms penetrate further the shape becomes less elongated until at u=0 a composite spherical
atom is formed.
Verified !!!
Drude Oscillator Model
Drude oscillators to carbon and oxygen
atoms (qDC and qDO) via harmonic springs
LP: Lone pairs on the oxygen atom
α11, α22, and α33 are the tensor components
of the anisotropic atomic polarizability
The isotropic polarizability, α, of a given atom is described by the distribution
of the total charge on the atom, q, between the core atom and its Drude
oscillator, which is attached to the core atom via a harmonic spring
The charges assigned to the Drude
oscillator (qD) and to the core atom (qA)
are such that total charge q = qA + qD
kDis the force constant of the
Drude−atom harmonic bond
displacement, d, between the Drude oscillator and
the core atom in response to an electric field, E,
the induced atomic dipole
harmonic self-polarization term
anisotropic form is
achieved by
Lemkul, Justin A., et al. Chemical reviews 116.9 (2016): 4983-5013.
1−2
These interactions are treated explicitly
via Thole screening in the Drude-2013
force field, contributing to anisotropic
molecular polarizability.
Explicit inclusion of dipole−dipole interactions for atoms within three bonds
1−3
a is a damping constant
Thole smeared charge distribution
Subsequently, atom-specific values were
implemented, as they may be assigned to
fine-tune the molecular polarizability tensor
a = ai + aj
Directional response of the Drude oscillators to external electric fields (E)
due to the 1−2 dipole−dipole interactions when
the electric field is perpendicular to the bond,
the 1−2 dipoles damp each other, leading to a
smaller polarization response compared to the
electric field being parallel to the bond leading
to the 1−2 dipoles enhancing each other.
Molecular Dynamics Algorithm
In the induced point dipole representation, each
atomic site is assigned a dipole that is calculated using
an iterative self consistent field (SCF) procedure.
the positions of the Drude oscillators are solved
via energy minimization at each integration step
computationally demanding
Solution: Lagrangian approachwhere each Drude oscillator is
ascribed a small mass (mD) that is subtracted from the parent
atom such that the total mass of the Drude−atom pair
remained equal to the atomic mass.
Velocity−Verlet algorithm is used to numerically
integrate the equations of motion calculating
the forces on the Drude−atom centers of
mass (Ri) and Drude−atom displacements (di)
Polarization catastrophe can occur when the
forces acting on the Drude oscillator cause
large displacement from its core atom
The constraint reflects the Drude oscillator back
toward its core atom, along the Drude−atom
bond, if it reaches a predefined limit
prevented by using “hard-wall constraint”
Divalent ions are prone to “polarization catastrophe” due to a
Coulombic singularity that can arise when the Drude particle
on the ion interacts with the core charge on the water oxygen.
circumvented by using a through-space Thole screening same
as used for bonded dipoles
1
2