Deep neural network for the dielectric response of insulators

TL;DR

This paper presents a deep neural network model that accurately predicts the dielectric response of insulators by learning from ab-initio data, enabling simulations of complex environments and phase transitions.

Contribution

The authors develop a symmetry-preserving neural network integrated with the Deep Potential model to predict dielectric responses from electronic structure data.

Findings

Successfully predicts infrared spectra of liquid water and ice under pressure.

Captures dielectric response in mutating chemical environments.

Enables simulations of phase transitions beyond direct ab-initio methods.

Abstract

We introduce a deep neural network to model in a symmetry preserving way the environmental dependence of the centers of the electronic charge. The model learns from ab-initio density functional theory, wherein the electronic centers are uniquely assigned by the maximally localized Wannier functions. When combined with the Deep Potential model of the atomic potential energy surface, the scheme predicts the dielectric response of insulators for trajectories inaccessible to direct ab-initio simulation. The scheme is non-perturbative and can capture the response of a mutating chemical environment. We demonstrate the approach by calculating the infrared spectra of liquid water at standard conditions, and of ice under extreme pressure, when it transforms from a molecular to an ionic crystal.

Equations

$${\mathbf{w}}_i=\frac{1}{4}\sum _{l_i=1}^4{\mathbf{w}}_{l_i}$$

$${\mathbf{w}}_i=\mathbf{f}(\{{\mathbf{r}}_k,k\in {\mathcal{N}}_i\}).$$

$${\mathbf{r}}_{ki}^{\prime }\equiv {\mathbf{r}}_k-{\mathbf{r}}_i,{\mathbf{w}}_i^{\prime }\equiv {\mathbf{w}}_i-{\mathbf{r}}_i.$$

$$T_{i\lambda }=\sum _{k=1}^{N_i}{E}_i(s(r_{ki}^{\prime }))Q_{k\lambda }.$$

$$\mathbf{D}=\mathbf{T}{\mathbf{S}}^T.$$

$$w^{\prime }{}_i^{\lambda }=\sum _{j=1}^M{F}_j(\mathbf{D})T_{j,\lambda +1}.$$

$$\alpha (\omega )n(\omega )=\frac{2\pi \beta }{3cV}{\int }_{-\infty }^{+\infty }dt{e}^{-i\omega t}⟨\dot{\text{M}}(0)\cdot \dot{\text{M}}(t)⟩,$$

Full text

Machine learning (ML) schemes introduced in the last decade Behler and Parrinello (2007); Bartók et al. (2010); Rupp et al. (2012); Montavon et al. (2013); Botu et al. (2016); Chmiela et al. (2017); Schütt et al. (2017); Smith et al. (2017); Han et al. (2018); Zhang et al. (2018a, b) can model accurately the potential energy surface (PES) of a multi-atomic system upon training with first-principle electronic density functional theory (DFT) data Kohn and Sham (1965). These approaches extend the size and time range of ab initio molecular dynamics (AIMD) Car and Parrinello (1985); Marx and Hutter (2009), and make possible studies of rare events, such as crystal nucleation, with enhanced sampling methodologies in simulations of ab initio quality Bonati and Parrinello (2018).

Most approaches, so far, focused on representing the dependence of the PES, a scalar quantity, on the atomic coordinates. Recently, methods to fit the environmental dependence of electronic properties have been proposed Gastegger et al. (2017); Grisafi et al. (2018a, b); Wilkins et al. (2019); Chandrasekaran et al. (2019); Zepeda-Núñez et al. (2019); Raimbault et al. (2019); Kapil et al. (2019). In particular, kernel based methods have been used to represent the polarization and its time derivatives Grisafi et al. (2018a); Raimbault et al. (2019); Kapil et al. (2019), which are needed in many studies of materials, including calculations of infrared (IR) Sharma et al. (2005), Raman Putrino et al. (2000); Wan et al. (2013), and sum frequency generation (SFG) Wan and Galli (2015) spectra, transport calculations in ionic liquids Rozsa et al. (2018) and superionic crystals Sun et al. (2015); Wood and Marzari (2006); Schwegler et al. (2008), and simulations of ferroelectric phase transformations Srinivasan et al. (2003); Fluri et al. (2017). These schemes learn from AIMD trajectories, but, so far, the calculated spectra of liquid water Kapil et al. (2019) do not match the quality of direct many-body expansions of the dipole moment Liu et al. (2015).

Here we propose an alternative approach based on deep neural networks (DNNs) and maximally localized Wannier functions (MLWFs) Marzari and Vanderbilt (1997); Marzari et al. (2012), i.e. electronic orbitals with minimal spatial spread obtained from a unitary transformation of the occupied orbitals in insulators. In spin saturated systems the MLWFs describe electron pairs. Thus, upon assigning charges of 2$e^{-}$ to the Wannier centers (WCs), i.e. the centers of the MLWF distributions, the electric polarization is the dipole moment $\bm{M}$ of the neutral system of point charges made by the WCs and the atomic nuclei. In extended periodic systems, $\bm{M}$ includes an arbitrary quantum Resta (1994), but its derivatives and correlation functions are well defined and describe observable properties.

In systems with different atomic species, the WCs fluctuate near the most electronegative atoms during molecular evolution. As a consequence of the nearsightedness of the electronic matter Kohn (1996); Prodan and Kohn (2005), the WCs only depend on the atoms in their local environment, and their positions can be accurately represented by the DNN model, called Deep Wannier (DW), which is introduced here. It is an end-to-end scheme that does not use any ad hoc construction, in addition to the coordinate information and the network itself, to map input atomic coordinates into output WC coordinates. The model is size extensive, preserves translational, rotational, and permutational symmetry, and yields a polarization that varies continuously with the atomic coordinates. As such it can describe how the polarization responds to chemical bond changes in electronic insulators.

To demonstrate the methodology, we compute (a) the IR absorption spectrum of liquid water at standard temperature and pressure, and (b) the frequency dependent imaginary dielectric function in the high pressure transformation of ice VII to ice VII’, and from ice VII’ to ice X Hernandez and Caracas (2018). In ice VII there are 2 donor and 2 acceptor H-bonds per oxygen, as stipulated by the ice rule Bernal and Fowler (1933); Pauling (1935), and the water molecules are always well defined. In ice VII’ the hydrogens hop between two equivalent sites near each oxygen of an O-O bond, occasionally violating the ice rule. In ice X the hydrogens sit in the mid O-O bonds, the water molecules cannot be identified, and the system is better described in terms of O2- and H+ ions. This behavior is reflected in the O-H pair correlation functions in Fig. 1. For more sophisticated analyses and a phase diagram see e.g. Ref Hernandez and Caracas (2018). The polarization change accompanying the above transformations is seamlessly described by DW. The scheme should also work at higher temperature when ice X becomes superionic Millot et al. (2019), or at even higher temperature when the superionic crystal melts into an ionic liquid Rozsa et al. (2018). Similarly, the scheme should work for proton transfer events in water at standard conditions when neutral molecules interconvert with hydronium and hydroxide ion complexes (see e.g. Ref. Chen et al. (2018)).

In our implementation, DW is trained with valence only pseudopotential electronic structure calculations. In the laboratory frame, the positions of the atoms and of the WCs are $\bm{r}_{1},\bm{r}_{2},...,\bm{r}_{i},...,\bm{r}_{N}$ and $\bm{w}_{1},\bm{w}_{2},...,\bm{w}_{i},...,\bm{w}_{N_{w}}$, respectively. We assume, for simplicity, that the WCs are only associated to one atomic species and consider water, in which there are 4 WCs per oxygen, as a concrete example 111Generalizations to more than one reference atom would be feasible to deal with situations in which the WC are shared among covalently bonded atoms of same electronegativity.. It is easy to select with a cutoff distance the 4 WCs with coordinates $\bm{w}_{l_{i}=1,..,4}$ that are closer to the oxygen $i$ located at $\bm{r}_{i}$. Their centroid

[TABLE]

is well defined even when water molecules cannot be identified. Our aim is to construct a vector function $\bm{f}$ that gives $\bm{w}_{i}$ if the positions of the $N_{i}$ atoms $\bm{r}_{k}$ in the neighborhood $\mathcal{N}_{i}$ of $\bm{r}_{i}$, defined by the radius $r_{c}$, are known, i.e.:

[TABLE]

Here $i$ is an oxygen atom but the atoms $k$ include oxygens and hydrogens. To ensure that $\bm{f}$ preserves the translational, rotational, and permutational symmetry, we generalize the scheme for the Deep Potential (DP), the DNN representing the PES, introduced in Refs Zhang et al. (2018a, b).

First, we make a $local$ frame transformation to the primed coordinates, which preserve translational symmetry:

[TABLE]

Then, we introduce a weight function $s(r_{ki}^{\prime})$ equal to $1/r_{ki}^{\prime}$ at short distance $r^{\prime}_{ki}=(\bm{r}_{ki}^{\prime}\cdot\bm{r}_{ki}^{\prime})^{1/2}$, and decaying smoothly to zero as $r_{ki}^{\prime}$ approaches $r_{c}$. Using $s(r_{ki}^{\prime})$, we describe the atomic coordinates $\bm{r}_{ki}^{\prime}=(x_{ki}^{\prime},y_{ki}^{\prime},z_{ki}^{\prime})$ with the 4-vector \bm{q}_{ki}=(q_{ki}^{1},q_{ki}^{2},q_{ki}^{3},q_{ki}^{4})$$=(s(r_{ki}^{\prime}),s(r_{ki}^{\prime})x_{ki}^{\prime},s(r_{ki}^{\prime})y_{ki}^{\prime},s(r_{ki}^{\prime})z_{ki}^{\prime}) to enforce continuous evolution when atoms enter/exit the neighborhood.

Next, we enforce permutational invariance and rotational covariance by introducing two DNNs, an $embedding$ DNN and a $fitting$ DNN. The number of hidden layers and outputs is refined in the training procedure. The embedding DNN is the matrix $\bm{E}=(E_{ik})=(E_{i}(s(r_{ki}^{\prime})))$ with $M$ rows and $N_{i}$ columns optimized by training, which maps the set $\{s(r_{ki}^{\prime}),~{}k\in\mathcal{N}_{i}\}$ onto $M$ outputs.

The set of generalized coordinates $\{\bm{q}_{ki}\}$ in a neighborhood is represented by the matrix $\bm{Q}=(Q_{k\lambda})=(q_{ki}^{\lambda})$ with $N_{i}$ rows and 4 columns. Multiplication of $\bm{E}$ by $\bm{Q}$ gives the matrix $\bm{T}=\bm{E}\bm{Q}$ with $M$ rows and 4 columns, whose generic element is:

[TABLE]

Let $\bm{S}$ be the matrix formed by the first $M^{\prime}$ (<$M$) rows of $\bm{T}$. Multiplication of $\bm{T}$ by $\bm{S}^{T}$, the transpose of $\bm{S}$, gives the matrix $\bm{D}$ of dimension $M\times{M^{\prime}}$, called the feature matrix:

[TABLE]

$\bm{D}$ is the argument of the fitting DNN, a row matrix $\bm{F}(\bm{D})=\{F_{j}(\bm{D}),j=1,...,M\}$ that converts the atomic coordinate information encoded in $\bm{D}$ onto $M$ outputs, which are mapped onto the centroid $\bm{w}_{i}^{\prime}=(w^{\prime}{}_{i}^{1},w^{\prime}{}_{i}^{2},w^{\prime}{}_{i}^{3})$ upon multiplication with with the last three columns of $\bm{T}$:

[TABLE]

Finally, $\bm{w}_{i}$ is retrieved from $\bm{w}_{i}^{\prime}$ using Eq. 3 and one obtains the desired representation of Eq. 2.

We notice that $\bm{w}_{i}$ constructed in the way introduced above naturally preserves all the symmetry requirements. Translational symmetry is preserved by the adoption of a local frame and relative positions in Eq. 3. Permutational symmetry is preserved by the smooth sum over the neighboring atoms in Eq. 4. Finally, as shown in Eq. 4, the last three columns of $\bm{T}$ ($\lambda=2,3,4$) transform covariantly under rotation because $(Q_{k2},Q_{k3},Q_{k4})$ transforms like $\bm{r}_{ki}^{\prime}$. Then it is straightforward to verify that the elements of $\bm{D}$, and hence $\bm{F}(\bm{D})$, are invariant under rotation. Therefore, $w^{\prime}{}_{i}^{\lambda}$ in Eq. 6 transforms like $T_{j,\lambda+1}$, and is hence rotationally covariant. The values of $M$, $M^{\prime}$, and of the number of layers of the DNNs are chosen empirically based on performance. In the applications discussed in this paper we adopt $M$=100 (of the same order of the number of atoms in a neighborhood), $M^{\prime}$=6, and use a 3-layer representation for all the DNNs. The parameters $\bm{\gamma}$ of the embedding and fitting networks are determined by training, i.e., an optimization process that minimizes a loss function, here the mean square difference between the DW prediction and the training data. The Adam stochastic gradient descent method Kingma and Ba (2015) is adopted for the optimization. Generalization of this formalism to tensor properties (like the polarizability) is introduced in Ref. Sommers et al. (2020).

DW should be combined with a DNN for the PES to study the evolution of the polarization along MD trajectories. For consistency, the two networks should be trained with electronic structure data at the same level of theory, as in the applications below, which used DW and the DP representation of the PES Zhang et al. (2018a, b). Ab initio electronic structure data are expensive and efficient learning strategies are crucial. We used the iterative learning scheme of Ref. Zhang et al. (2019). In this approach, a DNN, initially trained with a limited pool of ab initio data, is used to explore inexpensively the configuration space. A small subset of the visited configurations is selected with a suitable error indicator and single shot ab initio calculations are performed at these configurations. Training with the new data improves the model for further exploration and selection, followed by new data acquisition and learning. The protocol is repeated until all the explored configurations are described with satisfactory accuracy. The error indicator exploits the highly non-linear dependence of the DNN models on the network parameters. As a consequence, different initializations of the parameters lead to different local minima in the landscape of the loss function, originating an ensemble of minimizing DNNs. The variance of the predictions within this ensemble is an intrinsic property of a DNN model and is often a reliable indicator of its accuracy Podryabinkin and Shapeev (2017); Zhang et al. (2018c). In our experience, good DNN models constructed with the above procedure require significantly less ab initio data in the target thermodynamic range than learning approaches based on independent AIMD sampling data.

In the following we report calculations on liquid water at STP and on ice undergoing pressure induced structural phase transitions. For liquid water at STP, we did not use the incremental data generation scheme, since the training data were available and accessible online from previous work Ko et al. (2019). For high-pressure ice, electronic structure data were not available, and we constructed DP and DW from scratch using the incremental learning procedure outlined above. Full details on the implementation, training, and validation of the models (DP and DW) are given in the Supplemental Material [SM]. The code for this work has been integrated into the open-source software package DeePMD-kit Wang et al. (2018) and we used the DP-GEN package Zhang et al. (2020) for the iterative scheme.

We use DFT at the hybrid functional level (PBE0 Adamo and Barone (1999)) with dispersion corrections Tkatchenko and Scheffler (2009) for STP water. Using DW and DP Zhang et al. (2018b) we calculate the IR absorption spectra of liquid H2O and D2O for a cell with 512 molecules under periodic boundary conditions. We use two microcanonical trajectories lasting 0.5 ns each, for H2O and D2O, at an average temperature of $\sim$300 K at the equilibrium density of the simulation. The frequency dependent absorption coefficient per unit length, $\alpha(\omega)$, times the refractive index, $n(\omega)$, is given by the Fourier transform of the time correlation function of the time derivative of the cell polarization $\dot{\textbf{{M}}}=+6e\sum_{l}\dot{\bm{r}}_{O_{l}}+e\sum_{m}\dot{\bm{r}}_{H_{m}}-8e\sum_{n}\dot{\bm{w}}_{n}$ according to:

[TABLE]

where $V$ is the volume, $\beta=1/k_{B}T$ is the inverse temperature, and $k_{B}$ is Boltzmann’s constant. Fig.  2 shows that the calculated spectra are in good agreement with the corresponding experimental observations.

Similarly accurate IR spectra of liquid water can be obtained from representations of the PES Babin et al. (2013); Shank et al. (2009) and the dipole moment Liu et al. (2015) based on many-body molecular expansions. These powerful approaches are limited to molecular liquids and crystals. By contrast, our non-perturbative method works also for non-molecular systems, as we demonstrate by considering ice at T=300 K in the pressure range from 20 to 110 GPa, wherein structural phase transitions from ice VII to ice VII’ and to ice X occur. We adopted the PBE functional approximation of DFT as in Refs. Schwegler et al. (2008); Hernandez and Caracas (2018). We constructed DP/DW networks for ice in the temperature range from 240 to 330 K and pressure range from 20 to 120 GPa with the iterative learning approach set forth earlier in this paper. This procedure required a total of 2248 single shot DFT calculations with a 16 molecule cell and a total of 2400 single shot DFT calculations with a 128 molecule cell. For each cell size, the corresponding computational effort was less than the cost of a short ($\sim$5 ps) AIMD trajectory.

We sampled the ice configurations at 300 K in the pressure range 20 - 110 GPa with constant pressure DPMD on a variable periodic cell with 128 water molecules, using a mild Nosé-Hoover thermostat Nosé (1984); Hoover (1985) with damping time of 5 ps, much longer than the vibrational periods, to control the temperature. Using DP+DW, we calculated $\alpha(\omega)n(\omega)$ according to Eq. 7 with a set of 0.5 ns long trajectories at various pressures. A direct AIMD study of the spectral changes in the closely related transformations from ice VIII (the proton ordered form of ice VII) to ice VII’ and to ice X, reported in a pioneering paper by Bernasconi et al. Bernasconi et al. (1998), used a small 16 molecule cell and $\sim$10 ps long trajectories.

Our results in Fig. 3 show a dramatic change with pressure of the product $\alpha(\omega)n(\omega)$. We also report the same quantity obtained from reflectance measurements of pressurized ice in a diamond anvil cell Goncharov et al. (1996). Experimental data are not available for $\omega\lesssim 800$ cm*-1*, and for 1800 cm${}^{-1}<\omega<$2400 cm*-1*. The theoretical curves are displayed at approximately 10 GPa higher pressure than the experimental curves to empirically correct for the missing nuclear quantum effects in the simulation and the inaccuracy of the adopted functional approximation.

Taking into account the limitations of theory and experiment, the experimental trend is reproduced nicely by the DP+DW simulations. At the lowest pressures, molecular vibration features can be discerned, such as the H stretching band at about 3000 cm*-1*. These modes soften dramatically and broaden with pressure, indicating a progressive weakening of the covalent O-H bonds. In the simulation, the spectral changes are quite abrupt at $\sim$70 GPa, suggesting that this should be approximately the transition pressure to symmetric ice X. The same behavior is observed in experiment at $\sim$60 GPa. By further increasing the pressure, two strong features emerge, characteristic of ice X, that harden and sharpen with pressure. The higher frequency feature has dominant weight on H while the lower frequency feature has more O weight. Interestingly, there is no close correspondence between the IR spectrum and the power spectrum of atomic dynamics reported in the [SM], suggesting that extreme anharmonicity affects the IR spectrum, as pointed out in Ref. Bernasconi et al. (1998).

Quantum effects in the dynamics were ignored in our calculations. These effects are small but non-negligible in liquid water and should be more pronounced near the ice VII to X transition in view of the relatively large isotope effect on the transition pressure, which is $\sim$10 GPa higher for D2O than for H2O Goncharov et al. (1996), suggesting that the transition is facilitated by proton tunneling. The calculation of dynamic quantum correlation functions is a major challenge for statistical simulations. Quantum IR correlations have been calculated recently for liquid water using approximate methods like ring polymer Marsalek and Markland (2017) and centroid Marsalek and Markland (2017); Medders and Paesani (2015) MD, indicating that quantum corrections tend to red shift the classical spectral features. It would be extremely interesting to study how quantum corrections affect the dielectric properties in the ice VII, VII’, and X transition sequence. Further studies of these issues will be facilitated by methods like DPMD and DW, which improve significantly the statistical quality of ab initio simulations, as they are orders of magnitude faster than DFT methods and scale linearly with system size. Quantitative cost comparisons between direct AIMD and DP/DW simulations are reported in Figs. S2 and S3 in the [SM].

In summary, DW is a useful tool to parametrize the dependence of the polarization on the atomic environment. The approach can be naturally extended Sommers et al. (2020) to the environmental dependence of the polarizability tensor $P_{\delta\sigma}=\frac{\partial M_{\delta}}{\partial E_{\sigma}}$ ($M_{\delta}$ and $E_{\sigma}$ are Cartesian components of the polarization and of the electric field), allowing us to compute Raman Wan et al. (2013); Sommers et al. (2020) and sum frequency generation spectra Mukamel (1998); Nagata and Mukamel (2010). Access to the concerted evolution of atomic coordinates and polarization in simulations of large systems over long time scales should also open the way to studies of ferroelectric phase transitions with MD simulations of ab initio quality, rather than relying on effective Hamiltonian models Resta and Vanderbilt (2007); Srinivasan et al. (2003); Vanderbilt and Zhong (1998). Finally, related symmetry preserving DNN schemes have been considered in Refs. Zaheer et al. (2017); Battaglia et al. (2018); Esteves et al. (2018); Han et al. (2019) and we defer to future work a discussion of the mathematical and machine learning aspects of the DP/DW models.