De glaciēbus or deductive molecular mechanics of ice polymorphs

Abstract

Yet in 1960s Del Re with coworkers considered the electronic structure of organic molecules using hybrid orbitals. They studied the simplest molecules: CH₄, NH₃, H₂O, and more complex ones: cyclopropane, cyclobutane, cubane, using either the optimal overlap or maximal localization principles to determine the hybrids. Later Malrieu with coworkers used hybrid orbitals in the PCILO method. Later, we determined either the form and orientation of the hybrid orbitals or two-electron functions of the two-center bonds constructed on the basis of these hybrids from the minimum condition for total electronic energy as implemented in the SLG method. This gave us significant improvement in the efficiency: the dependence of the required computational resources on the molecule size reduces down to O(N). The paradigm based on the usage of the variation principle for determination of either the hybrid orbitals or the elements of the reduced density matrices in their basis allows one to formulate and prove exact statements about electronic structure. We start from establishing the energy expression for highly symmetric non-molecular ice X and prove mathematically the stability of this polymorph above a critical pressure. Below it, we derive the pressure dependence of the interaction energy of the effective dipoles emerging in the system when the symmetric layout of the hydrogen atoms, specific for ice X, breaks down. This reproduces semiquantitatively the characteristic and unusual (as compared to the others—practically vertical) form of the boundary between the areas of the ordered and disordered ice VIII and VII. We also discuss the possibility of describing the differences between the ice phases existing at lower pressures (down to normal) by including the long-range electrostatic contributions: charge-charge and dipole-dipole in the crystal energy.

Appendix: Parameterization

The most ridiculous part of this story is the parameterization. It used to be a common place to complain about numerous parameters characteristic for semiempirical theories. Our days it is considered to be acceptable to upload a ca. 500 kB file of parameters not having whatever physical sense per atom within a PAW/DFT procedure and to call this ab initio.

Nevertheless, we find even the number of parameters required by the standard CNDO setting to be excessive for analytical treatment. Thus, we extended the procedure Ref. [39] and express required parameters through the Slater orbital exponents characteristic for atoms O and H: $\zeta _{\mathrm{H}}=1;\,\zeta _{\mathrm{O}}=2.275$ [13]. For the atomic parameters, the following exact expressions are used

$$\begin{aligned}&h=\frac{5}{8}\zeta _{\mathrm{H}};\,H=\frac{1}{2}\zeta _{\mathrm{H}};\,G=\frac{93}{256}\zeta _{\mathrm{O}};\\&U=3\zeta _{\mathrm{O}}-\frac{5}{12}\zeta _{\mathrm{O}}^{2};\\&W=U_{ss}-U_{pp}=\frac{1}{3}\zeta _{\mathrm{O}}^{2};\\&\beta =\frac{1}{2}\left( G+h\right) . \end{aligned}$$

which produce the “theoretical” values given in Table 2. Corresponding “experimental” values are extracted from different parameterization schemes previously developed for the CNDO approximation [50,51,52,53,54,55,56,57,58] to comply with experiment. In case of W, some scaling is obviously necessary (factor of ca. 0.3), but it is understandable due to the absence of the radial node in the 2s-Slater AO. The overestimate of U with respect to “ experimental” values is compensated by the respective overestimate of G and h which non-trivially enter only together.

Table 2 Comparison of “theoretical” and “experimental” CNDO parameters used (eV)

One-electron resonance integrals $t_{\sigma \sigma },t_{\zeta \sigma }$ are expressed as $t_{\mu \nu }=\beta S_{\mu \nu }$, where $S_{\mu \nu }$ is the corresponding overlap integral between the relevant Slater AOs.

For the two-center Coulomb repulsion integrals $\gamma _{\mathrm {AB}}$ the exact formulae [59] are employed which guarantees the “Yukawa-like” exponential decay of the core-core repulsion $Y_{\mathrm {AB}}$.

Finally, the value of G is scaled by the factor of 0.9151 to get the critical pressure exactly at 61 GPa. This, of course, does not affect the qualitative picture.