Communication: Dielectric properties of condensed systems composed of fragments

The polarizability of molecules and nano-structures is an important property determining the assembly and behavior of solids and liquids composed of well-defined building blocks (BBs). In addition, molecular polarizabilities play a key role in vibrational spectroscopy, e.g., in the calculation of Raman¹ and sum-frequency generation² spectra, in the determination of van der Waals interactions in solids and liquids,^3,4 and in the development of polarizable force fields.

While several methods are available to compute and predict polarizabilities of isolated molecules and nanostructures, their definition and calculation in condensed phases have been challenging and various levels of approximations have been adopted in the literature (e.g., Refs. 5–9). For example, the Clausius-Mossotti (CM) equation relates the average atomic or molecular polarizability α of a material building block to its electronic dielectric constant ϵ_∞,^5,6

where N is the number density of atoms or molecules. In Eq. (1), if we substitute ϵ_∞ with the refractive index of the material n using ϵ_∞ = n², the Lorentz-Lorentz equation is recovered. The validity of the CM relation in condensed phases, such as molecular liquids or assembly of nanostructured solids, depends on the system and general rules to establish its regime of applicability are not available. We note that often times, the variation of polarizabilities from the gas to the condensed phase are neglected. For example, many molecular dynamics simulations of aqueous solutions using force fields assume a fixed molecular polarizability of water,¹⁰ though first-principles electronic structure studies of water at ambient conditions have shown that molecular polarizabilities have rather broad distributions and are not isotropic.^7–9 Most electronic structure methods^7–9 consider only the dipole-dipole interaction when calculating the variation of polarizabilities upon assembly of molecular fluids or solids, and in many cases, such approximations have remained untested. Recently, substantial progress has been reported in computing polarizabilities of building blocks using maximally localized Wannier functions (MLWFs).¹¹ However, the method of Ref. 11 requires separate calculations of the dielectric properties of each constituent self-consistently, and it is not based on global induced fields within the condensed system.^12,13

In this paper, we propose a first-principles method to compute the polarizabilities of building blocks in condensed phases. The method, based solely on electronic structure calculations for the condensed phase, is applicable to semiconductors and insulators (i.e., to any systems with a gap between occupied and virtual states such that maximally localized Wannier functions can be defined¹¹). We present results for the molecular polarizabilities of water in a wide pressure-temperature (P-T) range, and we validate the CM relation for water at ambient conditions and the dipole-induced-dipole approximation (DID). We found that the DID becomes increasingly less accurate under pressure and breaks down when covalent bonds are present and oxygen ions are formed within the solid.

We start by summarizing our formulation. A building block (BB) composing of a condensed system (e.g., a molecule in a molecular crystal) is defined by its ionic coordinates and by electronic wave functions spatially localized at the BB site, for example, maximally localized Wannier functions $wi(r→)$¹¹ constructed from the Bloch orbitals of the condensed phase. The linearly induced electron polarization density of the BB, in response to a macroscopic field $E→$ is

where N_orb is the number of localized electronic orbitals (e.g., four doubly-degenerate orbitals for a water molecule with 8 valence electrons) and Δ$wi$ is the variation of the ith Wannier function under the field $E→$⁠. The local field (⁠$E→loc$⁠) acting on the BB is given by two contributions,

where $E→env$ denotes the field produced by the environment surrounding the BB, that is, by all the electrons that are not associated, through the definition of the $wi$⁠, to the BB. In most previous studies, $E→env$ was approximated by dipole-induced-dipole (DID) electrostatic interactions.^7–9 

The polarizability tensor α_BB of the BB is defined by the equation

where $μ→BB$ is the dipole moment computed as

In Eq. (5), q_e is the elementary charge and Δρ_BB is from Eq. (2). To compute α_BB, we need to calculate $E→loc$⁠. Since $E→$ is fixed, we only need to determine $E→env$⁠, which is simply

where $e′→(r→)$ is the microscopic electric field induced by all the electrons outside the BB.

The microscopic electric field $e′→$ can be evaluated within the random phase approximation (RPA) or by including the variation of the exchange and correlation potential [we denote the latter with density functional theory (DFT)]. Within RPA, $e′→RPA$ is obtained using Gauss’ law

where Δρ′ = Δρ − Δρ_BB and Δρ is the electron polarization density of the whole system. Historically, it is from the wing part of the inverse dielectric matrix.^27,28 At the DFT level, the exchange-correlation (xc) potential V_xc also contributes to the microscopic local field,

Once E_loc and μ are computed from Eqs. (3) and (5), respectively, α_BB is known. Therefore the procedure outlined here to obtain polarizabilities of BB within a condensed system is rather simple. Once the density and single particle wave functions are computed, e.g., by solving the Kohn-Sham equations, the electron polarization density is obtained by performing a single self-consistent calculation for the whole system. Equation (7) or (8) is then solved non-self-consistently, including multipole interactions at all orders. Solvers to obtain linear variations of the electron density exist in most DFT codes, using either density-functional perturbation theory (DFPT)¹⁴ or finite fields.¹⁵ 

We now turn to applying the method outlined above to the study of the molecular polarizabilities of water in a broad P-T range from ambient to supercritical conditions. The electron polarization density was obtained by DFPT, as implemented in the plane-wave pseudopotential code Qbox (http://qboxcode.org/).^9,16 We used Hamann-Schluter-Chiang-Vanderbilt norm-conserving pseudopotentials^29,30 with a plane-wave kinetic energy cutoff of 85 Ry. The MD trajectories of water at ambient conditions were taken from the water PBE400 dataset:³¹ http://www.quantum-simulation.org, obtained using 64 water molecules in the simulation box. The supercritical water trajectories were from our previous simulations,^18,19 with 128 water molecules in the supercell. The MD trajectories of the Na⁺-water solution are from Ref. 20. At least 60 snapshots from each MD trajectory were employed in our electronic structure calculations. We calculated the polarizability of each water molecule for each of the MD snapshots and obtained averages and standard deviations. For ice VIII and ice X, we used a 96 and 128-molecule supercells, respectively; the accuracy of the cell size was validated using a Monkhorst-Pack k-point mesh of 8 × 8 × 8 with the primitive cells.³² The Perdew-Burke-Ernzerhof (PBE) exchange-correlation (xc) functional¹⁷ was used. Although PBE overestimates the molecular polarizability of an isolated water molecule by ∼10%, for water under pressure PBE gives both the static and the electronic dielectric constants in better agreement with experimental values than at ambient conditions, as shown in previous studies.^18,19 Here we used one xc functional to analyze trends of polarizabilities as a function of P and T; however, the method is general and can be used with any functional. In particular, we note that using finite field methods to compute polarizabilities (http://qboxcode.org/),⁹ calculations with hybrid functionals are readily carried out.

Figure 1 shows that at ambient conditions, the molecular polarizabilities of water given by DFT [Eq. (8)] are anisotropic. The out of plane polarizability is the largest, and the ones in-plane and perpendicular to the water dipole direction are smaller, consistent with the reports of other authors using just DID interactions to compute E_env.^7–9 We found that at high pressures and high temperatures, the anisotropy substantially decreases as shown in Fig. 1. Note that the polarizability of an isolated water molecule is less anisotropic than in the liquid at ambient conditions,^9,13 suggesting that the anisotropy is critically related to the formation of hydrogen bonds. Indeed, also in supercritical water, the polarizability components are less dissimilar than at ambient conditions (see Fig. 1).

In Table I, four different methods to compute polarizabilities are compared from conditions varying from ambient to 11 GPa and from 0 to 2000 K. All methods show that with increasing pressure along an isotherm, the average molecular polarizability of water (⁠$α¯mol=13Tr{αmol}$⁠) decreases, while with increasing temperature along an isobar, it increases. Our previous study showed that the average dipole moment of water molecules increases with pressure, but decreases with temperature, so the present results indicate that varying the molecular dipole moments of water becomes more difficult when the values of the moments increase.

The polarizabilities obtained by DFT [Eq. (8)] are slightly larger than those from RPA [Eq. (7)] by ∼0.02 ų. When applying the two methods, we used the same electron polarization density Δρ, which is obtained when the exchange-correlation functional is included. The local electric field mainly comes from the electrostatic interactions, so the DFT and RPA values turn out to be similar.

In order to test the validity of the CM relation, we substituted the electronic dielectric constant ϵ_∞, obtained by DFPT into the CM relation to calculate the average molecular polarizability. It is interesting to see that the CM relation yields nearly identical results as the DID approximation. The standard deviations obtained for the CM relation are smaller than those from the DID approximation by one order of magnitude, as they only arise from the thermal fluctuation of ϵ_∞, not from molecular distributions as shown in Fig. 1. The CM relation holds when the Lorentz relation holds in an isotropic material, i.e., when the field at the center of a fictitious spherical cavity created by molecules inside the cavity vanishes.⁶ A well-known example for the Lorentz relation is the lattice with cubic symmetry, where only dipole-dipole interactions are considered.⁶ The agreement between the results obtained using the CM relation and the DID approximation suggests that the Lorentz relation is accurate when we consider only the dipole-dipole interaction for the water systems studied in Table I.

We now turn to comparing molecular polarizabilities of water in various phases. If we substitute the refractive index of 1.333, the experimental value for water at 293 K, and ambient pressure, into the CM relation, we get a molecular polarizability of 1.47 ų, which is the same as the experimental value for water vapor. At the PBE level of theory, the polarizability of an isolated water molecule is 1.60 ų, the same as the values obtained at ambient conditions using the CM relation and DID approximation (see Table I), consistent with previous studies.⁹ Hence, within the DID approximation, the average molecular polarizability of water does not change from gas to liquid phase. However, both the RPA and DFT methods give slightly larger values (∼3% larger than that obtained with CM and DID methods). We note that recently, Ge and Lu reported the molecular polarizabilities of water and ice at ambient conditions calculated using the local dielectric response of orbitals,¹³ where the electron polarization density of each molecule Δρ_mol is evaluated individually. In general, the sum of Δρ_mol does not equate the total Δρ. In Ref. 13, the $α¯mol$ of water was found to increase by ∼10% (instead of ∼3%) from gas to the liquid at ambient conditions.

Using the DFT, we also calculated the molecular polarizabilities of water in the first solvation shell of the Na⁺ ion at ambient conditions:²⁰ $α¯mol$ is 1.62 ų, which is again slightly larger than that obtained by the DID approximation by ∼3%. For water molecules with dangling bonds in the basal surface layer of ice Ih,^21,22 the difference in $α¯mol$ given by the DFT and DID approaches is even smaller: 1.62 ų vs 1.61 ų, only ∼1%. Our results suggest that at ambient conditions, the DID approximation works remarkably well.

Table I shows that the polarizabilities obtained by our method and the CM relation or the DID approximation differ with increasing pressure at a fixed temperature. For ice VIII, a high pressure phase consisting of two interpenetrating cubic ice sublattices,²³ when increasing pressure from 0 to 30 GPa, the difference between DFT values and the DID approximation increases from 6% to 13%, as shown in Fig. 2. These results indicate that interactions higher than dipole-dipole play a bigger role for denser water.

We close by considering the case of extremely dense water, ice X, the highest pressure phase ever determined experimentally.²⁴ In ice X, the oxygen atoms are in a body-centered cubic lattice and the hydrogen atoms sit right between the two nearest O atoms (see Fig. 3). Because the H atom is equidistant to two O atoms, it is no longer possible to define H₂O molecules; however, since the four maximally localized Wannier orbitals are still closely localized around O atoms, a new BB can be defined, and the molecular polarizability discussed below refers to the polarizability of the O²⁻ anion.

Figure 3 shows that the CM relation gives the same results as the DID approximation, whereas the molecular polarizabilities given by the DFT and RPA methods are about 20% larger. The reason is that in ice X covalent bonds are present, and indeed the BB identified by our calculation is no longer a water molecule, but rather an anion, for which higher-order interactions play an important role.

For ice X, another interesting finding of our calculation is that the electronic dielectric constant ϵ_∞ has a minimum at around 250 GPa, and accordingly the bandgap increases up to 150 GPa and then decreases slowly, as shown in Fig. 4. The inverse correlation between the electronic dielectric constant and the bandgap of ice X is consistent with the Penn model,^25,26 and differs from what we found in ice VII/VIII and water up to 30 GPa in our previous study.¹⁹ Generally, ϵ_∞ increases when both the molecular polarizability and material density become larger. With increasing pressure, the molecular polarizability of ice X decreases as shown in Fig. 3, whereas the material density increases due to volume shrinking, hence the molecular polarizability and the material density of ice X are two competing factors determining ϵ_∞; this is also the reason why the variation of ϵ_∞ is not very pronounced (see Fig. 4). From 50 to 250 GPa, the molecular polarizability dominates the change of ϵ_∞, but above 250 GPa, the rate of its decrease becomes slower and thus the material density becomes a more important factor. As a result, ϵ_∞ increases slowly as shown in Fig. 4.

In order to predict the properties of solids and liquids composed of well-defined building blocks, it is important to determine the variation of the dielectric properties of the isolated molecular or nano-scale constituents upon assembly. Hence the ability to compute dipole moments and polarizabilities of building blocks in condensed phases is critical. In this paper, we proposed a first-principles method to compute polarizabilities of sub-entities in condensed phases, which includes multipole interactions at all orders and is applicable to semiconductors and insulators. The method only requires a single self-consistent calculation for the entire condensed system, as opposed to multiple calculations for each building block, and it is readily applicable within and beyond the RPA.

As an example, we presented results for the molecular polarizabilities of liquid water in a wide pressure and temperature range. We found that at ambient conditions, the dipole-induced-dipole approximation is sufficiently accurate and the Clausius-Mossotti relation may be used, e.g., to obtain molecular polarizabilities from experimental refractive indexes. However with increasing pressure, this approximation becomes unreliable. For example, in ice VIII, the contribution of multipoles beyond the dipole is ∼13% at 30 GPa, and in ice X, the difference between all multipoles and the DID contribution is about ∼20% at 350 GPa, indicating that when hydrogen bonds are replaced by covalent bonds, higher-order interactions cannot be ignored. In the case of ice X, the CM relation is not valid, though the Lorentz relation still holds under the DID approximation. We also found that the bandgap of ice X has a maximum, while the electronic dielectric constant of ice X has a minimum, as a function of pressure. Finally we note that the knowledge of the polarizabilities of sub-entities under pressure may help to design polarizable force fields suitable for extreme P-T conditions. The method presented here can be used to study the local dielectric response of a wide range of semiconductors and insulators.

We thank Deyu Lu, He Ma, and Ikutaro Hamada for their helpful discussions. D.P. acknowledges support from Hong Kong Research Grants Council (Project No. ECS-26305017), the National Natural Science Foundation of China (Project No. 11774072), the Alfred P. Sloan Foundation through the Deep Carbon Observatory, and the Croucher Foundation through the Croucher Innovation Grant. M.G. and G.G. were supported by MICCoM, as part of the Computational Materials Sciences Program funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division. This research used resources of the Research Computing Center at the University of Chicago, the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract No. DE-AC02-06CH11357.