Shell-model and first-principles calculations of vibrational, structural and ferroelectric properties of KH$_2$PO$_4$

We develop a shell model (SM) for potassium dihydrogen phosphate (KDP) which is fitted to ab initio (AI) results that include nonlocal van der Waals corrections. The SM is comprehensively tested by comparing results of structural, vibrational and ferroelectric properties with AI and experimental data. The relaxed structural parameters are in very good agreement with the AI results and the available experimental data. The $\Gamma$-point phonons and the total phonon densities of states (DOSs) in the ferroelectric and paraelectric phases calculated with the developed SM are in good overall agreement with the corresponding AI and experimental data. We also compute the effective Debye temperature as a function of $T$ which shows good accordance with the corresponding AI and experimental results. Classical molecular dynamics (MD) simulations obtained with the developed SM show a FE-PE phase transition at $\approx 360$ K in remarkable agreement with ab initio MD calculations.


Introduction
The prototype of the H-bonded ferroelectric (FE) compounds, KH 2 PO 4 or KDP, was extensively studied in the past due to its important technological applications as well as for fundamental interest in its phenomenology [1]. In subsequent years after its discovery, KDP found extensive applications as electrooptical devices, as well as filter modulators in the wide field of laser spectroscopy. Besides the research in technological aspects, many efforts were also directed in the past to unveil the puzzle in the origin of the huge isotope effect that KDP manifests: the critical temperature for the paraelectric-ferroelectric (PE-FE) phase transition nearly doubles when the system is deuterated [2,3]. In spite of these efforts, the controversy about whether tunneling [4,5] or geometrical effects [2,3] are at the root of the giant isotope effect still remains.
This phenomenology is widely observed in the whole family of H-bonded ferroelectric compounds [6][7][8][9][10][11][12][13] including organic ferroelectrics that were recently discovered which attract much attention because of their potential for ecofriendly technological applications [14,15]. Ab initio (AI) calculations have shown that tunneling and geometrical effects are involved in a self-consistent phenomenon that greatly amplifies the effect of their correlation leading to the huge isotope effect observed [7,16].
With the aim of studying the phase transition and the isotope effects in KDP, many models were developed in the past [4,6,[17][18][19][20], e.g., the tunneling model and subsequent modifications [4,6,17,18]. The coupled proton-phonon model, which accounts for the interaction between the proton collective mode exhibiting tunneling and the optical lattice phonon of B 2 symmetry polarized along the -axis, was invoked to explain the emergence of polarization in the -direction [17,18]. The lower-frequency mode arising from this interaction is the FE soft mode that softens with temperature when the critical temperature is approached. Infrared measurements have also shown strong correlations between these used to study different ferroelectric, structural and vibrational properties. The SM calculations included full structural relaxations, zone-center phonons, phonon densities of states, the effective Debye temperature and a "classical-nuclei" SM molecular dynamics study of the phase transition which were compared to AI molecular dynamics simulations also performed here. In this paper we have also determined new full structural AI relaxations in the PE and FE phases, and computed the AI total phonon density of states in the PE phase. These data together with previous AI data from reference [22] and available experimental results were used to extensively test the results of the new model developed. The SM results for the effective Debye temperature show good agreement compared to AI [22] and experimental data. Moreover, we have found a very good agreement in the results of the SM and AI molecular dynamics simulations performed for the study of the phase transition in deuterated KDP (DKDP). This paper is organized as follows: in section 2 we give details about the SM developed and the AI and SM calculations performed. The results are presented in section 3. In the first subsection 3.1 we analyze the structures obtained with the present SM for the PE and FE phases and compare them with AI (vdW-DF) results as well as with the available experimental data. The SM calculations of the phonons and related properties are presented in section 3.2, which are in turn divided into two subsections. In the first subsection 3.2.1 we report some of the SM phonons at the Brillouin zone center and compare them with AI and experimental data. Besides, we also present a comparison of the total DOS in the FE and PE phases with the corresponding AI results. In the second subsection 3.2.2, we analyze the SM results for the effective Debye temperature and compare them with the corresponding AI and experimental data. In section 3.3 we present molecular dynamics simulations for the developed SM and compare them with AI results. Finally, a summary is presented in section 4.

Shell model and ab initio methods. Calculation details
The PE phase of KDP has a body-centered tetragonal structure with shorter lattice constant along the tetragonal axis. The space group is I42d or D 12 2 . In the FE phase, the crystal is 0.8% shear distorted along the [110] direction and becomes orthorhombic with space group Fdd2 or C 19 2 . The primitive cell in both phases contains 16 atoms (two formula units).
In the PE phase of KDP, the hydrogens have two equilibrium positions equidistant to the middle of the O-H-O bond and separated by a distance . Both positions are occupied by these protons with equal probability in this phase, and hence the averaged proton position is = /2 = 0. Therefore, in the AI and SM simulations for this phase, we performed full structural optimizations with the H atoms exactly at the middle of the O-H-O bonds [13], which remained in their positions after relaxations due to symmetry reasons. This hypothetical phase at 0 K is used to keep the macroscopic center inversion symmetry of I42d phase. For this phase, the phonon calculations give three unstable modes of imaginary frequency, one of which is shown in table 3. In the case of the FE phase, the protons are originally slightly displaced from the middle of the O-H-O bonds following the pattern of the FE mode. After that, the full structural SM and AI relaxations lead the system in each case to the ordered FE phase. In the SM and AI optimized-cell structural relaxations for both phases, the cell plus atomic structural parameters were allowed to relax at zero pressure.
Our starting point for the SM is the one developed in reference [29]. It contains ionic polarizabilities for the P and O ions which are described by an electronic shell with charge harmonically coupled with a spring to the core. In addition, we consider short-range shell-shell repulsive interactions arising from the wavefunction overlap between neighboring ions and long-range Coulomb interactions between cores and shells of every ion. The short-range interactions are of the Born-Mayer type, e − / , for the P-O and O-H bonds as well as of the Buckingham type, e − / − / 6 , for the K-O bonds [32,33,36]. Three-body angular potentials of the form (1/2) ( − 0 ) 2 for the covalent O-P-O and P-O-H bonds [36] and a three-body potential represented by the term ( / 12 − / 10 ) cos 4 ( − 180 • ) for the O-H-O bond are also included [29]. Considering the total charge of the ions and using the condition of charge neutrality, the model finally contains 20 adjustable parameters.
The SM simulations were carried out using the GULP code [39] which can perform structure optimizations at zero temperature as well as phonon calculations (phonon frequencies and eigenvectors at the Brillouin zone center, phonon dispersion curves and density of states, etc.) and classical molecular 43709-3  [40,41]. The plane-wave basis was expanded to an energy cutoff of 750 eV. In the calculations we use an automatic Monkhorst-Pack 5 × 5 × 5 grid sampling of the electronic Brillouin zone, which proved sufficient to achieve converged results. The calculations were carried out using the AI scheme vdW-DF that includes nonlocal van der Waals dispersion corrections [42][43][44]. The geometry optimizations were performed until the forces on every atom were smaller than 5 meV/Å. Shell model molecular dynamics (SMMD) and ab initio molecular dynamics (AIMD) simulations were performed in the NVT ensemble with the GULP and VASP programs, respectively, using a Nosé-Hoover thermostat to control the temperature. The Newton's equations of motion were integrated using the Verlet's algorithm with a typical time step of 0.3 fs (0.2 fs) in the AIMD (SMMD) simulations for an accurate integration of the electronic degrees of freedom in H-bonded systems. The AIMD (SMMD) calculations were carried out for a series of temperatures from 200 to 500 K using supercells of 128 (256) atoms subjected to periodic boundary conditions. At each temperature, a well equilibrated configuration is achieved after running at least 5 ps. Then, the system evolves further at least 5 and 15 ps to calculate the thermodynamic averages in the AIMD and SMMD calculations, respectively.
The AIMD calculations were performed with a high-accuracy cutoff energy of 600 eV and Γ-point Brillouin zone sampling, and include nonlocal dispersion corrections at the vdW-DF level. It is worth mentioning that in order to ensure that the energy is well converged we have verified that the drift produced by the Nosé-Hoover dynamics is less than 1 meV/atom within a time step of 1 ps [45].
The lattice parameters were fixed in the simulations for DKDP to those corresponding to the experimental values in the PE phase measured at + 5 K ≈ 234 K [46]. Accordingly, the ordered phase at low temperature arises in the calculations at fixed cell (NVT ensemble) as a FE phase with tetragonal lattice, which is slightly distorted with respect to the orthorhombic lattice [47]. This strategy enables the study of the phase transition using a supercell that conserves volume. We verified that all AIMD and SMMD simulations yield to average equilibrium structures which are stable up to the largest temperature considered = 500 K.
The total ab initio phonon density of states (PDOS) was derived from the phonon calculations using the finite difference (FD) method as implemented in the PHONOPY code [22,48,49]. Here, in order to compute the atomic forces for the different configurations generated by the PHONOPY FD method we used the VASP program with the vdW-DF method to account for nonlocal dispersion effects. In this case, the energy cutoff for the plane-wave basis was set to 450 eV, and a 4 × 4 × 4 grid for the Brillouin zone sampling was used. The atomic distortions were performed in a 2 × 2 × 2 supercell with 128 atoms. In this case we used a tighter tolerance in the forces of 0.5 meV/Å in order to achieve convergence in the phonon dispersion results. Once the phonon dispersion relations were obtained, the PHONOPY code allowed us to calculate the total PDOS by a Brillouin zone integration.

Structural optimizations
The model parameters of the SM were initially taken from reference [29] and were further adjusted in order to reproduce the AI (vdW-DF) results for the relaxed internal structure and lattice constants in the FE and PE phases. Special emphasis was made on reproducing the O-O distance and the parameter for the H-bonds, which were mainly controlled by adjusting the three-body potential for the O-H-O bond and the Born-Mayer O-H potential. Once the set of parameters satisfactorily reproduced the structure of both phases, we further refined the model to adjust better the Brillouin zone-center (Γ-point) phonons to the AI data. In doing this, due to the complexity of the system, special care was taken to adjust various properties at the same time without spoiling the whole fit.
In table 2, we present results of the structural parameters obtained with full AI and SM relaxations (atomic plus cell optimizations) in the PE and FE phases. These results are compared to the corresponding AI values from reference [22]. A good overall agreement between the AI (vdW-DF) and SM structural parameters results for both phases is observed in table 2. Moreover, the comparison of these results with the available experimental data is also satisfactory, with relative percentage differences of the order or less than 3% in most of the cases. The underestimation in the O-O distance in the PE phase can be mainly attributed to the static optimization for centered protons in the H-bonds [7,16], while in the actual PE phase observed at finite temperature, the protons are delocalized over two symmetric, off-center positions along the bond [2,3]. Notice that the theoretical calculations correspond to classical (infinite mass) nuclei. Actually, H-bond parameters such as the O-O, O-H and distances are especially affected by the quantum dynamics of the proton [10,22]. The inclusion of these dynamics approximately by means of PIMC calculations, leads to NQC that improve the correspondence of the AI results for these parameters [22] with the experiment. Notice that the present SM calculations overestimate the relaxed H-bond geometry in the FE phase in comparison with the experiment, e.g., the O-O distance and the parameter (see table 2). Therefore, due to the fact that NQC produce a contraction in the H-bonds [22], we speculate that if NQC were applied to the SM calculations, the agreement with the experiment for the H-bond geometry could be improved. Table 2. SM and AI (vdW-DF) results of the lattice and internal structure parameters for the FE and PE phases of KDP. Also shown is the experimental data of reference [46]. Distances in angstroms and angles in degrees. We also show in parenthesis the relative percentage differences with available experimental data.

-point phonons and density of states
We have checked the capability of the developed SM to reproduce vibrational properties of the system obtained from first principles calculations. To this end, we compare some of the SM Brillouin zone-center phonons with the corresponding AI (vdW-DF) and experimental data in table 3. The table is organized in view of the correspondence between representations in the PE and FE phases. [22] Therefore, frequencies in the same row correspond to qualitatively similar patterns of motion for both phases, although the correspondence should be taken with caution because in some cases there were ambiguities in the assignement. In table 3, we show the SM results for the zone-center phonons of symmetry A 1 , A 2 , B 1 , and B 2 in the PE phase, which are divided into the two subspaces, (A 1 + B 2 ) and (A 2 + B 1 ). In the right-hand panel we report the corresponding SM frequencies of the A 1 and A 2 modes in the FE phase. Table 3. SM results of the Γ-point phonons of symmetries 1 , 2 , 2 , 1 (PE phase) and 1 , 2 (FE phase) for KDP. Also shown are the AI (VASP) results obtained with the exchange-correlation functional vdW-DF from reference [22] and the experimental results of references [50][51][52] for the PE phase and references [21,50,53] for the FE phase. The phonon frequencies are shown in cm −1 . According to the calculated eigenvectors, the following classification is shown in the This classification can be applied to the Γ-point phonon results as well as to the phonon bands obtained from the phonon density of states (DOS) results which are shown below.
The SM results for the chosen Brillouin zone-center phonons are compared in table 3 with the AI (vdW-DF) results and with the available experimental data. A good overall agreement is observed between the SM frequencies and the corresponding AI [22] and experimental data, as well as with other AI phonon calculations for KDP [23]. Despite the fact that we are analysing a hypothetical PE phase at 0 K where the protons are fixed centered in the H-bonds, and that the low frequency phonons and modes involving hydrogen atoms may be affected by the three unstable modes arising from the used approximation, it is instructive to qualitatively compare the changes observed in the phonon spectrum between both PE and FE phases. For instance, one of the unstable modes in the PE phase of B 2 symmetry, which is identified by its imaginary frequency in table 3, corresponds to a similar one found with the AI calculations. Such B 2 mode is the FE soft mode in the PE phase, which becomes stable in the FE phase as a high-frequency HS phonon (see table 3).
When the phase changes from PE to FE some modes soften while others stiffen and, with a few exceptions, the SM phonons follow the same trend as the AI phonons. In particular, one of the most significant frequency changes between both phases is observed for the PO 4 -rotational A 1 mode obtained at 648 cm −1 in the PE phase for the SM calculation, which softens to 274 cm −1 in the FE phase (see  table 3). Similarly, this mode softens strongly in the AI calculation. It is demostrated that this phonon couples significantly with the FE soft mode and that this interaction becomes stronger as the amplitude of the latter is increased [22]. The FE soft-mode picture for the phase transition in H-bonded ferroelectrics is consistent with recent NMR experiments which show the importance of a displacive component near the critical temperature [54,55].
The total SM-DOSs in the FE and PE phases are plotted in figure 1 and compared to the corresponding AI (vdW-DF) results, showing a general good agreement between both calculations for both phases. Notice that the SM and AI (vdW-DF) unstable phonon bands observed in the PE phase at imaginary frequencies (plotted as negative frequencies in the right-hand panel of figure 1), are related to the unstable FE zone-center phonon reported in table 3. From these plots and the classification of the phononic displacement patterns made in table 3, we can identify the different types of bands present in the spectra for the FE phase. The highest-frequency band corresponds to the HS band centered at ≈ 2550 cm −1 for the SM calculation which is in very good agreement with the AI result for the HS band centered at ≈ 2600 cm −1 as shown in the left-hand panel of figure 1. Particularly, the bimodal shape of this band and the spectral weight are correctly described by the SM calculation, although the position of the highest-frequency peak of the SM bimodal distribution is a bit shifted towards lower frequencies with respect to the AI result.
The SM result for the HB band is also in remarkable agreement with the corresponding AI result, displaying similar peak positions (≈ 1300 cm −1 ) and band extensions in both calculations (see figure 1). The large gap between the HB and HS bands observed in the FE phase in the AI calculation is also qualitatively reproduced by the SM results as shown in figure 1. We also find that the SM DOS result is in qualitative agreement with another AI calculation for the FE phase [38].
The SM lowest-frequency band in the FE phase centered at ≈ 150 cm −1 is also in very good agreement with the corresponding AI data (see figure 1). This is a lattice external (ET+ER) band that extends up to ≈ 300 cm −1 , which is of K and O character mainly.
The intermediate frequency region of the DOSs spectra for the FE phase consists of two bands of internal molecular (IM) phonons involving primarily the PO 4 tetrahedron. The high-frequency band extends from ≈ 700 to ≈ 1200 cm −1 in the SM calculation in close agreement with the AI result [22] (see figure 1), and consists of the modes involving mainly the stretching of the P-O bonds (IMS band). The other IM band is observed in the SM calculation in the region ≈ 400-700 cm −1 , and involves the O-P-O bending modes (IMB band) as shown in figure 1. The SM IMB band is associated with two peaks similarly to the AI band, but appears shifted ≈ 150 cm −1 towards higher frequencies than the corresponding AI band. Unfortunately, any attempt to improve the fitting of these bands by modifying different parameters led to a worsening of the overall fit.
Infrared measurements of the imaginary dielectric function at 80 K in KDP with the electric field polarized along the FE axis show well-defined peaks up to ≈ 1400 cm −1 [21]. The main resonances  are: three in the region (0-700) cm −1 , centered at ≈ 200, 300 and 500 cm −1 , and the other three in the region (700-1400) cm −1 , centered at ≈ 900, 1000 and 1300 cm −1 . These resonances also appear in the infrared spectrum obtained with the electric field polarized in the plane perpendicular to the FE axis [21]. Our SM phononic band structure in the FE phase is in close agreement with the infrared spectra and enables us to identify the main experimental peaks. For instance, the experimental peaks centered at ≈ 200 and 1300 cm −1 are in remarkable agreement with the SM lattice external (ET+ER) and HB bands centered at ≈ 150 and 1300 cm −1 , respectively. On the other hand, the bands observed experimentally at ≈ 900 and 1000 cm −1 can be assigned to our SM bands centered at ≈ 800 and 1000 cm −1 which are of IMS+HB character (see the left-hand panel of figure 1 and table 3). The experimental peak observed at ≈ 500 cm −1 should be related to the SM IMB+HB band centered at nearly the same frecuency, as shown in the left-hand panel of figure 1. These assignments coincide with those made by the AI calculation of reference [22] (see also the similarities in the band distributions for the SM and AI calculations in the left-hand panel of figure 1). On the other hand, the experimental peak at 300 cm −1 can be ascribed to a mixed IMB+lattice band centered at ≈ 250 cm −1 observed in the SM and AI results which mainly involve O, K and P displacements (see the left-hand panel of figure 1), or it can be related to the welldefined AI IMB band centered at ≈ 350 cm −1 . Notice that in the latter case, the SM IMB band cannot account for it because it is shifted towards higher frequencies.

Effective Debye temperature
With the aid of the developed model, we have calculated the effective (or equivalent) Debye temperature in the FE phase, Θ D ( ), as an integral property of the total DOS [22,35,56,57]. This temperature is associated to the value of the specific heat at each temperature [56]. The SM results for Θ D ( ) are plotted in figure 2 as a function of and compared to the corresponding AI data of reference [22] and the available experimental data [58][59][60]. A very good overall agreement of the SM results with the AI and experimental data is obtained, especially at temperatures 10 K, as can be judged from figure 2. The low-temperature Debye region, where Θ D ( ) is temperature independent and the Debye model is strictly valid, is observed at 3 K and 4 K for both SM and AI calculations, respectively. This is in qualitative agreement with the experimental Debye region observed at 5 K considering the results of references [59,60], as shown in figure 2. As → 0, Θ D ( ) takes the value of ≈ 450 K for the SM which is in qualitative agreement with the corresponding experimental and AI values, ≈ 380 K and ≈ 350 K, respectively. The SM curve for Θ D ( ) reaches a minimum value at ≈ 13 K in very good agreement with the corresponding minima observed in the AI and experimental results, as shown in 43709-8 figure 2. However, the large drop observed in the SM result for Θ D ( ) from the Debye region to the minimum value is overestimated in comparison with the corresponding falls in the AI and experimental curves. The observed drop is a typical feature in the effective Debye-temperature profile which is produced by strong hybridizations of the acoustic and low-frequency optical branches towards the Brillouin zone boundary [57].

Molecular dynamics simulations
Finally we have performed classical-nuclei SMMD and AIMD simulations to study the phase transition in DKDP, which is expected to have less quantum effects than KDP. To track the transition we must first define the order parameter. To accomplish this, we assign to each proton a positive (+) or negative (−) displacement from the instantaneous O-H-O bond center. The values of are considered positive if the proton displacement direction in the particular O-H-O bond coincides with that observed for the global FE mode which causes polarization in the + direction, or negative otherwise [22]. It is worth noting here that the coordinated proton displacements around each phosphate following the FE mode pattern are strongly correlated with the local ionic plus electronic polarization along the direction [7,16]. The order parameter is then computed as the spatial and time average of all proton displacements . We define the critical temperature as the temperature at which the order parameter falls to half of its maximum value in the ordered phase. Figure 3 depicts vs. for the AI and SM simulations. Simulations at different temperatures show that the system has a FE ordering up to a critical temperature . For larger temperatures, the average total order parameter vanishes and the PE phase arises. We observe an excellent agreement for the values of obtained with the SMMD and AIMD simulations for DKDP, SM (DKDP) ≈ 360 K and AI (DKDP) ≈ 365 K, respectively. Notice that the present molecular dynamics calculations do not consider the quantum nature of the protons (and hence the possibility of proton tunneling), and for this reason the theoretical values of are expected to be larger than the corresponding experimental value for DKDP, which is exp (DKDP) ≈ 229 K [2,3].

Summary
We have developed a SM fitted to AI results that include non-local vdW corrections which are important to properly describe H-bond geometries and global proton-transfer energy barriers. The development  of this model is a first step for a future dynamical study in large-size systems of the elusive nature of the ferroelectric phase transition in KDP and other compounds of the H-bonded FE family. The SM was extensively tested by comparing structural and vibrational results to the corresponding first principles and experimental data. The relaxed lattice and atomic structural parameters in the FE and PE phases are in very good correspondence with the respective vdW-corrected AI data as well as with the available experimental data.
Regarding the vibrational properties, the SM Γ-point phonons obtained for the FE and PE phases are in good overall agreement with the corresponding AI results and the available experimental data. The FE mode at the Brillouin zone center calculated with the developed SM is unstable in the PE phase, in agreement with the AI results. The total SM-DOSs in the FE and PE phases are also in very good accordance with the corresponding AI data suggesting that the model represents well the vibrational properties of KDP throughout the whole Brillouin zone. On the other hand, the obtained SM IMB band is in qualitative agreement with the corresponding AI band, but is shifted ≈ 150 cm −1 towards higher frequencies. The SM and AI bands in the FE phase are also in very good agreement with the resonances observed in infrared measurements of the imaginary dielectric function at 80 K in KDP, which enables us to associate the measured peaks with the corresponding phononic displacement patterns. Using the phonon DOS obtained in the FE phase with the developed SM, we have computed the effective Debye temperature which is in very good agreement with the corresponding AI and experimental data. Finally, we have carried out SMMD simulations for classical nuclei and found a FE-PE phase transition for DKDP at ≈ 360 K which is in very good agreement with the corresponding AIMD simulations that include nonlocal dispersion effects via the vdW-DF approach.