Introduction

KDP (KH2PO4) is the prototypical member of a class of isomorphous hydrogen-bonded molecular crystals that display unusually large isotope effects. Crystals in this class have the formula unit MH2XO4 (M = K, Rb, Cs, NH4; X = P, As)1. The XO4 groups are held together by strong hydrogen bonds (H-bonds) (Fig. 1a). The H-bonds are correlated with the polar distortion of the XO4 groups, giving rise to ferroelectricity when M = K, Rb, Cs, and anti-ferroelectricity when M = NH4. In these materials, Tc, the paraelectric to the ferroelectric (antiferroelectric) phase transition temperature, decreases with increasing pressure. Moreover, Tc can drop to zero at some critical pressure1. When hydrogen is substituted with its heavier isotopes, changes in Tc of the order of 100 K are observed1. At ambient pressure, KDP and DKDP (KD2PO4) have Tc equal to 123 K and 229 K, respectively. The isotope effects modulate piezoelectric, electro-optical, and nonlinear optical properties of KDP/DKDP2,3, leading to different applications such as electro-optic modulators, frequency converters, and piezoelectric devices4,5. A detailed quantum mechanical model is still missing to explain these isotope effects. Generally speaking, isotope effects are a manifestation of nuclear quantum effects (NQE) in atomic dynamics. Empirical observations indicate that NQE make strong H-bonds stronger6, an effect that should favor more symmetric bond configurations in KDP than in DKDP. However, the mechanism by which the strength of the H-bonds affects drastically the paraelectric–ferroelectric phase transition is not well understood.

Fig. 1: Structure of KDP.
figure 1

a Schematic representation71 of the Fdd2 equilibrium structure of KDP. The z-axis is parallel to the c crystallographic axis. b Sketch of a ferroelectric KDP unit. K (shared), H (shared), O, and P ions are represented by purple, light blue, red, and yellow spheres. The gray arrow shows the direction of the molecular dipole moment d. c Left: schematics of an O–H…O bridge. Colored shadows suggest the quantum delocalization of the particles. The geometric center of the oxygen pair is the origin of the coordinates. δxy is the projection on the xy plane of the position vector of the proton. Right: contour plot of the probability density function n(δxy), below and above Tc, inferred from neutron diffraction experiments12,70. Darker contour lines indicate higher probability density. Δxy(T), the distance of the maximum of n(δxy) from the origin, quantifies proton off-centering. d 16 types of KDP units are divided into four classes with decreasing energetic stability from left to right. Units in the same class have similar energetic stability. K ions and H-bonded H ions are not shown.

Full size image

Earlier theories attributed the emergence of spontaneous polarization to a uniform structural distortion combining irreversible off-centering of the H/D ions with polar distortion of the PO4 groups associated with a displacement of the P ions as in Fig. 1b. In the paraelectric phase, the H-bonds are uniformly symmetrized, and the PO4 groups are non-polar. These assumptions underlie the tunneling model7 and its variants8,9, which are defined by effective Hamiltonians that treat the protons as two-level systems and ignore or substantially reduce the degrees of freedom of the heavy ions (K, P, O). These models attribute the phase transition to the softening of a polar phonon mode combining collective proton tunneling with a polar distortion of the PO4 groups. The displacive mechanism was challenged by experiments9,10,11, which found that, in the paraelectric phase, the O–H…O bridges are asymmetric but disordered, so that the PO4 groups are polarized in random directions with distortions not substantially weaker than in the ferroelectric phase. This scenario attributes spontaneous polarization to the establishment of long-range order, consistent with neutron diffraction showing a lack of displacive-type structural distortions across the phase transition12. Indeed, the spatial distribution of the H/D ion relative to the neighboring O ions shows off-centered peaks in both ferroelectric and paraelectric phases, as illustrated in Fig. 1c.

It is crucial to understand precisely how NQE affects the disordering of the H-bond network. The insight gained on KDP should help in understanding the phase transition in the entire KDP family of crystals1. It may also elucidate the isotope shifts larger than 100 K observed in hydrogen-bonded H2O4C4 and some PbHPO4-type crystals13,14.

Key aspects of the order–disorder mechanism were elucidated by Slater15 and Takagi16,17, who pointed out that the local H-bond structure could be associated with the 16 phosphate configurations depicted in Fig. 1d, which groups the phosphates into four classes, i.e., Ferro, Slater, Takagi, and Nonpolar units, in order of decreasing energy stability from left to right. The more stable Ferro and Slater units satisfy the “ice rules”, because each phosphate has two donor and two acceptor H-bonds, in analogy with H2O ice. The less stable Takagi and Nonpolar units break the “ice rules” as they have unequal numbers of donor and acceptor H-bonds. Attributing nominal charges to the ions and taking K+ into account, Ferro and Slater units (\({{\rm{H}}}_{2}{{\rm{PO}}}_{4}^{-}\)) satisfy local charge neutrality, while Takagi and Nonpolar units violate it. All these units, with the exception of the Nonpolar units, have non-zero dipole moments. For T < Tc, broken symmetry Ferro units dominate, and the other units act as dipolar defects, destabilizing the ordered phase. For T > Tc, extensive excitation of dipolar defects occurs, and the symmetry of the two types of Ferro units is restored. Based on this phenomenology, early order-disorder models18 explained semi-quantitatively the phase transition by reducing the degrees of freedom of each phosphate unit to the 16 discrete states associated with the 16 unit types just described. Potentially, these models oversimplify the physics by ignoring the fluctuations that connect distorted configurations. At low temperatures, these fluctuations should be negligible in classical but not in quantum ferroelectrics.

We overcome the limitation of tunneling and discrete order-disorder models with all-atom path integral molecular dynamics (PIMD) simulations that sample all possible distortions of the phosphate groups according to a rigorously defined quantum partition function. Thermal and quantal fluctuations, including tunneling and anharmonicity19,20, are treated with controllable accuracy, revealing the microscopic mechanism that differentiates KDP from DKDP. We find that a substantial fraction of quantum dipolar defects is present in KDP at low temperatures. These intrinsic defects are almost eliminated by deuteration, leaving low-temperature DKDP in a nearly classical ferroelectric state. The quantum dipolar defects in KDP are mainly associated with Takagi units, despite their energy being higher than that of Slater units in classical order-disorder theory. Quantum dipolar defects explain the 25% smaller spontaneous polarization of KDP relative to DKDP measured in experiments at low temperatures. Quantum dipolar defects, as a source of zero-point entropy, also provide a natural explanation for the 30% smaller transition entropy of KDP relative to DKDP.

Results

PIMD models of KDP and DKDP

In PIMD N distinguishable quantum ions are mapped onto N classical ring polymers made of an equal number, L, of beads, subject to intra- and inter-polymer interactions. The intra-polymer interaction is harmonic and describes the quantum delocalization of a particle, while the bead inter-polymer interaction, given by the potential energy surface of the system rescaled by L−1, describes the physical interaction of the ions. For large enough L, the equilibrium values of position-dependent observables estimated by sampling the ring polymer distribution reproduce the corresponding quantum values (see “Methods”).

In our PIMD models of KDP and DKDP, the inter-polymer interaction does not rely on a phenomenological Hamiltonian or empirical force field. Instead, we derive the potential energy surface from quantum mechanics as in ab initio PIMD21, an approach in which the potential energy and polarization surfaces are computed on the fly within Kohn–Sham density functional theory (DFT)22. Due to the large computational cost of ab initio PIMD, only a small periodic system of four KDP units was studied in this way without finding a stable ferroelectric phase23. We overcome the limitations of computational cost with neural network representations of the potential energy24,25 and polarization surfaces26 (see Methods). The neural networks are trained on DFT data generated with the strongly constrained and appropriately normed (SCAN) functional approximation27, which, for KDP, predicts a good ferroelectric classical ground state28.

In the simulations, we adopt periodic boundary conditions. Then, the total dipole of the supercell, given by the product of the polarization \({\boldsymbol{{\mathcal{P}}}}\) and the volume V of the supercell, equals the sum of the molecular dipoles d associated with the phosphate groups sharing H and K ions with their neighbors, as depicted in Fig. 1b, modulo an immaterial constant (“the quantum of polarization”29) that does not affect the changes of \({\boldsymbol{{\mathcal{P}}}}\)29,30. Each dipole d is defined in terms of the positions and charges of the ions (nuclei plus core electrons) and of the centroids of the valence electrons (see “Methods”). The shared ions (K, H) carry half charge. The non-shared ions (P, O) and electronic centroids carry full charge. Similar definitions have been adopted in the literature31. For a periodic structure, various equivalent local dipole decompositions are allowed when they lead to uniquely defined physical observables such as a polarization change30. Given its connection with the macroscopic polarization, the average d is a proper local order parameter for the phase transition. We find that Nonpolar units have nearly vanishing d, Takagi and Slater units have non-zero d along directions approximately orthogonal to z, while Ferro units have non-zero d along z.

When modeling NQE on hydrogen bonds, the accuracy of the DFT approximation is crucial ref. 20 showed that common functional approximations predict double-well potential energy barriers of varying magnitude for the H-bonds in KDP. Since tunneling depends exponentially on the barrier height, a small error in the latter may significantly affect Tc. We find that SCAN underestimates the barrier height as it predicts a paraelectric KDP down to the lowest temperature of the simulations, while DKDP, which has reduced NQE, becomes ferroelectric at a temperature of ≈ 80 K, amounting to a downward shift of Tc by more than 100 K. Other semilocal functional approximations similarly destabilize the ferroelectric phase20, likely due to the self-interaction or delocalization error of DFT32,33. On the other hand, the nonlocal functional vdW-DF34,35 predicts a significantly larger tunneling barrier20. To check how it affects the ferroelectric transition we trained a DP model based on vdW-DF with the same protocol adopted for the SCAN-based models. PIMD simulations using this model yield Tc ≈ 215 K for KDP and Tc ≈ 285 K for DKDP (see Supplementary Note for details), amounting to an upward shift of Tc of ≈100 K and a reduced isotope effect. While SCAN errs by delocalizing the H/D ions too much, vdW-DF errs in the opposite direction, albeit to a slightly minor extent. This suggests that empirical fixing is unavoidable to accurately reproduce the experimental isotope effects with existing functional approximations. We pursue this goal within SCAN by correcting the excessive quantum delocalization of the H/D ions by stiffening the springs of the associated ring polymers, while leaving the other polymers unaffected. We do so by multiplying the spring constant by a stiffening factor μ chosen so that the off-centering in the xy plane of H and D, i.e., Δxy in Fig. 1c, match the respective neutron diffraction estimates12,36,37,38 (see “Methods”). This procedure leads to μ = 2.5 and μ = 16, respectively, for H and D ions. This empirical fix leaves the DFT potential energy surface unaffected but lowers the zero point energy of H/D to compensate for the underestimation of the double-well barrier in the H-bonds, enhancing the weight of the asymmetric O–H…O configurations. This procedure compromises momentum-dependent observables, like the momentum distribution measured in deep inelastic neutron scattering experiments39,40, but we find that it predicts static position-dependent observables with remarkable accuracy.

Using the above scheme we estimate the phase transition temperature Tc from the temperature evolution of the spontaneous polarization \(| {{\mathcal{P}}}_{z}|\), reported with the corresponding experimental data in Fig. 2a. The estimated Tc is ≈120 K for KDP, and ≈230 K for DKDP, with an error relative to experiments of less than 10 K in both cases. The predicted saturated polarizations, for T < Tc − 20 K, closely match their experimental counterparts. Notably, the simulations reproduce the experimental observation that \(| {{\mathcal{P}}}_{z}|\) in KDP is ≈25% smaller than in DKDP, further supporting the validity of the model. Due to the finite size of the simulations (see “Methods”), which include 512 KDP/DKDP units, the observed onset of spontaneous polarization is smoother than in experiments. In Supplementary Note, we also report the calculated temperature dependence of the electric susceptibility, which shows approximate Curie–Weiss behavior near Tc for both KDP and DKDP.

Fig. 2: Spontaneous polarization and lattice constants.
figure 2

a The longitudinal spontaneous polarization \(| {{\mathcal{P}}}_{z}|\) as a function of T. The experimental data are for pure KDP and KD1.96H0.04PO472. b Predicted and experimental12 lattice constants a and b as functions of T.

Full size image

Experimentally, the phase transition in KDP/DKDP is first-order41, structurally associated with a shear distortion of the lattice in the ferroelectric phase. The corresponding tetragonal to orthorhombic transformation is illustrated by the changes of the a and b lattice constants in Fig. 2b. The agreement of theory and experiment is excellent. The models reproduce the measured thermal expansion for T > Tc. For T < Tc, they slightly underestimate the orthorhombic distortion, while predicting a less sharp transition than in experiments due to finite-size effects. The shift of the a and b curves upon deuteration is well described. In Supplementary Note, we also analyze the isotope effects on the lattice constant c, the unit cell volume Ω, and the H-bond length dO…H. The simulation reproduces the Ubbelohde effect42,43, i.e., the elongation of the H-bonds following deuteration, a fine effect in the order of the hundredths of an angstrom.

Quantum dipolar defects

The remarkable isotope effects displayed in Fig. 2 are due to the presence of dipolar defects induced by quantum fluctuations. In PIMD, each configuration of the beads at equal imaginary time defines a distinct replica of the system (see “Methods”). Each KDP/DKDP unit in a replica can be assigned to one of the four classes displayed in Fig. 1d, depending on the donor/acceptor character of its four H-bonds. The character of each bond is linearly correlated, to excellent approximation, with the length and character of the OP bond and the position of the electronic centroid about the oxygen, as illustrated in Fig. 3a. In view of this correlation, the two types of P–O–H(D) structure in Fig. 3a can be distinguished by the length of rO→H. In the following, we choose the value r* = 1.25 Å for the divider between donor and acceptor type structures for both KDP and DKDP below and above Tc. This choice is supported by the calculated joint PDF of the angle α and of rO→H, i.e., n(α, rO→H), reported in Fig. 3b for DKDP and KDP at T = Tc − 10 K, and at T = Tc + 20 K in the insets. In all cases, n(α, rO→H) is strongly bimodal, with the peaks corresponding to typical donor and acceptor type P–O–H(D) structures, respectively. Moreover, n(α, rO→H) changes little across the phase transition. In DKDP, the divider r*, indicated by the gray dashed line in the figure, separates well donor and acceptor configurations. In KDP, quantum fluctuations make donor and acceptor distributions overlap, but symmetric H-bond configurations, in which H lies approximately in the middle of the O–O bond, have a substantially smaller statistical weight than asymmetric configurations, even at the lowest temperature of the simulations (T = 50 K), indicating that our choice of r* is meaningful.

Fig. 3: Structure and population of dipolar defects.
figure 3

a Schematic representation of the rules for the classification of each KDP/DKDP replica. Each P–O–H (yellow–red–blue) structure can be donor (left) or acceptor (right) type. The gray sphere indicates the ___location of the electronic centroid. b The joint probability distribution of α and rO→H for KDP and DKDP, respectively, in the ferroelectric (T = Tc − 10 K) and the paraelectric phase (insets, T = Tc + 20 K). The vertical gray dashed line shows r* = 1.25 Å. c Percent population of the different classes of KDP/DKDP units vs temperature. The vertical gray dashed line shows Tc. d RDF of the Takagi defects in KDP at 50 K. The distance r between the defects is the spatial separation of the respective P ions.

Full size image

The temperature variation of the population of phosphate units, classified using r*, is plotted in Fig. 3c. For T < Tc, broken symmetry Ferro units dominate both DKDP and KDP. Long-range order is lost in the paraelectric phase where the Ferro units, while still present, are not dominant, and the majority of the units are Slater and/or Takagi. In this phase, the frequent interconversions of Ferro, Slater, and Takagi units, restore the z → −z symmetry of the Ferro units. The population of Nonpolar units is negligible at all temperatures. The most important message of Fig. 3c is that Slater and Takagi populations are markedly different in DKDP and KDP. In DKDP, when T < 150 K, almost all units are Ferro. Dipolar defects, mostly associated with Slater units, appear only for T > 150 K, causing a gradual reduction of spontaneous polarization. By contrast, in KDP, a substantial fraction of Takagi, rather than the more stable Slater units, persists down to T = 50 K, whereby Ferro units only account for ≈75% of the phosphates. On the time scale of the PIMD simulations, the Takagi units are dynamic and fluctuate with H-bond switches every several picoseconds. By contrast, H-bond switches require tens of picoseconds at low temperatures in DKDP.

The Takagi defects are spatially correlated as shown by their radial distribution function (RDF) in KDP at T = 50 K in Fig. 3d. Due to averaging over the replicas in imaginary time, this RDF is already well converged in a single PIMD snapshot. The first peak corresponds to pairs of adjacent defects created mostly by switches of the H-bonds between neighboring Ferro units. The weak temperature dependence of the Takagi fraction in Fig. 3c suggests that the effect should originate from the quantum tunneling of the protons observed in deep inelastic neutron scattering39. Interestingly, similar experiments did not provide tunneling evidence in DKDP40. The prominent second peak of the RDF indicates that an additional defect is likely to be formed near a defect pair, while the non-monotonic decay of the RDF with the Euclidean distance actually reflects a monotonic decay with the graph distance defined by the minimal number of O–H…O bridges connecting pairs of Takagi defects. The defects are also dynamically correlated as suggested by analyses of their joint PDF reported in Supplementary Note. Recombination and migration to a neighboring site are two sources of dynamical correlation. The Takagi defects are nominally charged, and correlation acts to reduce local charge fluctuations.

Because of quantum disorder, the ferroelectric state of KDP is fundamentally different from that of DKDP. Linear extrapolation in Fig. 3c suggests that a substantial fraction (≈15%) of fluctuating Takagi defects should be present near absolute zero, conferring considerable zero-point entropy to KDP. Configurational disorder in the ferroelectric phase, which is significantly larger in KDP than in DKDP, is likely the main reason why the experimental transition entropy is ≈30% smaller in the former than in the latter44. The zero-point entropy here is an equilibrium property originating from quantum fluctuations. It is different from the residual entropy of ice, which is a non-equilibrium property due to quenched proton disorder45. Interestingly, the linear temperature dependence of the heat capacity has been observed near absolute zero in KDP46,47,48. The fluctuating Takagi defects might contribute to this effect.

Distribution of the local dipoles

To explain the roughly 25% smaller polarization of KDP relative to DKDP at low temperature, we study the probability density function (PDF) of the local dipole moments d = (dx, dy, dz) in KDP and DKDP.

The distributions n(dx, dz) ≡ n(dx, 0, dz) and n(dx, dy) ≡ n(dx, dy, 0), above and below Tc, are displayed in Fig. 4 for DKDP and KDP, respectively. In DKDP, breaking of the z-reflection symmetry of n(dx, dz) for T < Tc is evident in the upper (a) panel. This symmetry is restored for T > Tc (upper (b) panel). The magnitude of dz in the peak of the distribution in panel (a) closely matches the dipole of the Ferro units (Fig. 1d). The peak magnitude of dz does not change appreciably across the transition ((b) panel), as expected for a phase transition dominated by order-disorder effects. In the lower (a) panel, n(dx, dy) is uniformly close to zero, reflecting the z alignment of the dipoles in the ferroelectric phase. In the paraelectric phase (lower (b) panel), n(dx, dy) shows eight off-center spots with octagonal symmetry, associated with the Slater units of Fig. 1d. Each one of the four Slater units in that figure contributes to two off-center spots due to the glide reflection symmetry of the crystal. The central peak in the lower (b) panel originates from the diffusive part of the two dominant peaks of the upper (b) panel. The dipole distributions of KDP are reported in panels (c) and (d). Spontaneous symmetry breaking is again evident in n(dx, dz), with a dipole magnitude that does not change much across the transition, consistent with the order-disorder character. Relative to DKDP, the KDP distributions have a more strongly diffusive nature, signaling the presence of significantly more disorder above and below Tc. This is particularly evident in the paraelectric n(dx, dy) distribution, which, in KDP, shows roughly circular symmetry (lower (d) panel). This is due to the prominent role played by Takagi units generated by quantum fluctuations. While in DKDP Slater units dominate the paraelectric phase near Tc, in KDP Takagi and Slater units are equally important, blurring the octagonal pattern of the lower (b) panel.

Fig. 4: The joint PDFs n(dx, 0, dz) (upper) and n(dx, dy, 0) (lower) for KDP and DKDP.
figure 4

Panels a–d are associated respectively with T = 200 K, 250 K, 90 K, and 130 K.

Full size image

Comparing Fig. 4a and c, we find that the distribution of the local dipole moments associated with the Ferro units differs only marginally in KDP and DKDP. The value of dz at the peak of the distribution is equal to ≈1.7D in DKPD and to ≈1.6D in KDP, amounting to a dipole reduction of approximately 6% when deuterium is substituted by hydrogen. It suggests that the observed isotope effect on the polarization should be considered as strong experimental evidence for the existence of quantum dipolar defects in KDP at low temperatures.

Discussion

All-atom PIMD simulations of KDP and DKDP that agree well with experiments on the transition temperatures, the structural changes at the transition, and the magnitudes of the polarization, reveal the microscopic mechanisms that underlie the huge isotope effects. The ferroelectric phase is dominated by phosphate configurations with dipoles aligned in the polarization direction, but fluctuations introduce dynamic bond switching coupled with significant distortions of the phosphate groups, as shown in Fig. 3a. Ferroelectrically ordered phosphate groups have two donor and two acceptor H-bonds like in water ice, with zero nominal charge. A bond switch between two neighboring units creates a pair of phosphates with nominal charges of +1 and −1, which break the ice rules and act as dipolar defects having molecular dipoles approximately orthogonal to the polarization direction. The large effect of deuteration stems from the different roles played by these defects in KDP and DKDP. In the latter, they are only present near the transition temperature, while, in the former, they are present in sizeable amounts at all temperatures. Quantum disorder explains the 25% smaller spontaneous polarization of KDP relative to DKDP and also rationalizes the 30% smaller transition entropy of KDP relative to DKDP.

The huge effect of quantum fluctuations in KDP requires a fine balance between the tunneling barrier and the zero-point energy. This balance, which is broken by deuteration, is realized in the KDP family of crystals due to the presence of unusually strong H-bonds. H-bonds can be made stronger by pressure. In water ice, where the H-bonds at ambient conditions are relatively strong but significantly weaker than in KDP, isotope effects are relatively small, but they become large at pressures in the hundreds of kbar range49. In KDP, experiments show that the temperature of the ferroelectric transition is reduced with increasing pressure until, at a pressure in excess of 17 kbar, the phase transition disappears and the system becomes a quantum paraelectric50,51. We understand this phenomenology as the result of an increasing Takagi population originating from increasing H-bond symmetrization with pressure. We speculate that the phase point at which the phase transition disappears (T = 0 K, P = 17 kbar) may be associated with quantum criticality52,53, a hypothesis consistent with the large dielectric susceptibility measured in the vicinity of this point51. In future studies, it would be of interest to look for possible non-classical temperature dependence of the inverse dielectric susceptibility at very low temperature53.

Methods

Density functional theory and Deep Potential model

We adopt SCAN, the strongly constrained and appropriately normed functional approximation for exchange and correlation27. SCAN predicts a classical ferroelectric ground state whose structural properties are in better agreement with experiments than in the classical ferroelectric ground state predicted by other popular functional approximations28. We perform electronic structure calculations with Quantum Espresso54, using norm-conserving pseudopotentials55 that include in the valence the 2s and 2p states of O, the 3s and 3p states of P, and the 3s, 3p, and 4s states of K. We find a stable classical ferroelectric structure of space group Fdd2 at T = 0 K (Fig. 1a), with lattice constants a = 10.62 Å, b = 10.67 Å, and c = 6.86 Å. The length of the O-H covalent bond is 1.042 Å, that of the O…H H-bond is 1.438 Å, and the H off-centering is equal to \({\Delta }_{xy}^{\,\text{Fdd2}\,}=0.198\) Å. The above results agree well with ref. 28 with minor differences due to the different numerical methods adopted in the two calculations.

Then, we use the active learning agent DPGEN56 to collect SCAN-DFT data of KDP in the temperature range [50, 400] K and the pressure range [−20, 100] kbar. The parameters for the active learning procedure are publicly available at ref. 57. The active-learning exploration of atomic configurations is performed for 23 iterations. The final training dataset contains 5800 atomic configurations. We use DeePMD-Kit25 to train the neural network model, called the Deep Potential (DP) model24, for representing the potential energy surface. The performance of the DP model on the training set is shown in Fig. 5. On the training set, the root mean squared error of the model relative to DFT is around 0.33 meV/atom in energy prediction and around 0.07 eV/Å in Cartesian components of the force prediction. Using the DP model, we perform PIMD simulations of KDP with 16 beads, without stiffening spring constants, at respectively T = 100 K and T = 150 K. From the simulations, we sample 256 atomic configurations as the testing dataset, where SCAN-DFT calculations are performed to generate the data labels. On the testing set, the root mean squared error of the model relative to DFT is around 0.45 meV/atom in energy prediction and around 0.08 eV/Å in Cartesian components of the force prediction.

Fig. 5: Comparison of energies (zero of energy is shifted by an immaterial constant) and forces predicted by the DP model with the reference SCAN-DFT results in the training dataset.
figure 5

The insets show the error distribution.

Full size image

The error levels in the training and testing datasets are consistent, showing the DP model is an accurate representation of SCAN-DFT for PIMD simulations.

Deep Wannier model and molecular dipoles

In addition, we train a Deep Wannier (DW) model26 for the polarization surface in terms of the Wannier decomposition of the valence electron density29. Maximally localized Wannier functions are computed with Wannier9058 for a subset of the atomic configurations in the DP training set. Each Wannier function accommodates two spin-degenerate electrons. In all the atomic configurations four Wannier centers (geometric centers of Wannier probability distributions) are uniquely associated to each O and four others to each K atom, forming approximate atom-centered tetrahedral structures. The Wannier centroids (WCs), i.e., the geometric centers of the four Wannier centers associated with O or K atoms provide all the necessary information to compute the polarization26. The DW model gives the environmental dependence of the WCs. While the oxygen WCs depend crucially on the chemical environment, as illustrated in Fig. 3a, the potassium WCs, associated with the 3s and 3p semi-core states, are essentially independent of the environment and their contribution to polarization changes is negligible. We do not include them in the DW model, i.e., to calculate the polarization we treat the semi-core states as frozen core electrons of the K+ ion. The details of training the DW model are as follows.

From the atomic configurations generated through the first 10 iterations of active-learning exploration, we randomly collect 100 atomic configurations from each iteration. The 1000 collected configurations form a reduced dataset for the DW model. The reduction of the dataset is due to the high computational cost of generating training labels for the DW model, and the fact that the full dataset generated by active learning is typically redundant. The WCs in these atomic configurations are calculated with Wannier9058. The training labels for the DW model are the displacement vectors from the oxygens to the associated WCs. The training of the DW model is done with DeePMD-Kit. As shown in Fig. 6, the standard error of the DW model is 0.002 Å for the Cartesian components of the displacement vectors, much smaller than the typical distance from the oxygen of the associated WC.

Fig. 6: Comparison of the displacement vectors of the Wannier centroids from their associated O atom predicted using the DW model with the reference SCAN-DFT results in the training dataset.
figure 6

The insets show the error distribution.

Full size image

Following the theory of polarization29,30, we assign to each elementary KDP/DKDP unit (Fig. 1b) a molecular dipole d given by:

$$\begin{array}{rcl}{\bf{d}}&=&\mathop{\sum }\limits_{i=1}^{2}{q}_{{\rm{K}}}\left({{\bf{r}}}_{{\rm{K}}}^{(i)}-{{\bf{r}}}_{{\rm{P}}}\right)+\mathop{\sum }\limits_{i=1}^{4}{q}_{{\rm{O}}}\left({{\bf{r}}}_{{\rm{O}}}^{(i)}-{{\bf{r}}}_{{\rm{P}}}\right)\\ &&+\mathop{\sum }\limits_{i=1}^{4}{q}_{{\rm{H}}}\left({{\bf{r}}}_{{\rm{H}}}^{(i)}-{{\bf{r}}}_{{\rm{P}}}\right)+\mathop{\sum }\limits_{i=1}^{4}{q}_{{\rm{WC}}}\left({{\bf{r}}}_{{\rm{WC}}}^{(i)}-{{\bf{r}}}_{{\rm{P}}}\right).\end{array}$$
(1)

Here, \({{\bf{r}}}_{{\rm{K}}}^{(i)},{{\bf{r}}}_{{\rm{O}}}^{(i)},{{\bf{r}}}_{{\rm{H}}}^{(i)}\), and rP are the position vectors of the K, O, H/D, and P ions belonging to the elementary unit, and qK = 0.5e, qO = 6e, and qH = 0.5e are the charges carried by the K, O, and H/D ions, respectively. Since the K and H/D ions are shared with a neighboring unit, they carry half charge. The P ion, carrying a charge of 5e, is taken as the local reference and does not contribute to the dipole. Finally, \({{\bf{r}}}_{{\rm{WC}}}^{(i)}\) are the position vectors of the WCs associated with the O ions and qWC = −8e are their charges. When all the ionic and electronic charges are taken into account, each elementary unit is electrically neutral. The positions of the ions and of the WCs corresponding to each atomic configuration R = {r1, ⋯, rN} visited by PIMD in a periodic simulation box of volume V and N physical particles are provided by the DP and DW models. The global dipole of the configuration R is D(R) = ∑jdj, where the sum extends to all the local molecular dipoles dj contained in the simulation box. The polarization surface \({\boldsymbol{{\mathcal{P}}}}({\bf{R}})\), i.e., the polarization of the configuration R, is given, modulo an immaterial constant, by \({\boldsymbol{{\mathcal{P}}}}({\bf{R}})={\bf{D}}({\bf{R}})/V\)29.

Path integral molecular dynamics

Within Feynman’s path integral formulation of quantum statistical mechanics59, the canonical partition function of N distinguishable quantum particles can be approximated with a Trotter factorization in imaginary time:

$${Z}_{L}\approx \int{\mathcal{D}}{\bf{R}}{e}^{-\beta {{\mathcal{H}}}_{L}^{{\rm{eff}}}},$$
(2)

where

$${\mathcal{D}}{\bf{R}}\equiv \mathop{\prod }\limits_{i=1}^{N}{\left(\frac{{m}_{i}L}{2\pi \beta {\hslash }^{2}}\right)}^{\frac{3L}{2}}\,\text{d}\,{{\bf{r}}}_{i}^{(1)}\cdots \,\text{d}\,{{\bf{r}}}_{i}^{(L)}.$$
(3)

The effective Hamiltonian \({{\mathcal{H}}}_{L}^{\,\text{eff}\,}\) is

$${{\mathcal{H}}}_{L}^{\,\text{eff}\,}=\mathop{\sum }\limits_{s=1}^{L}\left(\mathop{\sum }\limits_{i=1}^{N}\frac{{k}_{i}}{2}{\left({{\bf{r}}}_{i}^{(s)}-{{\bf{r}}}_{i}^{(s+1)}\right)}^{2}+\frac{1}{L}U({{\bf{R}}}^{(s)})\right),$$
(4)

where \({k}_{i}=\frac{{m}_{i}L}{{\beta }^{2}{\hslash }^{2}}\) and \({{\bf{R}}}^{(s)}=\{{{\bf{r}}}_{1}^{(s)},\cdots \,,{{\bf{r}}}_{N}^{(s)}\}\). The index s labels imaginary times, the positive integer L is the total number of imaginary time slices, and the condition L + 1 = 1 applies in the argument of the sum over the beads so that \({{\bf{r}}}_{i}^{(L+1)}={{\bf{r}}}_{i}^{(1)}\). Thus, a discretized Feynman path in imaginary time, \(\{{{\bf{r}}}_{i}^{(1)},\cdots \,,{{\bf{r}}}_{i}^{(L)}\}\), is equivalent to a ring polymer. U is the potential energy surface provided by the DP model. ZL recovers the exact quantum partition function in the limit L → ∞. For finite L, Eq. (2) can be seen as the configurational partition function of a classical system of N ring polymers each made of L beads and described by the effective Hamiltonian \({{\mathcal{H}}}_{L}^{\,\text{eff}\,}\), in which the constant ki defines the spring that links neighboring beads in imaginary time, and \(\frac{1}{L}U({{\bf{R}}}^{(s)})\) is the potential of interaction between different ring polymers at equal imaginary time. The ensemble of all the beads at equal imaginary time constitutes a replica of the system. PIMD60 is a molecular dynamics scheme for sampling the configurations of the ring polymers with Boltzmann weights proportional to \({e}^{-\beta {{\mathcal{H}}}_{L}^{{\rm{eff}}}}\). Then, a canonical average of a position-dependent observable of the quantum system is estimated from a classical canonical average of the ring-polymer system. For example, the average macroscopic polarization is calculated as follows:

$$\begin{array}{rcl}{\boldsymbol{{\mathcal{P}}}}&=&\frac{1}{L}{\left\langle \mathop{\sum }\limits_{s = 1}^{L}{\boldsymbol{{\mathcal{P}}}}({{\bf{R}}}^{(s)})\right\rangle }_{{Z}_{L}}\\ &=&\frac{\int{\mathcal{D}}{\bf{R}}\mathop{\sum }\nolimits_{s = 1}^{L}{\boldsymbol{{\mathcal{P}}}}({{\bf{R}}}^{(s)}){e}^{-\beta {{\mathcal{H}}}_{L}^{{\rm{eff}}}}}{L\int{\mathcal{D}}{\bf{R}}{e}^{-\beta {{\mathcal{H}}}_{L}^{{\rm{eff}}}}},\end{array}$$
(5)

where \({\boldsymbol{{\mathcal{P}}}}({{\bf{R}}}^{(s)})\), the polarization of the s-th replica, is provided by the DW model.

PIMD simulations

We use I-PI61, LAMMPS62, and DeePMD-Kit to perform NPT-PIMD simulations on a supercell with 512 KDP/DKDP units in the temperature range [50, 300] K. The temperature is controlled with the PILE thermostat63 and the pressure is maintained with an anisotropic barostat64. The spring constants ki (Eq. (4)) associated with K, P, and O ions take their physical values, but the ki associated with H is multiplied by a stiffening factor μ = 2.5, while the one corresponding to D is multiplied by μ = 16. These values are chosen to fit the experimental neutron data on the off-center displacement of H/D, as explained in the next section.

To control convergence with L, we use the estimator KL of the thermodynamic quantum kinetic energy per particle suggested in ref. 65:

$${K}_{L}=\frac{3}{2\beta }+\frac{1}{2LN}{\left\langle \mathop{\sum }\limits_{s = 1}^{L}\mathop{\sum }\limits_{i = 1}^{N}\left({{\bf{r}}}_{i}^{(s)}-\overline{{{\bf{r}}}_{i}}\right)\cdot \frac{\partial U}{\partial {{\bf{r}}}_{i}^{(s)}}\right\rangle }_{{Z}_{L}}.$$
(6)

Here, \(\overline{{{\bf{r}}}_{i}}\) is the centroid of all the beads of particle i. In the classical limit, L = 1 and KL = 3kBT/2, the classical kinetic energy per particle. By monitoring the convergence of KL with L, we estimate ∣KL − K∞∣ at different temperatures as follows.

We perform PIMD simulations of systems with 512 KDP/DKDP units at different temperatures (T = 50, 100, 150, 200, 250, 300 K) with the adopted stiffening factor for H and D ions, respectively. The calculated KL as a function of L is shown in the upper panel of Fig. 7. KL converges more slowly at lower temperatures. About 128 beads are necessary at T = 50 K, the lowest temperature probed in our simulations. To estimate the error of a finite L calculation, we estimate \({K}_{\infty }=\mathop{\lim }\nolimits_{L\to \infty }{K}_{L}\), by fitting the KL data with the function KL = K∞ − h/L2 suggested by the L dependence of the Trotter factorization error for a path of length L. We employ a non-linear least square curve fitting using the Python package SciPy66, treating K∞ and h as parameters. At T = 50 K, KL data with L ≥ 32 are used for fitting. At all the other temperatures, KL data with L ≥ 16 are used. The fitted curves for KL − K∞ are displayed in the lower panel of Fig. 7. We consider a calculation converged with L when the error is less then 1 meV/atom, which is close to the typical error of our DP model. With this criterion we find that in KDP at T ≥ 150 K one needs L = 32. When 150K > T ≥ 100 K, L is set to 64, while for 100 K > T ≥ 50 K, L is set to 128. In DKDP, when T > 200 K we set L ≥ 16, and we set L to L = 64 when T ≤ 200 K.

Fig. 7
figure 7

(Upper) The quantum kinetic energy KL as a function of the number of beads. (Lower) KL − K∞ as a function of the number of beads. The diamonds are the numerical results from PIMD. The solid lines are the curves from optimal fitting. The two horizontal dotted guiding lines are at KL − K∞ = 0 and KL − K∞ = −1 meV, respectively.

Full size image

Classical finite size effects

A study of the finite size effect is impractical with PIMD simulations as L can be as large as 128. For that reason, We only study finite size effects on the ferroelectric phase transition of KDP with classical MD, i.e., by setting L = 1 in PIMD.

Calling \({\mathcal{N}}\) the number of KDP units in the simulation supercell, we compute 500 ps long classical NPT trajectories for \({\mathcal{N}}=64\), 512, and 1728 at different temperatures and standard pressure conditions. The finite size effect on the phase transition is illustrated by the temperature dependence of the polarization in Fig. 8a and by the temperature dependence of the average nearest neighbor O-O distance in Fig. 8b. The figure suggests that the size effect is small when \({\mathcal{N}}\ge 512\). A minor sharpening of the transition onset can be seen as \({\mathcal{N}}\) increases from 512 to 1728, suggesting that significantly larger sizes would be needed to reproduce more closely the sharp onset observed in experiment. From the transition onset in Fig. 8 we estimate that Tc should be in the range [245, 255] K. The estimated change of lOO upon the phase transition is approximately 0.005 Å, close to experimental findings12. We expect that similar size effects should affect the PIMD simulations reported in this work.

Fig. 8: Size dependence of the spontaneous polarization and O-O distance.
figure 8

a \(| {{\mathcal{P}}}_{z}|\) as a function of temperature from classical MD. b Average lOO as a function of temperature from classical MD.

Full size image

The spring stiffening factor μ for H and D ions

In KDP, SCAN-DFT predicts a ferroelectrically ordered classical ground state67, but this structure is destabilized by NQE. The absence of ferroelectricity was also found in ab initio PIMD simulations of KDP with PBE-DFT23. This occurs because common semilocal functionals tend to underestimate the double-well potential energy barrier in the H-bonds, wiping out the classical ferroelectric ground state20.

In the presence of NQE, a correct description of the ferroelectric phase requires a delicate balance of tunneling barrier and zero-point energy of the protons, which are both strongly correlated with the geometric distortion of the phosphate groups, as shown in numerous experimental12,13,36,37,38 and theoretical19,68,69 studies. Experiments find that, in KDP and DKDP, Tc at different pressure conditions is to very good approximation linearly correlated with \({\tilde{\Delta }}_{xy}^{{\rm{EXP}}}\), the off-center displacement of the H/D ions in the paraelectric phase just above the phase transition12,36,37,38. This correlation, shown in Fig. 9, strongly supports the notion that H/D off-centering is the key geometric parameter controlling Tc, and hints at a possible way of improving the experimental agreement of DFT-based PIMD simulations by fixing the predicted off-centering. In KDP, SCAN-based PIMD always yields unimodal proton distributions with Δxy(T) = 0, but the excessive quantum delocalization of H/D, responsible for the effect, can be reduced as desired by multiplying the H/D spring constant \(k=\frac{mL}{{\beta }^{2}{\hslash }^{2}}\) by a stiffening factor μ greater than 1, a procedure equivalent to assigning an effective mass m* = μm to H/D. The displacements \({\Delta }_{xy}(T\to {T}_{c}^{+})\) and the corresponding Tc calculated with PIMD for three different values of m* at ambient pressure (AP) are reported in Fig. 9, showing a linear behavior of Δxy with Tc that closely mimics the experimental one. For m* = 2.5 Da, i.e., μ = 2.5, and for m* = 32 Da, i.e., \(\mu =16,{\Delta }_{xy}(T\to {T}_{c}^{+})\) and Tc match the experimental displacements and transition temperatures of KDP and DKDP, respectively. This is the choice of μ that we adopt in our PIMD simulations of KDP and DKDP. In addition, we report the lattice constants and the H-bond length of KDP/DKDP simulated with μ = 1, KDP simulated with μ = 2.5, and the experimental KDP, in Supplementary Fig. 1.

Fig. 9
figure 9

(Upper) The relation between Tc and H/D off-centering. The black stars mark the experimental measurements of Tc and \({\tilde{\Delta }}_{xy}^{{\rm{EXP}}}\) for KDP and DKDP at different external pressures12,36,37,38. Blue spheres mark the numerical Tc and \({\Delta }_{xy}(T\to {T}_{c}^{+})\) obtained with PIMD with different m* at atmospheric pressure. The error bar of Tc is ±10 K, reflecting the T spacing adopted in the calculations of the order parameter near Tc. The error bar of \({\Delta }_{xy}(T\to {T}_{c}^{+})\) is ±0.005 Å, reflecting the space resolution of the histogram representing n(δxy) from which we extract \({\Delta }_{xy}(T\to {T}_{c}^{+})\). Blue shadows are a guide for the eye, suggesting an almost linear Tc-\({\Delta }_{xy}(T\to {T}_{c}^{+})\) relation. (Lower) The heatmap of n(δxy) computed at T > Tc. The left figure is for m* = 2.5 Da and T = 150 K. The right figure is for m* = 32 Da and T = 250 K.

Full size image

The displacement \({\Delta }_{xy}(T\to {T}_{c}^{+})\) is estimated from thermal averages at T ≈ Tc + 20 K to reduce the effect of the fluctuations that are large near Tc. The probability distributions n(δxy) for H and D are also reported in Fig. 9, showing a bimodal character in qualitative agreement with neutron diffraction experiments70. Relative to KDP, DKDP has a more localized n(δxy) and a larger Δxy, consistent with the more classical behavior of the deuterated H-bonds. In fact, in the classical limit, when all the atoms in PIMD have infinite mass, Tc is approximately equal to 250 K, slightly above the predicted transition temperature of DKDP, and consistent with a predominantly classical character of the transition in DKDP.

An argument based on a one-dimensional double-well potential model can explain qualitatively the significantly larger μ needed for D (μD = 16) than for H (μH = 2.5). When T = Tc, KQ, the excess kinetic energy of the tunneling mode, should be of the order of the barrier height (γ) minus the classical thermal energy of the mode (\(\frac{1}{2}{k}_{B}{T}_{c}\)). To good approximation, KQ should be equal to \(\frac{1}{2}K{\langle {x}^{2}\rangle }_{{T}_{c}}\) in terms of K, the spring constant of the potential, and of the spread at T = Tc. Approximating the latter with the spread at T = 0, leads to simple scaling of KQ with \(\frac{1}{\sqrt{{m}^{* }}}\). Then, using the SCAN-DFT value of γ (≈12 meV) and the values of 5meV and 10meV, respectively, for the classical thermal energies of the H and D modes, we get \({R}_{\mu }=\frac{{\mu }_{D}}{{\mu }_{H}}\approx 6.1\), almost equal to the value of 6.4 adopted in the manuscript. However, this close agreement is, to large extent, fortuitous because the scaling law of \({\langle {x}^{2}\rangle }_{{T}_{c}}\) is more complex than that of \({\langle {x}^{2}\rangle }_{0}\), and we have ignored the anharmonic coupling of the tunneling mode with the environment, which is essential to explain the physics of KDP/DKDP.